Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
RooDataHist.cxx
Go to the documentation of this file.
1/*****************************************************************************
2 * Project: RooFit *
3 * Package: RooFitCore *
4 * @(#)root/roofitcore:$Id$
5 * Authors: *
6 * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu *
7 * DK, David Kirkby, UC Irvine, dkirkby@uci.edu *
8 * *
9 * Copyright (c) 2000-2005, Regents of the University of California *
10 * and Stanford University. All rights reserved. *
11 * *
12 * Redistribution and use in source and binary forms, *
13 * with or without modification, are permitted according to the terms *
14 * listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
15 *****************************************************************************/
16
17/**
18\file RooDataHist.cxx
19\class RooDataHist
20\ingroup Roofitcore
21
22Container class to hold N-dimensional binned data. Each bin's central
23coordinates in N-dimensional space are represented by a RooArgSet containing RooRealVar, RooCategory
24or RooStringVar objects, thus data can be binned in real and/or discrete dimensions.
25
26There is an unbinned equivalent, RooDataSet.
27
28### Inspecting a datahist
29Inspect a datahist using Print() to get the coordinates and `weight()` to get the bin contents:
30```
31datahist->Print("V");
32datahist->get(0)->Print("V"); std::cout << "w=" << datahist->weight(0) << std::endl;
33datahist->get(1)->Print("V"); std::cout << "w=" << datahist->weight(1) << std::endl;
34...
35```
36
37### Plotting data.
38See RooAbsData::plotOn().
39
40### Creating a datahist using RDataFrame
41\see RooAbsDataHelper, rf408_RDataFrameToRooFit.C
42
43**/
44
45#include "RooDataHist.h"
46
47#include "Riostream.h"
48#include "RooMsgService.h"
50#include "RooAbsLValue.h"
51#include "RooArgList.h"
52#include "RooRealVar.h"
53#include "RooMath.h"
54#include "RooBinning.h"
55#include "RooPlot.h"
56#include "RooHistError.h"
57#include "RooCategory.h"
58#include "RooCmdConfig.h"
59#include "RooLinkedListIter.h"
60#include "RooTreeDataStore.h"
61#include "RooVectorDataStore.h"
62#include "RooTrace.h"
63#include "RooFormulaVar.h"
64#include "RooFormula.h"
65#include "RooUniformBinning.h"
66
67#include "RooFitImplHelpers.h"
68
69#include <ROOT/RSpan.hxx>
70#include <ROOT/StringUtils.hxx>
71
72#include "TAxis.h"
73#include "TH1.h"
74#include "TTree.h"
75#include "TBuffer.h"
76#include "TMath.h"
77#include "Math/Util.h"
78
79using std::cout, std::endl, std::string, std::ostream;
80
82
83
84////////////////////////////////////////////////////////////////////////////////
85/// Default constructor
86
88{
90}
91
92
93std::unique_ptr<RooAbsDataStore>
95{
97 ? static_cast<std::unique_ptr<RooAbsDataStore>>(std::make_unique<RooTreeDataStore>(name, title, vars))
98 : static_cast<std::unique_ptr<RooAbsDataStore>>(std::make_unique<RooVectorDataStore>(name, title, vars));
99}
100
101
102////////////////////////////////////////////////////////////////////////////////
103/// Constructor of an empty data hist from a RooArgSet defining the dimensions
104/// of the data space. The range and number of bins in each dimensions are taken
105/// from getMin()getMax(),getBins() of each RooAbsArg representing that
106/// dimension.
107///
108/// For real dimensions, the fit range and number of bins can be set independently
109/// of the plot range and number of bins, but it is advisable to keep the
110/// ratio of the plot bin width and the fit bin width an integer value.
111/// For category dimensions, the fit ranges always comprises all defined states
112/// and each state is always has its individual bin
113///
114/// To effectively bin real dimensions with variable bin sizes,
115/// construct a RooThresholdCategory of the real dimension to be binned variably.
116/// Set the thresholds at the desired bin boundaries, and construct the
117/// data hist as a function of the threshold category instead of the real variable.
118RooDataHist::RooDataHist(RooStringView name, RooStringView title, const RooArgSet& vars, const char* binningName) :
119 RooAbsData(name,title,vars)
120{
121 // Initialize datastore
123
124 initialize(binningName) ;
125
127
128 appendToDir(this,true) ;
130}
131
132
133
134////////////////////////////////////////////////////////////////////////////////
135/// Constructor of a data hist from an existing data collection (binned or unbinned)
136/// The RooArgSet 'vars' defines the dimensions of the histogram.
137/// The range and number of bins in each dimensions are taken
138/// from getMin(), getMax(), getBins() of each argument passed.
139///
140/// For real dimensions, the fit range and number of bins can be set independently
141/// of the plot range and number of bins, but it is advisable to keep the
142/// ratio of the plot bin width and the fit bin width an integer value.
143/// For category dimensions, the fit ranges always comprises all defined states
144/// and each state is always has its individual bin
145///
146/// To effectively bin real dimensions with variable bin sizes,
147/// construct a RooThresholdCategory of the real dimension to be binned variably.
148/// Set the thresholds at the desired bin boundaries, and construct the
149/// data hist as a function of the threshold category instead of the real variable.
150///
151/// If the constructed data hist has less dimensions that in source data collection,
152/// all missing dimensions will be projected.
153
155 RooDataHist(name,title,vars)
156{
157 add(data,static_cast<const RooFormulaVar*>(nullptr),wgt);
158}
159
160
161
162////////////////////////////////////////////////////////////////////////////////
163/// Constructor of a data hist from a map of TH1,TH2 or TH3 that are collated into a x+1 dimensional
164/// RooDataHist where the added dimension is a category that labels the input source as defined
165/// in the histMap argument. The state names used in histMap must correspond to predefined states
166/// 'indexCat'
167///
168/// The RooArgList 'vars' defines the dimensions of the histogram.
169/// The ranges and number of bins are taken from the input histogram and must be the same in all histograms
170
172 std::map<string,TH1*> histMap, double wgt) :
173 RooAbsData(name,title,RooArgSet(vars,&indexCat))
174{
175 // Initialize datastore
177
178 importTH1Set(vars, indexCat, histMap, wgt, false) ;
179
182}
183
184
185
186////////////////////////////////////////////////////////////////////////////////
187/// Constructor of a data hist from a map of RooDataHists that are collated into a x+1 dimensional
188/// RooDataHist where the added dimension is a category that labels the input source as defined
189/// in the histMap argument. The state names used in histMap must correspond to predefined states
190/// 'indexCat'
191///
192/// The RooArgList 'vars' defines the dimensions of the histogram.
193/// The ranges and number of bins are taken from the input histogram and must be the same in all histograms
194
196 std::map<string,RooDataHist*> dhistMap, double wgt) :
197 RooAbsData(name,title,RooArgSet(vars,&indexCat))
198{
199 // Initialize datastore
201
202 importDHistSet(vars, indexCat, dhistMap, wgt) ;
203
206}
207
208
209
210////////////////////////////////////////////////////////////////////////////////
211/// Constructor of a data hist from an TH1,TH2 or TH3
212/// The RooArgSet 'vars' defines the dimensions of the histogram. The ranges
213/// and number of bins are taken from the input histogram, and the corresponding
214/// values are set accordingly on the arguments in 'vars'
215
216RooDataHist::RooDataHist(RooStringView name, RooStringView title, const RooArgList& vars, const TH1* hist, double wgt) :
217 RooAbsData(name,title,vars)
218{
219 // Initialize datastore
221
222 // Check consistency in number of dimensions
223 if (int(vars.size()) != hist->GetDimension()) {
224 std::stringstream errorMsgStream;
225 errorMsgStream << "RooDataHist::ctor(" << GetName() << ") ERROR: dimension of input histogram must match "
226 << "number of dimension variables";
227 const std::string errorMsg = errorMsgStream.str();
228 coutE(InputArguments) << errorMsg << std::endl;
229 throw std::invalid_argument(errorMsg);
230 }
231
232 importTH1(vars,*hist,wgt, false) ;
233
236}
237
238
239
240////////////////////////////////////////////////////////////////////////////////
241/// Constructor of a binned dataset from a RooArgSet defining the dimensions
242/// of the data space. The range and number of bins in each dimensions are taken
243/// from getMin() getMax(),getBins() of each RooAbsArg representing that
244/// dimension.
245///
246/// <table>
247/// <tr><th> Optional Argument <th> Effect
248/// <tr><td> Import(TH1&, bool impDens) <td> Import contents of the given TH1/2/3 into this binned dataset. The
249/// ranges and binning of the binned dataset are automatically adjusted to
250/// match those of the imported histogram.
251///
252/// Please note: for TH1& with unequal binning _only_,
253/// you should decide if you want to import the absolute bin content,
254/// or the bin content expressed as density. The latter is default and will
255/// result in the same histogram as the original TH1. For certain types of
256/// bin contents (containing efficiencies, asymmetries, or ratio is general)
257/// you should import the absolute value and set impDens to false
258///
259///
260/// <tr><td> Weight(double) <td> Apply given weight factor when importing histograms
261///
262/// <tr><td> Index(RooCategory&) <td> Prepare import of multiple TH1/1/2/3 into a N+1 dimensional RooDataHist
263/// where the extra discrete dimension labels the source of the imported histogram
264/// If the index category defines states for which no histogram is be imported
265/// the corresponding bins will be left empty.
266///
267/// <tr><td> Import(const char*, TH1&) <td> Import a THx to be associated with the given state name of the index category
268/// specified in Index(). If the given state name is not yet defined in the index
269/// category it will be added on the fly. The import command can be specified
270/// multiple times.
271/// <tr><td> Import(map<string,TH1*>&) <td> As above, but allows specification of many imports in a single operation
272/// <tr><td> `GlobalObservables(const RooArgSet&)` <td> Define the set of global observables to be stored in this RooDataHist.
273/// A snapshot of the passed RooArgSet is stored, meaning the values wont't change unexpectedly.
274/// </table>
275///
276
277RooDataHist::RooDataHist(RooStringView name, RooStringView title, const RooArgList& vars, const RooCmdArg& arg1, const RooCmdArg& arg2, const RooCmdArg& arg3,
278 const RooCmdArg& arg4,const RooCmdArg& arg5,const RooCmdArg& arg6,const RooCmdArg& arg7,const RooCmdArg& arg8) :
279 RooAbsData(name,title,RooArgSet(vars,static_cast<RooAbsArg*>(RooCmdConfig::decodeObjOnTheFly("RooDataHist::RooDataHist", "IndexCat",0,nullptr,arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8))))
280{
281 // Initialize datastore
283
284 // Define configuration for this method
285 RooCmdConfig pc("RooDataHist::ctor(" + std::string(GetName()) + ")");
286 pc.defineObject("impHist","ImportHisto",0) ;
287 pc.defineInt("impDens","ImportHisto",0) ;
288 pc.defineObject("indexCat","IndexCat",0) ;
289 pc.defineObject("impSliceData","ImportDataSlice",0,nullptr,true) ; // array
290 pc.defineString("impSliceState","ImportDataSlice",0,"",true) ; // array
291 pc.defineDouble("weight","Weight",0,1) ;
292 pc.defineObject("dummy1","ImportDataSliceMany",0) ;
293 pc.defineSet("glObs","GlobalObservables",0,nullptr) ;
294 pc.defineMutex("ImportHisto","ImportDataSlice");
295 pc.defineDependency("ImportDataSlice","IndexCat") ;
296
298 l.Add((TObject*)&arg1) ; l.Add((TObject*)&arg2) ;
299 l.Add((TObject*)&arg3) ; l.Add((TObject*)&arg4) ;
300 l.Add((TObject*)&arg5) ; l.Add((TObject*)&arg6) ;
301 l.Add((TObject*)&arg7) ; l.Add((TObject*)&arg8) ;
302
303 // Process & check varargs
304 pc.process(l) ;
305 if (!pc.ok(true)) {
306 throw std::invalid_argument("Invalid command arguments passed to RooDataHist constructor!");
307 }
308
309 if(pc.getSet("glObs")) setGlobalObservables(*pc.getSet("glObs"));
310
311 TH1* impHist = static_cast<TH1*>(pc.getObject("impHist")) ;
312 bool impDens = pc.getInt("impDens") ;
313 double initWgt = pc.getDouble("weight") ;
314 RooCategory* indexCat = static_cast<RooCategory*>(pc.getObject("indexCat")) ;
315 const char* impSliceNames = pc.getString("impSliceState","",true) ;
316 const RooLinkedList& impSliceHistos = pc.getObjectList("impSliceData") ;
317
318
319 if (impHist) {
320
321 // Initialize importing contents from TH1
322 importTH1(vars,*impHist,initWgt, impDens) ;
323
324 } else if (indexCat) {
325
326
327 // Initialize importing mapped set of RooDataHists and TH1s
328 std::map<std::string,RooDataHist*> dmap ;
329 std::map<std::string,TH1*> hmap ;
330 auto hiter = impSliceHistos.begin() ;
331 for (const auto& token : ROOT::Split(impSliceNames, ",", /*skipEmpty=*/true)) {
332 if(auto dHist = dynamic_cast<RooDataHist*>(*hiter)) {
333 dmap[token] = dHist;
334 }
335 if(auto hHist = dynamic_cast<TH1*>(*hiter)) {
336 hmap[token] = hHist;
337 }
338 ++hiter;
339 }
340 if(!dmap.empty() && !hmap.empty()) {
341 std::stringstream errorMsgStream;
342 errorMsgStream << "RooDataHist::ctor(" << GetName() << ") ERROR: you can't import mix of TH1 and RooDataHist";
343 const std::string errorMsg = errorMsgStream.str();
344 coutE(InputArguments) << errorMsg << std::endl;
345 throw std::invalid_argument(errorMsg);
346 }
347 if (!dmap.empty()) {
348 importDHistSet(vars,*indexCat,dmap,initWgt);
349 }
350 if (!hmap.empty()) {
351 importTH1Set(vars,*indexCat,hmap,initWgt,false);
352 }
353
354
355 } else {
356
357 // Initialize empty
358 initialize() ;
359 appendToDir(this,true) ;
360
361 }
362
365
366}
367
368
369
370
371////////////////////////////////////////////////////////////////////////////////
372/// Import data from given TH1/2/3 into this RooDataHist
373
374void RooDataHist::importTH1(const RooArgList& vars, const TH1& histo, double wgt, bool doDensityCorrection)
375{
376 // Adjust binning of internal observables to match that of input THx
377 Int_t offset[3]{0, 0, 0};
378 adjustBinning(vars, histo, offset) ;
379
380 // Initialize internal data structure
381 initialize() ;
382 appendToDir(this,true) ;
383
384 // Define x,y,z as 1st, 2nd and 3rd observable
385 RooRealVar* xvar = static_cast<RooRealVar*>(_vars.find(vars.at(0)->GetName())) ;
386 RooRealVar* yvar = static_cast<RooRealVar*>(vars.at(1) ? _vars.find(vars.at(1)->GetName()) : nullptr ) ;
387 RooRealVar* zvar = static_cast<RooRealVar*>(vars.at(2) ? _vars.find(vars.at(2)->GetName()) : nullptr ) ;
388
389 // Transfer contents
390 Int_t xmin(0);
391 Int_t ymin(0);
392 Int_t zmin(0);
393 RooArgSet vset(*xvar) ;
394 xmin = offset[0] ;
395 if (yvar) {
396 vset.add(*yvar) ;
397 ymin = offset[1] ;
398 }
399 if (zvar) {
400 vset.add(*zvar) ;
401 zmin = offset[2] ;
402 }
403
404 Int_t iX(0);
405 Int_t iY(0);
406 Int_t iz(0);
407 for (iX=0 ; iX < xvar->getBins() ; iX++) {
408 xvar->setBin(iX) ;
409 if (yvar) {
410 for (iY=0 ; iY < yvar->getBins() ; iY++) {
411 yvar->setBin(iY) ;
412 if (zvar) {
413 for (iz=0 ; iz < zvar->getBins() ; iz++) {
414 zvar->setBin(iz) ;
415 double bv = doDensityCorrection ? binVolume(vset) : 1;
416 add(vset,bv*histo.GetBinContent(iX+1+xmin,iY+1+ymin,iz+1+zmin)*wgt,bv*std::pow(histo.GetBinError(iX+1+xmin,iY+1+ymin,iz+1+zmin)*wgt,2)) ;
417 }
418 } else {
419 double bv = doDensityCorrection ? binVolume(vset) : 1;
420 add(vset,bv*histo.GetBinContent(iX+1+xmin,iY+1+ymin)*wgt,bv*std::pow(histo.GetBinError(iX+1+xmin,iY+1+ymin)*wgt,2)) ;
421 }
422 }
423 } else {
424 double bv = doDensityCorrection ? binVolume(vset) : 1 ;
425 add(vset,bv*histo.GetBinContent(iX+1+xmin)*wgt,bv*std::pow(histo.GetBinError(iX+1+xmin)*wgt,2)) ;
426 }
427 }
428
429}
430
431namespace {
432bool checkConsistentAxes(const TH1* first, const TH1* second) {
433 return first->GetDimension() == second->GetDimension()
434 && first->GetNbinsX() == second->GetNbinsX()
435 && first->GetNbinsY() == second->GetNbinsY()
436 && first->GetNbinsZ() == second->GetNbinsZ()
437 && first->GetXaxis()->GetXmin() == second->GetXaxis()->GetXmin()
438 && first->GetXaxis()->GetXmax() == second->GetXaxis()->GetXmax()
439 && (first->GetNbinsY() == 1 || (first->GetYaxis()->GetXmin() == second->GetYaxis()->GetXmin()
440 && first->GetYaxis()->GetXmax() == second->GetYaxis()->GetXmax() ) )
441 && (first->GetNbinsZ() == 1 || (first->GetZaxis()->GetXmin() == second->GetZaxis()->GetXmin()
442 && first->GetZaxis()->GetXmax() == second->GetZaxis()->GetXmax() ) );
443}
444
445bool hasConsistentLayoutAndBinning(RooDataHist const& h1, RooDataHist const& h2) {
446
447 // Relative tolerance for bin boundary comparison
448 constexpr double tolerance = 1e-6;
449
450 auto const& vars1 = *h1.get();
451 auto const& vars2 = *h2.get();
452
453 // Check if number of variables and names is consistent
454 if(!vars1.hasSameLayout(vars2)) {
455 return false;
456 }
457
458 for(std::size_t iVar = 0; iVar < vars1.size(); ++iVar) {
459 auto * var1 = dynamic_cast<RooRealVar*>(vars1[iVar]);
460 auto * var2 = dynamic_cast<RooRealVar*>(vars2[iVar]);
461
462 // Check if variables are consistently real-valued
463 if((!var1 && var2) || (var1 && !var2)) return false;
464
465 // Not a real-valued variable
466 if(!var1) continue;
467
468 // Now check the binning
469 auto const& bng1 = var1->getBinning();
470 auto const& bng2 = var2->getBinning();
471
472 // Compare bin numbers
473 if(bng1.numBins() != bng2.numBins()) return false;
474
475 std::size_t nBins = bng1.numBins();
476
477 // Compare bin boundaries
478 for(std::size_t iBin = 0; iBin < nBins; ++iBin) {
479 double v1 = bng1.binLow(iBin);
480 double v2 = bng2.binLow(iBin);
481 if(std::abs((v1 - v2) / v1) > tolerance) return false;
482 }
483 double v1 = bng1.binHigh(nBins - 1);
484 double v2 = bng2.binHigh(nBins - 1);
485 if(std::abs((v1 - v2) / v1) > tolerance) return false;
486 }
487 return true;
488}
489}
490
491
492////////////////////////////////////////////////////////////////////////////////
493/// Import data from given set of TH1/2/3 into this RooDataHist. The category indexCat labels the sources
494/// in the constructed RooDataHist. The stl map provides the mapping between the indexCat state labels
495/// and the import source
496
497void RooDataHist::importTH1Set(const RooArgList& vars, RooCategory& indexCat, std::map<string,TH1*> hmap, double wgt, bool doDensityCorrection)
498{
499 RooCategory* icat = static_cast<RooCategory*>(_vars.find(indexCat.GetName())) ;
500
501 TH1* histo(nullptr) ;
502 bool init(false) ;
503 for (const auto& hiter : hmap) {
504 // Store pointer to first histogram from which binning specification will be taken
505 if (!histo) {
506 histo = hiter.second;
507 } else {
508 if (!checkConsistentAxes(histo, hiter.second)) {
509 coutE(InputArguments) << "Axes of histogram " << hiter.second->GetName() << " are not consistent with first processed "
510 << "histogram " << histo->GetName() << std::endl;
511 throw std::invalid_argument("Axes of inputs for RooDataHist are inconsistent");
512 }
513 }
514 // Define state labels in index category (both in provided indexCat and in internal copy in dataset)
515 if (!indexCat.hasLabel(hiter.first)) {
516 indexCat.defineType(hiter.first) ;
517 coutI(InputArguments) << "RooDataHist::importTH1Set(" << GetName() << ") defining state \"" << hiter.first << "\" in index category " << indexCat.GetName() << endl ;
518 }
519 if (!icat->hasLabel(hiter.first)) {
520 icat->defineType(hiter.first) ;
521 }
522 }
523
524 // Check consistency in number of dimensions
525 if (histo && int(vars.size()) != histo->GetDimension()) {
526 coutE(InputArguments) << "RooDataHist::importTH1Set(" << GetName() << "): dimension of input histogram must match "
527 << "number of continuous variables" << endl ;
528 throw std::invalid_argument("Inputs histograms for RooDataHist are not compatible with dimensions of variables.");
529 }
530
531 // Copy bins and ranges from THx to dimension observables
532 Int_t offset[3] ;
533 adjustBinning(vars,*histo,offset) ;
534
535 // Initialize internal data structure
536 if (!init) {
537 initialize() ;
538 appendToDir(this,true) ;
539 init = true ;
540 }
541
542 // Define x,y,z as 1st, 2nd and 3rd observable
543 RooRealVar* xvar = static_cast<RooRealVar*>(_vars.find(vars.at(0)->GetName())) ;
544 RooRealVar* yvar = static_cast<RooRealVar*>(vars.at(1) ? _vars.find(vars.at(1)->GetName()) : nullptr ) ;
545 RooRealVar* zvar = static_cast<RooRealVar*>(vars.at(2) ? _vars.find(vars.at(2)->GetName()) : nullptr ) ;
546
547 // Transfer contents
548 Int_t xmin(0);
549 Int_t ymin(0);
550 Int_t zmin(0);
551 RooArgSet vset(*xvar) ;
552 double volume = xvar->getMax()-xvar->getMin() ;
553 xmin = offset[0] ;
554 if (yvar) {
555 vset.add(*yvar) ;
556 ymin = offset[1] ;
557 volume *= (yvar->getMax()-yvar->getMin()) ;
558 }
559 if (zvar) {
560 vset.add(*zvar) ;
561 zmin = offset[2] ;
562 volume *= (zvar->getMax()-zvar->getMin()) ;
563 }
564 double avgBV = volume / numEntries() ;
565
566 Int_t ic(0);
567 Int_t iX(0);
568 Int_t iY(0);
569 Int_t iz(0);
570 for (ic=0 ; ic < icat->numBins(nullptr) ; ic++) {
571 icat->setBin(ic) ;
572 histo = hmap[icat->getCurrentLabel()] ;
573 for (iX=0 ; iX < xvar->getBins() ; iX++) {
574 xvar->setBin(iX) ;
575 if (yvar) {
576 for (iY=0 ; iY < yvar->getBins() ; iY++) {
577 yvar->setBin(iY) ;
578 if (zvar) {
579 for (iz=0 ; iz < zvar->getBins() ; iz++) {
580 zvar->setBin(iz) ;
581 double bv = doDensityCorrection ? binVolume(vset)/avgBV : 1;
582 add(vset,bv*histo->GetBinContent(iX+1+xmin,iY+1+ymin,iz+1+zmin)*wgt,bv*std::pow(histo->GetBinError(iX+1+xmin,iY+1+ymin,iz+1+zmin)*wgt,2)) ;
583 }
584 } else {
585 double bv = doDensityCorrection ? binVolume(vset)/avgBV : 1;
586 add(vset,bv*histo->GetBinContent(iX+1+xmin,iY+1+ymin)*wgt,bv*std::pow(histo->GetBinError(iX+1+xmin,iY+1+ymin)*wgt,2)) ;
587 }
588 }
589 } else {
590 double bv = doDensityCorrection ? binVolume(vset)/avgBV : 1;
591 add(vset,bv*histo->GetBinContent(iX+1+xmin)*wgt,bv*std::pow(histo->GetBinError(iX+1+xmin)*wgt,2)) ;
592 }
593 }
594 }
595
596}
597
598
599
600////////////////////////////////////////////////////////////////////////////////
601/// Import data from given set of TH1/2/3 into this RooDataHist. The category indexCat labels the sources
602/// in the constructed RooDataHist. The stl map provides the mapping between the indexCat state labels
603/// and the import source
604
605void RooDataHist::importDHistSet(const RooArgList & /*vars*/, RooCategory &indexCat,
606 std::map<std::string, RooDataHist *> dmap, double initWgt)
607{
608 auto *icat = static_cast<RooCategory *>(_vars.find(indexCat.GetName()));
609
610 RooDataHist *dhistForBinning = nullptr;
611
612 for (const auto &diter : dmap) {
613
614 std::string const &label = diter.first;
615 RooDataHist *dhist = diter.second;
616
617 if (!dhistForBinning) {
618 dhistForBinning = dhist;
619 } else {
620 if (!hasConsistentLayoutAndBinning(*dhistForBinning, *dhist)) {
621 coutE(InputArguments) << "Layout or binning of histogram " << dhist->GetName()
622 << " is not consistent with first processed "
623 << "histogram " << dhistForBinning->GetName() << std::endl;
624 throw std::invalid_argument("Layout or binning of inputs for RooDataHist is inconsistent");
625 }
626 }
627
628 // Define state labels in index category (both in provided indexCat and in internal copy in dataset)
629 if (!indexCat.hasLabel(label)) {
630 indexCat.defineType(label);
631 coutI(InputArguments) << "RooDataHist::importDHistSet(" << GetName() << ") defining state \"" << label
632 << "\" in index category " << indexCat.GetName() << endl;
633 }
634 if (!icat->hasLabel(label)) {
635 icat->defineType(label);
636 }
637 }
638
639 // adjust the binning of the created histogram
640 for (auto *theirVar : dynamic_range_cast<RooRealVar *>(dhistForBinning->_vars)) {
641 auto *ourVar = dynamic_cast<RooRealVar *>(_vars.find(theirVar->GetName()));
642 if (!theirVar || !ourVar)
643 continue;
644 ourVar->setBinning(theirVar->getBinning());
645 }
646
647 initialize();
648 appendToDir(this, true);
649
650 for (const auto &diter : dmap) {
651 std::string const &label = diter.first;
652 RooDataHist *dhist = diter.second;
653
654 icat->setLabel(label.c_str());
655
656 // Transfer contents
657 for (Int_t i = 0; i < dhist->numEntries(); i++) {
658 _vars.assign(*dhist->get(i));
659 add(_vars, dhist->weight() * initWgt, pow(dhist->weightError(SumW2), 2));
660 }
661 }
662}
663
664////////////////////////////////////////////////////////////////////////////////
665/// Helper doing the actual work of adjustBinning().
666
667void RooDataHist::_adjustBinning(RooRealVar &theirVar, const TAxis &axis,
668 RooRealVar *ourVar, Int_t *offset)
669{
670 const std::string ourVarName(ourVar->GetName() ? ourVar->GetName() : "");
671 const std::string ownName(GetName() ? GetName() : "");
672 // RooRealVar is derived from RooAbsRealLValue which is itself
673 // derived from RooAbsReal and a virtual class RooAbsLValue
674 // supplying setter functions, check if ourVar is indeed derived
675 // as real
676 if (!dynamic_cast<RooAbsReal *>(ourVar)) {
677 coutE(InputArguments) << "RooDataHist::adjustBinning(" << ownName << ") ERROR: dimension " << ourVarName
678 << " must be real\n";
679 throw std::logic_error("Incorrect type object (" + ourVarName +
680 ") passed as argument to RooDataHist::_adjustBinning. Please report this issue.");
681 }
682
683 const double xlo = theirVar.getMin();
684 const double xhi = theirVar.getMax();
685
686 const bool isUniform = !axis.GetXbins()->GetArray();
687 std::unique_ptr<RooAbsBinning> xbins;
688
689 if (!isUniform) {
690 xbins = std::make_unique<RooBinning>(axis.GetNbins(), axis.GetXbins()->GetArray());
691 } else {
692 xbins = std::make_unique<RooUniformBinning>(axis.GetXmin(), axis.GetXmax(), axis.GetNbins());
693 }
694
695 const double tolerance = 1e-6 * xbins->averageBinWidth();
696
697 // Adjust xlo/xhi to nearest boundary
698 const int iBinLo = xbins->binNumber(xlo + tolerance);
699 const int iBinHi = xbins->binNumber(xhi - tolerance);
700 const int nBinsAdj = iBinHi - iBinLo + 1;
701 const double xloAdj = xbins->binLow(iBinLo);
702 const double xhiAdj = xbins->binHigh(iBinHi);
703
704 if (isUniform) {
705 xbins = std::make_unique<RooUniformBinning>(xloAdj, xhiAdj, nBinsAdj);
706 theirVar.setRange(xloAdj, xhiAdj);
707 } else {
708 xbins->setRange(xloAdj, xhiAdj);
709 theirVar.setBinning(*xbins);
710 }
711
712 if (std::abs(xloAdj - xlo) > tolerance || std::abs(xhiAdj - xhi) > tolerance) {
713 coutI(DataHandling) << "RooDataHist::adjustBinning(" << ownName << "): fit range of variable " << ourVarName
714 << " expanded to nearest bin boundaries: [" << xlo << "," << xhi << "] --> [" << xloAdj << ","
715 << xhiAdj << "]"
716 << "\n";
717 }
718
719 ourVar->setBinning(*xbins);
720
721 // The offset is the bin number of the adjusted lower limit of the RooFit
722 // variable in the original TH1 histogram, starting from zero.
723 if (offset) {
724 *offset = axis.FindFixBin(xloAdj + tolerance) - 1;
725 }
726}
727
728////////////////////////////////////////////////////////////////////////////////
729/// Adjust binning specification on first and optionally second and third
730/// observable to binning in given reference TH1. Used by constructors
731/// that import data from an external TH1.
732/// Both the variables in vars and in this RooDataHist are adjusted.
733/// @param vars List with variables that are supposed to have their binning adjusted.
734/// @param href Reference histogram that dictates the binning
735/// @param offset If not nullptr, a possible bin count offset for the axes x,y,z is saved here as Int_t[3]
736
737void RooDataHist::adjustBinning(const RooArgList& vars, const TH1& href, Int_t* offset)
738{
739 auto xvar = static_cast<RooRealVar*>(_vars.find(*vars.at(0)) );
740 _adjustBinning(*static_cast<RooRealVar*>(vars.at(0)), *href.GetXaxis(), xvar, offset ? &offset[0] : nullptr);
741
742 if (vars.at(1)) {
743 auto yvar = static_cast<RooRealVar*>(_vars.find(*vars.at(1)));
744 if (yvar)
745 _adjustBinning(*static_cast<RooRealVar*>(vars.at(1)), *href.GetYaxis(), yvar, offset ? &offset[1] : nullptr);
746 }
747
748 if (vars.at(2)) {
749 auto zvar = static_cast<RooRealVar*>(_vars.find(*vars.at(2)));
750 if (zvar)
751 _adjustBinning(*static_cast<RooRealVar*>(vars.at(2)), *href.GetZaxis(), zvar, offset ? &offset[2] : nullptr);
752 }
753
754}
755
756
757namespace {
758/// Clone external weight arrays, unless the external array is nullptr.
759void cloneArray(double*& ours, const double* theirs, std::size_t n) {
760 if (ours) delete[] ours;
761 ours = nullptr;
762 if (!theirs) return;
763 ours = new double[n];
764 std::copy(theirs, theirs+n, ours);
765}
766
767/// Allocate and initialise an array with desired size and values.
768void initArray(double*& arr, std::size_t n, double val) {
769 if (arr) delete[] arr;
770 arr = nullptr;
771 if (n == 0) return;
772 arr = new double[n];
773 std::fill(arr, arr+n, val);
774}
775}
776
777
778////////////////////////////////////////////////////////////////////////////////
779/// Initialization procedure: allocate weights array, calculate
780/// multipliers needed for N-space to 1-dim array jump table,
781/// and fill the internal tree with all bin center coordinates
782
783void RooDataHist::initialize(const char* binningName, bool fillTree)
784{
785 _lvvars.clear();
786 _lvbins.clear();
787
788 // Fill array of LValue pointers to variables
789 for (unsigned int i = 0; i < _vars.size(); ++i) {
790 if (binningName) {
791 RooRealVar* rrv = dynamic_cast<RooRealVar*>(_vars[i]);
792 if (rrv) {
793 rrv->setBinning(rrv->getBinning(binningName));
794 }
795 }
796
797 auto lvarg = dynamic_cast<RooAbsLValue*>(_vars[i]);
798 assert(lvarg);
799 _lvvars.push_back(lvarg);
800
801 const RooAbsBinning* binning = lvarg->getBinningPtr(nullptr);
802 _lvbins.emplace_back(binning ? binning->clone() : nullptr);
803 }
804
805
806 // Allocate coefficients array
807 _idxMult.resize(_vars.size()) ;
808
809 _arrSize = 1 ;
810 unsigned int n = 0u;
811 for (const auto var : _vars) {
812 auto arg = dynamic_cast<const RooAbsLValue*>(var);
813 assert(arg);
814
815 // Calculate sub-index multipliers for master index
816 for (unsigned int i = 0u; i<n; i++) {
817 _idxMult[i] *= arg->numBins() ;
818 }
819 _idxMult[n++] = 1 ;
820
821 // Calculate dimension of weight array
822 _arrSize *= arg->numBins() ;
823 }
824
825 // Allocate and initialize weight array if necessary
826 if (!_wgt) {
827 initArray(_wgt, _arrSize, 0.);
828 delete[] _errLo; _errLo = nullptr;
829 delete[] _errHi; _errHi = nullptr;
830 delete[] _sumw2; _sumw2 = nullptr;
831 initArray(_binv, _arrSize, 0.);
832
833 // Refill array pointers in data store when reading
834 // from Streamer
835 if (!fillTree) {
837 }
838 }
839
840 if (!fillTree) return ;
841
842 // Fill TTree with bin center coordinates
843 // Calculate plot bins of components from master index
844
845 for (Int_t ibin=0 ; ibin < _arrSize ; ibin++) {
846 Int_t j(0);
847 Int_t idx(0);
848 Int_t tmp(ibin);
849 double theBinVolume(1) ;
850 for (auto arg2 : _lvvars) {
851 idx = tmp / _idxMult[j] ;
852 tmp -= idx*_idxMult[j++] ;
853 arg2->setBin(idx) ;
854 theBinVolume *= arg2->getBinWidth(idx) ;
855 }
856 _binv[ibin] = theBinVolume ;
857
858 fill() ;
859 }
860
861
862}
863
864
865////////////////////////////////////////////////////////////////////////////////
866
868{
869 if (!_binbounds.empty()) return;
870 for (auto& it : _lvbins) {
871 _binbounds.push_back(std::vector<double>());
872 if (it) {
873 std::vector<double>& bounds = _binbounds.back();
874 bounds.reserve(2 * it->numBins());
875 for (Int_t i = 0; i < it->numBins(); ++i) {
876 bounds.push_back(it->binLow(i));
877 bounds.push_back(it->binHigh(i));
878 }
879 }
880 }
881}
882
883
884////////////////////////////////////////////////////////////////////////////////
885/// Copy constructor
886
887RooDataHist::RooDataHist(const RooDataHist& other, const char* newname) :
888 RooAbsData(other,newname), RooDirItem(), _arrSize(other._arrSize), _idxMult(other._idxMult), _pbinvCache(other._pbinvCache)
889{
890 // Allocate and initialize weight array
891 assert(_arrSize == other._arrSize);
892 cloneArray(_wgt, other._wgt, other._arrSize);
893 cloneArray(_errLo, other._errLo, other._arrSize);
894 cloneArray(_errHi, other._errHi, other._arrSize);
895 cloneArray(_binv, other._binv, other._arrSize);
896 cloneArray(_sumw2, other._sumw2, other._arrSize);
897
898 // Fill array of LValue pointers to variables
899 for (const auto rvarg : _vars) {
900 auto lvarg = dynamic_cast<RooAbsLValue*>(rvarg);
901 assert(lvarg);
902 _lvvars.push_back(lvarg);
903 const RooAbsBinning* binning = lvarg->getBinningPtr(nullptr);
904 _lvbins.emplace_back(binning ? binning->clone() : nullptr) ;
905 }
906
908
909 appendToDir(this,true) ;
910}
911
912
913////////////////////////////////////////////////////////////////////////////////
914/// Implementation of RooAbsData virtual method that drives the RooAbsData::reduce() methods
915
916std::unique_ptr<RooAbsData> RooDataHist::reduceEng(const RooArgSet& varSubset, const RooFormulaVar* cutVar, const char* cutRange,
917 std::size_t nStart, std::size_t nStop) const
918{
919 checkInit() ;
920 RooArgSet myVarSubset;
921 _vars.selectCommon(varSubset, myVarSubset);
922 auto rdh = std::make_unique<RooDataHist>(GetName(), GetTitle(), myVarSubset);
923
924 RooFormulaVar* cloneVar = nullptr;
925 std::unique_ptr<RooArgSet> tmp;
926 if (cutVar) {
927 tmp = std::make_unique<RooArgSet>();
928 // Deep clone cutVar and attach clone to this dataset
929 if (RooArgSet(*cutVar).snapshot(*tmp)) {
930 coutE(DataHandling) << "RooDataHist::reduceEng(" << GetName() << ") Couldn't deep-clone cut variable, abort," << endl ;
931 return nullptr;
932 }
933 cloneVar = static_cast<RooFormulaVar*>(tmp->find(*cutVar));
934 cloneVar->attachDataSet(*this) ;
935 }
936
937 double lo;
938 double hi;
939 const std::size_t nevt = nStop < static_cast<std::size_t>(numEntries()) ? nStop : static_cast<std::size_t>(numEntries());
940 for (auto i=nStart; i<nevt ; i++) {
941 const RooArgSet* row = get(i) ;
942
943 bool doSelect(true) ;
944 if (cutRange) {
945 for (const auto arg : *row) {
946 if (!arg->inRange(cutRange)) {
947 doSelect = false ;
948 break ;
949 }
950 }
951 }
952 if (!doSelect) continue ;
953
954 if (!cloneVar || cloneVar->getVal()) {
955 weightError(lo,hi,SumW2) ;
956 rdh->add(*row,weight(),lo*lo) ;
957 }
958 }
959
960 return rdh ;
961}
962
963
964
965////////////////////////////////////////////////////////////////////////////////
966/// Destructor
967
969{
970 delete[] _wgt;
971 delete[] _errLo;
972 delete[] _errHi;
973 delete[] _sumw2;
974 delete[] _binv;
975
976 removeFromDir(this) ;
978}
979
980
981
982
983////////////////////////////////////////////////////////////////////////////////
984/// Calculate bin number of the given coordinates. If only a subset of the internal
985/// coordinates are passed, the missing coordinates are taken at their current value.
986/// \param[in] coord Variables that are representing the coordinates.
987/// \param[in] fast If the variables in `coord` and the ones of the data hist have the
988/// same size and layout, `fast` can be set to skip checking that all variables are
989/// present in `coord`.
990Int_t RooDataHist::getIndex(const RooAbsCollection& coord, bool fast) const {
991 checkInit() ;
992 return calcTreeIndex(coord, fast);
993}
994
996 bool correctForBinSize) const
997{
998 std::vector<double> vals(_arrSize);
999 if (correctForBinSize) {
1000 for (std::size_t i = 0; i < vals.size(); ++i) {
1001 vals[i] = _wgt[i] / _binv[i];
1002 }
1003 } else {
1004 for (std::size_t i = 0; i < vals.size(); ++i) {
1005 vals[i] = _wgt[i];
1006 }
1007 }
1008 return ctx.buildArg(vals);
1009}
1010
1012 const RooAbsCollection &coords, bool reverse) const
1013{
1014 assert(coords.size() == _vars.size());
1015
1016 std::string code;
1017 int idxMult = 1;
1018
1019 for (std::size_t i = 0; i < _vars.size(); ++i) {
1020
1021 std::size_t iVar = reverse ? _vars.size() - 1 - i : i;
1022 const RooAbsArg *internalVar = _vars[iVar];
1023 const RooAbsArg *theVar = coords[iVar];
1024
1025 const RooAbsBinning *binning = _lvbins[iVar].get();
1026 if (!binning) {
1027 coutE(InputArguments) << "RooHistPdf::weight(" << GetName()
1028 << ") ERROR: Code Squashing currently does not support category values." << std::endl;
1029 return "";
1030 } else if (!dynamic_cast<RooUniformBinning const *>(binning)) {
1031 coutE(InputArguments) << "RooHistPdf::weight(" << GetName()
1032 << ") ERROR: Code Squashing currently only supports uniformly binned cases."
1033 << std::endl;
1034 return "";
1035 }
1036
1037 std::string const &bin = ctx.buildCall("RooFit::Detail::MathFuncs::getUniformBinning", binning->lowBound(),
1038 binning->highBound(), *theVar, binning->numBins());
1039 code += " + " + std::to_string(idxMult) + " * " + bin;
1040
1041 // Use RooAbsLValue here because it also generalized to categories, which
1042 // is useful in the future. dynamic_cast because it's a cross-cast.
1043 idxMult *= dynamic_cast<RooAbsLValue const *>(internalVar)->numBins();
1044 }
1045
1046 return "(" + code + ")";
1047}
1048
1049////////////////////////////////////////////////////////////////////////////////
1050/// Calculate the bin index corresponding to the coordinates passed as argument.
1051/// \param[in] coords Coordinates. If `fast == false`, these can be partial.
1052/// \param[in] fast Promise that the coordinates in `coords` have the same order
1053/// as the internal coordinates. In this case, values are looked up only by index.
1054std::size_t RooDataHist::calcTreeIndex(const RooAbsCollection& coords, bool fast) const
1055{
1056 // With fast, caller promises that layout of `coords` is identical to our internal `vars`.
1057 // Previously, this was verified with an assert in debug mode like this:
1058 //
1059 // assert(!fast || coords.hasSameLayout(_vars));
1060 //
1061 // However, there are usecases where the externally provided `coords` have
1062 // different names than the internal variables, even though they correspond
1063 // to each other. For example, if the observables in the computation graph
1064 // are renamed with `redirectServers`. Hence, we can't do a meaningful assert
1065 // here.
1066
1067 if (&_vars == &coords)
1068 fast = true;
1069
1070 std::size_t masterIdx = 0;
1071
1072 for (unsigned int i=0; i < _vars.size(); ++i) {
1073 const RooAbsArg* internalVar = _vars[i];
1074 const RooAbsBinning* binning = _lvbins[i].get();
1075
1076 // Find the variable that we need values from.
1077 // That's either the variable directly from the external coordinates
1078 // or we find the external one that has the same name as "internalVar".
1079 const RooAbsArg* theVar = fast ? coords[i] : coords.find(*internalVar);
1080 if (!theVar) {
1081 // Variable is not in external coordinates. Use current internal value.
1082 theVar = internalVar;
1083 }
1084 // If fast is on, users promise that the sets have the same layout:
1085 //
1086 // assert(!fast || strcmp(internalVar->GetName(), theVar->GetName()) == 0);
1087 //
1088 // This assert is commented out for the same reasons that applied to the
1089 // other assert explained above.
1090
1091 if (binning) {
1092 assert(dynamic_cast<const RooAbsReal*>(theVar));
1093 const double val = static_cast<const RooAbsReal*>(theVar)->getVal();
1094 masterIdx += _idxMult[i] * binning->binNumber(val);
1095 } else {
1096 // We are a category. No binning.
1097 assert(dynamic_cast<const RooAbsCategoryLValue*>(theVar));
1098 auto cat = static_cast<const RooAbsCategoryLValue*>(theVar);
1099 masterIdx += _idxMult[i] * cat->getBin(static_cast<const char*>(nullptr));
1100 }
1101 }
1102
1103 return masterIdx ;
1104}
1105
1106
1107////////////////////////////////////////////////////////////////////////////////
1108/// Back end function to plotting functionality. Plot RooDataHist on given
1109/// frame in mode specified by plot options 'o'. The main purpose of
1110/// this function is to match the specified binning on 'o' to the
1111/// internal binning of the plot observable in this RooDataHist.
1112/// \see RooAbsData::plotOn() for plotting options.
1114{
1115 checkInit() ;
1116 if (o.bins) return RooAbsData::plotOn(frame,o) ;
1117
1118 if(!frame) {
1119 coutE(InputArguments) << ClassName() << "::" << GetName() << ":plotOn: frame is null" << endl;
1120 return nullptr;
1121 }
1122 auto var= static_cast<RooAbsRealLValue*>(frame->getPlotVar());
1123 if(!var) {
1124 coutE(InputArguments) << ClassName() << "::" << GetName()
1125 << ":plotOn: frame does not specify a plot variable" << endl;
1126 return nullptr;
1127 }
1128
1129 auto dataVar = static_cast<RooRealVar*>(_vars.find(*var));
1130 if (!dataVar) {
1131 coutE(InputArguments) << ClassName() << "::" << GetName()
1132 << ":plotOn: dataset doesn't contain plot frame variable" << endl;
1133 return nullptr;
1134 }
1135
1136 o.bins = &dataVar->getBinning() ;
1137 return RooAbsData::plotOn(frame,o) ;
1138}
1139
1140
1141////////////////////////////////////////////////////////////////////////////////
1142/// A vectorized version of interpolateDim for boundary safe quadratic
1143/// interpolation of one dimensional histograms.
1144///
1145/// \param[out] output An array of interpolated weights corresponding to the
1146/// values in xVals.
1147/// \param[in] xVals An array of event coordinates for which the weights should be
1148/// calculated.
1149/// \param[in] correctForBinSize Enable the inverse bin volume correction factor.
1150/// \param[in] cdfBoundaries Enable the special boundary condition for a cdf:
1151/// Underflow bins are assumed to have weight zero and
1152/// overflow bins have weight one. Otherwise, the
1153/// histogram is mirrored at the boundaries for the
1154/// interpolation.
1155
1156void RooDataHist::interpolateQuadratic(double* output, std::span<const double> xVals,
1157 bool correctForBinSize, bool cdfBoundaries)
1158{
1159 const std::size_t nBins = numEntries();
1160 const std::size_t nEvents = xVals.size();
1161
1162 RooAbsBinning const& binning = *_lvbins[0];
1163 // Reuse the output buffer for bin indices and zero-initialize it
1164 auto binIndices = reinterpret_cast<int*>(output + nEvents) - nEvents;
1165 std::fill(binIndices, binIndices + nEvents, 0);
1166 binning.binNumbers(xVals.data(), binIndices, nEvents);
1167
1168 // Extend coordinates and weights with one extra point before the first bin
1169 // and one extra point after the last bin. This means the original histogram
1170 // bins span elements 1 to nBins in coordsExt and weightsExt
1171 std::vector<double> coordsExt(nBins+3);
1172 double* binCoords = coordsExt.data() + 2;
1173 binCoords[0] = binning.lowBound() + 0.5*_binv[0];
1174 for (std::size_t binIdx = 1; binIdx < nBins ; ++binIdx) {
1175 if (binning.isUniform()) {
1176 double binWidth = _binv[0];
1177 binCoords[binIdx] = binIdx*binWidth + binCoords[0];
1178 }
1179 else {
1180 double binCentDiff = 0.5*_binv[binIdx-1] + 0.5*_binv[binIdx];
1181 binCoords[binIdx] = binCoords[binIdx-1] + binCentDiff;
1182 }
1183 }
1184
1185 std::vector<double> weightsExt(nBins+3);
1186 // Fill weights for bins that are inside histogram boundaries
1187 for (std::size_t binIdx = 0; binIdx < nBins; ++binIdx) {
1188 weightsExt[binIdx+2] = correctForBinSize ? _wgt[binIdx] / _binv[binIdx] : _wgt[binIdx];
1189 }
1190
1191 if (cdfBoundaries) {
1192 coordsExt[0] = - 1e-10 + binning.lowBound();
1193 weightsExt[0] = 0.;
1194
1195 coordsExt[1] = binning.lowBound();
1196 weightsExt[1] = 0.;
1197
1198 coordsExt[nBins+2] = binning.highBound();
1199 weightsExt[nBins+2] = 1.;
1200 }
1201 else {
1202 // Mirror first two bins and last bin
1203 coordsExt[0] = binCoords[1] - 2*_binv[0] - _binv[1];
1204 weightsExt[0] = weightsExt[3];
1205
1206 coordsExt[1] = binCoords[0] - _binv[0];
1207 weightsExt[1] = weightsExt[2];
1208
1209 coordsExt[nBins+2] = binCoords[nBins-1] + _binv[nBins-1];
1210 weightsExt[nBins+2] = weightsExt[nBins+1];
1211 }
1212
1213 // We use the current bin center and two bin centers on the left for
1214 // interpolation if xVal is to the left of the current bin center
1215 for (std::size_t i = 0; i < nEvents ; ++i) {
1216 double xVal = xVals[i];
1217 std::size_t binIdx = binIndices[i] + 2;
1218
1219 // If xVal is to the right of the current bin center, shift all bin
1220 // coordinates one step to the right and use that for the interpolation
1221 if (xVal > coordsExt[binIdx]) {
1222 binIdx += 1;
1223 }
1224
1225 double x1 = coordsExt[binIdx-2];
1226 double y1 = weightsExt[binIdx-2];
1227
1228 double x2 = coordsExt[binIdx-1];
1229 double y2 = weightsExt[binIdx-1];
1230
1231 double x3 = coordsExt[binIdx];
1232 double y3 = weightsExt[binIdx];
1233
1234 // Evaluate a few repeated factors
1235 double quotient = (x3-x1) / (x2-x1);
1236 double x1Sqrd = x1*x1;
1237 double x3Sqrd = x3*x3;
1238 // Solve coefficients in system of three quadratic equations!
1239 double secondCoeff = (y3 - y1 - (y2-y1) * quotient) / (x3Sqrd - x1Sqrd - (x2*x2 - x1Sqrd) * quotient);
1240 double firstCoeff = (y3 - y1 - secondCoeff*(x3Sqrd - x1Sqrd)) / (x3-x1);
1241 double zerothCoeff = y1 - secondCoeff*x1Sqrd - firstCoeff*x1;
1242 // Get the interpolated weight using the equation of a second degree polynomial
1243 output[i] = secondCoeff * xVal * xVal + firstCoeff * xVal + zerothCoeff;
1244 }
1245}
1246
1247
1248////////////////////////////////////////////////////////////////////////////////
1249/// A vectorized version of interpolateDim for boundary safe linear
1250/// interpolation of one dimensional histograms.
1251///
1252/// \param[out] output An array of interpolated weights corresponding to the
1253/// values in xVals.
1254/// \param[in] xVals An array of event coordinates for which the weights should be
1255/// calculated.
1256/// \param[in] correctForBinSize Enable the inverse bin volume correction factor.
1257/// \param[in] cdfBoundaries Enable the special boundary condition for a cdf:
1258/// Underflow bins are assumed to have weight zero and
1259/// overflow bins have weight one. Otherwise, the
1260/// histogram is mirrored at the boundaries for the
1261/// interpolation.
1262
1263void RooDataHist::interpolateLinear(double* output, std::span<const double> xVals,
1264 bool correctForBinSize, bool cdfBoundaries)
1265{
1266 const std::size_t nBins = numEntries();
1267 const std::size_t nEvents = xVals.size();
1268
1269 RooAbsBinning const& binning = *_lvbins[0];
1270 // Reuse the output buffer for bin indices and zero-initialize it
1271 auto binIndices = reinterpret_cast<int*>(output + nEvents) - nEvents;
1272 std::fill(binIndices, binIndices + nEvents, 0);
1273 binning.binNumbers(xVals.data(), binIndices, nEvents);
1274
1275 // Extend coordinates and weights with one extra point before the first bin
1276 // and one extra point after the last bin. This means the original histogram
1277 // bins span elements 1 to nBins in coordsExt and weightsExt
1278 std::vector<double> coordsExt(nBins+2);
1279 double* binCoords = coordsExt.data() + 1;
1280 binCoords[0] = binning.lowBound() + 0.5*_binv[0];
1281 for (std::size_t binIdx = 1; binIdx < nBins ; ++binIdx) {
1282 if (binning.isUniform()) {
1283 double binWidth = _binv[0];
1284 binCoords[binIdx] = binIdx*binWidth + binCoords[0];
1285 }
1286 else {
1287 double binCentDiff = 0.5*_binv[binIdx-1] + 0.5*_binv[binIdx];
1288 binCoords[binIdx] = binCoords[binIdx-1] + binCentDiff;
1289 }
1290 }
1291
1292 std::vector<double> weightsExt(nBins+2);
1293 // Fill weights for bins that are inside histogram boundaries
1294 for (std::size_t binIdx = 0; binIdx < nBins; ++binIdx) {
1295 weightsExt[binIdx+1] = correctForBinSize ? _wgt[binIdx] / _binv[binIdx] : _wgt[binIdx];
1296 }
1297
1298 // Fill weights for bins that are outside histogram boundaries
1299 if (cdfBoundaries) {
1300 coordsExt[0] = binning.lowBound();
1301 weightsExt[0] = 0.;
1302 coordsExt[nBins+1] = binning.highBound();
1303 weightsExt[nBins+1] = 1.;
1304 }
1305 else {
1306 // Mirror first and last bins
1307 coordsExt[0] = binCoords[0] - _binv[0];
1308 weightsExt[0] = weightsExt[1];
1309 coordsExt[nBins+1] = binCoords[nBins-1] + _binv[nBins-1];
1310 weightsExt[nBins+1] = weightsExt[nBins];
1311 }
1312
1313 // Interpolate between current bin center and one bin center to the left
1314 // if xVal is to the left of the current bin center
1315 for (std::size_t i = 0; i < nEvents ; ++i) {
1316 double xVal = xVals[i];
1317 std::size_t binIdx = binIndices[i] + 1;
1318
1319 // If xVal is to the right of the current bin center, interpolate between
1320 // current bin center and one bin center to the right instead
1321 if (xVal > coordsExt[binIdx]) { binIdx += 1; }
1322
1323 double x1 = coordsExt[binIdx-1];
1324 double y1 = weightsExt[binIdx-1];
1325 double x2 = coordsExt[binIdx];
1326 double y2 = weightsExt[binIdx];
1327
1328 // Find coefficients by solving a system of two linear equations
1329 double firstCoeff = (y2-y1) / (x2-x1);
1330 double zerothCoeff = y1 - firstCoeff * x1;
1331 // Get the interpolated weight using the equation of a straight line
1332 output[i] = firstCoeff * xVal + zerothCoeff;
1333 }
1334}
1335
1336
1337////////////////////////////////////////////////////////////////////////////////
1338/// A vectorized version of RooDataHist::weight() for one dimensional histograms
1339/// with up to one dimensional interpolation.
1340/// \param[out] output An array of weights corresponding the values in xVals.
1341/// \param[in] xVals An array of coordinates for which the weights should be
1342/// calculated.
1343/// \param[in] intOrder Interpolation order; 0th and 1st order are supported.
1344/// \param[in] correctForBinSize Enable the inverse bin volume correction factor.
1345/// \param[in] cdfBoundaries Enable the special boundary condition for a cdf:
1346/// Underflow bins are assumed to have weight zero and
1347/// overflow bins have weight one. Otherwise, the
1348/// histogram is mirrored at the boundaries for the
1349/// interpolation.
1350
1351void RooDataHist::weights(double* output, std::span<double const> xVals, int intOrder, bool correctForBinSize, bool cdfBoundaries)
1352{
1353 auto const nEvents = xVals.size();
1354
1355 if (intOrder == 0) {
1356 RooAbsBinning const& binning = *_lvbins[0];
1357
1358 // Reuse the output buffer for bin indices and zero-initialize it
1359 auto binIndices = reinterpret_cast<int*>(output + nEvents) - nEvents;
1360 std::fill(binIndices, binIndices + nEvents, 0);
1361 binning.binNumbers(xVals.data(), binIndices, nEvents);
1362
1363 for (std::size_t i=0; i < nEvents; ++i) {
1364 auto binIdx = binIndices[i];
1365 output[i] = correctForBinSize ? _wgt[binIdx] / _binv[binIdx] : _wgt[binIdx];
1366 }
1367 }
1368 else if (intOrder == 1) {
1369 interpolateLinear(output, xVals, correctForBinSize, cdfBoundaries);
1370 }
1371 else if (intOrder == 2) {
1372 interpolateQuadratic(output, xVals, correctForBinSize, cdfBoundaries);
1373 }
1374 else {
1375 // Higher dimensional scenarios not yet implemented
1376 coutE(InputArguments) << "RooDataHist::weights(" << GetName() << ") interpolation in "
1377 << intOrder << " dimensions not yet implemented" << std::endl ;
1378 // Fall back to 1st order interpolation
1379 weights(output, xVals, 1, correctForBinSize, cdfBoundaries);
1380 }
1381}
1382
1383
1384////////////////////////////////////////////////////////////////////////////////
1385/// A faster version of RooDataHist::weight that assumes the passed arguments
1386/// are aligned with the histogram variables.
1387/// \param[in] bin Coordinates for which the weight should be calculated.
1388/// Has to be aligned with the internal histogram variables.
1389/// \param[in] intOrder Interpolation order, i.e. how many neighbouring bins are
1390/// used for the interpolation. If zero, the bare weight for
1391/// the bin enclosing the coordinatesis returned.
1392/// \param[in] correctForBinSize Enable the inverse bin volume correction factor.
1393/// \param[in] cdfBoundaries Enable the special boundary condition for a cdf:
1394/// underflow bins are assumed to have weight zero and
1395/// overflow bins have weight one. Otherwise, the
1396/// histogram is mirrored at the boundaries for the
1397/// interpolation.
1398
1399double RooDataHist::weightFast(const RooArgSet& bin, Int_t intOrder, bool correctForBinSize, bool cdfBoundaries)
1400{
1401 checkInit() ;
1402
1403 // Handle illegal intOrder values
1404 if (intOrder<0) {
1405 coutE(InputArguments) << "RooDataHist::weight(" << GetName() << ") ERROR: interpolation order must be positive" << endl ;
1406 return 0 ;
1407 }
1408
1409 // Handle no-interpolation case
1410 if (intOrder==0) {
1411 const auto idx = calcTreeIndex(bin, true);
1412 return correctForBinSize ? _wgt[idx] / _binv[idx] : _wgt[idx];
1413 }
1414
1415 // Handle all interpolation cases
1416 return weightInterpolated(bin, intOrder, correctForBinSize, cdfBoundaries);
1417}
1418
1419
1420////////////////////////////////////////////////////////////////////////////////
1421/// Return the weight at given coordinates with optional interpolation.
1422/// \param[in] bin Coordinates for which the weight should be calculated.
1423/// \param[in] intOrder Interpolation order, i.e. how many neighbouring bins are
1424/// used for the interpolation. If zero, the bare weight for
1425/// the bin enclosing the coordinatesis returned.
1426/// \param[in] correctForBinSize Enable the inverse bin volume correction factor.
1427/// \param[in] cdfBoundaries Enable the special boundary condition for a cdf:
1428/// underflow bins are assumed to have weight zero and
1429/// overflow bins have weight one. Otherwise, the
1430/// histogram is mirrored at the boundaries for the
1431/// interpolation.
1432/// \param[in] oneSafe Ignored.
1433
1434double RooDataHist::weight(const RooArgSet& bin, Int_t intOrder, bool correctForBinSize, bool cdfBoundaries, bool /*oneSafe*/)
1435{
1436 checkInit() ;
1437
1438 // Handle illegal intOrder values
1439 if (intOrder<0) {
1440 coutE(InputArguments) << "RooDataHist::weight(" << GetName() << ") ERROR: interpolation order must be positive" << endl ;
1441 return 0 ;
1442 }
1443
1444 // Handle no-interpolation case
1445 if (intOrder==0) {
1446 const auto idx = calcTreeIndex(bin, false);
1447 return correctForBinSize ? _wgt[idx] / _binv[idx] : _wgt[idx];
1448 }
1449
1450 // Handle all interpolation cases
1451 _vars.assignValueOnly(bin) ;
1452
1453 return weightInterpolated(_vars, intOrder, correctForBinSize, cdfBoundaries);
1454}
1455
1456
1457////////////////////////////////////////////////////////////////////////////////
1458/// Return the weight at given coordinates with interpolation.
1459/// \param[in] bin Coordinates for which the weight should be calculated.
1460/// Has to be aligned with the internal histogram variables.
1461/// \param[in] intOrder Interpolation order, i.e. how many neighbouring bins are
1462/// used for the interpolation.
1463/// \param[in] correctForBinSize Enable the inverse bin volume correction factor.
1464/// \param[in] cdfBoundaries Enable the special boundary condition for a cdf:
1465/// underflow bins are assumed to have weight zero and
1466/// overflow bins have weight one. Otherwise, the
1467/// histogram is mirrored at the boundaries for the
1468/// interpolation.
1469
1470double RooDataHist::weightInterpolated(const RooArgSet& bin, int intOrder, bool correctForBinSize, bool cdfBoundaries) {
1471 VarInfo const& varInfo = getVarInfo();
1472
1473 const auto centralIdx = calcTreeIndex(bin, true);
1474
1475 double wInt{0} ;
1476 if (varInfo.nRealVars == 1) {
1477
1478 // buffer needs to be 2 x (interpolation order + 1), with the factor 2 for x and y.
1479 _interpolationBuffer.resize(2 * intOrder + 2);
1480
1481 // 1-dimensional interpolation
1482 auto const& realX = static_cast<RooRealVar const&>(*bin[varInfo.realVarIdx1]);
1483 wInt = interpolateDim(varInfo.realVarIdx1, realX.getVal(), centralIdx, intOrder, correctForBinSize, cdfBoundaries) ;
1484
1485 } else if (varInfo.nRealVars == 2) {
1486
1487 // buffer needs to be 2 x 2 x (interpolation order + 1), with one factor 2
1488 // for x and y, and the other for the number of dimensions.
1489 _interpolationBuffer.resize(4 * intOrder + 4);
1490
1491 // 2-dimensional interpolation
1492 auto const& realX = static_cast<RooRealVar const&>(*bin[varInfo.realVarIdx1]);
1493 auto const& realY = static_cast<RooRealVar const&>(*bin[varInfo.realVarIdx2]);
1494 double xval = realX.getVal() ;
1495 double yval = realY.getVal() ;
1496
1497 RooAbsBinning const& binningY = realY.getBinning();
1498
1499 int ybinC = binningY.binNumber(yval) ;
1500 int ybinLo = ybinC-intOrder/2 - ((yval<binningY.binCenter(ybinC))?1:0) ;
1501 int ybinM = binningY.numBins() ;
1502
1503 auto idxMultY = _idxMult[varInfo.realVarIdx2];
1504 auto offsetIdx = centralIdx - idxMultY * ybinC;
1505
1506 // Use a class-member buffer to avoid repeated heap allocations.
1507 double * yarr = _interpolationBuffer.data() + 2 * intOrder + 2; // add offset to skip part reserved for other dim
1508 double * xarr = yarr + intOrder + 1;
1509 for (int i=ybinLo ; i<=intOrder+ybinLo ; i++) {
1510 int ibin ;
1511 if (i>=0 && i<ybinM) {
1512 // In range
1513 ibin = i ;
1514 xarr[i-ybinLo] = binningY.binCenter(ibin) ;
1515 } else if (i>=ybinM) {
1516 // Overflow: mirror
1517 ibin = 2*ybinM-i-1 ;
1518 xarr[i-ybinLo] = 2*binningY.highBound()-binningY.binCenter(ibin) ;
1519 } else {
1520 // Underflow: mirror
1521 ibin = -i -1;
1522 xarr[i-ybinLo] = 2*binningY.lowBound()-binningY.binCenter(ibin) ;
1523 }
1524 auto centralIdxX = offsetIdx + idxMultY * ibin;
1525 yarr[i-ybinLo] = interpolateDim(varInfo.realVarIdx1,xval,centralIdxX,intOrder,correctForBinSize,false) ;
1526 }
1527
1528 if (gDebug>7) {
1529 cout << "RooDataHist interpolating data is" << endl ;
1530 cout << "xarr = " ;
1531 for (int q=0; q<=intOrder ; q++) cout << xarr[q] << " " ;
1532 cout << " yarr = " ;
1533 for (int q=0; q<=intOrder ; q++) cout << yarr[q] << " " ;
1534 cout << endl ;
1535 }
1536 wInt = RooMath::interpolate(xarr,yarr,intOrder+1,yval) ;
1537
1538 } else {
1539
1540 // Higher dimensional scenarios not yet implemented
1541 coutE(InputArguments) << "RooDataHist::weight(" << GetName() << ") interpolation in "
1542 << varInfo.nRealVars << " dimensions not yet implemented" << endl ;
1543 return weightFast(bin,0,correctForBinSize,cdfBoundaries) ;
1544
1545 }
1546
1547 return wInt ;
1548}
1549
1550
1552 if (!_errLo || !_errHi) {
1553 initArray(_errLo, _arrSize, -1.);
1554 initArray(_errHi, _arrSize, -1.);
1556 }
1557}
1558
1559
1560////////////////////////////////////////////////////////////////////////////////
1561/// Return the asymmetric errors on the current weight.
1562/// \see weightError(ErrorType) const for symmetric error.
1563/// \param[out] lo Low error.
1564/// \param[out] hi High error.
1565/// \param[in] etype Type of error to compute. May throw if not supported.
1566/// Supported errors are
1567/// - `Poisson` Default. Asymmetric Poisson errors (68% CL).
1568/// - `SumW2` The square root of the sum of weights. (Symmetric).
1569/// - `None` Return zero.
1570void RooDataHist::weightError(double& lo, double& hi, ErrorType etype) const
1571{
1572 checkInit() ;
1573
1574 switch (etype) {
1575
1576 case Auto:
1577 throw std::invalid_argument("RooDataHist::weightError(" + std::string(GetName()) + ") error type Auto not allowed here");
1578 break ;
1579
1580 case Expected:
1581 throw std::invalid_argument("RooDataHist::weightError(" + std::string(GetName()) + ") error type Expected not allowed here");
1582 break ;
1583
1584 case Poisson:
1585 if (_errLo && _errLo[_curIndex] >= 0.0) {
1586 // Weight is preset or precalculated
1587 lo = _errLo[_curIndex];
1588 hi = _errHi[_curIndex];
1589 return ;
1590 }
1591
1592 // We didn't track asymmetric errors so far, so now we need to allocate
1594
1595 // Calculate poisson errors
1596 double ym;
1597 double yp;
1599 _errLo[_curIndex] = weight()-ym;
1600 _errHi[_curIndex] = yp-weight();
1601 lo = _errLo[_curIndex];
1602 hi = _errHi[_curIndex];
1603 return ;
1604
1605 case SumW2:
1606 lo = std::sqrt(weightSquared(_curIndex));
1607 hi = lo;
1608 return ;
1609
1610 case None:
1611 lo = 0 ;
1612 hi = 0 ;
1613 return ;
1614 }
1615}
1616
1617
1618// wve adjust for variable bin sizes
1619
1620////////////////////////////////////////////////////////////////////////////////
1621/// Perform boundary safe 'intOrder'-th interpolation of weights in dimension 'dim'
1622/// at current value 'xval'
1623
1624/// \param[in] iDim Index of the histogram dimension along which to interpolate.
1625/// \param[in] xval Value of histogram variable at dimension `iDim` for which
1626/// we want to interpolate the histogram weight.
1627/// \param[in] centralIdx Index of the bin that the point at which we
1628/// interpolate the histogram weight falls into
1629/// (can be obtained with `RooDataHist::calcTreeIndex`).
1630/// \param[in] intOrder Interpolation order, i.e. how many neighbouring bins are
1631/// used for the interpolation.
1632/// \param[in] correctForBinSize Enable the inverse bin volume correction factor.
1633/// \param[in] cdfBoundaries Enable the special boundary condition for a cdf:
1634/// underflow bins are assumed to have weight zero and
1635/// overflow bins have weight one. Otherwise, the
1636/// histogram is mirrored at the boundaries for the
1637/// interpolation.
1638double RooDataHist::interpolateDim(int iDim, double xval, size_t centralIdx, int intOrder, bool correctForBinSize, bool cdfBoundaries)
1639{
1640 auto const& binning = static_cast<RooRealVar&>(*_vars[iDim]).getBinning();
1641
1642 // Fill workspace arrays spanning interpolation area
1643 int fbinC = binning.binNumber(xval) ;
1644 int fbinLo = fbinC-intOrder/2 - ((xval<binning.binCenter(fbinC))?1:0) ;
1645 int fbinM = binning.numBins() ;
1646
1647 auto idxMult = _idxMult[iDim];
1648 auto offsetIdx = centralIdx - idxMult * fbinC;
1649
1650 // Use a class-member buffer to avoid repeated heap allocations.
1651 double * yarr = _interpolationBuffer.data();
1652 double * xarr = yarr + intOrder + 1;
1653
1654 for (int i=fbinLo ; i<=intOrder+fbinLo ; i++) {
1655 int ibin ;
1656 if (i>=0 && i<fbinM) {
1657 // In range
1658 ibin = i ;
1659 xarr[i-fbinLo] = binning.binCenter(ibin) ;
1660 auto idx = offsetIdx + idxMult * ibin;
1661 yarr[i - fbinLo] = _wgt[idx];
1662 if (correctForBinSize) yarr[i-fbinLo] /= _binv[idx] ;
1663 } else if (i>=fbinM) {
1664 // Overflow: mirror
1665 ibin = 2*fbinM-i-1 ;
1666 if (cdfBoundaries) {
1667 xarr[i-fbinLo] = binning.highBound()+1e-10*(i-fbinM+1) ;
1668 yarr[i-fbinLo] = 1.0 ;
1669 } else {
1670 auto idx = offsetIdx + idxMult * ibin;
1671 xarr[i-fbinLo] = 2*binning.highBound()-binning.binCenter(ibin) ;
1672 yarr[i - fbinLo] = _wgt[idx];
1673 if (correctForBinSize)
1674 yarr[i - fbinLo] /= _binv[idx];
1675 }
1676 } else {
1677 // Underflow: mirror
1678 ibin = -i - 1 ;
1679 if (cdfBoundaries) {
1680 xarr[i-fbinLo] = binning.lowBound()-ibin*(1e-10) ;
1681 yarr[i-fbinLo] = 0.0 ;
1682 } else {
1683 auto idx = offsetIdx + idxMult * ibin;
1684 xarr[i-fbinLo] = 2*binning.lowBound()-binning.binCenter(ibin) ;
1685 yarr[i - fbinLo] = _wgt[idx];
1686 if (correctForBinSize)
1687 yarr[i - fbinLo] /= _binv[idx];
1688 }
1689 }
1690 }
1691 return RooMath::interpolate(xarr,yarr,intOrder+1,xval) ;
1692}
1693
1694
1695
1696
1697////////////////////////////////////////////////////////////////////////////////
1698/// Increment the bin content of the bin enclosing the given coordinates.
1699///
1700/// \param[in] row Coordinates of the bin.
1701/// \param[in] wgt Increment by this weight.
1702/// \param[in] sumw2 Optionally, track the sum of squared weights. If a value > 0 or
1703/// a weight != 1. is passed for the first time, a vector for the squared weights will be allocated.
1704void RooDataHist::add(const RooArgSet& row, double wgt, double sumw2)
1705{
1706 checkInit() ;
1707
1708 if ((sumw2 > 0. || wgt != 1.) && !_sumw2) {
1709 // Receiving a weighted entry. SumW2 != sumw from now on.
1710 _sumw2 = new double[_arrSize];
1711 std::copy(_wgt, _wgt+_arrSize, _sumw2);
1712
1714 }
1715
1716 const auto idx = calcTreeIndex(row, false);
1717
1718 _wgt[idx] += wgt ;
1719 if (_sumw2) _sumw2[idx] += (sumw2 > 0 ? sumw2 : wgt*wgt);
1720
1721 _cache_sum_valid = false;
1722}
1723
1724
1725
1726////////////////////////////////////////////////////////////////////////////////
1727/// Set a bin content.
1728/// \param[in] row Coordinates of the bin to be set.
1729/// \param[in] wgt New bin content.
1730/// \param[in] wgtErrLo Low error of the bin content.
1731/// \param[in] wgtErrHi High error of the bin content.
1732void RooDataHist::set(const RooArgSet& row, double wgt, double wgtErrLo, double wgtErrHi)
1733{
1734 checkInit() ;
1735
1737
1738 const auto idx = calcTreeIndex(row, false);
1739
1740 _wgt[idx] = wgt ;
1741 _errLo[idx] = wgtErrLo ;
1742 _errHi[idx] = wgtErrHi ;
1743
1744 _cache_sum_valid = false;
1745}
1746
1747
1748
1749////////////////////////////////////////////////////////////////////////////////
1750/// Set bin content of bin that was last loaded with get(std::size_t).
1751/// \param[in] binNumber Optional bin number to set. If empty, currently active bin is set.
1752/// \param[in] wgt New bin content.
1753/// \param[in] wgtErr Error of the new bin content. If the weight need not have an error, use 0. or a negative number.
1754void RooDataHist::set(std::size_t binNumber, double wgt, double wgtErr) {
1755 checkInit() ;
1756
1757 if (wgtErr > 0. && !_sumw2) {
1758 // Receiving a weighted entry. Need to track sumw2 from now on:
1759 cloneArray(_sumw2, _wgt, _arrSize);
1760
1762 }
1763
1764 _wgt[binNumber] = wgt ;
1765 if (_errLo) _errLo[binNumber] = wgtErr;
1766 if (_errHi) _errHi[binNumber] = wgtErr;
1767 if (_sumw2) _sumw2[binNumber] = wgtErr*wgtErr;
1768
1769 _cache_sum_valid = false ;
1770}
1771
1772
1773////////////////////////////////////////////////////////////////////////////////
1774/// Set bin content of bin that was last loaded with get(std::size_t).
1775/// \param[in] wgt New bin content.
1776/// \param[in] wgtErr Optional error of the bin content.
1777void RooDataHist::set(double wgt, double wgtErr) {
1778 if (_curIndex == std::numeric_limits<std::size_t>::max()) {
1779 _curIndex = calcTreeIndex(_vars, true) ;
1780 }
1781
1782 set(_curIndex, wgt, wgtErr);
1783}
1784
1785
1786////////////////////////////////////////////////////////////////////////////////
1787/// Set a bin content.
1788/// \param[in] row Coordinates to compute the bin from.
1789/// \param[in] wgt New bin content.
1790/// \param[in] wgtErr Optional error of the bin content.
1791void RooDataHist::set(const RooArgSet& row, double wgt, double wgtErr) {
1792 set(calcTreeIndex(row, false), wgt, wgtErr);
1793}
1794
1795
1796
1797////////////////////////////////////////////////////////////////////////////////
1798/// Add all data points contained in 'dset' to this data set with given weight.
1799/// Optional cut string expression selects the data points to be added and can
1800/// reference any variable contained in this data set
1801
1802void RooDataHist::add(const RooAbsData& dset, const char* cut, double wgt)
1803{
1804 RooFormulaVar cutVar("select",cut,*dset.get()) ;
1805 add(dset,&cutVar,wgt) ;
1806}
1807
1808
1809
1810////////////////////////////////////////////////////////////////////////////////
1811/// Add all data points contained in 'dset' to this data set with given weight.
1812/// Optional RooFormulaVar pointer selects the data points to be added.
1813
1814void RooDataHist::add(const RooAbsData& dset, const RooFormulaVar* cutVar, double wgt)
1815{
1816 checkInit() ;
1817
1818 RooFormulaVar* cloneVar = nullptr;
1819 std::unique_ptr<RooArgSet> tmp;
1820 if (cutVar) {
1821 // Deep clone cutVar and attach clone to this dataset
1822 tmp = std::make_unique<RooArgSet>();
1823 if(RooArgSet(*cutVar).snapshot(*tmp)) {
1824 coutE(DataHandling) << "RooDataHist::add(" << GetName() << ") Couldn't deep-clone cut variable, abort," << endl ;
1825 return ;
1826 }
1827
1828 cloneVar = static_cast<RooFormulaVar*>(tmp->find(*cutVar)) ;
1829 cloneVar->attachDataSet(dset) ;
1830 }
1831
1832
1833 Int_t i ;
1834 for (i=0 ; i<dset.numEntries() ; i++) {
1835 const RooArgSet* row = dset.get(i) ;
1836 if (!cloneVar || cloneVar->getVal()) {
1837 add(*row,wgt*dset.weight(), wgt*wgt*dset.weightSquared()) ;
1838 }
1839 }
1840
1841 _cache_sum_valid = false ;
1842}
1843
1844
1845
1846////////////////////////////////////////////////////////////////////////////////
1847/// Return the sum of the weights of all bins in the histogram.
1848///
1849/// \param[in] correctForBinSize Multiply the sum of weights in each bin
1850/// with the N-dimensional bin volume, making the return value
1851/// the integral over the function represented by this histogram.
1852/// \param[in] inverseBinCor Divide by the N-dimensional bin volume.
1853double RooDataHist::sum(bool correctForBinSize, bool inverseBinCor) const
1854{
1855 checkInit() ;
1856
1857 // Check if result was cached
1858 const CacheSumState_t cache_code = !correctForBinSize ? kNoBinCorrection : (inverseBinCor ? kInverseBinCorr : kCorrectForBinSize);
1859 if (_cache_sum_valid == static_cast<Int_t>(cache_code)) {
1860 return _cache_sum ;
1861 }
1862
1864 for (Int_t i=0; i < _arrSize; i++) {
1865 const double theBinVolume = correctForBinSize ? (inverseBinCor ? 1/_binv[i] : _binv[i]) : 1.0 ;
1866 kahanSum += _wgt[i] * theBinVolume;
1867 }
1868
1869 // Store result in cache
1870 _cache_sum_valid = cache_code;
1871 _cache_sum = kahanSum.Sum();
1872
1873 return kahanSum.Sum();
1874}
1875
1876
1877
1878////////////////////////////////////////////////////////////////////////////////
1879/// Return the sum of the weights of a multi-dimensional slice of the histogram
1880/// by summing only over the dimensions specified in sumSet.
1881///
1882/// The coordinates of all other dimensions are fixed to those given in sliceSet
1883///
1884/// If correctForBinSize is specified, the sum of weights
1885/// is multiplied by the M-dimensional bin volume, (M = N(sumSet)),
1886/// making the return value the integral over the function
1887/// represented by this histogram
1888
1889double RooDataHist::sum(const RooArgSet& sumSet, const RooArgSet& sliceSet, bool correctForBinSize, bool inverseBinCor)
1890{
1891 checkInit() ;
1892
1893 RooArgSet varSave ;
1894 varSave.addClone(_vars) ;
1895
1896 RooArgSet sliceOnlySet(sliceSet);
1897 sliceOnlySet.remove(sumSet,true,true) ;
1898
1899 _vars.assign(sliceOnlySet);
1900 std::vector<double> const * pbinv = nullptr;
1901
1902 if(correctForBinSize && inverseBinCor) {
1903 pbinv = &calculatePartialBinVolume(sliceOnlySet);
1904 } else if(correctForBinSize && !inverseBinCor) {
1905 pbinv = &calculatePartialBinVolume(sumSet);
1906 }
1907
1908 // Calculate mask and reference plot bins for non-iterating variables
1909 std::vector<bool> mask(_vars.size());
1910 std::vector<int> refBin(_vars.size());
1911
1912 for (unsigned int i = 0; i < _vars.size(); ++i) {
1913 const RooAbsArg* arg = _vars[i];
1914 const RooAbsLValue* argLv = _lvvars[i]; // Same as above, but cross-cast
1915
1916 if (sumSet.find(*arg)) {
1917 mask[i] = false ;
1918 } else {
1919 mask[i] = true ;
1920 refBin[i] = argLv->getBin();
1921 }
1922 }
1923
1924 // Loop over entire data set, skipping masked entries
1926 for (Int_t ibin=0; ibin < _arrSize; ++ibin) {
1927
1928 std::size_t tmpibin = ibin;
1929 bool skip(false) ;
1930
1931 // Check if this bin belongs in selected slice
1932 for (unsigned int ivar = 0; !skip && ivar < _vars.size(); ++ivar) {
1933 const Int_t idx = tmpibin / _idxMult[ivar] ;
1934 tmpibin -= idx*_idxMult[ivar] ;
1935 if (mask[ivar] && idx!=refBin[ivar])
1936 skip = true ;
1937 }
1938
1939 if (!skip) {
1940 const double theBinVolume = correctForBinSize ? (inverseBinCor ? 1/(*pbinv)[ibin] : (*pbinv)[ibin] ) : 1.0 ;
1941 total += _wgt[ibin] * theBinVolume;
1942 }
1943 }
1944
1945 _vars.assign(varSave) ;
1946
1947 return total.Sum();
1948}
1949
1950////////////////////////////////////////////////////////////////////////////////
1951/// Return the sum of the weights of a multi-dimensional slice of the histogram
1952/// by summing only over the dimensions specified in sumSet.
1953///
1954/// The coordinates of all other dimensions are fixed to those given in sliceSet
1955///
1956/// If correctForBinSize is specified, the sum of weights
1957/// is multiplied by the M-dimensional bin volume, (M = N(sumSet)),
1958/// or the fraction of it that falls inside the range rangeName,
1959/// making the return value the integral over the function
1960/// represented by this histogram.
1961///
1962/// If correctForBinSize is not specified, the weights are multiplied by the
1963/// fraction of the bin volume that falls inside the range, i.e. a factor of
1964/// binVolumeInRange/totalBinVolume.
1965
1966double RooDataHist::sum(const RooArgSet& sumSet, const RooArgSet& sliceSet,
1967 bool correctForBinSize, bool inverseBinCor,
1968 const std::map<const RooAbsArg*, std::pair<double, double> >& ranges,
1969 std::function<double(int)> getBinScale)
1970{
1971 checkInit();
1973 RooArgSet varSave;
1974 varSave.addClone(_vars);
1975 {
1976 RooArgSet sliceOnlySet(sliceSet);
1977 sliceOnlySet.remove(sumSet, true, true);
1978 _vars.assign(sliceOnlySet);
1979 }
1980
1981 // Calculate mask and reference plot bins for non-iterating variables,
1982 // and get ranges for iterating variables
1983 std::vector<bool> mask(_vars.size());
1984 std::vector<int> refBin(_vars.size());
1985 std::vector<double> rangeLo(_vars.size(), -std::numeric_limits<double>::infinity());
1986 std::vector<double> rangeHi(_vars.size(), +std::numeric_limits<double>::infinity());
1987
1988 for (std::size_t i = 0; i < _vars.size(); ++i) {
1989 const RooAbsArg* arg = _vars[i];
1990 const RooAbsLValue* argLV = _lvvars[i]; // Same object as above, but cross cast
1991
1992 RooAbsArg* sumsetv = sumSet.find(*arg);
1993 RooAbsArg* slicesetv = sliceSet.find(*arg);
1994 mask[i] = !sumsetv;
1995 if (mask[i]) {
1996 assert(argLV);
1997 refBin[i] = argLV->getBin();
1998 }
1999
2000 auto it = ranges.find(sumsetv ? sumsetv : slicesetv);
2001 if (ranges.end() != it) {
2002 rangeLo[i] = it->second.first;
2003 rangeHi[i] = it->second.second;
2004 }
2005 }
2006
2007 // Loop over entire data set, skipping masked entries
2009 for (Int_t ibin = 0; ibin < _arrSize; ++ibin) {
2010 // Check if this bin belongs in selected slice
2011 bool skip{false};
2012 for (int ivar = 0, tmp = ibin; !skip && ivar < int(_vars.size()); ++ivar) {
2013 const Int_t idx = tmp / _idxMult[ivar];
2014 tmp -= idx*_idxMult[ivar];
2015 if (mask[ivar] && idx!=refBin[ivar]) skip = true;
2016 }
2017
2018 if (skip) continue;
2019
2020 // Work out bin volume
2021 // It's not necessary to figure out the bin volume for the slice-only set explicitly here.
2022 // We need to loop over the sumSet anyway to get the partial bin containment correction,
2023 // so we can get the slice-only set volume later by dividing _binv[ibin] / binVolumeSumSetFull.
2024 double binVolumeSumSetFull = 1.;
2025 double binVolumeSumSetInRange = 1.;
2026 for (Int_t ivar = 0, tmp = ibin; ivar < (int)_vars.size(); ++ivar) {
2027 const Int_t idx = tmp / _idxMult[ivar];
2028 tmp -= idx*_idxMult[ivar];
2029
2030 // If the current variable is not in the sumSet, it should not be considered for the bin volume
2031 const auto arg = _vars[ivar];
2032 if (!sumSet.find(*arg)) {
2033 continue;
2034 }
2035
2036 if (_binbounds[ivar].empty()) continue;
2037 const double binLo = _binbounds[ivar][2 * idx];
2038 const double binHi = _binbounds[ivar][2 * idx + 1];
2039 if (binHi < rangeLo[ivar] || binLo > rangeHi[ivar]) {
2040 // bin is outside of allowed range - effective bin volume is zero
2041 binVolumeSumSetInRange = 0.;
2042 break;
2043 }
2044
2045 binVolumeSumSetFull *= binHi - binLo;
2046 binVolumeSumSetInRange *= std::min(rangeHi[ivar], binHi) - std::max(rangeLo[ivar], binLo);
2047 }
2048 const double corrPartial = binVolumeSumSetInRange / binVolumeSumSetFull;
2049 if (0. == corrPartial) continue;
2050 const double corr = correctForBinSize ? (inverseBinCor ? binVolumeSumSetFull / _binv[ibin] : binVolumeSumSetFull ) : 1.0;
2051 total += getBinScale(ibin)*(_wgt[ibin] * corr * corrPartial);
2052 }
2053
2054 _vars.assign(varSave);
2055
2056 return total.Sum();
2057}
2058
2059
2060
2061////////////////////////////////////////////////////////////////////////////////
2062/// Fill the transient cache with partial bin volumes with up-to-date
2063/// values for the partial volume specified by observables 'dimSet'
2064
2065const std::vector<double>& RooDataHist::calculatePartialBinVolume(const RooArgSet& dimSet) const
2066{
2067 // The code bitset has all bits set to one whose position corresponds to arguments in dimSet.
2068 // It is used as the key for the bin volume caching hash map.
2069 int code{0};
2070 {
2071 int i{0} ;
2072 for (auto const& v : _vars) {
2073 code += ((dimSet.find(*v) ? 1 : 0) << i) ;
2074 ++i;
2075 }
2076 }
2077
2078 auto& pbinv = _pbinvCache[code];
2079 if(!pbinv.empty()) {
2080 return pbinv;
2081 }
2082 pbinv.resize(_arrSize);
2083
2084 // Calculate plot bins of components from master index
2085 std::vector<bool> selDim(_vars.size());
2086 for (std::size_t i = 0; i < selDim.size(); ++i) {
2087 selDim[i] = (code >> i) & 1 ;
2088 }
2089
2090 // Recalculate partial bin volume cache
2091 for (Int_t ibin=0; ibin < _arrSize ;ibin++) {
2092 Int_t idx(0);
2093 Int_t tmp(ibin);
2094 double theBinVolume(1) ;
2095 for (unsigned int j=0; j < _lvvars.size(); ++j) {
2096 const RooAbsLValue* arg = _lvvars[j];
2097 assert(arg);
2098
2099 idx = tmp / _idxMult[j];
2100 tmp -= idx*_idxMult[j];
2101 if (selDim[j]) {
2102 theBinVolume *= arg->getBinWidth(idx) ;
2103 }
2104 }
2105 pbinv[ibin] = theBinVolume ;
2106 }
2107
2108 return pbinv;
2109}
2110
2111
2112////////////////////////////////////////////////////////////////////////////////
2113/// Sum the weights of all bins.
2116}
2117
2118
2119
2120////////////////////////////////////////////////////////////////////////////////
2121/// Return the sum of weights in all entries matching cutSpec (if specified)
2122/// and in named range cutRange (if specified)
2123/// Return the
2124
2125double RooDataHist::sumEntries(const char* cutSpec, const char* cutRange) const
2126{
2127 checkInit() ;
2128
2129 if (cutSpec==nullptr && cutRange==nullptr) {
2130 return sumEntries();
2131 } else {
2132
2133 // Setup RooFormulaVar for cutSpec if it is present
2134 std::unique_ptr<RooFormula> select;
2135 if (cutSpec) {
2136 select = std::make_unique<RooFormula>("select",cutSpec,*get());
2137 }
2138
2139 // Otherwise sum the weights in the event
2140 ROOT::Math::KahanSum<> kahanSum;
2141 for (Int_t i=0; i < _arrSize; i++) {
2142 get(i) ;
2143 if ((select && select->eval() == 0.) || (cutRange && !_vars.allInRange(cutRange)))
2144 continue;
2145
2146 kahanSum += weight(i);
2147 }
2148
2149 return kahanSum.Sum();
2150 }
2151}
2152
2153
2154
2155////////////////////////////////////////////////////////////////////////////////
2156/// Reset all bin weights to zero
2157
2159{
2160 // WVE DO NOT CALL RooTreeData::reset() for binned
2161 // datasets as this will delete the bin definitions
2162
2163 std::fill(_wgt, _wgt + _arrSize, 0.);
2164 delete[] _errLo; _errLo = nullptr;
2165 delete[] _errHi; _errHi = nullptr;
2166 delete[] _sumw2; _sumw2 = nullptr;
2167
2169
2170 _cache_sum_valid = false;
2171}
2172
2173
2174
2175////////////////////////////////////////////////////////////////////////////////
2176/// Load bin `binNumber`, and return an argset with the coordinates of the bin centre.
2177/// \note The argset is owned by this data hist, and this function has a side effect, because
2178/// it alters the currently active bin.
2179const RooArgSet* RooDataHist::get(Int_t binNumber) const
2180{
2181 checkInit() ;
2182 _curIndex = binNumber;
2183
2184 return RooAbsData::get(_curIndex);
2185}
2186
2187
2188
2189////////////////////////////////////////////////////////////////////////////////
2190/// Return a RooArgSet with whose coordinates denote the bin centre of the bin
2191/// enclosing the point in `coord`.
2192/// \note The argset is owned by this data hist, and this function has a side effect, because
2193/// it alters the currently active bin.
2194const RooArgSet* RooDataHist::get(const RooArgSet& coord) const {
2195 return get(calcTreeIndex(coord, false));
2196}
2197
2198
2199
2200////////////////////////////////////////////////////////////////////////////////
2201/// Return the volume of the bin enclosing coordinates 'coord'.
2202double RooDataHist::binVolume(const RooArgSet& coord) const {
2203 checkInit() ;
2204 return _binv[calcTreeIndex(coord, false)] ;
2205}
2206
2207
2208////////////////////////////////////////////////////////////////////////////////
2209/// Create an iterator over all bins in a slice defined by the subset of observables
2210/// listed in sliceArg. The position of the slice is given by otherArgs
2211
2213{
2214 // Update to current position
2215 _vars.assign(otherArgs) ;
2216 _curIndex = calcTreeIndex(_vars, true);
2217
2218 RooAbsArg* intArg = _vars.find(sliceArg) ;
2219 if (!intArg) {
2220 coutE(InputArguments) << "RooDataHist::sliceIterator() variable " << sliceArg.GetName() << " is not part of this RooDataHist" << endl ;
2221 return nullptr ;
2222 }
2223 return new RooDataHistSliceIter(*this,*intArg) ;
2224}
2225
2226
2227////////////////////////////////////////////////////////////////////////////////
2228/// Change the name of the RooDataHist
2229
2231{
2232 if (_dir) _dir->GetList()->Remove(this);
2233 // We need to use the function from RooAbsData, because it already overrides TNamed::SetName
2235 if (_dir) _dir->GetList()->Add(this);
2236}
2237
2238
2239////////////////////////////////////////////////////////////////////////////////
2240/// Change the title of this RooDataHist
2241
2242void RooDataHist::SetNameTitle(const char *name, const char* title)
2243{
2244 SetName(name);
2245 SetTitle(title);
2246}
2247
2248
2249////////////////////////////////////////////////////////////////////////////////
2250/// Print value of the dataset, i.e. the sum of weights contained in the dataset
2251
2252void RooDataHist::printValue(ostream& os) const
2253{
2254 os << numEntries() << " bins (" << sumEntries() << " weights)" ;
2255}
2256
2257
2258
2259
2260////////////////////////////////////////////////////////////////////////////////
2261/// Print argument of dataset, i.e. the observable names
2262
2263void RooDataHist::printArgs(ostream& os) const
2264{
2265 os << "[" ;
2266 bool first(true) ;
2267 for (const auto arg : _vars) {
2268 if (first) {
2269 first=false ;
2270 } else {
2271 os << "," ;
2272 }
2273 os << arg->GetName() ;
2274 }
2275 os << "]" ;
2276}
2277
2278
2279
2280////////////////////////////////////////////////////////////////////////////////
2281/// Returns true if dataset contains entries with a non-integer weight.
2282
2284{
2285 for (Int_t i=0; i < _arrSize; ++i) {
2286 const double wgt = _wgt[i];
2287 double intpart;
2288 if (std::abs(std::modf(wgt, &intpart)) > 1.E-10)
2289 return true;
2290 }
2291
2292 return false;
2293}
2294
2295
2296////////////////////////////////////////////////////////////////////////////////
2297/// Print the details on the dataset contents
2298
2299void RooDataHist::printMultiline(ostream& os, Int_t content, bool verbose, TString indent) const
2300{
2301 RooAbsData::printMultiline(os,content,verbose,indent) ;
2302
2303 os << indent << "Binned Dataset " << GetName() << " (" << GetTitle() << ")" << endl ;
2304 os << indent << " Contains " << numEntries() << " bins with a total weight of " << sumEntries() << endl;
2305
2306 if (!verbose) {
2307 os << indent << " Observables " << _vars << endl ;
2308 } else {
2309 os << indent << " Observables: " ;
2311 }
2312
2313 if(verbose) {
2314 if (!_cachedVars.empty()) {
2315 os << indent << " Caches " << _cachedVars << endl ;
2316 }
2317 }
2318}
2319
2320void RooDataHist::printDataHistogram(ostream& os, RooRealVar* obs) const
2321{
2322 for(Int_t i=0; i<obs->getBins(); ++i){
2323 this->get(i);
2324 obs->setBin(i);
2325 os << this->weight() << " +/- " << this->weightSquared() << endl;
2326 }
2327}
2328
2329
2330////////////////////////////////////////////////////////////////////////////////
2331/// Stream an object of class RooDataHist.
2333 if (R__b.IsReading()) {
2334
2335 UInt_t R__s;
2336 UInt_t R__c;
2337 Version_t R__v = R__b.ReadVersion(&R__s, &R__c);
2338
2339 if (R__v > 2) {
2340 R__b.ReadClassBuffer(RooDataHist::Class(),this,R__v,R__s,R__c);
2341 R__b.CheckByteCount(R__s, R__c, RooDataHist::IsA());
2342 initialize(nullptr, false);
2343 } else {
2344
2345 // Legacy dataset conversion happens here. Legacy RooDataHist inherits from RooTreeData
2346 // which in turn inherits from RooAbsData. Manually stream RooTreeData contents on
2347 // file here and convert it into a RooTreeDataStore which is installed in the
2348 // new-style RooAbsData base class
2349
2350 // --- This is the contents of the streamer code of RooTreeData version 2 ---
2351 UInt_t R__s1;
2352 UInt_t R__c1;
2353 Version_t R__v1 = R__b.ReadVersion(&R__s1, &R__c1); if (R__v1) { }
2354
2356 TTree* X_tree(nullptr) ; R__b >> X_tree;
2357 RooArgSet X_truth ; X_truth.Streamer(R__b);
2358 TString X_blindString ; X_blindString.Streamer(R__b);
2359 R__b.CheckByteCount(R__s1, R__c1, TClass::GetClass("RooTreeData"));
2360 // --- End of RooTreeData-v1 streamer
2361
2362 // Construct RooTreeDataStore from X_tree and complete initialization of new-style RooAbsData
2363 _dstore = std::make_unique<RooTreeDataStore>(X_tree,_vars);
2364 _dstore->SetName(GetName()) ;
2365 _dstore->SetTitle(GetTitle()) ;
2366 _dstore->checkInit() ;
2367
2369 R__b >> _arrSize;
2370 delete [] _wgt;
2371 _wgt = new double[_arrSize];
2373 delete [] _errLo;
2374 _errLo = new double[_arrSize];
2376 delete [] _errHi;
2377 _errHi = new double[_arrSize];
2379 delete [] _sumw2;
2380 _sumw2 = new double[_arrSize];
2382 delete [] _binv;
2383 _binv = new double[_arrSize];
2384 RooArgSet tmpSet;
2385 tmpSet.Streamer(R__b);
2386 double tmp;
2387 R__b >> tmp; //_curWeight;
2388 R__b >> tmp; //_curWgtErrLo;
2389 R__b >> tmp; //_curWgtErrHi;
2390 R__b >> tmp; //_curSumW2;
2391 R__b >> tmp; //_curVolume;
2392 R__b >> _curIndex;
2393 R__b.CheckByteCount(R__s, R__c, RooDataHist::IsA());
2394 }
2395
2396 } else {
2397
2399 }
2400}
2401
2402
2403////////////////////////////////////////////////////////////////////////////////
2404/// Return event weights of all events in range [first, first+len).
2405/// If cacheValidEntries() has been called, out-of-range events will have a weight of 0.
2406std::span<const double> RooDataHist::getWeightBatch(std::size_t first, std::size_t len, bool sumW2 /*=false*/) const {
2407 return {(sumW2 && _sumw2 ? _sumw2 : _wgt) + first, len};
2408}
2409
2410
2411////////////////////////////////////////////////////////////////////////////////
2412/// Hand over pointers to our weight arrays to the data store implementation.
2414 _dstore->setExternalWeightArray(_wgt, _errLo, _errHi, _sumw2);
2415}
2416
2417
2418////////////////////////////////////////////////////////////////////////////////
2419/// Return reference to VarInfo struct with cached histogram variable
2420/// information that is frequently used for histogram weights retrieval.
2421///
2422/// If the `_varInfo` struct was not initialized yet, it will be initialized in
2423/// this function.
2425
2426 if(_varInfo.initialized) return _varInfo;
2427
2428 auto& info = _varInfo;
2429
2430 {
2431 // count the number of real vars and get their indices
2432 info.nRealVars = 0;
2433 size_t iVar = 0;
2434 for (const auto real : _vars) {
2435 if (dynamic_cast<RooRealVar*>(real)) {
2436 if(info.nRealVars == 0) info.realVarIdx1 = iVar;
2437 if(info.nRealVars == 1) info.realVarIdx2 = iVar;
2438 ++info.nRealVars;
2439 }
2440 ++iVar;
2441 }
2442 }
2443
2444 {
2445 // assert that the variables are either real values or categories
2446 for (unsigned int i=0; i < _vars.size(); ++i) {
2447 if (_lvbins[i].get()) {
2448 assert(dynamic_cast<const RooAbsReal*>(_vars[i]));
2449 } else {
2450 assert(dynamic_cast<const RooAbsCategoryLValue*>(_vars[i]));
2451 }
2452 }
2453 }
2454
2455 info.initialized = true;
2456
2457 return info;
2458}
#define e(i)
Definition RSha256.hxx:103
#define coutI(a)
#define coutE(a)
#define TRACE_DESTROY
Definition RooTrace.h:24
#define TRACE_CREATE
Definition RooTrace.h:23
int Int_t
Definition RtypesCore.h:45
short Version_t
Definition RtypesCore.h:65
#define ClassImp(name)
Definition Rtypes.h:382
static void indent(ostringstream &buf, int indent_level)
static unsigned int total
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t mask
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char Pixmap_t Pixmap_t PictureAttributes_t attr const char char ret_data h unsigned char height h offset
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char Pixmap_t Pixmap_t PictureAttributes_t attr const char char ret_data h unsigned char height h Atom_t Int_t ULong_t ULong_t unsigned char prop_list Atom_t Atom_t Atom_t Time_t UChar_t len
Option_t Option_t TPoint TPoint const char x2
Option_t Option_t TPoint TPoint const char x1
Option_t Option_t TPoint TPoint const char y2
Option_t Option_t TPoint TPoint const char y1
char name[80]
Definition TGX11.cxx:110
float xmin
#define hi
float * q
float ymin
Int_t gDebug
Definition TROOT.cxx:597
The Kahan summation is a compensated summation algorithm, which significantly reduces numerical error...
Definition Util.h:122
static KahanSum< T, N > Accumulate(Iterator begin, Iterator end, T initialValue=T{})
Iterate over a range and return an instance of a KahanSum.
Definition Util.h:211
Common abstract base class for objects that represent a value and a "shape" in RooFit.
Definition RooAbsArg.h:79
void attachDataSet(const RooAbsData &set)
Replace server nodes with names matching the dataset variable names with those data set variables,...
Abstract base class for RooRealVar binning definitions.
int binNumber(double x) const
Returns the bin number corresponding to the value x.
virtual void binNumbers(double const *x, int *bins, std::size_t n, int coef=1) const =0
Compute the bin indices for multiple values of x.
virtual double binCenter(Int_t bin) const =0
Int_t numBins() const
Return number of bins.
virtual bool isUniform() const
virtual double highBound() const =0
virtual double lowBound() const =0
virtual RooAbsBinning * clone(const char *name=nullptr) const =0
Abstract base class for objects that represent a discrete value that can be set from the outside,...
Int_t numBins(const char *rangeName=nullptr) const override
Return the number of fit bins ( = number of types )
void setBin(Int_t ibin, const char *rangeName=nullptr) override
Set category to i-th fit bin, which is the i-th registered state.
Int_t getBin(const char *=nullptr) const override
Get the index of the plot bin for the current value of this category.
bool hasLabel(const std::string &label) const
Check if a state with name label exists.
virtual const char * getCurrentLabel() const
Return label string of current state.
Abstract container object that can hold multiple RooAbsArg objects.
virtual bool remove(const RooAbsArg &var, bool silent=false, bool matchByNameOnly=false)
Remove the specified argument from our list.
RooAbsCollection & assignValueOnly(const RooAbsCollection &other, bool forceIfSizeOne=false)
Sets the value of any argument in our set that also appears in the other set.
bool allInRange(const char *rangeSpec) const
Return true if all contained object report to have their value inside the specified range.
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to list.
void assign(const RooAbsCollection &other) const
Sets the value, cache and constant attribute of any argument in our set that also appears in the othe...
Storage_t::size_type size() const
virtual RooAbsArg * addClone(const RooAbsArg &var, bool silent=false)
Add a clone of the specified argument to list.
RooAbsArg * find(const char *name) const
Find object with given name in list.
Abstract base class for binned and unbinned datasets.
Definition RooAbsData.h:57
virtual double weight() const =0
virtual const RooArgSet * get() const
Definition RooAbsData.h:101
void printMultiline(std::ostream &os, Int_t contents, bool verbose=false, TString indent="") const override
Interface for detailed printing of object.
void SetName(const char *name) override
Set the name of the TNamed.
void setGlobalObservables(RooArgSet const &globalObservables)
Sets the global observables stored in this data.
void checkInit() const
virtual RooPlot * plotOn(RooPlot *frame, const RooCmdArg &arg1={}, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={}, const RooCmdArg &arg6={}, const RooCmdArg &arg7={}, const RooCmdArg &arg8={}) const
static StorageType defaultStorageType
Definition RooAbsData.h:312
std::unique_ptr< RooAbsDataStore > _dstore
Data storage implementation.
Definition RooAbsData.h:351
virtual void fill()
RooArgSet _vars
Dimensions of this data set.
Definition RooAbsData.h:348
RooArgSet _cachedVars
! External variables cached with this data set
Definition RooAbsData.h:349
virtual Int_t numEntries() const
Return number of entries in dataset, i.e., count unweighted entries.
virtual double weightSquared() const =0
void Streamer(TBuffer &) override
Stream an object of class RooAbsData.
Abstract base class for objects that are lvalues, i.e.
virtual Int_t getBin(const char *rangeName=nullptr) const =0
virtual double getBinWidth(Int_t i, const char *rangeName=nullptr) const =0
Abstract base class for objects that represent a real value that may appear on the left hand side of ...
virtual Int_t getBins(const char *name=nullptr) const
Get number of bins of currently defined range.
virtual double getMax(const char *name=nullptr) const
Get maximum of currently defined range.
void setBin(Int_t ibin, const char *rangeName=nullptr) override
Set value to center of bin 'ibin' of binning 'rangeName' (or of default binning if no range is specif...
virtual double getMin(const char *name=nullptr) const
Get minimum of currently defined range.
Abstract base class for objects that represent a real value and implements functionality common to al...
Definition RooAbsReal.h:59
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition RooAbsReal.h:103
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition RooArgList.h:22
RooAbsArg * at(Int_t idx) const
Return object at given index, or nullptr if index is out of range.
Definition RooArgList.h:110
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:24
RooArgSet * snapshot(bool deepCopy=true) const
Use RooAbsCollection::snapshot(), but return as RooArgSet.
Definition RooArgSet.h:154
void Streamer(TBuffer &) override
Stream an object of class TObject.
RooArgSet * selectCommon(const RooAbsCollection &refColl) const
Use RooAbsCollection::selecCommon(), but return as RooArgSet.
Definition RooArgSet.h:149
Object to represent discrete states.
Definition RooCategory.h:28
bool defineType(const std::string &label)
Define a state with given name.
Named container for two doubles, two integers two object points and three string pointers that can be...
Definition RooCmdArg.h:26
Configurable parser for RooCmdArg named arguments.
void defineMutex(const char *head, Args_t &&... tail)
Define arguments where any pair is mutually exclusive.
bool process(const RooCmdArg &arg)
Process given RooCmdArg.
double getDouble(const char *name, double defaultValue=0.0) const
Return double property registered with name 'name'.
void defineDependency(const char *refArgName, const char *neededArgName)
Define that processing argument name refArgName requires processing of argument named neededArgName t...
bool defineDouble(const char *name, const char *argName, int doubleNum, double defValue=0.0)
Define double property name 'name' mapped to double in slot 'doubleNum' in RooCmdArg with name argNam...
RooArgSet * getSet(const char *name, RooArgSet *set=nullptr) const
Return RooArgSet property registered with name 'name'.
bool defineSet(const char *name, const char *argName, int setNum, const RooArgSet *set=nullptr)
Define TObject property name 'name' mapped to object in slot 'setNum' in RooCmdArg with name argName ...
bool ok(bool verbose) const
Return true of parsing was successful.
bool defineObject(const char *name, const char *argName, int setNum, const TObject *obj=nullptr, bool isArray=false)
Define TObject property name 'name' mapped to object in slot 'setNum' in RooCmdArg with name argName ...
const char * getString(const char *name, const char *defaultValue="", bool convEmptyToNull=false) const
Return string property registered with name 'name'.
bool defineString(const char *name, const char *argName, int stringNum, const char *defValue="", bool appendMode=false)
Define double property name 'name' mapped to double in slot 'stringNum' in RooCmdArg with name argNam...
const RooLinkedList & getObjectList(const char *name) const
Return list of objects registered with name 'name'.
bool defineInt(const char *name, const char *argName, int intNum, int defValue=0)
Define integer property name 'name' mapped to integer in slot 'intNum' in RooCmdArg with name argName...
int getInt(const char *name, int defaultValue=0) const
Return integer property registered with name 'name'.
TObject * getObject(const char *name, TObject *obj=nullptr) const
Return TObject property registered with name 'name'.
Container class to hold N-dimensional binned data.
Definition RooDataHist.h:40
std::span< const double > getWeightBatch(std::size_t first, std::size_t len, bool sumW2=false) const override
Return event weights of all events in range [first, first+len).
void interpolateQuadratic(double *output, std::span< const double > xVals, bool correctForBinSize, bool cdfBoundaries)
A vectorized version of interpolateDim for boundary safe quadratic interpolation of one dimensional h...
double sum(bool correctForBinSize, bool inverseCorr=false) const
Return the sum of the weights of all bins in the histogram.
void weights(double *output, std::span< double const > xVals, int intOrder, bool correctForBinSize, bool cdfBoundaries)
A vectorized version of RooDataHist::weight() for one dimensional histograms with up to one dimension...
Int_t _cache_sum_valid
! Is cache sum valid? Needs to be Int_t instead of CacheSumState_t for subclasses.
double interpolateDim(int iDim, double xval, size_t centralIdx, int intOrder, bool correctForBinSize, bool cdfBoundaries)
Perform boundary safe 'intOrder'-th interpolation of weights in dimension 'dim' at current value 'xva...
double weightSquared() const override
Return squared weight of last bin that was requested with get().
friend class RooDataHistSliceIter
void importTH1(const RooArgList &vars, const TH1 &histo, double initWgt, bool doDensityCorrection)
Import data from given TH1/2/3 into this RooDataHist.
void printDataHistogram(std::ostream &os, RooRealVar *obs) const
static TClass * Class()
TClass * IsA() const override
void SetNameTitle(const char *name, const char *title) override
Change the title of this RooDataHist.
double _cache_sum
! Cache for sum of entries ;
void initialize(const char *binningName=nullptr, bool fillTree=true)
Initialization procedure: allocate weights array, calculate multipliers needed for N-space to 1-dim a...
std::string calculateTreeIndexForCodeSquash(RooAbsArg const *klass, RooFit::Detail::CodeSquashContext &ctx, const RooAbsCollection &coords, bool reverse=false) const
VarInfo _varInfo
!
Int_t getIndex(const RooAbsCollection &coord, bool fast=false) const
Calculate bin number of the given coordinates.
void add(const RooArgSet &row, double wgt=1.0) override
Add wgt to the bin content enclosed by the coordinates passed in row.
Definition RooDataHist.h:72
const std::vector< double > & calculatePartialBinVolume(const RooArgSet &dimSet) const
Fill the transient cache with partial bin volumes with up-to-date values for the partial volume speci...
static std::unique_ptr< RooAbsDataStore > makeDefaultDataStore(RooStringView name, RooStringView title, RooArgSet const &vars)
double weightInterpolated(const RooArgSet &bin, int intOrder, bool correctForBinSize, bool cdfBoundaries)
Return the weight at given coordinates with interpolation.
std::unordered_map< int, std::vector< double > > _pbinvCache
! Cache for arrays of partial bin volumes
void checkBinBounds() const
void initializeAsymErrArrays() const
double weight(std::size_t i) const
Return weight of i-th bin.
void set(std::size_t binNumber, double weight, double wgtErr)
Set bin content of bin that was last loaded with get(std::size_t).
void weightError(double &lo, double &hi, ErrorType etype=Poisson) const override
Return the asymmetric errors on the current weight.
double * _errHi
[_arrSize] High-side error on weight array
void importTH1Set(const RooArgList &vars, RooCategory &indexCat, std::map< std::string, TH1 * > hmap, double initWgt, bool doDensityCorrection)
Import data from given set of TH1/2/3 into this RooDataHist.
void adjustBinning(const RooArgList &vars, const TH1 &href, Int_t *offset=nullptr)
Adjust binning specification on first and optionally second and third observable to binning in given ...
double * _binv
[_arrSize] Bin volume array
RooDataHist()
Default constructor.
ULong64_t _curIndex
Current index.
double weightFast(const RooArgSet &bin, int intOrder, bool correctForBinSize, bool cdfBoundaries)
A faster version of RooDataHist::weight that assumes the passed arguments are aligned with the histog...
RooPlot * plotOn(RooPlot *frame, PlotOpt o) const override
Back end function to plotting functionality.
double weight() const override
Return weight of last bin that was requested with get().
std::vector< std::vector< double > > _binbounds
! list of bin bounds per dimension
void printArgs(std::ostream &os) const override
Print argument of dataset, i.e. the observable names.
void importDHistSet(const RooArgList &vars, RooCategory &indexCat, std::map< std::string, RooDataHist * > dmap, double initWgt)
Import data from given set of TH1/2/3 into this RooDataHist.
void _adjustBinning(RooRealVar &theirVar, const TAxis &axis, RooRealVar *ourVar, Int_t *offset)
Helper doing the actual work of adjustBinning().
void printMultiline(std::ostream &os, Int_t content, bool verbose=false, TString indent="") const override
Print the details on the dataset contents.
double * _sumw2
[_arrSize] Sum of weights^2
TIterator * sliceIterator(RooAbsArg &sliceArg, const RooArgSet &otherArgs)
Create an iterator over all bins in a slice defined by the subset of observables listed in sliceArg.
Int_t calcTreeIndex() const
Legacy overload to calculate the tree index from the current value of _vars.
~RooDataHist() override
Destructor.
std::string declWeightArrayForCodeSquash(RooFit::Detail::CodeSquashContext &ctx, bool correctForBinSize) const
bool isNonPoissonWeighted() const override
Returns true if dataset contains entries with a non-integer weight.
std::vector< RooAbsLValue * > _lvvars
! List of observables casted as RooAbsLValue
void SetName(const char *name) override
Change the name of the RooDataHist.
std::vector< std::unique_ptr< const RooAbsBinning > > _lvbins
! List of used binnings associated with lvalues
void Streamer(TBuffer &) override
Stream an object of class RooDataHist.
std::vector< double > _interpolationBuffer
! Buffer to contain values used for weight interpolation
std::vector< Int_t > _idxMult
void registerWeightArraysToDataStore() const
Hand over pointers to our weight arrays to the data store implementation.
void reset() override
Reset all bin weights to zero.
double * _errLo
[_arrSize] Low-side error on weight array
double * _wgt
[_arrSize] Weight array
void printValue(std::ostream &os) const override
Print value of the dataset, i.e. the sum of weights contained in the dataset.
VarInfo const & getVarInfo()
Return reference to VarInfo struct with cached histogram variable information that is frequently used...
std::unique_ptr< RooAbsData > reduceEng(const RooArgSet &varSubset, const RooFormulaVar *cutVar, const char *cutRange=nullptr, std::size_t nStart=0, std::size_t nStop=std::numeric_limits< std::size_t >::max()) const override
Implementation of RooAbsData virtual method that drives the RooAbsData::reduce() methods.
const RooArgSet * get() const override
Get bin centre of current bin.
Definition RooDataHist.h:82
void interpolateLinear(double *output, std::span< const double > xVals, bool correctForBinSize, bool cdfBoundaries)
A vectorized version of interpolateDim for boundary safe linear interpolation of one dimensional hist...
double binVolume() const
Return volume of current bin.
double sumEntries() const override
Sum the weights of all bins.
Utility base class for RooFit objects that are to be attached to ROOT directories.
Definition RooDirItem.h:22
virtual void Streamer(TBuffer &)
void removeFromDir(TObject *obj)
Remove object from directory it was added to.
TDirectory * _dir
! Associated directory
Definition RooDirItem.h:33
void appendToDir(TObject *obj, bool forceMemoryResident=false)
Append object to directory.
A class to maintain the context for squashing of RooFit models into code.
std::string buildCall(std::string const &funcname, Args_t const &...args)
Build the code to call the function with name funcname, passing some arguments.
std::string buildArg(RooAbsCollection const &x)
Function to save a RooListProxy as an array in the squashed code.
A RooFormulaVar is a generic implementation of a real-valued object, which takes a RooArgList of serv...
static const RooHistError & instance()
Return a reference to a singleton object that is created the first time this method is called.
bool getPoissonInterval(Int_t n, double &mu1, double &mu2, double nSigma=1) const
Return a confidence interval for the expected number of events given n observed (unweighted) events.
Collection class for internal use, storing a collection of RooAbsArg pointers in a doubly linked list...
RooLinkedListIterImpl begin() const
static double interpolate(double yArr[], Int_t nOrder, double x)
Definition RooMath.cxx:78
Plot frame and a container for graphics objects within that frame.
Definition RooPlot.h:45
RooAbsRealLValue * getPlotVar() const
Definition RooPlot.h:143
virtual void printStream(std::ostream &os, Int_t contents, StyleOption style, TString indent="") const
Print description of object on ostream, printing contents set by contents integer,...
Variable that can be changed from the outside.
Definition RooRealVar.h:37
const RooAbsBinning & getBinning(const char *name=nullptr, bool verbose=true, bool createOnTheFly=false) const override
Return binning definition with name.
void setBinning(const RooAbsBinning &binning, const char *name=nullptr)
Add given binning under name 'name' with this variable.
void setRange(const char *name, double min, double max)
Set a fit or plotting range.
The RooStringView is a wrapper around a C-style string that can also be constructed from a std::strin...
Implementation of RooAbsBinning that provides a uniform binning in 'n' bins between the range end poi...
const Double_t * GetArray() const
Definition TArrayD.h:43
Class to manage histogram axis.
Definition TAxis.h:31
const TArrayD * GetXbins() const
Definition TAxis.h:136
Double_t GetXmax() const
Definition TAxis.h:140
virtual Int_t FindFixBin(Double_t x) const
Find bin number corresponding to abscissa x.
Definition TAxis.cxx:419
Double_t GetXmin() const
Definition TAxis.h:139
Int_t GetNbins() const
Definition TAxis.h:125
Buffer base class used for serializing objects.
Definition TBuffer.h:43
virtual Version_t ReadVersion(UInt_t *start=nullptr, UInt_t *bcnt=nullptr, const TClass *cl=nullptr)=0
virtual Int_t CheckByteCount(UInt_t startpos, UInt_t bcnt, const TClass *clss)=0
virtual void ReadFastArray(Bool_t *b, Int_t n)=0
virtual Int_t ReadClassBuffer(const TClass *cl, void *pointer, const TClass *onfile_class=nullptr)=0
Bool_t IsReading() const
Definition TBuffer.h:86
virtual Int_t WriteClassBuffer(const TClass *cl, void *pointer)=0
static TClass * GetClass(const char *name, Bool_t load=kTRUE, Bool_t silent=kFALSE)
Static method returning pointer to TClass of the specified class name.
Definition TClass.cxx:2968
virtual TList * GetList() const
Definition TDirectory.h:222
TH1 is the base class of all histogram classes in ROOT.
Definition TH1.h:59
TAxis * GetZaxis()
Definition TH1.h:326
virtual Int_t GetNbinsY() const
Definition TH1.h:298
virtual Double_t GetBinError(Int_t bin) const
Return value of error associated to bin number bin.
Definition TH1.cxx:9054
virtual Int_t GetNbinsZ() const
Definition TH1.h:299
virtual Int_t GetDimension() const
Definition TH1.h:283
TAxis * GetXaxis()
Definition TH1.h:324
virtual Int_t GetNbinsX() const
Definition TH1.h:297
TAxis * GetYaxis()
Definition TH1.h:325
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition TH1.cxx:5052
Iterator abstract base class.
Definition TIterator.h:30
void Add(TObject *obj) override
Definition TList.h:83
TObject * Remove(TObject *obj) override
Remove object from the list.
Definition TList.cxx:820
virtual void SetTitle(const char *title="")
Set the title of the TNamed.
Definition TNamed.cxx:164
const char * GetName() const override
Returns name of object.
Definition TNamed.h:47
const char * GetTitle() const override
Returns title of object.
Definition TNamed.h:48
Mother of all ROOT objects.
Definition TObject.h:41
virtual const char * ClassName() const
Returns name of class to which the object belongs.
Definition TObject.cxx:213
Basic string class.
Definition TString.h:139
virtual void Streamer(TBuffer &)
Stream a string object.
Definition TString.cxx:1412
A TTree represents a columnar dataset.
Definition TTree.h:79
const Int_t n
Definition legend1.C:16
TH1F * h1
Definition legend1.C:5
std::vector< std::string > Split(std::string_view str, std::string_view delims, bool skipEmpty=false)
Splits a string at each character in delims.
RooAbsBinning * bins
Definition RooAbsData.h:181
Structure to cache information on the histogram variable that is frequently used for histogram weight...
TLine l
Definition textangle.C:4
static void output()