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