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)
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 if (_arrSize == 1)
999 return std::to_string(_wgt[0]) + (correctForBinSize ? " / " + std::to_string(_binv[0]) : "");
1000 std::string weightName = +"_HW" + RooFit::Detail::makeValidVarName(klass->GetName());
1001
1002 std::string arrayDecl = "double " + weightName + "[" + std::to_string(_arrSize) + "] = {";
1003 for (Int_t i = 0; i < _arrSize; i++) {
1004 arrayDecl += " " + std::to_string(_wgt[i]) + (correctForBinSize ? " / " + std::to_string(_binv[i]) : "") + ",";
1005 }
1006 arrayDecl.back() = ' ';
1007 arrayDecl += "};\n";
1008
1009 ctx.addToCodeBody(arrayDecl, true);
1010
1011 return weightName;
1012}
1013
1015 const RooAbsCollection &coords, bool reverse) const
1016{
1017 assert(coords.size() == _vars.size());
1018
1019 std::string idxName = "idx_" + ctx.getTmpVarName();
1020 ctx.addToCodeBody(klass, "unsigned int " + idxName + " = 0;\n");
1021 int idxMult = 1;
1022
1023 for (std::size_t i = 0; i < _vars.size(); ++i) {
1024
1025 std::size_t iVar = reverse ? _vars.size() - 1 - i : i;
1026 const RooAbsArg *internalVar = _vars[iVar];
1027 const RooAbsArg *theVar = coords[iVar];
1028
1029 const RooAbsBinning *binning = _lvbins[iVar].get();
1030 if (!binning) {
1031 coutE(InputArguments) << "RooHistPdf::weight(" << GetName()
1032 << ") ERROR: Code Squashing currently does not support category values." << std::endl;
1033 return "";
1034 } else if (!dynamic_cast<RooUniformBinning const *>(binning)) {
1035 coutE(InputArguments) << "RooHistPdf::weight(" << GetName()
1036 << ") ERROR: Code Squashing currently only supports uniformly binned cases."
1037 << std::endl;
1038 return "";
1039 }
1040
1041 std::string const &bin = ctx.buildCall("RooFit::Detail::MathFuncs::getUniformBinning", binning->lowBound(),
1042 binning->highBound(), *theVar, binning->numBins());
1043 ctx.addToCodeBody(klass, idxName + " += " + std::to_string(idxMult) + " * " + bin + ";\n");
1044
1045 // Use RooAbsLValue here because it also generalized to categories, which
1046 // is useful in the future. dynamic_cast because it's a cross-cast.
1047 idxMult *= dynamic_cast<RooAbsLValue const *>(internalVar)->numBins();
1048 }
1049
1050 return idxName;
1051}
1052
1053////////////////////////////////////////////////////////////////////////////////
1054/// Calculate the bin index corresponding to the coordinates passed as argument.
1055/// \param[in] coords Coordinates. If `fast == false`, these can be partial.
1056/// \param[in] fast Promise that the coordinates in `coords` have the same order
1057/// as the internal coordinates. In this case, values are looked up only by index.
1058std::size_t RooDataHist::calcTreeIndex(const RooAbsCollection& coords, bool fast) const
1059{
1060 // With fast, caller promises that layout of `coords` is identical to our internal `vars`.
1061 // Previously, this was verified with an assert in debug mode like this:
1062 //
1063 // assert(!fast || coords.hasSameLayout(_vars));
1064 //
1065 // However, there are usecases where the externally provided `coords` have
1066 // different names than the internal variables, even though they correspond
1067 // to each other. For example, if the observables in the computation graph
1068 // are renamed with `redirectServers`. Hence, we can't do a meaningful assert
1069 // here.
1070
1071 if (&_vars == &coords)
1072 fast = true;
1073
1074 std::size_t masterIdx = 0;
1075
1076 for (unsigned int i=0; i < _vars.size(); ++i) {
1077 const RooAbsArg* internalVar = _vars[i];
1078 const RooAbsBinning* binning = _lvbins[i].get();
1079
1080 // Find the variable that we need values from.
1081 // That's either the variable directly from the external coordinates
1082 // or we find the external one that has the same name as "internalVar".
1083 const RooAbsArg* theVar = fast ? coords[i] : coords.find(*internalVar);
1084 if (!theVar) {
1085 // Variable is not in external coordinates. Use current internal value.
1086 theVar = internalVar;
1087 }
1088 // If fast is on, users promise that the sets have the same layout:
1089 //
1090 // assert(!fast || strcmp(internalVar->GetName(), theVar->GetName()) == 0);
1091 //
1092 // This assert is commented out for the same reasons that applied to the
1093 // other assert explained above.
1094
1095 if (binning) {
1096 assert(dynamic_cast<const RooAbsReal*>(theVar));
1097 const double val = static_cast<const RooAbsReal*>(theVar)->getVal();
1098 masterIdx += _idxMult[i] * binning->binNumber(val);
1099 } else {
1100 // We are a category. No binning.
1101 assert(dynamic_cast<const RooAbsCategoryLValue*>(theVar));
1102 auto cat = static_cast<const RooAbsCategoryLValue*>(theVar);
1103 masterIdx += _idxMult[i] * cat->getBin(static_cast<const char*>(nullptr));
1104 }
1105 }
1106
1107 return masterIdx ;
1108}
1109
1110
1111////////////////////////////////////////////////////////////////////////////////
1112/// Back end function to plotting functionality. Plot RooDataHist on given
1113/// frame in mode specified by plot options 'o'. The main purpose of
1114/// this function is to match the specified binning on 'o' to the
1115/// internal binning of the plot observable in this RooDataHist.
1116/// \see RooAbsData::plotOn() for plotting options.
1118{
1119 checkInit() ;
1120 if (o.bins) return RooAbsData::plotOn(frame,o) ;
1121
1122 if(!frame) {
1123 coutE(InputArguments) << ClassName() << "::" << GetName() << ":plotOn: frame is null" << endl;
1124 return nullptr;
1125 }
1126 auto var= static_cast<RooAbsRealLValue*>(frame->getPlotVar());
1127 if(!var) {
1128 coutE(InputArguments) << ClassName() << "::" << GetName()
1129 << ":plotOn: frame does not specify a plot variable" << endl;
1130 return nullptr;
1131 }
1132
1133 auto dataVar = static_cast<RooRealVar*>(_vars.find(*var));
1134 if (!dataVar) {
1135 coutE(InputArguments) << ClassName() << "::" << GetName()
1136 << ":plotOn: dataset doesn't contain plot frame variable" << endl;
1137 return nullptr;
1138 }
1139
1140 o.bins = &dataVar->getBinning() ;
1141 return RooAbsData::plotOn(frame,o) ;
1142}
1143
1144
1145////////////////////////////////////////////////////////////////////////////////
1146/// A vectorized version of interpolateDim for boundary safe quadratic
1147/// interpolation of one dimensional histograms.
1148///
1149/// \param[out] output An array of interpolated weights corresponding to the
1150/// values in xVals.
1151/// \param[in] xVals An array of event coordinates for which the weights should be
1152/// calculated.
1153/// \param[in] correctForBinSize Enable the inverse bin volume correction factor.
1154/// \param[in] cdfBoundaries Enable the special boundary condition for a cdf:
1155/// Underflow bins are assumed to have weight zero and
1156/// overflow bins have weight one. Otherwise, the
1157/// histogram is mirrored at the boundaries for the
1158/// interpolation.
1159
1160void RooDataHist::interpolateQuadratic(double* output, std::span<const double> xVals,
1161 bool correctForBinSize, bool cdfBoundaries)
1162{
1163 const std::size_t nBins = numEntries();
1164 const std::size_t nEvents = xVals.size();
1165
1166 RooAbsBinning const& binning = *_lvbins[0];
1167 // Reuse the output buffer for bin indices and zero-initialize it
1168 auto binIndices = reinterpret_cast<int*>(output + nEvents) - nEvents;
1169 std::fill(binIndices, binIndices + nEvents, 0);
1170 binning.binNumbers(xVals.data(), binIndices, nEvents);
1171
1172 // Extend coordinates and weights with one extra point before the first bin
1173 // and one extra point after the last bin. This means the original histogram
1174 // bins span elements 1 to nBins in coordsExt and weightsExt
1175 std::vector<double> coordsExt(nBins+3);
1176 double* binCoords = coordsExt.data() + 2;
1177 binCoords[0] = binning.lowBound() + 0.5*_binv[0];
1178 for (std::size_t binIdx = 1; binIdx < nBins ; ++binIdx) {
1179 if (binning.isUniform()) {
1180 double binWidth = _binv[0];
1181 binCoords[binIdx] = binIdx*binWidth + binCoords[0];
1182 }
1183 else {
1184 double binCentDiff = 0.5*_binv[binIdx-1] + 0.5*_binv[binIdx];
1185 binCoords[binIdx] = binCoords[binIdx-1] + binCentDiff;
1186 }
1187 }
1188
1189 std::vector<double> weightsExt(nBins+3);
1190 // Fill weights for bins that are inside histogram boundaries
1191 for (std::size_t binIdx = 0; binIdx < nBins; ++binIdx) {
1192 weightsExt[binIdx+2] = correctForBinSize ? _wgt[binIdx] / _binv[binIdx] : _wgt[binIdx];
1193 }
1194
1195 if (cdfBoundaries) {
1196 coordsExt[0] = - 1e-10 + binning.lowBound();
1197 weightsExt[0] = 0.;
1198
1199 coordsExt[1] = binning.lowBound();
1200 weightsExt[1] = 0.;
1201
1202 coordsExt[nBins+2] = binning.highBound();
1203 weightsExt[nBins+2] = 1.;
1204 }
1205 else {
1206 // Mirror first two bins and last bin
1207 coordsExt[0] = binCoords[1] - 2*_binv[0] - _binv[1];
1208 weightsExt[0] = weightsExt[3];
1209
1210 coordsExt[1] = binCoords[0] - _binv[0];
1211 weightsExt[1] = weightsExt[2];
1212
1213 coordsExt[nBins+2] = binCoords[nBins-1] + _binv[nBins-1];
1214 weightsExt[nBins+2] = weightsExt[nBins+1];
1215 }
1216
1217 // We use the current bin center and two bin centers on the left for
1218 // interpolation if xVal is to the left of the current bin center
1219 for (std::size_t i = 0; i < nEvents ; ++i) {
1220 double xVal = xVals[i];
1221 std::size_t binIdx = binIndices[i] + 2;
1222
1223 // If xVal is to the right of the current bin center, shift all bin
1224 // coordinates one step to the right and use that for the interpolation
1225 if (xVal > coordsExt[binIdx]) {
1226 binIdx += 1;
1227 }
1228
1229 double x1 = coordsExt[binIdx-2];
1230 double y1 = weightsExt[binIdx-2];
1231
1232 double x2 = coordsExt[binIdx-1];
1233 double y2 = weightsExt[binIdx-1];
1234
1235 double x3 = coordsExt[binIdx];
1236 double y3 = weightsExt[binIdx];
1237
1238 // Evaluate a few repeated factors
1239 double quotient = (x3-x1) / (x2-x1);
1240 double x1Sqrd = x1*x1;
1241 double x3Sqrd = x3*x3;
1242 // Solve coefficients in system of three quadratic equations!
1243 double secondCoeff = (y3 - y1 - (y2-y1) * quotient) / (x3Sqrd - x1Sqrd - (x2*x2 - x1Sqrd) * quotient);
1244 double firstCoeff = (y3 - y1 - secondCoeff*(x3Sqrd - x1Sqrd)) / (x3-x1);
1245 double zerothCoeff = y1 - secondCoeff*x1Sqrd - firstCoeff*x1;
1246 // Get the interpolated weight using the equation of a second degree polynomial
1247 output[i] = secondCoeff * xVal * xVal + firstCoeff * xVal + zerothCoeff;
1248 }
1249}
1250
1251
1252////////////////////////////////////////////////////////////////////////////////
1253/// A vectorized version of interpolateDim for boundary safe linear
1254/// interpolation of one dimensional histograms.
1255///
1256/// \param[out] output An array of interpolated weights corresponding to the
1257/// values in xVals.
1258/// \param[in] xVals An array of event coordinates for which the weights should be
1259/// calculated.
1260/// \param[in] correctForBinSize Enable the inverse bin volume correction factor.
1261/// \param[in] cdfBoundaries Enable the special boundary condition for a cdf:
1262/// Underflow bins are assumed to have weight zero and
1263/// overflow bins have weight one. Otherwise, the
1264/// histogram is mirrored at the boundaries for the
1265/// interpolation.
1266
1267void RooDataHist::interpolateLinear(double* output, std::span<const double> xVals,
1268 bool correctForBinSize, bool cdfBoundaries)
1269{
1270 const std::size_t nBins = numEntries();
1271 const std::size_t nEvents = xVals.size();
1272
1273 RooAbsBinning const& binning = *_lvbins[0];
1274 // Reuse the output buffer for bin indices and zero-initialize it
1275 auto binIndices = reinterpret_cast<int*>(output + nEvents) - nEvents;
1276 std::fill(binIndices, binIndices + nEvents, 0);
1277 binning.binNumbers(xVals.data(), binIndices, nEvents);
1278
1279 // Extend coordinates and weights with one extra point before the first bin
1280 // and one extra point after the last bin. This means the original histogram
1281 // bins span elements 1 to nBins in coordsExt and weightsExt
1282 std::vector<double> coordsExt(nBins+2);
1283 double* binCoords = coordsExt.data() + 1;
1284 binCoords[0] = binning.lowBound() + 0.5*_binv[0];
1285 for (std::size_t binIdx = 1; binIdx < nBins ; ++binIdx) {
1286 if (binning.isUniform()) {
1287 double binWidth = _binv[0];
1288 binCoords[binIdx] = binIdx*binWidth + binCoords[0];
1289 }
1290 else {
1291 double binCentDiff = 0.5*_binv[binIdx-1] + 0.5*_binv[binIdx];
1292 binCoords[binIdx] = binCoords[binIdx-1] + binCentDiff;
1293 }
1294 }
1295
1296 std::vector<double> weightsExt(nBins+2);
1297 // Fill weights for bins that are inside histogram boundaries
1298 for (std::size_t binIdx = 0; binIdx < nBins; ++binIdx) {
1299 weightsExt[binIdx+1] = correctForBinSize ? _wgt[binIdx] / _binv[binIdx] : _wgt[binIdx];
1300 }
1301
1302 // Fill weights for bins that are outside histogram boundaries
1303 if (cdfBoundaries) {
1304 coordsExt[0] = binning.lowBound();
1305 weightsExt[0] = 0.;
1306 coordsExt[nBins+1] = binning.highBound();
1307 weightsExt[nBins+1] = 1.;
1308 }
1309 else {
1310 // Mirror first and last bins
1311 coordsExt[0] = binCoords[0] - _binv[0];
1312 weightsExt[0] = weightsExt[1];
1313 coordsExt[nBins+1] = binCoords[nBins-1] + _binv[nBins-1];
1314 weightsExt[nBins+1] = weightsExt[nBins];
1315 }
1316
1317 // Interpolate between current bin center and one bin center to the left
1318 // if xVal is to the left of the current bin center
1319 for (std::size_t i = 0; i < nEvents ; ++i) {
1320 double xVal = xVals[i];
1321 std::size_t binIdx = binIndices[i] + 1;
1322
1323 // If xVal is to the right of the current bin center, interpolate between
1324 // current bin center and one bin center to the right instead
1325 if (xVal > coordsExt[binIdx]) { binIdx += 1; }
1326
1327 double x1 = coordsExt[binIdx-1];
1328 double y1 = weightsExt[binIdx-1];
1329 double x2 = coordsExt[binIdx];
1330 double y2 = weightsExt[binIdx];
1331
1332 // Find coefficients by solving a system of two linear equations
1333 double firstCoeff = (y2-y1) / (x2-x1);
1334 double zerothCoeff = y1 - firstCoeff * x1;
1335 // Get the interpolated weight using the equation of a straight line
1336 output[i] = firstCoeff * xVal + zerothCoeff;
1337 }
1338}
1339
1340
1341////////////////////////////////////////////////////////////////////////////////
1342/// A vectorized version of RooDataHist::weight() for one dimensional histograms
1343/// with up to one dimensional interpolation.
1344/// \param[out] output An array of weights corresponding the values in xVals.
1345/// \param[in] xVals An array of coordinates for which the weights should be
1346/// calculated.
1347/// \param[in] intOrder Interpolation order; 0th and 1st order are supported.
1348/// \param[in] correctForBinSize Enable the inverse bin volume correction factor.
1349/// \param[in] cdfBoundaries Enable the special boundary condition for a cdf:
1350/// Underflow bins are assumed to have weight zero and
1351/// overflow bins have weight one. Otherwise, the
1352/// histogram is mirrored at the boundaries for the
1353/// interpolation.
1354
1355void RooDataHist::weights(double* output, std::span<double const> xVals, int intOrder, bool correctForBinSize, bool cdfBoundaries)
1356{
1357 auto const nEvents = xVals.size();
1358
1359 if (intOrder == 0) {
1360 RooAbsBinning const& binning = *_lvbins[0];
1361
1362 // Reuse the output buffer for bin indices and zero-initialize it
1363 auto binIndices = reinterpret_cast<int*>(output + nEvents) - nEvents;
1364 std::fill(binIndices, binIndices + nEvents, 0);
1365 binning.binNumbers(xVals.data(), binIndices, nEvents);
1366
1367 for (std::size_t i=0; i < nEvents; ++i) {
1368 auto binIdx = binIndices[i];
1369 output[i] = correctForBinSize ? _wgt[binIdx] / _binv[binIdx] : _wgt[binIdx];
1370 }
1371 }
1372 else if (intOrder == 1) {
1373 interpolateLinear(output, xVals, correctForBinSize, cdfBoundaries);
1374 }
1375 else if (intOrder == 2) {
1376 interpolateQuadratic(output, xVals, correctForBinSize, cdfBoundaries);
1377 }
1378 else {
1379 // Higher dimensional scenarios not yet implemented
1380 coutE(InputArguments) << "RooDataHist::weights(" << GetName() << ") interpolation in "
1381 << intOrder << " dimensions not yet implemented" << std::endl ;
1382 // Fall back to 1st order interpolation
1383 weights(output, xVals, 1, correctForBinSize, cdfBoundaries);
1384 }
1385}
1386
1387
1388////////////////////////////////////////////////////////////////////////////////
1389/// A faster version of RooDataHist::weight that assumes the passed arguments
1390/// are aligned with the histogram variables.
1391/// \param[in] bin Coordinates for which the weight should be calculated.
1392/// Has to be aligned with the internal histogram variables.
1393/// \param[in] intOrder Interpolation order, i.e. how many neighbouring bins are
1394/// used for the interpolation. If zero, the bare weight for
1395/// the bin enclosing the coordinatesis returned.
1396/// \param[in] correctForBinSize Enable the inverse bin volume correction factor.
1397/// \param[in] cdfBoundaries Enable the special boundary condition for a cdf:
1398/// underflow bins are assumed to have weight zero and
1399/// overflow bins have weight one. Otherwise, the
1400/// histogram is mirrored at the boundaries for the
1401/// interpolation.
1402
1403double RooDataHist::weightFast(const RooArgSet& bin, Int_t intOrder, bool correctForBinSize, bool cdfBoundaries)
1404{
1405 checkInit() ;
1406
1407 // Handle illegal intOrder values
1408 if (intOrder<0) {
1409 coutE(InputArguments) << "RooDataHist::weight(" << GetName() << ") ERROR: interpolation order must be positive" << endl ;
1410 return 0 ;
1411 }
1412
1413 // Handle no-interpolation case
1414 if (intOrder==0) {
1415 const auto idx = calcTreeIndex(bin, true);
1416 return correctForBinSize ? _wgt[idx] / _binv[idx] : _wgt[idx];
1417 }
1418
1419 // Handle all interpolation cases
1420 return weightInterpolated(bin, intOrder, correctForBinSize, cdfBoundaries);
1421}
1422
1423
1424////////////////////////////////////////////////////////////////////////////////
1425/// Return the weight at given coordinates with optional interpolation.
1426/// \param[in] bin Coordinates for which the weight should be calculated.
1427/// \param[in] intOrder Interpolation order, i.e. how many neighbouring bins are
1428/// used for the interpolation. If zero, the bare weight for
1429/// the bin enclosing the coordinatesis returned.
1430/// \param[in] correctForBinSize Enable the inverse bin volume correction factor.
1431/// \param[in] cdfBoundaries Enable the special boundary condition for a cdf:
1432/// underflow bins are assumed to have weight zero and
1433/// overflow bins have weight one. Otherwise, the
1434/// histogram is mirrored at the boundaries for the
1435/// interpolation.
1436/// \param[in] oneSafe Ignored.
1437
1438double RooDataHist::weight(const RooArgSet& bin, Int_t intOrder, bool correctForBinSize, bool cdfBoundaries, bool /*oneSafe*/)
1439{
1440 checkInit() ;
1441
1442 // Handle illegal intOrder values
1443 if (intOrder<0) {
1444 coutE(InputArguments) << "RooDataHist::weight(" << GetName() << ") ERROR: interpolation order must be positive" << endl ;
1445 return 0 ;
1446 }
1447
1448 // Handle no-interpolation case
1449 if (intOrder==0) {
1450 const auto idx = calcTreeIndex(bin, false);
1451 return correctForBinSize ? _wgt[idx] / _binv[idx] : _wgt[idx];
1452 }
1453
1454 // Handle all interpolation cases
1455 _vars.assignValueOnly(bin) ;
1456
1457 return weightInterpolated(_vars, intOrder, correctForBinSize, cdfBoundaries);
1458}
1459
1460
1461////////////////////////////////////////////////////////////////////////////////
1462/// Return the weight at given coordinates with interpolation.
1463/// \param[in] bin Coordinates for which the weight should be calculated.
1464/// Has to be aligned with the internal histogram variables.
1465/// \param[in] intOrder Interpolation order, i.e. how many neighbouring bins are
1466/// used for the interpolation.
1467/// \param[in] correctForBinSize Enable the inverse bin volume correction factor.
1468/// \param[in] cdfBoundaries Enable the special boundary condition for a cdf:
1469/// underflow bins are assumed to have weight zero and
1470/// overflow bins have weight one. Otherwise, the
1471/// histogram is mirrored at the boundaries for the
1472/// interpolation.
1473
1474double RooDataHist::weightInterpolated(const RooArgSet& bin, int intOrder, bool correctForBinSize, bool cdfBoundaries) {
1475 VarInfo const& varInfo = getVarInfo();
1476
1477 const auto centralIdx = calcTreeIndex(bin, true);
1478
1479 double wInt{0} ;
1480 if (varInfo.nRealVars == 1) {
1481
1482 // buffer needs to be 2 x (interpolation order + 1), with the factor 2 for x and y.
1483 _interpolationBuffer.resize(2 * intOrder + 2);
1484
1485 // 1-dimensional interpolation
1486 auto const& realX = static_cast<RooRealVar const&>(*bin[varInfo.realVarIdx1]);
1487 wInt = interpolateDim(varInfo.realVarIdx1, realX.getVal(), centralIdx, intOrder, correctForBinSize, cdfBoundaries) ;
1488
1489 } else if (varInfo.nRealVars == 2) {
1490
1491 // buffer needs to be 2 x 2 x (interpolation order + 1), with one factor 2
1492 // for x and y, and the other for the number of dimensions.
1493 _interpolationBuffer.resize(4 * intOrder + 4);
1494
1495 // 2-dimensional interpolation
1496 auto const& realX = static_cast<RooRealVar const&>(*bin[varInfo.realVarIdx1]);
1497 auto const& realY = static_cast<RooRealVar const&>(*bin[varInfo.realVarIdx2]);
1498 double xval = realX.getVal() ;
1499 double yval = realY.getVal() ;
1500
1501 RooAbsBinning const& binningY = realY.getBinning();
1502
1503 int ybinC = binningY.binNumber(yval) ;
1504 int ybinLo = ybinC-intOrder/2 - ((yval<binningY.binCenter(ybinC))?1:0) ;
1505 int ybinM = binningY.numBins() ;
1506
1507 auto idxMultY = _idxMult[varInfo.realVarIdx2];
1508 auto offsetIdx = centralIdx - idxMultY * ybinC;
1509
1510 // Use a class-member buffer to avoid repeated heap allocations.
1511 double * yarr = _interpolationBuffer.data() + 2 * intOrder + 2; // add offset to skip part reserved for other dim
1512 double * xarr = yarr + intOrder + 1;
1513 for (int i=ybinLo ; i<=intOrder+ybinLo ; i++) {
1514 int ibin ;
1515 if (i>=0 && i<ybinM) {
1516 // In range
1517 ibin = i ;
1518 xarr[i-ybinLo] = binningY.binCenter(ibin) ;
1519 } else if (i>=ybinM) {
1520 // Overflow: mirror
1521 ibin = 2*ybinM-i-1 ;
1522 xarr[i-ybinLo] = 2*binningY.highBound()-binningY.binCenter(ibin) ;
1523 } else {
1524 // Underflow: mirror
1525 ibin = -i -1;
1526 xarr[i-ybinLo] = 2*binningY.lowBound()-binningY.binCenter(ibin) ;
1527 }
1528 auto centralIdxX = offsetIdx + idxMultY * ibin;
1529 yarr[i-ybinLo] = interpolateDim(varInfo.realVarIdx1,xval,centralIdxX,intOrder,correctForBinSize,false) ;
1530 }
1531
1532 if (gDebug>7) {
1533 cout << "RooDataHist interpolating data is" << endl ;
1534 cout << "xarr = " ;
1535 for (int q=0; q<=intOrder ; q++) cout << xarr[q] << " " ;
1536 cout << " yarr = " ;
1537 for (int q=0; q<=intOrder ; q++) cout << yarr[q] << " " ;
1538 cout << endl ;
1539 }
1540 wInt = RooMath::interpolate(xarr,yarr,intOrder+1,yval) ;
1541
1542 } else {
1543
1544 // Higher dimensional scenarios not yet implemented
1545 coutE(InputArguments) << "RooDataHist::weight(" << GetName() << ") interpolation in "
1546 << varInfo.nRealVars << " dimensions not yet implemented" << endl ;
1547 return weightFast(bin,0,correctForBinSize,cdfBoundaries) ;
1548
1549 }
1550
1551 return wInt ;
1552}
1553
1554
1556 if (!_errLo || !_errHi) {
1557 initArray(_errLo, _arrSize, -1.);
1558 initArray(_errHi, _arrSize, -1.);
1560 }
1561}
1562
1563
1564////////////////////////////////////////////////////////////////////////////////
1565/// Return the asymmetric errors on the current weight.
1566/// \see weightError(ErrorType) const for symmetric error.
1567/// \param[out] lo Low error.
1568/// \param[out] hi High error.
1569/// \param[in] etype Type of error to compute. May throw if not supported.
1570/// Supported errors are
1571/// - `Poisson` Default. Asymmetric Poisson errors (68% CL).
1572/// - `SumW2` The square root of the sum of weights. (Symmetric).
1573/// - `None` Return zero.
1574void RooDataHist::weightError(double& lo, double& hi, ErrorType etype) const
1575{
1576 checkInit() ;
1577
1578 switch (etype) {
1579
1580 case Auto:
1581 throw std::invalid_argument("RooDataHist::weightError(" + std::string(GetName()) + ") error type Auto not allowed here");
1582 break ;
1583
1584 case Expected:
1585 throw std::invalid_argument("RooDataHist::weightError(" + std::string(GetName()) + ") error type Expected not allowed here");
1586 break ;
1587
1588 case Poisson:
1589 if (_errLo && _errLo[_curIndex] >= 0.0) {
1590 // Weight is preset or precalculated
1591 lo = _errLo[_curIndex];
1592 hi = _errHi[_curIndex];
1593 return ;
1594 }
1595
1596 // We didn't track asymmetric errors so far, so now we need to allocate
1598
1599 // Calculate poisson errors
1600 double ym;
1601 double yp;
1603 _errLo[_curIndex] = weight()-ym;
1604 _errHi[_curIndex] = yp-weight();
1605 lo = _errLo[_curIndex];
1606 hi = _errHi[_curIndex];
1607 return ;
1608
1609 case SumW2:
1610 lo = std::sqrt(weightSquared(_curIndex));
1611 hi = lo;
1612 return ;
1613
1614 case None:
1615 lo = 0 ;
1616 hi = 0 ;
1617 return ;
1618 }
1619}
1620
1621
1622// wve adjust for variable bin sizes
1623
1624////////////////////////////////////////////////////////////////////////////////
1625/// Perform boundary safe 'intOrder'-th interpolation of weights in dimension 'dim'
1626/// at current value 'xval'
1627
1628/// \param[in] iDim Index of the histogram dimension along which to interpolate.
1629/// \param[in] xval Value of histogram variable at dimension `iDim` for which
1630/// we want to interpolate the histogram weight.
1631/// \param[in] centralIdx Index of the bin that the point at which we
1632/// interpolate the histogram weight falls into
1633/// (can be obtained with `RooDataHist::calcTreeIndex`).
1634/// \param[in] intOrder Interpolation order, i.e. how many neighbouring bins are
1635/// used for the interpolation.
1636/// \param[in] correctForBinSize Enable the inverse bin volume correction factor.
1637/// \param[in] cdfBoundaries Enable the special boundary condition for a cdf:
1638/// underflow bins are assumed to have weight zero and
1639/// overflow bins have weight one. Otherwise, the
1640/// histogram is mirrored at the boundaries for the
1641/// interpolation.
1642double RooDataHist::interpolateDim(int iDim, double xval, size_t centralIdx, int intOrder, bool correctForBinSize, bool cdfBoundaries)
1643{
1644 auto const& binning = static_cast<RooRealVar&>(*_vars[iDim]).getBinning();
1645
1646 // Fill workspace arrays spanning interpolation area
1647 int fbinC = binning.binNumber(xval) ;
1648 int fbinLo = fbinC-intOrder/2 - ((xval<binning.binCenter(fbinC))?1:0) ;
1649 int fbinM = binning.numBins() ;
1650
1651 auto idxMult = _idxMult[iDim];
1652 auto offsetIdx = centralIdx - idxMult * fbinC;
1653
1654 // Use a class-member buffer to avoid repeated heap allocations.
1655 double * yarr = _interpolationBuffer.data();
1656 double * xarr = yarr + intOrder + 1;
1657
1658 for (int i=fbinLo ; i<=intOrder+fbinLo ; i++) {
1659 int ibin ;
1660 if (i>=0 && i<fbinM) {
1661 // In range
1662 ibin = i ;
1663 xarr[i-fbinLo] = binning.binCenter(ibin) ;
1664 auto idx = offsetIdx + idxMult * ibin;
1665 yarr[i - fbinLo] = _wgt[idx];
1666 if (correctForBinSize) yarr[i-fbinLo] /= _binv[idx] ;
1667 } else if (i>=fbinM) {
1668 // Overflow: mirror
1669 ibin = 2*fbinM-i-1 ;
1670 if (cdfBoundaries) {
1671 xarr[i-fbinLo] = binning.highBound()+1e-10*(i-fbinM+1) ;
1672 yarr[i-fbinLo] = 1.0 ;
1673 } else {
1674 auto idx = offsetIdx + idxMult * ibin;
1675 xarr[i-fbinLo] = 2*binning.highBound()-binning.binCenter(ibin) ;
1676 yarr[i - fbinLo] = _wgt[idx];
1677 if (correctForBinSize)
1678 yarr[i - fbinLo] /= _binv[idx];
1679 }
1680 } else {
1681 // Underflow: mirror
1682 ibin = -i - 1 ;
1683 if (cdfBoundaries) {
1684 xarr[i-fbinLo] = binning.lowBound()-ibin*(1e-10) ;
1685 yarr[i-fbinLo] = 0.0 ;
1686 } else {
1687 auto idx = offsetIdx + idxMult * ibin;
1688 xarr[i-fbinLo] = 2*binning.lowBound()-binning.binCenter(ibin) ;
1689 yarr[i - fbinLo] = _wgt[idx];
1690 if (correctForBinSize)
1691 yarr[i - fbinLo] /= _binv[idx];
1692 }
1693 }
1694 }
1695 return RooMath::interpolate(xarr,yarr,intOrder+1,xval) ;
1696}
1697
1698
1699
1700
1701////////////////////////////////////////////////////////////////////////////////
1702/// Increment the bin content of the bin enclosing the given coordinates.
1703///
1704/// \param[in] row Coordinates of the bin.
1705/// \param[in] wgt Increment by this weight.
1706/// \param[in] sumw2 Optionally, track the sum of squared weights. If a value > 0 or
1707/// a weight != 1. is passed for the first time, a vector for the squared weights will be allocated.
1708void RooDataHist::add(const RooArgSet& row, double wgt, double sumw2)
1709{
1710 checkInit() ;
1711
1712 if ((sumw2 > 0. || wgt != 1.) && !_sumw2) {
1713 // Receiving a weighted entry. SumW2 != sumw from now on.
1714 _sumw2 = new double[_arrSize];
1715 std::copy(_wgt, _wgt+_arrSize, _sumw2);
1716
1718 }
1719
1720 const auto idx = calcTreeIndex(row, false);
1721
1722 _wgt[idx] += wgt ;
1723 if (_sumw2) _sumw2[idx] += (sumw2 > 0 ? sumw2 : wgt*wgt);
1724
1725 _cache_sum_valid = false;
1726}
1727
1728
1729
1730////////////////////////////////////////////////////////////////////////////////
1731/// Set a bin content.
1732/// \param[in] row Coordinates of the bin to be set.
1733/// \param[in] wgt New bin content.
1734/// \param[in] wgtErrLo Low error of the bin content.
1735/// \param[in] wgtErrHi High error of the bin content.
1736void RooDataHist::set(const RooArgSet& row, double wgt, double wgtErrLo, double wgtErrHi)
1737{
1738 checkInit() ;
1739
1741
1742 const auto idx = calcTreeIndex(row, false);
1743
1744 _wgt[idx] = wgt ;
1745 _errLo[idx] = wgtErrLo ;
1746 _errHi[idx] = wgtErrHi ;
1747
1748 _cache_sum_valid = false;
1749}
1750
1751
1752
1753////////////////////////////////////////////////////////////////////////////////
1754/// Set bin content of bin that was last loaded with get(std::size_t).
1755/// \param[in] binNumber Optional bin number to set. If empty, currently active bin is set.
1756/// \param[in] wgt New bin content.
1757/// \param[in] wgtErr Error of the new bin content. If the weight need not have an error, use 0. or a negative number.
1758void RooDataHist::set(std::size_t binNumber, double wgt, double wgtErr) {
1759 checkInit() ;
1760
1761 if (wgtErr > 0. && !_sumw2) {
1762 // Receiving a weighted entry. Need to track sumw2 from now on:
1763 cloneArray(_sumw2, _wgt, _arrSize);
1764
1766 }
1767
1768 _wgt[binNumber] = wgt ;
1769 if (_errLo) _errLo[binNumber] = wgtErr;
1770 if (_errHi) _errHi[binNumber] = wgtErr;
1771 if (_sumw2) _sumw2[binNumber] = wgtErr*wgtErr;
1772
1773 _cache_sum_valid = false ;
1774}
1775
1776
1777////////////////////////////////////////////////////////////////////////////////
1778/// Set bin content of bin that was last loaded with get(std::size_t).
1779/// \param[in] wgt New bin content.
1780/// \param[in] wgtErr Optional error of the bin content.
1781void RooDataHist::set(double wgt, double wgtErr) {
1782 if (_curIndex == std::numeric_limits<std::size_t>::max()) {
1783 _curIndex = calcTreeIndex(_vars, true) ;
1784 }
1785
1786 set(_curIndex, wgt, wgtErr);
1787}
1788
1789
1790////////////////////////////////////////////////////////////////////////////////
1791/// Set a bin content.
1792/// \param[in] row Coordinates to compute the bin from.
1793/// \param[in] wgt New bin content.
1794/// \param[in] wgtErr Optional error of the bin content.
1795void RooDataHist::set(const RooArgSet& row, double wgt, double wgtErr) {
1796 set(calcTreeIndex(row, false), wgt, wgtErr);
1797}
1798
1799
1800
1801////////////////////////////////////////////////////////////////////////////////
1802/// Add all data points contained in 'dset' to this data set with given weight.
1803/// Optional cut string expression selects the data points to be added and can
1804/// reference any variable contained in this data set
1805
1806void RooDataHist::add(const RooAbsData& dset, const char* cut, double wgt)
1807{
1808 RooFormulaVar cutVar("select",cut,*dset.get()) ;
1809 add(dset,&cutVar,wgt) ;
1810}
1811
1812
1813
1814////////////////////////////////////////////////////////////////////////////////
1815/// Add all data points contained in 'dset' to this data set with given weight.
1816/// Optional RooFormulaVar pointer selects the data points to be added.
1817
1818void RooDataHist::add(const RooAbsData& dset, const RooFormulaVar* cutVar, double wgt)
1819{
1820 checkInit() ;
1821
1822 RooFormulaVar* cloneVar = nullptr;
1823 std::unique_ptr<RooArgSet> tmp;
1824 if (cutVar) {
1825 // Deep clone cutVar and attach clone to this dataset
1826 tmp = std::make_unique<RooArgSet>();
1827 if(RooArgSet(*cutVar).snapshot(*tmp)) {
1828 coutE(DataHandling) << "RooDataHist::add(" << GetName() << ") Couldn't deep-clone cut variable, abort," << endl ;
1829 return ;
1830 }
1831
1832 cloneVar = static_cast<RooFormulaVar*>(tmp->find(*cutVar)) ;
1833 cloneVar->attachDataSet(dset) ;
1834 }
1835
1836
1837 Int_t i ;
1838 for (i=0 ; i<dset.numEntries() ; i++) {
1839 const RooArgSet* row = dset.get(i) ;
1840 if (!cloneVar || cloneVar->getVal()) {
1841 add(*row,wgt*dset.weight(), wgt*wgt*dset.weightSquared()) ;
1842 }
1843 }
1844
1845 _cache_sum_valid = false ;
1846}
1847
1848
1849
1850////////////////////////////////////////////////////////////////////////////////
1851/// Return the sum of the weights of all bins in the histogram.
1852///
1853/// \param[in] correctForBinSize Multiply the sum of weights in each bin
1854/// with the N-dimensional bin volume, making the return value
1855/// the integral over the function represented by this histogram.
1856/// \param[in] inverseBinCor Divide by the N-dimensional bin volume.
1857double RooDataHist::sum(bool correctForBinSize, bool inverseBinCor) const
1858{
1859 checkInit() ;
1860
1861 // Check if result was cached
1862 const CacheSumState_t cache_code = !correctForBinSize ? kNoBinCorrection : (inverseBinCor ? kInverseBinCorr : kCorrectForBinSize);
1863 if (_cache_sum_valid == static_cast<Int_t>(cache_code)) {
1864 return _cache_sum ;
1865 }
1866
1868 for (Int_t i=0; i < _arrSize; i++) {
1869 const double theBinVolume = correctForBinSize ? (inverseBinCor ? 1/_binv[i] : _binv[i]) : 1.0 ;
1870 kahanSum += _wgt[i] * theBinVolume;
1871 }
1872
1873 // Store result in cache
1874 _cache_sum_valid = cache_code;
1875 _cache_sum = kahanSum.Sum();
1876
1877 return kahanSum.Sum();
1878}
1879
1880
1881
1882////////////////////////////////////////////////////////////////////////////////
1883/// Return the sum of the weights of a multi-dimensional slice of the histogram
1884/// by summing only over the dimensions specified in sumSet.
1885///
1886/// The coordinates of all other dimensions are fixed to those given in sliceSet
1887///
1888/// If correctForBinSize is specified, the sum of weights
1889/// is multiplied by the M-dimensional bin volume, (M = N(sumSet)),
1890/// making the return value the integral over the function
1891/// represented by this histogram
1892
1893double RooDataHist::sum(const RooArgSet& sumSet, const RooArgSet& sliceSet, bool correctForBinSize, bool inverseBinCor)
1894{
1895 checkInit() ;
1896
1897 RooArgSet varSave ;
1898 varSave.addClone(_vars) ;
1899
1900 RooArgSet sliceOnlySet(sliceSet);
1901 sliceOnlySet.remove(sumSet,true,true) ;
1902
1903 _vars.assign(sliceOnlySet);
1904 std::vector<double> const * pbinv = nullptr;
1905
1906 if(correctForBinSize && inverseBinCor) {
1907 pbinv = &calculatePartialBinVolume(sliceOnlySet);
1908 } else if(correctForBinSize && !inverseBinCor) {
1909 pbinv = &calculatePartialBinVolume(sumSet);
1910 }
1911
1912 // Calculate mask and reference plot bins for non-iterating variables
1913 std::vector<bool> mask(_vars.size());
1914 std::vector<int> refBin(_vars.size());
1915
1916 for (unsigned int i = 0; i < _vars.size(); ++i) {
1917 const RooAbsArg* arg = _vars[i];
1918 const RooAbsLValue* argLv = _lvvars[i]; // Same as above, but cross-cast
1919
1920 if (sumSet.find(*arg)) {
1921 mask[i] = false ;
1922 } else {
1923 mask[i] = true ;
1924 refBin[i] = argLv->getBin();
1925 }
1926 }
1927
1928 // Loop over entire data set, skipping masked entries
1930 for (Int_t ibin=0; ibin < _arrSize; ++ibin) {
1931
1932 std::size_t tmpibin = ibin;
1933 bool skip(false) ;
1934
1935 // Check if this bin belongs in selected slice
1936 for (unsigned int ivar = 0; !skip && ivar < _vars.size(); ++ivar) {
1937 const Int_t idx = tmpibin / _idxMult[ivar] ;
1938 tmpibin -= idx*_idxMult[ivar] ;
1939 if (mask[ivar] && idx!=refBin[ivar])
1940 skip = true ;
1941 }
1942
1943 if (!skip) {
1944 const double theBinVolume = correctForBinSize ? (inverseBinCor ? 1/(*pbinv)[ibin] : (*pbinv)[ibin] ) : 1.0 ;
1945 total += _wgt[ibin] * theBinVolume;
1946 }
1947 }
1948
1949 _vars.assign(varSave) ;
1950
1951 return total.Sum();
1952}
1953
1954////////////////////////////////////////////////////////////////////////////////
1955/// Return the sum of the weights of a multi-dimensional slice of the histogram
1956/// by summing only over the dimensions specified in sumSet.
1957///
1958/// The coordinates of all other dimensions are fixed to those given in sliceSet
1959///
1960/// If correctForBinSize is specified, the sum of weights
1961/// is multiplied by the M-dimensional bin volume, (M = N(sumSet)),
1962/// or the fraction of it that falls inside the range rangeName,
1963/// making the return value the integral over the function
1964/// represented by this histogram.
1965///
1966/// If correctForBinSize is not specified, the weights are multiplied by the
1967/// fraction of the bin volume that falls inside the range, i.e. a factor of
1968/// binVolumeInRange/totalBinVolume.
1969
1970double RooDataHist::sum(const RooArgSet& sumSet, const RooArgSet& sliceSet,
1971 bool correctForBinSize, bool inverseBinCor,
1972 const std::map<const RooAbsArg*, std::pair<double, double> >& ranges,
1973 std::function<double(int)> getBinScale)
1974{
1975 checkInit();
1977 RooArgSet varSave;
1978 varSave.addClone(_vars);
1979 {
1980 RooArgSet sliceOnlySet(sliceSet);
1981 sliceOnlySet.remove(sumSet, true, true);
1982 _vars.assign(sliceOnlySet);
1983 }
1984
1985 // Calculate mask and reference plot bins for non-iterating variables,
1986 // and get ranges for iterating variables
1987 std::vector<bool> mask(_vars.size());
1988 std::vector<int> refBin(_vars.size());
1989 std::vector<double> rangeLo(_vars.size(), -std::numeric_limits<double>::infinity());
1990 std::vector<double> rangeHi(_vars.size(), +std::numeric_limits<double>::infinity());
1991
1992 for (std::size_t i = 0; i < _vars.size(); ++i) {
1993 const RooAbsArg* arg = _vars[i];
1994 const RooAbsLValue* argLV = _lvvars[i]; // Same object as above, but cross cast
1995
1996 RooAbsArg* sumsetv = sumSet.find(*arg);
1997 RooAbsArg* slicesetv = sliceSet.find(*arg);
1998 mask[i] = !sumsetv;
1999 if (mask[i]) {
2000 assert(argLV);
2001 refBin[i] = argLV->getBin();
2002 }
2003
2004 auto it = ranges.find(sumsetv ? sumsetv : slicesetv);
2005 if (ranges.end() != it) {
2006 rangeLo[i] = it->second.first;
2007 rangeHi[i] = it->second.second;
2008 }
2009 }
2010
2011 // Loop over entire data set, skipping masked entries
2013 for (Int_t ibin = 0; ibin < _arrSize; ++ibin) {
2014 // Check if this bin belongs in selected slice
2015 bool skip{false};
2016 for (int ivar = 0, tmp = ibin; !skip && ivar < int(_vars.size()); ++ivar) {
2017 const Int_t idx = tmp / _idxMult[ivar];
2018 tmp -= idx*_idxMult[ivar];
2019 if (mask[ivar] && idx!=refBin[ivar]) skip = true;
2020 }
2021
2022 if (skip) continue;
2023
2024 // Work out bin volume
2025 // It's not necessary to figure out the bin volume for the slice-only set explicitly here.
2026 // We need to loop over the sumSet anyway to get the partial bin containment correction,
2027 // so we can get the slice-only set volume later by dividing _binv[ibin] / binVolumeSumSetFull.
2028 double binVolumeSumSetFull = 1.;
2029 double binVolumeSumSetInRange = 1.;
2030 for (Int_t ivar = 0, tmp = ibin; ivar < (int)_vars.size(); ++ivar) {
2031 const Int_t idx = tmp / _idxMult[ivar];
2032 tmp -= idx*_idxMult[ivar];
2033
2034 // If the current variable is not in the sumSet, it should not be considered for the bin volume
2035 const auto arg = _vars[ivar];
2036 if (!sumSet.find(*arg)) {
2037 continue;
2038 }
2039
2040 if (_binbounds[ivar].empty()) continue;
2041 const double binLo = _binbounds[ivar][2 * idx];
2042 const double binHi = _binbounds[ivar][2 * idx + 1];
2043 if (binHi < rangeLo[ivar] || binLo > rangeHi[ivar]) {
2044 // bin is outside of allowed range - effective bin volume is zero
2045 binVolumeSumSetInRange = 0.;
2046 break;
2047 }
2048
2049 binVolumeSumSetFull *= binHi - binLo;
2050 binVolumeSumSetInRange *= std::min(rangeHi[ivar], binHi) - std::max(rangeLo[ivar], binLo);
2051 }
2052 const double corrPartial = binVolumeSumSetInRange / binVolumeSumSetFull;
2053 if (0. == corrPartial) continue;
2054 const double corr = correctForBinSize ? (inverseBinCor ? binVolumeSumSetFull / _binv[ibin] : binVolumeSumSetFull ) : 1.0;
2055 total += getBinScale(ibin)*(_wgt[ibin] * corr * corrPartial);
2056 }
2057
2058 _vars.assign(varSave);
2059
2060 return total.Sum();
2061}
2062
2063
2064
2065////////////////////////////////////////////////////////////////////////////////
2066/// Fill the transient cache with partial bin volumes with up-to-date
2067/// values for the partial volume specified by observables 'dimSet'
2068
2069const std::vector<double>& RooDataHist::calculatePartialBinVolume(const RooArgSet& dimSet) const
2070{
2071 // The code bitset has all bits set to one whose position corresponds to arguments in dimSet.
2072 // It is used as the key for the bin volume caching hash map.
2073 int code{0};
2074 {
2075 int i{0} ;
2076 for (auto const& v : _vars) {
2077 code += ((dimSet.find(*v) ? 1 : 0) << i) ;
2078 ++i;
2079 }
2080 }
2081
2082 auto& pbinv = _pbinvCache[code];
2083 if(!pbinv.empty()) {
2084 return pbinv;
2085 }
2086 pbinv.resize(_arrSize);
2087
2088 // Calculate plot bins of components from master index
2089 std::vector<bool> selDim(_vars.size());
2090 for (std::size_t i = 0; i < selDim.size(); ++i) {
2091 selDim[i] = (code >> i) & 1 ;
2092 }
2093
2094 // Recalculate partial bin volume cache
2095 for (Int_t ibin=0; ibin < _arrSize ;ibin++) {
2096 Int_t idx(0);
2097 Int_t tmp(ibin);
2098 double theBinVolume(1) ;
2099 for (unsigned int j=0; j < _lvvars.size(); ++j) {
2100 const RooAbsLValue* arg = _lvvars[j];
2101 assert(arg);
2102
2103 idx = tmp / _idxMult[j];
2104 tmp -= idx*_idxMult[j];
2105 if (selDim[j]) {
2106 theBinVolume *= arg->getBinWidth(idx) ;
2107 }
2108 }
2109 pbinv[ibin] = theBinVolume ;
2110 }
2111
2112 return pbinv;
2113}
2114
2115
2116////////////////////////////////////////////////////////////////////////////////
2117/// Sum the weights of all bins.
2120}
2121
2122
2123
2124////////////////////////////////////////////////////////////////////////////////
2125/// Return the sum of weights in all entries matching cutSpec (if specified)
2126/// and in named range cutRange (if specified)
2127/// Return the
2128
2129double RooDataHist::sumEntries(const char* cutSpec, const char* cutRange) const
2130{
2131 checkInit() ;
2132
2133 if (cutSpec==nullptr && cutRange==nullptr) {
2134 return sumEntries();
2135 } else {
2136
2137 // Setup RooFormulaVar for cutSpec if it is present
2138 std::unique_ptr<RooFormula> select;
2139 if (cutSpec) {
2140 select = std::make_unique<RooFormula>("select",cutSpec,*get());
2141 }
2142
2143 // Otherwise sum the weights in the event
2144 ROOT::Math::KahanSum<> kahanSum;
2145 for (Int_t i=0; i < _arrSize; i++) {
2146 get(i) ;
2147 if ((select && select->eval() == 0.) || (cutRange && !_vars.allInRange(cutRange)))
2148 continue;
2149
2150 kahanSum += weight(i);
2151 }
2152
2153 return kahanSum.Sum();
2154 }
2155}
2156
2157
2158
2159////////////////////////////////////////////////////////////////////////////////
2160/// Reset all bin weights to zero
2161
2163{
2164 // WVE DO NOT CALL RooTreeData::reset() for binned
2165 // datasets as this will delete the bin definitions
2166
2167 std::fill(_wgt, _wgt + _arrSize, 0.);
2168 delete[] _errLo; _errLo = nullptr;
2169 delete[] _errHi; _errHi = nullptr;
2170 delete[] _sumw2; _sumw2 = nullptr;
2171
2173
2174 _cache_sum_valid = false;
2175}
2176
2177
2178
2179////////////////////////////////////////////////////////////////////////////////
2180/// Load bin `binNumber`, and return an argset with the coordinates of the bin centre.
2181/// \note The argset is owned by this data hist, and this function has a side effect, because
2182/// it alters the currently active bin.
2183const RooArgSet* RooDataHist::get(Int_t binNumber) const
2184{
2185 checkInit() ;
2186 _curIndex = binNumber;
2187
2188 return RooAbsData::get(_curIndex);
2189}
2190
2191
2192
2193////////////////////////////////////////////////////////////////////////////////
2194/// Return a RooArgSet with whose coordinates denote the bin centre of the bin
2195/// enclosing the point in `coord`.
2196/// \note The argset is owned by this data hist, and this function has a side effect, because
2197/// it alters the currently active bin.
2198const RooArgSet* RooDataHist::get(const RooArgSet& coord) const {
2199 return get(calcTreeIndex(coord, false));
2200}
2201
2202
2203
2204////////////////////////////////////////////////////////////////////////////////
2205/// Return the volume of the bin enclosing coordinates 'coord'.
2206double RooDataHist::binVolume(const RooArgSet& coord) const {
2207 checkInit() ;
2208 return _binv[calcTreeIndex(coord, false)] ;
2209}
2210
2211
2212////////////////////////////////////////////////////////////////////////////////
2213/// Create an iterator over all bins in a slice defined by the subset of observables
2214/// listed in sliceArg. The position of the slice is given by otherArgs
2215
2217{
2218 // Update to current position
2219 _vars.assign(otherArgs) ;
2220 _curIndex = calcTreeIndex(_vars, true);
2221
2222 RooAbsArg* intArg = _vars.find(sliceArg) ;
2223 if (!intArg) {
2224 coutE(InputArguments) << "RooDataHist::sliceIterator() variable " << sliceArg.GetName() << " is not part of this RooDataHist" << endl ;
2225 return nullptr ;
2226 }
2227 return new RooDataHistSliceIter(*this,*intArg) ;
2228}
2229
2230
2231////////////////////////////////////////////////////////////////////////////////
2232/// Change the name of the RooDataHist
2233
2235{
2236 if (_dir) _dir->GetList()->Remove(this);
2237 // We need to use the function from RooAbsData, because it already overrides TNamed::SetName
2239 if (_dir) _dir->GetList()->Add(this);
2240}
2241
2242
2243////////////////////////////////////////////////////////////////////////////////
2244/// Change the title of this RooDataHist
2245
2246void RooDataHist::SetNameTitle(const char *name, const char* title)
2247{
2248 SetName(name);
2249 SetTitle(title);
2250}
2251
2252
2253////////////////////////////////////////////////////////////////////////////////
2254/// Print value of the dataset, i.e. the sum of weights contained in the dataset
2255
2256void RooDataHist::printValue(ostream& os) const
2257{
2258 os << numEntries() << " bins (" << sumEntries() << " weights)" ;
2259}
2260
2261
2262
2263
2264////////////////////////////////////////////////////////////////////////////////
2265/// Print argument of dataset, i.e. the observable names
2266
2267void RooDataHist::printArgs(ostream& os) const
2268{
2269 os << "[" ;
2270 bool first(true) ;
2271 for (const auto arg : _vars) {
2272 if (first) {
2273 first=false ;
2274 } else {
2275 os << "," ;
2276 }
2277 os << arg->GetName() ;
2278 }
2279 os << "]" ;
2280}
2281
2282
2283
2284////////////////////////////////////////////////////////////////////////////////
2285/// Returns true if dataset contains entries with a non-integer weight.
2286
2288{
2289 for (Int_t i=0; i < _arrSize; ++i) {
2290 const double wgt = _wgt[i];
2291 double intpart;
2292 if (std::abs(std::modf(wgt, &intpart)) > 1.E-10)
2293 return true;
2294 }
2295
2296 return false;
2297}
2298
2299
2300////////////////////////////////////////////////////////////////////////////////
2301/// Print the details on the dataset contents
2302
2303void RooDataHist::printMultiline(ostream& os, Int_t content, bool verbose, TString indent) const
2304{
2305 RooAbsData::printMultiline(os,content,verbose,indent) ;
2306
2307 os << indent << "Binned Dataset " << GetName() << " (" << GetTitle() << ")" << endl ;
2308 os << indent << " Contains " << numEntries() << " bins with a total weight of " << sumEntries() << endl;
2309
2310 if (!verbose) {
2311 os << indent << " Observables " << _vars << endl ;
2312 } else {
2313 os << indent << " Observables: " ;
2315 }
2316
2317 if(verbose) {
2318 if (!_cachedVars.empty()) {
2319 os << indent << " Caches " << _cachedVars << endl ;
2320 }
2321 }
2322}
2323
2324void RooDataHist::printDataHistogram(ostream& os, RooRealVar* obs) const
2325{
2326 for(Int_t i=0; i<obs->getBins(); ++i){
2327 this->get(i);
2328 obs->setBin(i);
2329 os << this->weight() << " +/- " << this->weightSquared() << endl;
2330 }
2331}
2332
2333
2334////////////////////////////////////////////////////////////////////////////////
2335/// Stream an object of class RooDataHist.
2337 if (R__b.IsReading()) {
2338
2339 UInt_t R__s;
2340 UInt_t R__c;
2341 Version_t R__v = R__b.ReadVersion(&R__s, &R__c);
2342
2343 if (R__v > 2) {
2344 R__b.ReadClassBuffer(RooDataHist::Class(),this,R__v,R__s,R__c);
2345 R__b.CheckByteCount(R__s, R__c, RooDataHist::IsA());
2346 initialize(nullptr, false);
2347 } else {
2348
2349 // Legacy dataset conversion happens here. Legacy RooDataHist inherits from RooTreeData
2350 // which in turn inherits from RooAbsData. Manually stream RooTreeData contents on
2351 // file here and convert it into a RooTreeDataStore which is installed in the
2352 // new-style RooAbsData base class
2353
2354 // --- This is the contents of the streamer code of RooTreeData version 2 ---
2355 UInt_t R__s1;
2356 UInt_t R__c1;
2357 Version_t R__v1 = R__b.ReadVersion(&R__s1, &R__c1); if (R__v1) { }
2358
2360 TTree* X_tree(nullptr) ; R__b >> X_tree;
2361 RooArgSet X_truth ; X_truth.Streamer(R__b);
2362 TString X_blindString ; X_blindString.Streamer(R__b);
2363 R__b.CheckByteCount(R__s1, R__c1, TClass::GetClass("RooTreeData"));
2364 // --- End of RooTreeData-v1 streamer
2365
2366 // Construct RooTreeDataStore from X_tree and complete initialization of new-style RooAbsData
2367 _dstore = std::make_unique<RooTreeDataStore>(X_tree,_vars);
2368 _dstore->SetName(GetName()) ;
2369 _dstore->SetTitle(GetTitle()) ;
2370 _dstore->checkInit() ;
2371
2373 R__b >> _arrSize;
2374 delete [] _wgt;
2375 _wgt = new double[_arrSize];
2377 delete [] _errLo;
2378 _errLo = new double[_arrSize];
2380 delete [] _errHi;
2381 _errHi = new double[_arrSize];
2383 delete [] _sumw2;
2384 _sumw2 = new double[_arrSize];
2386 delete [] _binv;
2387 _binv = new double[_arrSize];
2388 RooArgSet tmpSet;
2389 tmpSet.Streamer(R__b);
2390 double tmp;
2391 R__b >> tmp; //_curWeight;
2392 R__b >> tmp; //_curWgtErrLo;
2393 R__b >> tmp; //_curWgtErrHi;
2394 R__b >> tmp; //_curSumW2;
2395 R__b >> tmp; //_curVolume;
2396 R__b >> _curIndex;
2397 R__b.CheckByteCount(R__s, R__c, RooDataHist::IsA());
2398 }
2399
2400 } else {
2401
2403 }
2404}
2405
2406
2407////////////////////////////////////////////////////////////////////////////////
2408/// Return event weights of all events in range [first, first+len).
2409/// If cacheValidEntries() has been called, out-of-range events will have a weight of 0.
2410std::span<const double> RooDataHist::getWeightBatch(std::size_t first, std::size_t len, bool sumW2 /*=false*/) const {
2411 return {(sumW2 && _sumw2 ? _sumw2 : _wgt) + first, len};
2412}
2413
2414
2415////////////////////////////////////////////////////////////////////////////////
2416/// Hand over pointers to our weight arrays to the data store implementation.
2418 _dstore->setExternalWeightArray(_wgt, _errLo, _errHi, _sumw2);
2419}
2420
2421
2422////////////////////////////////////////////////////////////////////////////////
2423/// Return reference to VarInfo struct with cached histogram variable
2424/// information that is frequently used for histogram weights retrieval.
2425///
2426/// If the `_varInfo` struct was not initialized yet, it will be initialized in
2427/// this function.
2429
2430 if(_varInfo.initialized) return _varInfo;
2431
2432 auto& info = _varInfo;
2433
2434 {
2435 // count the number of real vars and get their indices
2436 info.nRealVars = 0;
2437 size_t iVar = 0;
2438 for (const auto real : _vars) {
2439 if (dynamic_cast<RooRealVar*>(real)) {
2440 if(info.nRealVars == 0) info.realVarIdx1 = iVar;
2441 if(info.nRealVars == 1) info.realVarIdx2 = iVar;
2442 ++info.nRealVars;
2443 }
2444 ++iVar;
2445 }
2446 }
2447
2448 {
2449 // assert that the variables are either real values or categories
2450 for (unsigned int i=0; i < _vars.size(); ++i) {
2451 if (_lvbins[i].get()) {
2452 assert(dynamic_cast<const RooAbsReal*>(_vars[i]));
2453 } else {
2454 assert(dynamic_cast<const RooAbsCategoryLValue*>(_vars[i]));
2455 }
2456 }
2457 }
2458
2459 info.initialized = true;
2460
2461 return info;
2462}
#define e(i)
Definition RSha256.hxx:103
bool initialize()
#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:377
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 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 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 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:77
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.
bool selectCommon(const RooAbsCollection &refColl, RooAbsCollection &outColl) const
Create a subset of the current collection, consisting only of those elements that are contained as we...
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:55
RooArgSet * snapshot(bool deepCopy=true) const
Use RooAbsCollection::snapshot(), but return as RooArgSet.
Definition RooArgSet.h:178
void Streamer(TBuffer &) override
Stream an object of class TObject.
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:39
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:66
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.
bool isNonPoissonWeighted() const override
Returns true if dataset contains entries with a non-integer weight.
std::string declWeightArrayForCodeSquash(RooAbsArg const *klass, RooFit::Detail::CodeSquashContext &ctx, bool correctForBinSize) const
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()) 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:76
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.
void addToCodeBody(RooAbsArg const *klass, std::string const &in)
Adds the input string to the squashed code body.
std::string getTmpVarName() const
Get a unique variable name to be used in the generated 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:43
RooAbsRealLValue * getPlotVar() const
Definition RooPlot.h:139
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:207
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.
std::string makeValidVarName(std::string const &in)
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()