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