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