Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
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 
22 The RooDataHist is a container class to hold N-dimensional binned data. Each bin's central
23 coordinates in N-dimensional space are represented by a RooArgSet containing RooRealVar, RooCategory
24 or RooStringVar objects, thus data can be binned in real and/or discrete dimensions.
25 **/
26 
27 #include "RooFit.h"
28 #include "Riostream.h"
29 
30 #include "TH1.h"
31 #include "TH1.h"
32 #include "TDirectory.h"
33 #include "TMath.h"
34 #include "RooMsgService.h"
35 #include "RooDataHist.h"
36 #include "RooDataHistSliceIter.h"
37 #include "RooAbsLValue.h"
38 #include "RooArgList.h"
39 #include "RooRealVar.h"
40 #include "RooMath.h"
41 #include "RooBinning.h"
42 #include "RooPlot.h"
43 #include "RooHistError.h"
44 #include "RooCategory.h"
45 #include "RooCmdConfig.h"
46 #include "RooLinkedListIter.h"
47 #include "RooTreeDataStore.h"
48 #include "RooVectorDataStore.h"
49 #include "TTree.h"
50 #include "RooTrace.h"
51 #include "RooTreeData.h"
52 #include "RooHelpers.h"
53 
54 using namespace std ;
55 
56 ClassImp(RooDataHist);
57 ;
58 
59 
60 
61 ////////////////////////////////////////////////////////////////////////////////
62 /// Default constructor
63 
64 RooDataHist::RooDataHist() : _pbinvCacheMgr(0,10)
65 {
66  _arrSize = 0 ;
67  _wgt = 0 ;
68  _errLo = 0 ;
69  _errHi = 0 ;
70  _sumw2 = 0 ;
71  _binv = 0 ;
72  _pbinv = 0 ;
73  _curWeight = 0 ;
74  _curIndex = -1 ;
75  _binValid = 0 ;
76  _curSumW2 = 0 ;
77  _curVolume = 1 ;
78  _curWgtErrHi = 0 ;
79  _curWgtErrLo = 0 ;
80  _cache_sum_valid = 0 ;
81  TRACE_CREATE
82 }
83 
84 
85 
86 ////////////////////////////////////////////////////////////////////////////////
87 /// Constructor of an empty data hist from a RooArgSet defining the dimensions
88 /// of the data space. The range and number of bins in each dimensions are taken
89 /// from getMin()getMax(),getBins() of each RooAbsArg representing that
90 /// dimension.
91 ///
92 /// For real dimensions, the fit range and number of bins can be set independently
93 /// of the plot range and number of bins, but it is advisable to keep the
94 /// ratio of the plot bin width and the fit bin width an integer value.
95 /// For category dimensions, the fit ranges always comprises all defined states
96 /// and each state is always has its individual bin
97 ///
98 /// To effective achive binning of real dimensions with variable bins sizes,
99 /// construct a RooThresholdCategory of the real dimension to be binned variably.
100 /// Set the thresholds at the desired bin boundaries, and construct the
101 /// data hist as function of the threshold category instead of the real variable.
102 
103 RooDataHist::RooDataHist(const char *name, const char *title, const RooArgSet& vars, const char* binningName) :
104  RooAbsData(name,title,vars), _wgt(0), _binValid(0), _curWeight(0), _curVolume(1), _pbinv(0), _pbinvCacheMgr(0,10), _cache_sum_valid(0)
105 {
106  // Initialize datastore
107  _dstore = (defaultStorageType==Tree) ? ((RooAbsDataStore*) new RooTreeDataStore(name,title,_vars)) :
108  ((RooAbsDataStore*) new RooVectorDataStore(name,title,_vars)) ;
109 
110  initialize(binningName) ;
111 
112  _dstore->setExternalWeightArray(_wgt,_errLo,_errHi,_sumw2) ;
113 
114  appendToDir(this,kTRUE) ;
115  TRACE_CREATE
116 }
117 
118 
119 
120 ////////////////////////////////////////////////////////////////////////////////
121 /// Constructor of a data hist from an existing data collection (binned or unbinned)
122 /// The RooArgSet 'vars' defines the dimensions of the histogram.
123 /// The range and number of bins in each dimensions are taken
124 /// from getMin()getMax(),getBins() of each RooAbsArg representing that
125 /// dimension.
126 ///
127 /// For real dimensions, the fit range and number of bins can be set independently
128 /// of the plot range and number of bins, but it is advisable to keep the
129 /// ratio of the plot bin width and the fit bin width an integer value.
130 /// For category dimensions, the fit ranges always comprises all defined states
131 /// and each state is always has its individual bin
132 ///
133 /// To effective achive binning of real dimensions with variable bins sizes,
134 /// construct a RooThresholdCategory of the real dimension to be binned variably.
135 /// Set the thresholds at the desired bin boundaries, and construct the
136 /// data hist as function of the threshold category instead of the real variable.
137 ///
138 /// If the constructed data hist has less dimensions that in source data collection,
139 /// all missing dimensions will be projected.
140 
141 RooDataHist::RooDataHist(const char *name, const char *title, const RooArgSet& vars, const RooAbsData& data, Double_t wgt) :
142  RooAbsData(name,title,vars), _wgt(0), _binValid(0), _curWeight(0), _curVolume(1), _pbinv(0), _pbinvCacheMgr(0,10), _cache_sum_valid(0)
143 {
144  // Initialize datastore
145  _dstore = (defaultStorageType==Tree) ? ((RooAbsDataStore*) new RooTreeDataStore(name,title,_vars)) :
146  ((RooAbsDataStore*) new RooVectorDataStore(name,title,_vars)) ;
147 
148  initialize() ;
149  _dstore->setExternalWeightArray(_wgt,_errLo,_errHi,_sumw2) ;
150 
151  add(data,(const RooFormulaVar*)0,wgt) ;
152  appendToDir(this,kTRUE) ;
153  TRACE_CREATE
154 }
155 
156 
157 
158 ////////////////////////////////////////////////////////////////////////////////
159 /// Constructor of a data hist from a map of TH1,TH2 or TH3 that are collated into a x+1 dimensional
160 /// RooDataHist where the added dimension is a category that labels the input source as defined
161 /// in the histMap argument. The state names used in histMap must correspond to predefined states
162 /// 'indexCat'
163 ///
164 /// The RooArgList 'vars' defines the dimensions of the histogram.
165 /// The ranges and number of bins are taken from the input histogram and must be the same in all histograms
166 
167 RooDataHist::RooDataHist(const char *name, const char *title, const RooArgList& vars, RooCategory& indexCat,
168  map<string,TH1*> histMap, Double_t wgt) :
169  RooAbsData(name,title,RooArgSet(vars,&indexCat)),
170  _wgt(0), _binValid(0), _curWeight(0), _curVolume(1), _pbinv(0), _pbinvCacheMgr(0,10), _cache_sum_valid(0)
171 {
172  // Initialize datastore
173  _dstore = (defaultStorageType==Tree) ? ((RooAbsDataStore*) new RooTreeDataStore(name,title,_vars)) :
174  ((RooAbsDataStore*) new RooVectorDataStore(name,title,_vars)) ;
175 
176  importTH1Set(vars, indexCat, histMap, wgt, kFALSE) ;
177 
178  _dstore->setExternalWeightArray(_wgt,_errLo,_errHi,_sumw2) ;
179  TRACE_CREATE
180 }
181 
182 
183 
184 ////////////////////////////////////////////////////////////////////////////////
185 /// Constructor of a data hist from a map of RooDataHists that are collated into a x+1 dimensional
186 /// RooDataHist where the added dimension is a category that labels the input source as defined
187 /// in the histMap argument. The state names used in histMap must correspond to predefined states
188 /// 'indexCat'
189 ///
190 /// The RooArgList 'vars' defines the dimensions of the histogram.
191 /// The ranges and number of bins are taken from the input histogram and must be the same in all histograms
192 
193 RooDataHist::RooDataHist(const char *name, const char *title, const RooArgList& vars, RooCategory& indexCat,
194  map<string,RooDataHist*> dhistMap, Double_t wgt) :
195  RooAbsData(name,title,RooArgSet(vars,&indexCat)),
196  _wgt(0), _binValid(0), _curWeight(0), _curVolume(1), _pbinv(0), _pbinvCacheMgr(0,10), _cache_sum_valid(0)
197 {
198  // Initialize datastore
199  _dstore = (defaultStorageType==Tree) ? ((RooAbsDataStore*) new RooTreeDataStore(name,title,_vars)) :
200  ((RooAbsDataStore*) new RooVectorDataStore(name,title,_vars)) ;
201 
202  importDHistSet(vars, indexCat, dhistMap, wgt) ;
203 
204  _dstore->setExternalWeightArray(_wgt,_errLo,_errHi,_sumw2) ;
205  TRACE_CREATE
206 }
207 
208 
209 
210 ////////////////////////////////////////////////////////////////////////////////
211 /// Constructor of a data hist from an TH1,TH2 or TH3
212 /// The RooArgSet 'vars' defines the dimensions of the histogram. The ranges
213 /// and number of bins are taken from the input histogram, and the corresponding
214 /// values are set accordingly on the arguments in 'vars'
215 
216 RooDataHist::RooDataHist(const char *name, const char *title, const RooArgList& vars, const TH1* hist, Double_t wgt) :
217  RooAbsData(name,title,vars), _wgt(0), _binValid(0), _curWeight(0), _curVolume(1), _pbinv(0), _pbinvCacheMgr(0,10), _cache_sum_valid(0)
218 {
219  // Initialize datastore
220  _dstore = (defaultStorageType==Tree) ? ((RooAbsDataStore*) new RooTreeDataStore(name,title,_vars)) :
221  ((RooAbsDataStore*) new RooVectorDataStore(name,title,_vars)) ;
222 
223  // Check consistency in number of dimensions
224  if (vars.getSize() != hist->GetDimension()) {
225  coutE(InputArguments) << "RooDataHist::ctor(" << GetName() << ") ERROR: dimension of input histogram must match "
226  << "number of dimension variables" << endl ;
227  assert(0) ;
228  }
229 
230  importTH1(vars,*hist,wgt, kFALSE) ;
231 
232  _dstore->setExternalWeightArray(_wgt,_errLo,_errHi,_sumw2) ;
233  TRACE_CREATE
234 }
235 
236 
237 
238 ////////////////////////////////////////////////////////////////////////////////
239 /// Constructor of a binned dataset from a RooArgSet defining the dimensions
240 /// of the data space. The range and number of bins in each dimensions are taken
241 /// from getMin() getMax(),getBins() of each RooAbsArg representing that
242 /// dimension.
243 ///
244 /// <table>
245 /// <tr><th> Optional Argument <th> Effect
246 /// <tr><td> Import(TH1&, Bool_t impDens) <td> Import contents of the given TH1/2/3 into this binned dataset. The
247 /// ranges and binning of the binned dataset are automatically adjusted to
248 /// match those of the imported histogram.
249 ///
250 /// Please note: for TH1& with unequal binning _only_,
251 /// you should decide if you want to import the absolute bin content,
252 /// or the bin content expressed as density. The latter is default and will
253 /// result in the same histogram as the original TH1. For certain types of
254 /// bin contents (containing efficiencies, asymmetries, or ratio is general)
255 /// you should import the absolute value and set impDens to kFALSE
256 ///
257 ///
258 /// <tr><td> Weight(Double_t) <td> Apply given weight factor when importing histograms
259 ///
260 /// <tr><td> Index(RooCategory&) <td> Prepare import of multiple TH1/1/2/3 into a N+1 dimensional RooDataHist
261 /// where the extra discrete dimension labels the source of the imported histogram
262 /// If the index category defines states for which no histogram is be imported
263 /// the corresponding bins will be left empty.
264 ///
265 /// <tr><td> Import(const char*, TH1&) <td> Import a THx to be associated with the given state name of the index category
266 /// specified in Index(). If the given state name is not yet defined in the index
267 /// category it will be added on the fly. The import command can be specified
268 /// multiple times.
269 /// <tr><td> Import(map<string,TH1*>&) <td> As above, but allows specification of many imports in a single operation
270 ///
271 
272 RooDataHist::RooDataHist(const char *name, const char *title, const RooArgList& vars, const RooCmdArg& arg1, const RooCmdArg& arg2, const RooCmdArg& arg3,
273  const RooCmdArg& arg4,const RooCmdArg& arg5,const RooCmdArg& arg6,const RooCmdArg& arg7,const RooCmdArg& arg8) :
274  RooAbsData(name,title,RooArgSet(vars,(RooAbsArg*)RooCmdConfig::decodeObjOnTheFly("RooDataHist::RooDataHist", "IndexCat",0,0,arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8))),
275  _wgt(0), _binValid(0), _curWeight(0), _curVolume(1), _pbinv(0), _pbinvCacheMgr(0,10), _cache_sum_valid(0)
276 {
277  // Initialize datastore
278  _dstore = (defaultStorageType==Tree) ? ((RooAbsDataStore*) new RooTreeDataStore(name,title,_vars)) :
279  ((RooAbsDataStore*) new RooVectorDataStore(name,title,_vars)) ;
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,kTRUE) ; // array
287  pc.defineString("impSliceState","ImportHistoSlice",0,"",kTRUE) ; // array
288  pc.defineObject("impSliceDHist","ImportDataHistSlice",0,0,kTRUE) ; // array
289  pc.defineString("impSliceDState","ImportDataHistSlice",0,"",kTRUE) ; // array
290  pc.defineDouble("weight","Weight",0,1) ;
291  pc.defineObject("dummy1","ImportDataHistSliceMany",0) ;
292  pc.defineObject("dummy2","ImportHistoSliceMany",0) ;
293  pc.defineMutex("ImportHisto","ImportHistoSlice","ImportDataHistSlice") ;
294  pc.defineDependency("ImportHistoSlice","IndexCat") ;
295  pc.defineDependency("ImportDataHistSlice","IndexCat") ;
296 
297  RooLinkedList l ;
298  l.Add((TObject*)&arg1) ; l.Add((TObject*)&arg2) ;
299  l.Add((TObject*)&arg3) ; l.Add((TObject*)&arg4) ;
300  l.Add((TObject*)&arg5) ; l.Add((TObject*)&arg6) ;
301  l.Add((TObject*)&arg7) ; l.Add((TObject*)&arg8) ;
302 
303  // Process & check varargs
304  pc.process(l) ;
305  if (!pc.ok(kTRUE)) {
306  assert(0) ;
307  return ;
308  }
309 
310  TH1* impHist = static_cast<TH1*>(pc.getObject("impHist")) ;
311  Bool_t impDens = pc.getInt("impDens") ;
312  Double_t initWgt = pc.getDouble("weight") ;
313  const char* impSliceNames = pc.getString("impSliceState","",kTRUE) ;
314  const RooLinkedList& impSliceHistos = pc.getObjectList("impSliceHist") ;
315  RooCategory* indexCat = static_cast<RooCategory*>(pc.getObject("indexCat")) ;
316  const char* impSliceDNames = pc.getString("impSliceDState","",kTRUE) ;
317  const RooLinkedList& impSliceDHistos = pc.getObjectList("impSliceDHist") ;
318 
319 
320  if (impHist) {
321 
322  // Initialize importing contents from TH1
323  importTH1(vars,*impHist,initWgt, impDens) ;
324 
325  } else if (indexCat) {
326 
327 
328  if (impSliceHistos.GetSize()>0) {
329 
330  // Initialize importing mapped set of TH1s
331  map<string,TH1*> hmap ;
332  TIterator* hiter = impSliceHistos.MakeIterator() ;
333  for (const auto& token : RooHelpers::tokenise(impSliceNames, ",")) {
334  auto histo = static_cast<TH1*>(hiter->Next());
335  assert(histo);
336  hmap[token] = histo;
337  }
338  importTH1Set(vars,*indexCat,hmap,initWgt,kFALSE) ;
339  } else {
340 
341  // Initialize importing mapped set of RooDataHists
342  map<string,RooDataHist*> dmap ;
343  TIterator* hiter = impSliceDHistos.MakeIterator() ;
344  for (const auto& token : RooHelpers::tokenise(impSliceDNames, ",")) {
345  dmap[token] = (RooDataHist*) hiter->Next() ;
346  }
347  importDHistSet(vars,*indexCat,dmap,initWgt) ;
348  }
349 
350 
351  } else {
352 
353  // Initialize empty
354  initialize() ;
355  appendToDir(this,kTRUE) ;
356 
357  }
358 
359  _dstore->setExternalWeightArray(_wgt,_errLo,_errHi,_sumw2) ;
360  TRACE_CREATE
361 
362 }
363 
364 
365 
366 
367 ////////////////////////////////////////////////////////////////////////////////
368 /// Import data from given TH1/2/3 into this RooDataHist
369 
370 void RooDataHist::importTH1(const RooArgList& vars, const TH1& histo, Double_t wgt, Bool_t doDensityCorrection)
371 {
372  // Adjust binning of internal observables to match that of input THx
373  Int_t offset[3]{0, 0, 0};
374  adjustBinning(vars, histo, offset) ;
375 
376  // Initialize internal data structure
377  initialize() ;
378  appendToDir(this,kTRUE) ;
379 
380  // Define x,y,z as 1st, 2nd and 3rd observable
381  RooRealVar* xvar = (RooRealVar*) _vars.find(vars.at(0)->GetName()) ;
382  RooRealVar* yvar = (RooRealVar*) (vars.at(1) ? _vars.find(vars.at(1)->GetName()) : 0 ) ;
383  RooRealVar* zvar = (RooRealVar*) (vars.at(2) ? _vars.find(vars.at(2)->GetName()) : 0 ) ;
384 
385  // Transfer contents
386  Int_t xmin(0),ymin(0),zmin(0) ;
387  RooArgSet vset(*xvar) ;
388  Double_t volume = xvar->getMax()-xvar->getMin() ;
389  xmin = offset[0] ;
390  if (yvar) {
391  vset.add(*yvar) ;
392  ymin = offset[1] ;
393  volume *= (yvar->getMax()-yvar->getMin()) ;
394  }
395  if (zvar) {
396  vset.add(*zvar) ;
397  zmin = offset[2] ;
398  volume *= (zvar->getMax()-zvar->getMin()) ;
399  }
400  //Double_t avgBV = volume / numEntries() ;
401 // cout << "average bin volume = " << avgBV << endl ;
402 
403  Int_t ix(0),iy(0),iz(0) ;
404  for (ix=0 ; ix < xvar->getBins() ; ix++) {
405  xvar->setBin(ix) ;
406  if (yvar) {
407  for (iy=0 ; iy < yvar->getBins() ; iy++) {
408  yvar->setBin(iy) ;
409  if (zvar) {
410  for (iz=0 ; iz < zvar->getBins() ; iz++) {
411  zvar->setBin(iz) ;
412  Double_t bv = doDensityCorrection ? binVolume(vset) : 1;
413  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)) ;
414  }
415  } else {
416  Double_t bv = doDensityCorrection ? binVolume(vset) : 1;
417  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)) ;
418  }
419  }
420  } else {
421  Double_t bv = doDensityCorrection ? binVolume(vset) : 1 ;
422  add(vset,bv*histo.GetBinContent(ix+1+xmin)*wgt,bv*TMath::Power(histo.GetBinError(ix+1+xmin)*wgt,2)) ;
423  }
424  }
425 
426 }
427 
428 
429 
430 
431 
432 ////////////////////////////////////////////////////////////////////////////////
433 /// Import data from given set of TH1/2/3 into this RooDataHist. The category indexCat labels the sources
434 /// in the constructed RooDataHist. The stl map provides the mapping between the indexCat state labels
435 /// and the import source
436 
437 void RooDataHist::importTH1Set(const RooArgList& vars, RooCategory& indexCat, map<string,TH1*> hmap, Double_t wgt, Bool_t doDensityCorrection)
438 {
439  RooCategory* icat = (RooCategory*) _vars.find(indexCat.GetName()) ;
440 
441  TH1* histo(0) ;
442  Bool_t init(kFALSE) ;
443  for (map<string,TH1*>::iterator hiter = hmap.begin() ; hiter!=hmap.end() ; ++hiter) {
444  // Store pointer to first histogram from which binning specification will be taken
445  if (!histo) {
446  histo = hiter->second ;
447  }
448  // Define state labels in index category (both in provided indexCat and in internal copy in dataset)
449  if (!indexCat.lookupType(hiter->first.c_str())) {
450  indexCat.defineType(hiter->first.c_str()) ;
451  coutI(InputArguments) << "RooDataHist::importTH1Set(" << GetName() << ") defining state \"" << hiter->first << "\" in index category " << indexCat.GetName() << endl ;
452  }
453  if (!icat->lookupType(hiter->first.c_str())) {
454  icat->defineType(hiter->first.c_str()) ;
455  }
456  }
457 
458  // Check consistency in number of dimensions
459  if (histo && (vars.getSize() != histo->GetDimension())) {
460  coutE(InputArguments) << "RooDataHist::ctor(" << GetName() << ") ERROR: dimension of input histogram must match "
461  << "number of continuous variables" << endl ;
462  assert(0) ;
463  }
464 
465  // Copy bins and ranges from THx to dimension observables
466  Int_t offset[3] ;
467  adjustBinning(vars,*histo,offset) ;
468 
469  // Initialize internal data structure
470  if (!init) {
471  initialize() ;
472  appendToDir(this,kTRUE) ;
473  init = kTRUE ;
474  }
475 
476  // Define x,y,z as 1st, 2nd and 3rd observable
477  RooRealVar* xvar = (RooRealVar*) _vars.find(vars.at(0)->GetName()) ;
478  RooRealVar* yvar = (RooRealVar*) (vars.at(1) ? _vars.find(vars.at(1)->GetName()) : 0 ) ;
479  RooRealVar* zvar = (RooRealVar*) (vars.at(2) ? _vars.find(vars.at(2)->GetName()) : 0 ) ;
480 
481  // Transfer contents
482  Int_t xmin(0),ymin(0),zmin(0) ;
483  RooArgSet vset(*xvar) ;
484  Double_t volume = xvar->getMax()-xvar->getMin() ;
485  xmin = offset[0] ;
486  if (yvar) {
487  vset.add(*yvar) ;
488  ymin = offset[1] ;
489  volume *= (yvar->getMax()-yvar->getMin()) ;
490  }
491  if (zvar) {
492  vset.add(*zvar) ;
493  zmin = offset[2] ;
494  volume *= (zvar->getMax()-zvar->getMin()) ;
495  }
496  Double_t avgBV = volume / numEntries() ;
497 
498  Int_t ic(0),ix(0),iy(0),iz(0) ;
499  for (ic=0 ; ic < icat->numBins(0) ; ic++) {
500  icat->setBin(ic) ;
501  histo = hmap[icat->getLabel()] ;
502  for (ix=0 ; ix < xvar->getBins() ; ix++) {
503  xvar->setBin(ix) ;
504  if (yvar) {
505  for (iy=0 ; iy < yvar->getBins() ; iy++) {
506  yvar->setBin(iy) ;
507  if (zvar) {
508  for (iz=0 ; iz < zvar->getBins() ; iz++) {
509  zvar->setBin(iz) ;
510  Double_t bv = doDensityCorrection ? binVolume(vset)/avgBV : 1;
511  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)) ;
512  }
513  } else {
514  Double_t bv = doDensityCorrection ? binVolume(vset)/avgBV : 1;
515  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)) ;
516  }
517  }
518  } else {
519  Double_t bv = doDensityCorrection ? binVolume(vset)/avgBV : 1;
520  add(vset,bv*histo->GetBinContent(ix+1+xmin)*wgt,bv*TMath::Power(histo->GetBinError(ix+1+xmin)*wgt,2)) ;
521  }
522  }
523  }
524 
525 }
526 
527 
528 
529 ////////////////////////////////////////////////////////////////////////////////
530 /// Import data from given set of TH1/2/3 into this RooDataHist. The category indexCat labels the sources
531 /// in the constructed RooDataHist. The stl map provides the mapping between the indexCat state labels
532 /// and the import source
533 
534 void RooDataHist::importDHistSet(const RooArgList& /*vars*/, RooCategory& indexCat, std::map<std::string,RooDataHist*> dmap, Double_t initWgt)
535 {
536  RooCategory* icat = (RooCategory*) _vars.find(indexCat.GetName()) ;
537 
538  for (map<string,RooDataHist*>::iterator diter = dmap.begin() ; diter!=dmap.end() ; ++diter) {
539 
540  // Define state labels in index category (both in provided indexCat and in internal copy in dataset)
541  if (!indexCat.lookupType(diter->first.c_str())) {
542  indexCat.defineType(diter->first.c_str()) ;
543  coutI(InputArguments) << "RooDataHist::importDHistSet(" << GetName() << ") defining state \"" << diter->first << "\" in index category " << indexCat.GetName() << endl ;
544  }
545  if (!icat->lookupType(diter->first.c_str())) {
546  icat->defineType(diter->first.c_str()) ;
547  }
548  }
549 
550  initialize() ;
551  appendToDir(this,kTRUE) ;
552 
553 
554  for (map<string,RooDataHist*>::iterator diter = dmap.begin() ; diter!=dmap.end() ; ++diter) {
555 
556  RooDataHist* dhist = diter->second ;
557 
558  icat->setLabel(diter->first.c_str()) ;
559 
560  // Transfer contents
561  for (Int_t i=0 ; i<dhist->numEntries() ; i++) {
562  _vars = *dhist->get(i) ;
563  add(_vars,dhist->weight()*initWgt, pow(dhist->weightError(SumW2),2) ) ;
564  }
565 
566  }
567 }
568 
569 ////////////////////////////////////////////////////////////////////////////////
570 /// Helper doing the actual work of adjustBinning().
571 
572 void RooDataHist::_adjustBinning(RooRealVar &theirVar, const TAxis &axis,
573  RooRealVar *ourVar, Int_t *offset)
574 {
575  if (!dynamic_cast<RooRealVar*>(ourVar)) {
576  coutE(InputArguments) << "RooDataHist::adjustBinning(" << GetName() << ") ERROR: dimension " << ourVar->GetName() << " must be real" << endl ;
577  assert(0) ;
578  }
579 
580  const double xlo = theirVar.getMin();
581  const double xhi = theirVar.getMax();
582 
583  if (axis.GetXbins()->GetArray()) {
584  RooBinning xbins(axis.GetNbins(), axis.GetXbins()->GetArray());
585 
586  const double tolerance = 1e-6 * xbins.averageBinWidth();
587 
588  // Adjust xlo/xhi to nearest boundary
589  const double xloAdj = xbins.binLow(xbins.binNumber(xlo + tolerance));
590  const double xhiAdj = xbins.binHigh(xbins.binNumber(xhi - tolerance));
591  xbins.setRange(xloAdj, xhiAdj);
592 
593  theirVar.setBinning(xbins);
594 
595  if (true || fabs(xloAdj - xlo) > tolerance || fabs(xhiAdj - xhi) > tolerance) {
596  coutI(DataHandling)<< "RooDataHist::adjustBinning(" << GetName() << "): fit range of variable " << ourVar->GetName() << " expanded to nearest bin boundaries: [" << xlo << "," << xhi << "] --> [" << xloAdj << "," << xhiAdj << "]" << endl;
597  }
598 
599  ourVar->setBinning(xbins);
600 
601  if (offset) {
602  *offset = xbins.rawBinNumber(xloAdj + tolerance);
603  }
604  } else {
605  RooBinning xbins(axis.GetXmin(), axis.GetXmax());
606  xbins.addUniform(axis.GetNbins(), axis.GetXmin(), axis.GetXmax());
607 
608  const double tolerance = 1e-6 * xbins.averageBinWidth();
609 
610  // Adjust xlo/xhi to nearest boundary
611  const double xloAdj = xbins.binLow(xbins.binNumber(xlo + tolerance));
612  const double xhiAdj = xbins.binHigh(xbins.binNumber(xhi - tolerance));
613  xbins.setRange(xloAdj, xhiAdj);
614  theirVar.setRange(xloAdj, xhiAdj);
615 
616  if (fabs(xloAdj - xlo) > tolerance || fabs(xhiAdj - xhi) > tolerance) {
617  coutI(DataHandling)<< "RooDataHist::adjustBinning(" << GetName() << "): fit range of variable " << ourVar->GetName() << " expanded to nearest bin boundaries: [" << xlo << "," << xhi << "] --> [" << xloAdj << "," << xhiAdj << "]" << endl;
618  }
619 
620  RooUniformBinning xbins2(xloAdj, xhiAdj, xbins.numBins());
621  ourVar->setBinning(xbins2);
622 
623  if (offset) {
624  *offset = xbins.rawBinNumber(xloAdj + tolerance);
625  }
626  }
627 }
628 
629 ////////////////////////////////////////////////////////////////////////////////
630 /// Adjust binning specification on first and optionally second and third
631 /// observable to binning in given reference TH1. Used by constructors
632 /// that import data from an external TH1.
633 /// Both the variables in vars and in this RooDataHist are adjusted.
634 /// @param List with variables that are supposed to have their binning adjusted.
635 /// @param Reference histogram that dictates the binning
636 /// @param offset If not nullptr, a possible bin count offset for the axes x,y,z is saved here as Int_t[3]
637 
638 void RooDataHist::adjustBinning(const RooArgList& vars, const TH1& href, Int_t* offset)
639 {
640  auto xvar = static_cast<RooRealVar*>( _vars.find(*vars.at(0)) );
641  _adjustBinning(*static_cast<RooRealVar*>(vars.at(0)), *href.GetXaxis(), xvar, offset ? &offset[0] : nullptr);
642 
643  if (vars.at(1)) {
644  auto yvar = static_cast<RooRealVar*>(_vars.find(*vars.at(1)));
645  if (yvar)
646  _adjustBinning(*static_cast<RooRealVar*>(vars.at(1)), *href.GetYaxis(), yvar, offset ? &offset[1] : nullptr);
647  }
648 
649  if (vars.at(2)) {
650  auto zvar = static_cast<RooRealVar*>(_vars.find(*vars.at(2)));
651  if (zvar)
652  _adjustBinning(*static_cast<RooRealVar*>(vars.at(2)), *href.GetZaxis(), zvar, offset ? &offset[2] : nullptr);
653  }
654 
655 }
656 
657 
658 
659 
660 
661 ////////////////////////////////////////////////////////////////////////////////
662 /// Initialization procedure: allocate weights array, calculate
663 /// multipliers needed for N-space to 1-dim array jump table,
664 /// and fill the internal tree with all bin center coordinates
665 
666 void RooDataHist::initialize(const char* binningName, Bool_t fillTree)
667 {
668 
669  // Save real dimensions of dataset separately
670  for (const auto real : _vars) {
671  if (dynamic_cast<RooAbsReal*>(real)) _realVars.add(*real);
672  }
673 
674  // Fill array of LValue pointers to variables
675  for (const auto rvarg : _vars) {
676  if (binningName) {
677  RooRealVar* rrv = dynamic_cast<RooRealVar*>(rvarg);
678  if (rrv) {
679  rrv->setBinning(rrv->getBinning(binningName));
680  }
681  }
682  // coverity[FORWARD_NULL]
683  _lvvars.push_back(dynamic_cast<RooAbsLValue*>(rvarg));
684  // coverity[FORWARD_NULL]
685  const RooAbsBinning* binning = dynamic_cast<RooAbsLValue*>(rvarg)->getBinningPtr(0);
686  _lvbins.push_back(binning ? binning->clone() : 0);
687  }
688 
689 
690  // Allocate coefficients array
691  _idxMult.resize(_vars.getSize()) ;
692 
693  _arrSize = 1 ;
694  Int_t n(0), i ;
695  for (const auto var : _vars) {
696  auto arg = dynamic_cast<const RooAbsLValue*>(var);
697 
698  // Calculate sub-index multipliers for master index
699  for (i=0 ; i<n ; i++) {
700  _idxMult[i] *= arg->numBins() ;
701  }
702  _idxMult[n++] = 1 ;
703 
704  // Calculate dimension of weight array
705  _arrSize *= arg->numBins() ;
706  }
707 
708  // Allocate and initialize weight array if necessary
709  if (!_wgt) {
710  _wgt = new Double_t[_arrSize] ;
711  _errLo = new Double_t[_arrSize] ;
712  _errHi = new Double_t[_arrSize] ;
713  _sumw2 = new Double_t[_arrSize] ;
714  _binv = new Double_t[_arrSize] ;
715 
716  // Refill array pointers in data store when reading
717  // from Streamer
718  if (fillTree==kFALSE) {
719  _dstore->setExternalWeightArray(_wgt,_errLo,_errHi,_sumw2) ;
720  }
721 
722  for (i=0 ; i<_arrSize ; i++) {
723  _wgt[i] = 0 ;
724  _errLo[i] = -1 ;
725  _errHi[i] = -1 ;
726  _sumw2[i] = 0 ;
727  }
728  }
729 
730  if (!fillTree) return ;
731 
732  // Fill TTree with bin center coordinates
733  // Calculate plot bins of components from master index
734 
735  Int_t ibin ;
736  for (ibin=0 ; ibin<_arrSize ; ibin++) {
737  Int_t j(0), idx(0), tmp(ibin) ;
738  Double_t theBinVolume(1) ;
739  for (auto arg2 : _vars) {
740  idx = tmp / _idxMult[j] ;
741  tmp -= idx*_idxMult[j++] ;
742  auto arglv = dynamic_cast<RooAbsLValue*>(arg2);
743  arglv->setBin(idx) ;
744  theBinVolume *= arglv->getBinWidth(idx) ;
745 // cout << "init: bin width at idx=" << idx << " = " << arglv->getBinWidth(idx) << " binv[" << idx << "] = " << theBinVolume << endl ;
746  }
747  _binv[ibin] = theBinVolume ;
748 // cout << "binv[" << ibin << "] = " << theBinVolume << endl ;
749  fill() ;
750  }
751 
752 
753 }
754 
755 
756 ////////////////////////////////////////////////////////////////////////////////
757 
758 void RooDataHist::checkBinBounds() const
759 {
760  if (!_binbounds.empty()) return;
761  for (std::vector<const RooAbsBinning*>::const_iterator it = _lvbins.begin();
762  _lvbins.end() != it; ++it) {
763  _binbounds.push_back(std::vector<Double_t>());
764  if (*it) {
765  std::vector<Double_t>& bounds = _binbounds.back();
766  bounds.reserve(2 * (*it)->numBins());
767  for (Int_t i = 0; i < (*it)->numBins(); ++i) {
768  bounds.push_back((*it)->binLow(i));
769  bounds.push_back((*it)->binHigh(i));
770  }
771  }
772  }
773 }
774 
775 ////////////////////////////////////////////////////////////////////////////////
776 /// Copy constructor
777 
778 RooDataHist::RooDataHist(const RooDataHist& other, const char* newname) :
779  RooAbsData(other,newname), RooDirItem(), _idxMult(other._idxMult), _binValid(0), _curWeight(0), _curVolume(1), _pbinv(0), _pbinvCacheMgr(other._pbinvCacheMgr,0), _cache_sum_valid(0)
780 {
781  Int_t i ;
782 
783  // Allocate and initialize weight array
784  _arrSize = other._arrSize ;
785  _wgt = new Double_t[_arrSize] ;
786  _errLo = new Double_t[_arrSize] ;
787  _errHi = new Double_t[_arrSize] ;
788  _binv = new Double_t[_arrSize] ;
789  _sumw2 = new Double_t[_arrSize] ;
790  for (i=0 ; i<_arrSize ; i++) {
791  _wgt[i] = other._wgt[i] ;
792  _errLo[i] = other._errLo[i] ;
793  _errHi[i] = other._errHi[i] ;
794  _sumw2[i] = other._sumw2[i] ;
795  _binv[i] = other._binv[i] ;
796  }
797 
798  // Save real dimensions of dataset separately
799  for (const auto arg : _vars) {
800  if (dynamic_cast<RooAbsReal*>(arg) != nullptr) _realVars.add(*arg) ;
801  }
802 
803  // Fill array of LValue pointers to variables
804  for (const auto rvarg : _vars) {
805  // coverity[FORWARD_NULL]
806  _lvvars.push_back(dynamic_cast<RooAbsLValue*>(rvarg)) ;
807  // coverity[FORWARD_NULL]
808  const RooAbsBinning* binning = dynamic_cast<RooAbsLValue*>(rvarg)->getBinningPtr(0) ;
809  _lvbins.push_back(binning ? binning->clone() : 0) ;
810  }
811 
812  _dstore->setExternalWeightArray(_wgt,_errLo,_errHi,_sumw2) ;
813 
814  appendToDir(this,kTRUE) ;
815 }
816 
817 
818 
819 ////////////////////////////////////////////////////////////////////////////////
820 /// Constructor of a data hist from (part of) an existing data hist. The dimensions
821 /// of the data set are defined by the 'vars' RooArgSet, which can be identical
822 /// to 'dset' dimensions, or a subset thereof. Reduced dimensions will be projected
823 /// in the output data hist. The optional 'cutVar' formula variable can used to
824 /// select the subset of bins to be copied.
825 ///
826 /// For most uses the RooAbsData::reduce() wrapper function, which uses this constructor,
827 /// is the most convenient way to create a subset of an existing data
828 
829 RooDataHist::RooDataHist(const char* name, const char* title, RooDataHist* h, const RooArgSet& varSubset,
830  const RooFormulaVar* cutVar, const char* cutRange, Int_t nStart, Int_t nStop, Bool_t copyCache) :
831  RooAbsData(name,title,varSubset),
832  _wgt(0), _binValid(0), _curWeight(0), _curVolume(1), _pbinv(0), _pbinvCacheMgr(0,10), _cache_sum_valid(0)
833 {
834  // Initialize datastore
835  _dstore = new RooTreeDataStore(name,title,*h->_dstore,_vars,cutVar,cutRange,nStart,nStop,copyCache) ;
836 
837  initialize(0,kFALSE) ;
838 
839  _dstore->setExternalWeightArray(_wgt,_errLo,_errHi,_sumw2) ;
840 
841  // Copy weight array etc
842  Int_t i ;
843  for (i=0 ; i<_arrSize ; i++) {
844  _wgt[i] = h->_wgt[i] ;
845  _errLo[i] = h->_errLo[i] ;
846  _errHi[i] = h->_errHi[i] ;
847  _sumw2[i] = h->_sumw2[i] ;
848  _binv[i] = h->_binv[i] ;
849  }
850 
851 
852  appendToDir(this,kTRUE) ;
853  TRACE_CREATE
854 }
855 
856 
857 ////////////////////////////////////////////////////////////////////////////////
858 /// Construct a clone of this dataset that contains only the cached variables
859 
860 RooAbsData* RooDataHist::cacheClone(const RooAbsArg* newCacheOwner, const RooArgSet* newCacheVars, const char* newName)
861 {
862  checkInit() ;
863 
864  RooDataHist* dhist = new RooDataHist(newName?newName:GetName(),GetTitle(),this,*get(),0,0,0,2000000000,kTRUE) ;
865 
866  RooArgSet* selCacheVars = (RooArgSet*) newCacheVars->selectCommon(dhist->_cachedVars) ;
867  dhist->attachCache(newCacheOwner, *selCacheVars) ;
868  delete selCacheVars ;
869 
870  return dhist ;
871 }
872 
873 
874 
875 ////////////////////////////////////////////////////////////////////////////////
876 /// Implementation of RooAbsData virtual method that drives the RooAbsData::reduce() methods
877 
878 RooAbsData* RooDataHist::reduceEng(const RooArgSet& varSubset, const RooFormulaVar* cutVar, const char* cutRange,
879  Int_t nStart, Int_t nStop, Bool_t /*copyCache*/)
880 {
881  checkInit() ;
882  RooArgSet* myVarSubset = (RooArgSet*) _vars.selectCommon(varSubset) ;
883  RooDataHist *rdh = new RooDataHist(GetName(), GetTitle(), *myVarSubset) ;
884  delete myVarSubset ;
885 
886  RooFormulaVar* cloneVar = 0;
887  RooArgSet* tmp(0) ;
888  if (cutVar) {
889  // Deep clone cutVar and attach clone to this dataset
890  tmp = (RooArgSet*) RooArgSet(*cutVar).snapshot() ;
891  if (!tmp) {
892  coutE(DataHandling) << "RooDataHist::reduceEng(" << GetName() << ") Couldn't deep-clone cut variable, abort," << endl ;
893  return 0 ;
894  }
895  cloneVar = (RooFormulaVar*) tmp->find(*cutVar) ;
896  cloneVar->attachDataSet(*this) ;
897  }
898 
899  Int_t i ;
900  Double_t lo,hi ;
901  Int_t nevt = nStop < numEntries() ? nStop : numEntries() ;
902  TIterator* vIter = get()->createIterator() ;
903  for (i=nStart ; i<nevt ; i++) {
904  const RooArgSet* row = get(i) ;
905 
906  Bool_t doSelect(kTRUE) ;
907  if (cutRange) {
908  RooAbsArg* arg ;
909  vIter->Reset() ;
910  while((arg=(RooAbsArg*)vIter->Next())) {
911  if (!arg->inRange(cutRange)) {
912  doSelect = kFALSE ;
913  break ;
914  }
915  }
916  }
917  if (!doSelect) continue ;
918 
919  if (!cloneVar || cloneVar->getVal()) {
920  weightError(lo,hi,SumW2) ;
921  rdh->add(*row,weight(),lo*lo) ;
922  }
923  }
924  delete vIter ;
925 
926  if (cloneVar) {
927  delete tmp ;
928  }
929 
930  return rdh ;
931  }
932 
933 
934 
935 ////////////////////////////////////////////////////////////////////////////////
936 /// Destructor
937 
938 RooDataHist::~RooDataHist()
939 {
940  if (_wgt) delete[] _wgt ;
941  if (_errLo) delete[] _errLo ;
942  if (_errHi) delete[] _errHi ;
943  if (_sumw2) delete[] _sumw2 ;
944  if (_binv) delete[] _binv ;
945  if (_binValid) delete[] _binValid ;
946  vector<const RooAbsBinning*>::iterator iter = _lvbins.begin() ;
947  while(iter!=_lvbins.end()) {
948  delete *iter ;
949  ++iter ;
950  }
951 
952  removeFromDir(this) ;
953  TRACE_DESTROY
954 }
955 
956 
957 
958 
959 ////////////////////////////////////////////////////////////////////////////////
960 
961 Int_t RooDataHist::getIndex(const RooArgSet& coord, Bool_t fast)
962 {
963  checkInit() ;
964  if (fast) {
965  _vars.assignFast(coord,kFALSE) ;
966  } else {
967  _vars.assignValueOnly(coord) ;
968  }
969  return calcTreeIndex() ;
970 }
971 
972 
973 
974 
975 ////////////////////////////////////////////////////////////////////////////////
976 /// Calculate the index for the weights array corresponding to
977 /// to the bin enclosing the current coordinates of the internal argset
978 
979 Int_t RooDataHist::calcTreeIndex() const
980 {
981  int masterIdx(0);
982  for (unsigned int i=0; i < _lvvars.size(); ++i) {
983  const RooAbsLValue* lvvar = _lvvars[i];
984  const RooAbsBinning* binning = _lvbins[i];
985  masterIdx += _idxMult[i] * lvvar->getBin(binning);
986  }
987 
988  return masterIdx ;
989 }
990 
991 
992 
993 ////////////////////////////////////////////////////////////////////////////////
994 /// Debug stuff, should go...
995 
996 void RooDataHist::dump2()
997 {
998  cout << "_arrSize = " << _arrSize << endl ;
999  for (Int_t i=0 ; i<_arrSize ; i++) {
1000  cout << "wgt[" << i << "] = " << _wgt[i] << "sumw2[" << i << "] = " << _sumw2[i] << " vol[" << i << "] = " << _binv[i] << endl ;
1001  }
1002 }
1003 
1004 
1005 
1006 ////////////////////////////////////////////////////////////////////////////////
1007 /// Back end function to plotting functionality. Plot RooDataHist on given
1008 /// frame in mode specified by plot options 'o'. The main purpose of
1009 /// this function is to match the specified binning on 'o' to the
1010 /// internal binning of the plot observable in this RooDataHist.
1011 
1012 RooPlot *RooDataHist::plotOn(RooPlot *frame, PlotOpt o) const
1013 {
1014  checkInit() ;
1015  if (o.bins) return RooAbsData::plotOn(frame,o) ;
1016 
1017  if(0 == frame) {
1018  coutE(InputArguments) << ClassName() << "::" << GetName() << ":plotOn: frame is null" << endl;
1019  return 0;
1020  }
1021  RooAbsRealLValue *var= (RooAbsRealLValue*) frame->getPlotVar();
1022  if(0 == var) {
1023  coutE(InputArguments) << ClassName() << "::" << GetName()
1024  << ":plotOn: frame does not specify a plot variable" << endl;
1025  return 0;
1026  }
1027 
1028  RooRealVar* dataVar = (RooRealVar*) _vars.find(*var) ;
1029  if (!dataVar) {
1030  coutE(InputArguments) << ClassName() << "::" << GetName()
1031  << ":plotOn: dataset doesn't contain plot frame variable" << endl;
1032  return 0;
1033  }
1034 
1035  o.bins = &dataVar->getBinning() ;
1036  o.correctForBinWidth = kFALSE ;
1037  return RooAbsData::plotOn(frame,o) ;
1038 }
1039 
1040 
1041 
1042 
1043 ////////////////////////////////////////////////////////////////////////////////
1044 
1045 Double_t RooDataHist::weightSquared() const {
1046  return _curSumW2 ;
1047 }
1048 
1049 
1050 
1051 ////////////////////////////////////////////////////////////////////////////////
1052 /// Return the weight at given coordinates with optional
1053 /// interpolation. If intOrder is zero, the weight
1054 /// for the bin enclosing the coordinates
1055 /// contained in 'bin' is returned. For higher values,
1056 /// the result is interpolated in the real dimensions
1057 /// of the dataset
1058 
1059 Double_t RooDataHist::weight(const RooArgSet& bin, Int_t intOrder, Bool_t correctForBinSize, Bool_t cdfBoundaries, Bool_t oneSafe)
1060 {
1061  //cout << "RooDataHist::weight(" << bin << "," << intOrder << "," << correctForBinSize << "," << cdfBoundaries << "," << oneSafe << ")" << endl ;
1062 
1063  checkInit() ;
1064 
1065  // Handle illegal intOrder values
1066  if (intOrder<0) {
1067  coutE(InputArguments) << "RooDataHist::weight(" << GetName() << ") ERROR: interpolation order must be positive" << endl ;
1068  return 0 ;
1069  }
1070 
1071  // Handle no-interpolation case
1072  if (intOrder==0) {
1073  _vars.assignValueOnly(bin,oneSafe) ;
1074  Int_t idx = calcTreeIndex() ;
1075  //cout << "intOrder 0, idx = " << idx << endl ;
1076  if (correctForBinSize) {
1077  //calculatePartialBinVolume(*get()) ;
1078  //cout << "binw[" << idx << "] = " << _wgt[idx] << " / " << _binv[idx] << endl ;
1079  return get_wgt(idx) / _binv[idx];
1080  } else {
1081  //cout << "binw[" << idx << "] = " << _wgt[idx] << endl ;
1082  return get_wgt(idx);
1083  }
1084  }
1085 
1086  // Handle all interpolation cases
1087  _vars.assignValueOnly(bin) ;
1088 
1089  Double_t wInt(0) ;
1090  if (_realVars.getSize()==1) {
1091 
1092  // 1-dimensional interpolation
1093  const auto real = static_cast<RooRealVar*>(_realVars[static_cast<std::size_t>(0)]);
1094  const RooAbsBinning* binning = real->getBinningPtr(0) ;
1095  wInt = interpolateDim(*real,binning,((RooAbsReal*)bin.find(*real))->getVal(), intOrder, correctForBinSize, cdfBoundaries) ;
1096 
1097  } else if (_realVars.getSize()==2) {
1098 
1099  // 2-dimensional interpolation
1100  const auto realX = static_cast<RooRealVar*>(_realVars[static_cast<std::size_t>(0)]);
1101  const auto realY = static_cast<RooRealVar*>(_realVars[static_cast<std::size_t>(1)]);
1102  Double_t xval = ((RooAbsReal*)bin.find(*realX))->getVal() ;
1103  Double_t yval = ((RooAbsReal*)bin.find(*realY))->getVal() ;
1104 
1105  Int_t ybinC = realY->getBin() ;
1106  Int_t ybinLo = ybinC-intOrder/2 - ((yval<realY->getBinning().binCenter(ybinC))?1:0) ;
1107  Int_t ybinM = realY->numBins() ;
1108 
1109  Int_t i ;
1110  Double_t yarr[10] ;
1111  Double_t xarr[10] ;
1112  const RooAbsBinning* binning = realX->getBinningPtr(0) ;
1113  for (i=ybinLo ; i<=intOrder+ybinLo ; i++) {
1114  Int_t ibin ;
1115  if (i>=0 && i<ybinM) {
1116  // In range
1117  ibin = i ;
1118  realY->setBin(ibin) ;
1119  xarr[i-ybinLo] = realY->getVal() ;
1120  } else if (i>=ybinM) {
1121  // Overflow: mirror
1122  ibin = 2*ybinM-i-1 ;
1123  realY->setBin(ibin) ;
1124  xarr[i-ybinLo] = 2*realY->getMax()-realY->getVal() ;
1125  } else {
1126  // Underflow: mirror
1127  ibin = -i -1;
1128  realY->setBin(ibin) ;
1129  xarr[i-ybinLo] = 2*realY->getMin()-realY->getVal() ;
1130  }
1131  yarr[i-ybinLo] = interpolateDim(*realX,binning,xval,intOrder,correctForBinSize,kFALSE) ;
1132  }
1133 
1134  if (gDebug>7) {
1135  cout << "RooDataHist interpolating data is" << endl ;
1136  cout << "xarr = " ;
1137  for (int q=0; q<=intOrder ; q++) cout << xarr[q] << " " ;
1138  cout << " yarr = " ;
1139  for (int q=0; q<=intOrder ; q++) cout << yarr[q] << " " ;
1140  cout << endl ;
1141  }
1142  wInt = RooMath::interpolate(xarr,yarr,intOrder+1,yval) ;
1143 
1144  } else {
1145 
1146  // Higher dimensional scenarios not yet implemented
1147  coutE(InputArguments) << "RooDataHist::weight(" << GetName() << ") interpolation in "
1148  << _realVars.getSize() << " dimensions not yet implemented" << endl ;
1149  return weight(bin,0) ;
1150 
1151  }
1152 
1153  // Cut off negative values
1154 // if (wInt<=0) {
1155 // wInt=0 ;
1156 // }
1157 
1158  //cout << "RooDataHist wInt = " << wInt << endl ;
1159  return wInt ;
1160 }
1161 
1162 
1163 
1164 
1165 ////////////////////////////////////////////////////////////////////////////////
1166 /// Return the error on current weight
1167 
1168 void RooDataHist::weightError(Double_t& lo, Double_t& hi, ErrorType etype) const
1169 {
1170  checkInit() ;
1171 
1172  switch (etype) {
1173 
1174  case Auto:
1175  throw string(Form("RooDataHist::weightError(%s) error type Auto not allowed here",GetName())) ;
1176  break ;
1177 
1178  case Expected:
1179  throw string(Form("RooDataHist::weightError(%s) error type Expected not allowed here",GetName())) ;
1180  break ;
1181 
1182  case Poisson:
1183  if (_curWgtErrLo>=0) {
1184  // Weight is preset or precalculated
1185  lo = _curWgtErrLo ;
1186  hi = _curWgtErrHi ;
1187  return ;
1188  }
1189 
1190  // Calculate poisson errors
1191  Double_t ym,yp ;
1192  RooHistError::instance().getPoissonInterval(Int_t(weight()+0.5),ym,yp,1) ;
1193  _curWgtErrLo = weight()-ym ;
1194  _curWgtErrHi = yp-weight() ;
1195  _errLo[_curIndex] = _curWgtErrLo ;
1196  _errHi[_curIndex] = _curWgtErrHi ;
1197  lo = _curWgtErrLo ;
1198  hi = _curWgtErrHi ;
1199  return ;
1200 
1201  case SumW2:
1202  lo = sqrt(_curSumW2) ;
1203  hi = sqrt(_curSumW2) ;
1204  return ;
1205 
1206  case None:
1207  lo = 0 ;
1208  hi = 0 ;
1209  return ;
1210  }
1211 }
1212 
1213 
1214 // wve adjust for variable bin sizes
1215 
1216 ////////////////////////////////////////////////////////////////////////////////
1217 /// Perform boundary safe 'intOrder'-th interpolation of weights in dimension 'dim'
1218 /// at current value 'xval'
1219 
1220 Double_t RooDataHist::interpolateDim(RooRealVar& dim, const RooAbsBinning* binning, Double_t xval, Int_t intOrder, Bool_t correctForBinSize, Bool_t cdfBoundaries)
1221 {
1222  // Fill workspace arrays spanning interpolation area
1223  Int_t fbinC = dim.getBin(*binning) ;
1224  Int_t fbinLo = fbinC-intOrder/2 - ((xval<binning->binCenter(fbinC))?1:0) ;
1225  Int_t fbinM = dim.numBins(*binning) ;
1226 
1227 
1228  Int_t i ;
1229  Double_t yarr[10] ;
1230  Double_t xarr[10] ;
1231  for (i=fbinLo ; i<=intOrder+fbinLo ; i++) {
1232  Int_t ibin ;
1233  if (i>=0 && i<fbinM) {
1234  // In range
1235  ibin = i ;
1236  dim.setBinFast(ibin,*binning) ;
1237  //cout << "INRANGE: dim.getVal(ibin=" << ibin << ") = " << dim.getVal() << endl ;
1238  xarr[i-fbinLo] = dim.getVal() ;
1239  Int_t idx = calcTreeIndex();
1240  yarr[i - fbinLo] = get_wgt(idx);
1241  if (correctForBinSize) yarr[i-fbinLo] /= _binv[idx] ;
1242  } else if (i>=fbinM) {
1243  // Overflow: mirror
1244  ibin = 2*fbinM-i-1 ;
1245  dim.setBinFast(ibin,*binning) ;
1246  //cout << "OVERFLOW: dim.getVal(ibin=" << ibin << ") = " << dim.getVal() << endl ;
1247  if (cdfBoundaries) {
1248  xarr[i-fbinLo] = dim.getMax()+1e-10*(i-fbinM+1) ;
1249  yarr[i-fbinLo] = 1.0 ;
1250  } else {
1251  Int_t idx = calcTreeIndex() ;
1252  xarr[i-fbinLo] = 2*dim.getMax()-dim.getVal() ;
1253  yarr[i - fbinLo] = get_wgt(idx);
1254  if (correctForBinSize)
1255  yarr[i - fbinLo] /= _binv[idx];
1256  }
1257  } else {
1258  // Underflow: mirror
1259  ibin = -i - 1 ;
1260  dim.setBinFast(ibin,*binning) ;
1261  //cout << "UNDERFLOW: dim.getVal(ibin=" << ibin << ") = " << dim.getVal() << endl ;
1262  if (cdfBoundaries) {
1263  xarr[i-fbinLo] = dim.getMin()-ibin*(1e-10) ; ;
1264  yarr[i-fbinLo] = 0.0 ;
1265  } else {
1266  Int_t idx = calcTreeIndex() ;
1267  xarr[i-fbinLo] = 2*dim.getMin()-dim.getVal() ;
1268  yarr[i - fbinLo] = get_wgt(idx);
1269  if (correctForBinSize)
1270  yarr[i - fbinLo] /= _binv[idx];
1271  }
1272  }
1273  //cout << "ibin = " << ibin << endl ;
1274  }
1275 // for (int k=0 ; k<=intOrder ; k++) {
1276 // cout << "k=" << k << " x = " << xarr[k] << " y = " << yarr[k] << endl ;
1277 // }
1278  dim.setBinFast(fbinC,*binning) ;
1279  Double_t ret = RooMath::interpolate(xarr,yarr,intOrder+1,xval) ;
1280  return ret ;
1281 }
1282 
1283 
1284 
1285 
1286 ////////////////////////////////////////////////////////////////////////////////
1287 /// Increment the weight of the bin enclosing the coordinates given
1288 /// by 'row' by the specified amount. Add the sum of weights squared
1289 /// for the bin by 'sumw2' rather than wgt^2
1290 
1291 void RooDataHist::add(const RooArgSet& row, Double_t wgt, Double_t sumw2)
1292 {
1293  checkInit() ;
1294 
1295 // cout << "RooDataHist::add() accepted coordinate is " << endl ;
1296 // _vars.Print("v") ;
1297 
1298  _vars = row ;
1299  Int_t idx = calcTreeIndex() ;
1300  _wgt[idx] += wgt ;
1301  _sumw2[idx] += (sumw2>0?sumw2:wgt*wgt) ;
1302  _errLo[idx] = -1 ;
1303  _errHi[idx] = -1 ;
1304 
1305  _cache_sum_valid = kFALSE ;
1306 }
1307 
1308 
1309 
1310 ////////////////////////////////////////////////////////////////////////////////
1311 /// Increment the weight of the bin enclosing the coordinates
1312 /// given by 'row' by the specified amount. Associate errors
1313 /// [wgtErrLo,wgtErrHi] with the event weight on this bin.
1314 
1315 void RooDataHist::set(const RooArgSet& row, Double_t wgt, Double_t wgtErrLo, Double_t wgtErrHi)
1316 {
1317  checkInit() ;
1318 
1319  _vars = row ;
1320  Int_t idx = calcTreeIndex() ;
1321  _wgt[idx] = wgt ;
1322  _errLo[idx] = wgtErrLo ;
1323  _errHi[idx] = wgtErrHi ;
1324 
1325  _cache_sum_valid = kFALSE ;
1326 }
1327 
1328 
1329 
1330 ////////////////////////////////////////////////////////////////////////////////
1331 /// Increment the weight of the bin enclosing the coordinates
1332 /// given by 'row' by the specified amount. Associate errors
1333 /// [wgtErrLo,wgtErrHi] with the event weight on this bin.
1334 
1335 void RooDataHist::set(Double_t wgt, Double_t wgtErr)
1336 {
1337  checkInit() ;
1338 
1339  if (_curIndex<0) {
1340  _curIndex = calcTreeIndex() ;
1341  }
1342 
1343  _wgt[_curIndex] = wgt ;
1344  _errLo[_curIndex] = wgtErr ;
1345  _errHi[_curIndex] = wgtErr ;
1346  _sumw2[_curIndex] = wgtErr*wgtErr ;
1347 
1348  _cache_sum_valid = kFALSE ;
1349 }
1350 
1351 
1352 
1353 ////////////////////////////////////////////////////////////////////////////////
1354 /// Increment the weight of the bin enclosing the coordinates
1355 /// given by 'row' by the specified amount. Associate errors
1356 /// [wgtErrLo,wgtErrHi] with the event weight on this bin.
1357 
1358 void RooDataHist::set(const RooArgSet& row, Double_t wgt, Double_t wgtErr)
1359 {
1360  checkInit() ;
1361 
1362  _vars = row ;
1363  Int_t idx = calcTreeIndex() ;
1364  _wgt[idx] = wgt ;
1365  _errLo[idx] = wgtErr ;
1366  _errHi[idx] = wgtErr ;
1367  _sumw2[idx] = wgtErr*wgtErr ;
1368 
1369  _cache_sum_valid = kFALSE ;
1370 }
1371 
1372 
1373 
1374 ////////////////////////////////////////////////////////////////////////////////
1375 /// Add all data points contained in 'dset' to this data set with given weight.
1376 /// Optional cut string expression selects the data points to be added and can
1377 /// reference any variable contained in this data set
1378 
1379 void RooDataHist::add(const RooAbsData& dset, const char* cut, Double_t wgt)
1380 {
1381  RooFormulaVar cutVar("select",cut,*dset.get()) ;
1382  add(dset,&cutVar,wgt) ;
1383 }
1384 
1385 
1386 
1387 ////////////////////////////////////////////////////////////////////////////////
1388 /// Add all data points contained in 'dset' to this data set with given weight.
1389 /// Optional RooFormulaVar pointer selects the data points to be added.
1390 
1391 void RooDataHist::add(const RooAbsData& dset, const RooFormulaVar* cutVar, Double_t wgt)
1392 {
1393  checkInit() ;
1394 
1395  RooFormulaVar* cloneVar = 0;
1396  RooArgSet* tmp(0) ;
1397  if (cutVar) {
1398  // Deep clone cutVar and attach clone to this dataset
1399  tmp = (RooArgSet*) RooArgSet(*cutVar).snapshot() ;
1400  if (!tmp) {
1401  coutE(DataHandling) << "RooDataHist::add(" << GetName() << ") Couldn't deep-clone cut variable, abort," << endl ;
1402  return ;
1403  }
1404 
1405  cloneVar = (RooFormulaVar*) tmp->find(*cutVar) ;
1406  cloneVar->attachDataSet(dset) ;
1407  }
1408 
1409 
1410  Int_t i ;
1411  for (i=0 ; i<dset.numEntries() ; i++) {
1412  const RooArgSet* row = dset.get(i) ;
1413  if (!cloneVar || cloneVar->getVal()) {
1414  add(*row,wgt*dset.weight(), wgt*wgt*dset.weightSquared()) ;
1415  }
1416  }
1417 
1418  if (cloneVar) {
1419  delete tmp ;
1420  }
1421 
1422  _cache_sum_valid = kFALSE ;
1423 }
1424 
1425 
1426 
1427 ////////////////////////////////////////////////////////////////////////////////
1428 /// Return the sum of the weights of all hist bins.
1429 ///
1430 /// If correctForBinSize is specified, the sum of weights
1431 /// is multiplied by the N-dimensional bin volume,
1432 /// making the return value the integral over the function
1433 /// represented by this histogram
1434 
1435 Double_t RooDataHist::sum(Bool_t correctForBinSize, Bool_t inverseBinCor) const
1436 {
1437  checkInit() ;
1438 
1439  // Check if result was cached
1440  Int_t cache_code = 1 + (correctForBinSize?1:0) + ((correctForBinSize&&inverseBinCor)?1:0) ;
1441  if (_cache_sum_valid==cache_code) {
1442  return _cache_sum ;
1443  }
1444 
1445  Int_t i ;
1446  Double_t total(0), carry(0);
1447  for (i=0 ; i<_arrSize ; i++) {
1448 
1449  Double_t theBinVolume = correctForBinSize ? (inverseBinCor ? 1/_binv[i] : _binv[i]) : 1.0 ;
1450  // cout << "total += " << _wgt[i] << "*" << theBinVolume << endl ;
1451  // Double_t y = _wgt[i]*theBinVolume - carry;
1452  Double_t y = get_wgt(i) * theBinVolume - carry;
1453  Double_t t = total + y;
1454  carry = (t - total) - y;
1455  total = t;
1456  }
1457 
1458  // Store result in cache
1459  _cache_sum_valid=cache_code ;
1460  _cache_sum = total ;
1461 
1462  return total ;
1463 }
1464 
1465 
1466 
1467 ////////////////////////////////////////////////////////////////////////////////
1468 /// Return the sum of the weights of a multi-dimensional slice of the histogram
1469 /// by summing only over the dimensions specified in sumSet.
1470 ///
1471 /// The coordinates of all other dimensions are fixed to those given in sliceSet
1472 ///
1473 /// If correctForBinSize is specified, the sum of weights
1474 /// is multiplied by the M-dimensional bin volume, (M = N(sumSet)),
1475 /// making the return value the integral over the function
1476 /// represented by this histogram
1477 
1478 Double_t RooDataHist::sum(const RooArgSet& sumSet, const RooArgSet& sliceSet, Bool_t correctForBinSize, Bool_t inverseBinCor)
1479 {
1480  checkInit() ;
1481 
1482  RooArgSet varSave ;
1483  varSave.addClone(_vars) ;
1484 
1485  RooArgSet* sliceOnlySet = new RooArgSet(sliceSet) ;
1486  sliceOnlySet->remove(sumSet,kTRUE,kTRUE) ;
1487 
1488  _vars = *sliceOnlySet ;
1489  calculatePartialBinVolume(*sliceOnlySet) ;
1490  delete sliceOnlySet ;
1491 
1492  // Calculate mask and refence plot bins for non-iterating variables
1493  Bool_t* mask = new Bool_t[_vars.getSize()] ;
1494  Int_t* refBin = new Int_t[_vars.getSize()] ;
1495 
1496  for (unsigned int i = 0; i < _vars.size(); ++i) {
1497  const auto arg = _vars[i];
1498 
1499  if (sumSet.find(*arg)) {
1500  mask[i] = kFALSE ;
1501  } else {
1502  mask[i] = kTRUE ;
1503  refBin[i] = dynamic_cast<RooAbsLValue*>(arg)->getBin();
1504  }
1505  }
1506 
1507  // Loop over entire data set, skipping masked entries
1508  Double_t total(0), carry(0);
1509  Int_t ibin ;
1510  for (ibin=0 ; ibin<_arrSize ; ibin++) {
1511 
1512  Int_t idx(0), tmp(ibin), ivar(0) ;
1513  Bool_t skip(kFALSE) ;
1514 
1515  // Check if this bin belongs in selected slice
1516  for (unsigned int i = 0; !skip && i < _vars.size(); ++i) {
1517  idx = tmp / _idxMult[ivar] ;
1518  tmp -= idx*_idxMult[ivar] ;
1519  if (mask[ivar] && idx!=refBin[ivar]) skip=kTRUE ;
1520  ivar++ ;
1521  }
1522 
1523  if (!skip) {
1524  Double_t theBinVolume = correctForBinSize ? (inverseBinCor ? 1/(*_pbinv)[_vars.size()] : (*_pbinv)[_vars.size()] ) : 1.0 ;
1525  // cout << "adding bin[" << ibin << "] to sum wgt = " << _wgt[ibin] << " binv = " << theBinVolume << endl ;
1526  // Double_t y = _wgt[ibin]*theBinVolume - carry;
1527  Double_t y = get_wgt(ibin) * theBinVolume - carry;
1528  Double_t t = total + y;
1529  carry = (t - total) - y;
1530  total = t;
1531  }
1532  }
1533 
1534  delete[] mask ;
1535  delete[] refBin ;
1536 
1537  _vars = varSave ;
1538 
1539  return total ;
1540 }
1541 
1542 ////////////////////////////////////////////////////////////////////////////////
1543 /// Return the sum of the weights of a multi-dimensional slice of the histogram
1544 /// by summing only over the dimensions specified in sumSet.
1545 ///
1546 /// The coordinates of all other dimensions are fixed to those given in sliceSet
1547 ///
1548 /// If correctForBinSize is specified, the sum of weights
1549 /// is multiplied by the M-dimensional bin volume, (M = N(sumSet)),
1550 /// or the fraction of it that falls inside the range rangeName,
1551 /// making the return value the integral over the function
1552 /// represented by this histogram.
1553 ///
1554 /// If correctForBinSize is not specified, the weights are multiplied by the
1555 /// fraction of the bin volume that falls inside the range, i.e. a factor of
1556 /// binVolumeInRange/totalBinVolume.
1557 
1558 Double_t RooDataHist::sum(const RooArgSet& sumSet, const RooArgSet& sliceSet,
1559  Bool_t correctForBinSize, Bool_t inverseBinCor,
1560  const std::map<const RooAbsArg*, std::pair<Double_t, Double_t> >& ranges)
1561 {
1562  checkInit();
1563  checkBinBounds();
1564  RooArgSet varSave;
1565  varSave.addClone(_vars);
1566  {
1567  RooArgSet sliceOnlySet(sliceSet);
1568  sliceOnlySet.remove(sumSet,kTRUE,kTRUE);
1569  _vars = sliceOnlySet;
1570  }
1571 
1572  // Calculate mask and reference plot bins for non-iterating variables,
1573  // and get ranges for iterating variables
1574  std::vector<bool> mask(_vars.getSize());
1575  std::vector<Int_t> refBin(_vars.getSize());
1576  std::vector<Double_t> rangeLo(_vars.getSize(), -std::numeric_limits<Double_t>::infinity());
1577  std::vector<Double_t> rangeHi(_vars.getSize(), +std::numeric_limits<Double_t>::infinity());
1578 
1579  for (std::size_t i = 0; i < _vars.size(); ++i) {
1580  const auto arg = _vars[i];
1581  RooAbsArg* sumsetv = sumSet.find(*arg);
1582  RooAbsArg* slicesetv = sliceSet.find(*arg);
1583  mask[i] = !sumsetv;
1584  if (mask[i]) {
1585  auto argLV = dynamic_cast<const RooAbsLValue*>(arg);
1586  assert(argLV);
1587  refBin[i] = argLV->getBin();
1588  }
1589 
1590  auto it = ranges.find(sumsetv ? sumsetv : slicesetv);
1591  if (ranges.end() != it) {
1592  rangeLo[i] = it->second.first;
1593  rangeHi[i] = it->second.second;
1594  }
1595  }
1596 
1597  // Loop over entire data set, skipping masked entries
1598  Double_t total(0), carry(0);
1599  for (Int_t ibin = 0; ibin < _arrSize; ++ibin) {
1600  // Check if this bin belongs in selected slice
1601  Bool_t skip(kFALSE);
1602  for (int ivar = 0, tmp = ibin; !skip && ivar < int(_vars.size()); ++ivar) {
1603  const Int_t idx = tmp / _idxMult[ivar];
1604  tmp -= idx*_idxMult[ivar];
1605  if (mask[ivar] && idx!=refBin[ivar]) skip=kTRUE;
1606  }
1607 
1608  if (skip) continue;
1609 
1610  // work out bin volume
1611  Double_t theBinVolume = 1.;
1612  for (Int_t ivar = 0, tmp = ibin; ivar < (int)_vars.size(); ++ivar) {
1613  const Int_t idx = tmp / _idxMult[ivar];
1614  tmp -= idx*_idxMult[ivar];
1615  if (_binbounds[ivar].empty()) continue;
1616  const Double_t binLo = _binbounds[ivar][2 * idx];
1617  const Double_t binHi = _binbounds[ivar][2 * idx + 1];
1618  if (binHi < rangeLo[ivar] || binLo > rangeHi[ivar]) {
1619  // bin is outside of allowed range - effective bin volume is zero
1620  theBinVolume = 0.;
1621  break;
1622  }
1623  theBinVolume *=
1624  (std::min(rangeHi[ivar], binHi) - std::max(rangeLo[ivar], binLo));
1625  }
1626  const Double_t corrPartial = theBinVolume / _binv[ibin];
1627  if (0. == corrPartial) continue;
1628  const Double_t corr = correctForBinSize ? (inverseBinCor ? 1. / _binv[ibin] : _binv[ibin] ) : 1.0;
1629  //cout << "adding bin[" << ibin << "] to sum wgt = " << _wgt[ibin] << " binv = " << theBinVolume << " _binv[" << ibin << "] " << _binv[ibin] << endl;
1630  // const Double_t y = _wgt[ibin] * corr * corrPartial - carry;
1631  const Double_t y = get_wgt(ibin) * corr * corrPartial - carry;
1632  const Double_t t = total + y;
1633  carry = (t - total) - y;
1634  total = t;
1635  }
1636 
1637  _vars = varSave;
1638 
1639  return total;
1640 }
1641 
1642 
1643 
1644 ////////////////////////////////////////////////////////////////////////////////
1645 /// Fill the transient cache with partial bin volumes with up-to-date
1646 /// values for the partial volume specified by observables 'dimSet'
1647 
1648 void RooDataHist::calculatePartialBinVolume(const RooArgSet& dimSet) const
1649 {
1650  // Allocate cache if not yet existing
1651  vector<Double_t> *pbinv = _pbinvCacheMgr.getObj(&dimSet) ;
1652  if (pbinv) {
1653  _pbinv = pbinv ;
1654  return ;
1655  }
1656 
1657  pbinv = new vector<Double_t>(_arrSize) ;
1658 
1659  // Calculate plot bins of components from master index
1660  Bool_t* selDim = new Bool_t[_vars.getSize()] ;
1661  Int_t i(0) ;
1662  for (const auto v : _vars) {
1663  selDim[i++] = dimSet.find(*v) ? kTRUE : kFALSE ;
1664  }
1665 
1666  // Recalculate partial bin volume cache
1667  Int_t ibin ;
1668  for (ibin=0 ; ibin<_arrSize ; ibin++) {
1669  Int_t j(0), idx(0), tmp(ibin) ;
1670  Double_t theBinVolume(1) ;
1671  for (const auto absArg : _vars) {
1672  auto arg = dynamic_cast<const RooAbsLValue*>(absArg);
1673  if (!arg)
1674  break;
1675 
1676  idx = tmp / _idxMult[j] ;
1677  tmp -= idx*_idxMult[j++] ;
1678  if (selDim[j-1]) {
1679  theBinVolume *= arg->getBinWidth(idx) ;
1680  }
1681  }
1682  (*pbinv)[ibin] = theBinVolume ;
1683  }
1684 
1685  delete[] selDim ;
1686 
1687  // Put in cache (which takes ownership)
1688  _pbinvCacheMgr.setObj(&dimSet,pbinv) ;
1689 
1690  // Publicize the array
1691  _pbinv = pbinv ;
1692 }
1693 
1694 
1695 
1696 ////////////////////////////////////////////////////////////////////////////////
1697 /// Return the number of bins
1698 
1699 Int_t RooDataHist::numEntries() const
1700 {
1701  return RooAbsData::numEntries() ;
1702 }
1703 
1704 
1705 
1706 ////////////////////////////////////////////////////////////////////////////////
1707 
1708 Double_t RooDataHist::sumEntries() const
1709 {
1710  Int_t i ;
1711  Double_t n(0), carry(0);
1712  for (i=0 ; i<_arrSize ; i++) {
1713  if (!_binValid || _binValid[i]) {
1714  // Double_t y = _wgt[i] - carry;
1715  Double_t y = get_wgt(i) - carry;
1716  Double_t t = n + y;
1717  carry = (t - n) - y;
1718  n = t;
1719  }
1720  }
1721  return n ;
1722 }
1723 
1724 
1725 
1726 ////////////////////////////////////////////////////////////////////////////////
1727 /// Return the sum of weights in all entries matching cutSpec (if specified)
1728 /// and in named range cutRange (if specified)
1729 /// Return the
1730 
1731 Double_t RooDataHist::sumEntries(const char* cutSpec, const char* cutRange) const
1732 {
1733  checkInit() ;
1734 
1735  if (cutSpec==0 && cutRange==0) {
1736  return sumEntries();
1737  } else {
1738 
1739  // Setup RooFormulaVar for cutSpec if it is present
1740  RooFormula* select = 0 ;
1741  if (cutSpec) {
1742  select = new RooFormula("select",cutSpec,*get()) ;
1743  }
1744 
1745  // Otherwise sum the weights in the event
1746  Double_t sumw(0), carry(0);
1747  Int_t i ;
1748  for (i=0 ; i<numEntries() ; i++) {
1749  get(i) ;
1750  if (select && select->eval()==0.) continue ;
1751  if (cutRange && !_vars.allInRange(cutRange)) continue ;
1752 
1753  if (!_binValid || _binValid[i]) {
1754  Double_t y = weight() - carry;
1755  Double_t t = sumw + y;
1756  carry = (t - sumw) - y;
1757  sumw = t;
1758  }
1759  }
1760 
1761  if (select) delete select ;
1762 
1763  return sumw ;
1764  }
1765 }
1766 
1767 
1768 
1769 ////////////////////////////////////////////////////////////////////////////////
1770 /// Reset all bin weights to zero
1771 
1772 void RooDataHist::reset()
1773 {
1774  // WVE DO NOT CALL RooTreeData::reset() for binned
1775  // datasets as this will delete the bin definitions
1776 
1777  Int_t i ;
1778  for (i=0 ; i<_arrSize ; i++) {
1779  _wgt[i] = 0. ;
1780  _errLo[i] = -1 ;
1781  _errHi[i] = -1 ;
1782  }
1783  _curWeight = 0 ;
1784  _curWgtErrLo = -1 ;
1785  _curWgtErrHi = -1 ;
1786  _curVolume = 1 ;
1787 
1788  _cache_sum_valid = kFALSE ;
1789 
1790 }
1791 
1792 
1793 
1794 ////////////////////////////////////////////////////////////////////////////////
1795 /// Return an argset with the bin center coordinates for
1796 /// bin sequential number 'masterIdx'. For iterative use.
1797 
1798 const RooArgSet* RooDataHist::get(Int_t masterIdx) const
1799 {
1800  checkInit() ;
1801  _curWeight = _wgt[masterIdx] ;
1802  _curWgtErrLo = _errLo[masterIdx] ;
1803  _curWgtErrHi = _errHi[masterIdx] ;
1804  _curSumW2 = _sumw2[masterIdx] ;
1805  _curVolume = _binv[masterIdx] ;
1806  _curIndex = masterIdx ;
1807  return RooAbsData::get(masterIdx) ;
1808 }
1809 
1810 
1811 
1812 ////////////////////////////////////////////////////////////////////////////////
1813 /// Return a RooArgSet with center coordinates of the bin
1814 /// enclosing the point 'coord'
1815 
1816 const RooArgSet* RooDataHist::get(const RooArgSet& coord) const
1817 {
1818  ((RooDataHist*)this)->_vars = coord ;
1819  return get(calcTreeIndex()) ;
1820 }
1821 
1822 
1823 
1824 ////////////////////////////////////////////////////////////////////////////////
1825 /// Return the volume of the bin enclosing coordinates 'coord'
1826 
1827 Double_t RooDataHist::binVolume(const RooArgSet& coord)
1828 {
1829  checkInit() ;
1830  ((RooDataHist*)this)->_vars = coord ;
1831  return _binv[calcTreeIndex()] ;
1832 }
1833 
1834 
1835 ////////////////////////////////////////////////////////////////////////////////
1836 /// Set all the event weight of all bins to the specified value
1837 
1838 void RooDataHist::setAllWeights(Double_t value)
1839 {
1840  for (Int_t i=0 ; i<_arrSize ; i++) {
1841  _wgt[i] = value ;
1842  }
1843 
1844  _cache_sum_valid = kFALSE ;
1845 }
1846 
1847 
1848 
1849 ////////////////////////////////////////////////////////////////////////////////
1850 /// Create an iterator over all bins in a slice defined by the subset of observables
1851 /// listed in sliceArg. The position of the slice is given by otherArgs
1852 
1853 TIterator* RooDataHist::sliceIterator(RooAbsArg& sliceArg, const RooArgSet& otherArgs)
1854 {
1855  // Update to current position
1856  _vars = otherArgs ;
1857  _curIndex = calcTreeIndex() ;
1858 
1859  RooAbsArg* intArg = _vars.find(sliceArg) ;
1860  if (!intArg) {
1861  coutE(InputArguments) << "RooDataHist::sliceIterator() variable " << sliceArg.GetName() << " is not part of this RooDataHist" << endl ;
1862  return 0 ;
1863  }
1864  return new RooDataHistSliceIter(*this,*intArg) ;
1865 }
1866 
1867 
1868 ////////////////////////////////////////////////////////////////////////////////
1869 /// Change the name of the RooDataHist
1870 
1871 void RooDataHist::SetName(const char *name)
1872 {
1873  if (_dir) _dir->GetList()->Remove(this);
1874  TNamed::SetName(name) ;
1875  if (_dir) _dir->GetList()->Add(this);
1876 }
1877 
1878 
1879 ////////////////////////////////////////////////////////////////////////////////
1880 /// Change the title of this RooDataHist
1881 
1882 void RooDataHist::SetNameTitle(const char *name, const char* title)
1883 {
1884  if (_dir) _dir->GetList()->Remove(this);
1885  TNamed::SetNameTitle(name,title) ;
1886  if (_dir) _dir->GetList()->Add(this);
1887 }
1888 
1889 
1890 ////////////////////////////////////////////////////////////////////////////////
1891 /// Print value of the dataset, i.e. the sum of weights contained in the dataset
1892 
1893 void RooDataHist::printValue(ostream& os) const
1894 {
1895  os << numEntries() << " bins (" << sumEntries() << " weights)" ;
1896 }
1897 
1898 
1899 
1900 
1901 ////////////////////////////////////////////////////////////////////////////////
1902 /// Print argument of dataset, i.e. the observable names
1903 
1904 void RooDataHist::printArgs(ostream& os) const
1905 {
1906  os << "[" ;
1907  Bool_t first(kTRUE) ;
1908  for (const auto arg : _vars) {
1909  if (first) {
1910  first=kFALSE ;
1911  } else {
1912  os << "," ;
1913  }
1914  os << arg->GetName() ;
1915  }
1916  os << "]" ;
1917 }
1918 
1919 
1920 
1921 ////////////////////////////////////////////////////////////////////////////////
1922 /// Cache the datahist entries with bin centers that are inside/outside the
1923 /// current observable definitio
1924 
1925 void RooDataHist::cacheValidEntries()
1926 {
1927  checkInit() ;
1928 
1929  if (!_binValid) {
1930  _binValid = new Bool_t[_arrSize] ;
1931  }
1932  TIterator* iter = _vars.createIterator() ;
1933  RooAbsArg* arg ;
1934  for (Int_t i=0 ; i<_arrSize ; i++) {
1935  get(i) ;
1936  _binValid[i] = kTRUE ;
1937  iter->Reset() ;
1938  while((arg=(RooAbsArg*)iter->Next())) {
1939  // coverity[CHECKED_RETURN]
1940  _binValid[i] &= arg->inRange(0) ;
1941  }
1942  }
1943  delete iter ;
1944 
1945 }
1946 
1947 
1948 ////////////////////////////////////////////////////////////////////////////////
1949 /// Return true if currently loaded coordinate is considered valid within
1950 /// the current range definitions of all observables
1951 
1952 Bool_t RooDataHist::valid() const
1953 {
1954  // If caching is enabled, use the precached result
1955  if (_binValid) {
1956  return _binValid[_curIndex] ;
1957  }
1958 
1959  return kTRUE ;
1960 }
1961 
1962 
1963 
1964 ////////////////////////////////////////////////////////////////////////////////
1965 /// Returns true if datasets contains entries with a non-integer weight
1966 
1967 Bool_t RooDataHist::isNonPoissonWeighted() const
1968 {
1969  for (int i=0 ; i<numEntries() ; i++) {
1970  if (fabs(_wgt[i]-Int_t(_wgt[i]))>1e-10) return kTRUE ;
1971  }
1972  return kFALSE ;
1973 }
1974 
1975 
1976 
1977 
1978 ////////////////////////////////////////////////////////////////////////////////
1979 /// Print the details on the dataset contents
1980 
1981 void RooDataHist::printMultiline(ostream& os, Int_t content, Bool_t verbose, TString indent) const
1982 {
1983  RooAbsData::printMultiline(os,content,verbose,indent) ;
1984 
1985  os << indent << "Binned Dataset " << GetName() << " (" << GetTitle() << ")" << endl ;
1986  os << indent << " Contains " << numEntries() << " bins with a total weight of " << sumEntries() << endl;
1987 
1988  if (!verbose) {
1989  os << indent << " Observables " << _vars << endl ;
1990  } else {
1991  os << indent << " Observables: " ;
1992  _vars.printStream(os,kName|kValue|kExtras|kTitle,kVerbose,indent+" ") ;
1993  }
1994 
1995  if(verbose) {
1996  if (_cachedVars.getSize()>0) {
1997  os << indent << " Caches " << _cachedVars << endl ;
1998  }
1999  }
2000 }
2001 
2002 
2003 
2004 ////////////////////////////////////////////////////////////////////////////////
2005 /// Stream an object of class RooDataHist.
2006 
2007 void RooDataHist::Streamer(TBuffer &R__b)
2008 {
2009  if (R__b.IsReading()) {
2010 
2011  UInt_t R__s, R__c;
2012  Version_t R__v = R__b.ReadVersion(&R__s, &R__c);
2013 
2014  if (R__v>2) {
2015 
2016  R__b.ReadClassBuffer(RooDataHist::Class(),this,R__v,R__s,R__c);
2017  initialize(0,kFALSE) ;
2018 
2019  } else {
2020 
2021  // Legacy dataset conversion happens here. Legacy RooDataHist inherits from RooTreeData
2022  // which in turn inherits from RooAbsData. Manually stream RooTreeData contents on
2023  // file here and convert it into a RooTreeDataStore which is installed in the
2024  // new-style RooAbsData base class
2025 
2026  // --- This is the contents of the streamer code of RooTreeData version 2 ---
2027  UInt_t R__s1, R__c1;
2028  Version_t R__v1 = R__b.ReadVersion(&R__s1, &R__c1); if (R__v1) { }
2029 
2030  RooAbsData::Streamer(R__b);
2031  TTree* X_tree(0) ; R__b >> X_tree;
2032  RooArgSet X_truth ; X_truth.Streamer(R__b);
2033  TString X_blindString ; X_blindString.Streamer(R__b);
2034  R__b.CheckByteCount(R__s1, R__c1, RooTreeData::Class());
2035  // --- End of RooTreeData-v1 streamer
2036 
2037  // Construct RooTreeDataStore from X_tree and complete initialization of new-style RooAbsData
2038  _dstore = new RooTreeDataStore(X_tree,_vars) ;
2039  _dstore->SetName(GetName()) ;
2040  _dstore->SetTitle(GetTitle()) ;
2041  _dstore->checkInit() ;
2042 
2043  RooDirItem::Streamer(R__b);
2044  R__b >> _arrSize;
2045  delete [] _wgt;
2046  _wgt = new Double_t[_arrSize];
2047  R__b.ReadFastArray(_wgt,_arrSize);
2048  delete [] _errLo;
2049  _errLo = new Double_t[_arrSize];
2050  R__b.ReadFastArray(_errLo,_arrSize);
2051  delete [] _errHi;
2052  _errHi = new Double_t[_arrSize];
2053  R__b.ReadFastArray(_errHi,_arrSize);
2054  delete [] _sumw2;
2055  _sumw2 = new Double_t[_arrSize];
2056  R__b.ReadFastArray(_sumw2,_arrSize);
2057  delete [] _binv;
2058  _binv = new Double_t[_arrSize];
2059  R__b.ReadFastArray(_binv,_arrSize);
2060  _realVars.Streamer(R__b);
2061  R__b >> _curWeight;
2062  R__b >> _curWgtErrLo;
2063  R__b >> _curWgtErrHi;
2064  R__b >> _curSumW2;
2065  R__b >> _curVolume;
2066  R__b >> _curIndex;
2067  R__b.CheckByteCount(R__s, R__c, RooDataHist::IsA());
2068 
2069  }
2070 
2071  } else {
2072 
2073  R__b.WriteClassBuffer(RooDataHist::Class(),this);
2074  }
2075 }
2076 
2077