Logo ROOT  
Reference Guide
RooNumRunningInt.cxx
Go to the documentation of this file.
1  /*****************************************************************************
2  * Project: RooFit *
3  * *
4  * Copyright (c) 2000-2005, Regents of the University of California *
5  * and Stanford University. All rights reserved. *
6  * *
7  * Redistribution and use in source and binary forms, *
8  * with or without modification, are permitted according to the terms *
9  * listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
10  *****************************************************************************/
11 
12 /**
13 \file RooNumRunningInt.cxx
14 \class RooNumRunningInt
15 \ingroup Roofitcore
16 
17 Class RooNumRunningInt is an implementation of RooAbsCachedReal that represents a running integral
18 \f[ RI(f(x)) = \int_{xlow}^{x} f(x') dx' \f]
19 that is calculated internally with a numeric technique: The input function
20 is first sampled into a histogram, which is then numerically integrated.
21 The output function is an interpolated version of the integrated histogram.
22 The sampling density is controlled by the binning named "cache" in the observable x.
23 The shape of the p.d.f is always calculated for the entire domain in x and
24 cached in a histogram. The cache histogram is automatically recalculated
25 when any of the parameters of the input p.d.f. has changed.
26 **/
27 
28 #include "Riostream.h"
29 
30 #include "RooAbsPdf.h"
31 #include "RooNumRunningInt.h"
32 #include "RooAbsReal.h"
33 #include "RooMsgService.h"
34 #include "RooDataHist.h"
35 #include "RooHistPdf.h"
36 #include "RooRealVar.h"
37 
38 using namespace std;
39 
41  ;
42 
43 
44 
45 ////////////////////////////////////////////////////////////////////////////////
46 /// Construct running integral of function '_func' over x_print from
47 /// the lower bound on _x to the present value of _x using a numeric
48 /// sampling technique. The sampling frequency is controlled by the
49 /// binning named 'bname' and a default second order interpolation
50 /// is applied to smooth the histogram-based c.d.f.
51 
52 RooNumRunningInt::RooNumRunningInt(const char *name, const char *title, RooAbsReal& _func, RooRealVar& _x, const char* bname) :
53  RooAbsCachedReal(name,title),
54  func("func","func",this,_func),
55  x("x","x",this,_x),
56  _binningName(bname?bname:"cache")
57  {
59  }
60 
61 
62 
63 
64 ////////////////////////////////////////////////////////////////////////////////
65 /// Copy constructor
66 
68  RooAbsCachedReal(other,name),
69  func("func",this,other.func),
70  x("x",this,other.x),
71  _binningName(other._binningName)
72  {
73  }
74 
75 
76 
77 ////////////////////////////////////////////////////////////////////////////////
78 /// Destructor
79 
81 {
82 }
83 
84 
85 ////////////////////////////////////////////////////////////////////////////////
86 /// Return unique name for RooAbsCachedPdf cache components
87 /// constructed from input function name
88 
90 {
91  static string ret ;
92  ret = func.arg().GetName() ;
93  ret += "_NUMRUNINT" ;
94  return ret.c_str() ;
95 } ;
96 
97 
98 
99 ////////////////////////////////////////////////////////////////////////////////
100 /// Construct RunningIntegral CacheElement
101 
103  FuncCacheElem(self,nset), _self(&const_cast<RooNumRunningInt&>(self))
104 {
105  // Instantiate temp arrays
106  _ax = new Double_t[hist()->numEntries()] ;
107  _ay = new Double_t[hist()->numEntries()] ;
108 
109  // Copy X values from histo
110  _xx = (RooRealVar*) hist()->get()->find(self.x.arg().GetName()) ;
111  for (int i=0 ; i<hist()->numEntries() ; i++) {
112  hist()->get(i) ;
113  _ax[i] = _xx->getVal() ;
114  _ay[i] = -1 ;
115  }
116 
117 }
118 
119 
120 ////////////////////////////////////////////////////////////////////////////////
121 /// Destructor
122 
124 {
125  // Delete temp arrays
126  delete[] _ax ;
127  delete[] _ay ;
128 }
129 
130 
131 ////////////////////////////////////////////////////////////////////////////////
132 /// Return all RooAbsArg components contained in cache element
133 
135 {
136  RooArgList ret ;
137  ret.add(FuncCacheElem::containedArgs(action)) ;
138  ret.add(*_self) ;
139  ret.add(*_xx) ;
140  return ret ;
141 }
142 
143 
144 
145 ////////////////////////////////////////////////////////////////////////////////
146 /// Calculate the numeric running integral and store
147 /// the result in the cache histogram provided
148 /// by RooAbsCachedPdf
149 
151 {
152  // Update contents of histogram
153  Int_t nbins = hist()->numEntries() ;
154 
155  Double_t xsave = _self->x ;
156 
157  Int_t lastHi=0 ;
158  Int_t nInitRange=32 ;
159  for (int i=1 ; i<=nInitRange ; i++) {
160  Int_t hi = (i*nbins)/nInitRange -1 ;
161  Int_t lo = lastHi ;
162  addRange(lo,hi,nbins) ;
163  lastHi=hi ;
164  }
165 
166  // Perform numeric integration
167  for (int i=1 ; i<nbins ; i++) {
168  _ay[i] += _ay[i-1] ;
169  }
170 
171  // Normalize and transfer to cache histogram
172  Double_t binv = (_self->x.max()-_self->x.min())/nbins ;
173  for (int i=0 ; i<nbins ; i++) {
174  hist()->get(i) ;
175  if (cdfmode) {
176  hist()->set(_ay[i]/_ay[nbins-1]) ;
177  } else {
178  hist()->set(_ay[i]*binv) ;
179  }
180  }
181 
182  if (cdfmode) {
183  func()->setCdfBoundaries(kTRUE) ;
184  }
185  _self->x = xsave ;
186 }
187 
188 
189 
190 ////////////////////////////////////////////////////////////////////////////////
191 /// Fill all empty histogram bins in the range [ixlo,ixhi] where nbins is the
192 /// total number of histogram bins. This method samples the mid-point of the
193 /// range and if the mid-point value is within small tolerance of the interpolated
194 /// mid-point value fills all remaining elements through linear interpolation.
195 /// If the tolerance is exceeded, the algorithm is recursed on the two subranges
196 /// [xlo,xmid] and [xmid,xhi]
197 
199 {
200  // Add first and last point, if not there already
201  if (_ay[ixlo]<0) {
202  addPoint(ixlo) ;
203  }
204  if (_ay[ixhi]<0) {
205  addPoint(ixhi) ;
206  }
207 
208  // Terminate here if there is no gap
209  if (ixhi-ixlo==1) {
210  return ;
211  }
212 
213  // If gap size is one, simply fill gap and return
214  if (ixhi-ixlo==2) {
215  addPoint(ixlo+1) ;
216  return ;
217  }
218 
219  // Add mid-point
220  Int_t ixmid = (ixlo+ixhi)/2 ;
221  addPoint(ixmid) ;
222 
223  // Calculate difference of mid-point w.r.t interpolated value
224  Double_t yInt = _ay[ixlo] + (_ay[ixhi]-_ay[ixlo])*(ixmid-ixlo)/(ixhi-ixlo) ;
225 
226  // If relative deviation is greater than tolerance divide and iterate
227  if (fabs(yInt-_ay[ixmid])*(_ax[nbins-1]-_ax[0])>1e-6) {
228  addRange(ixlo,ixmid,nbins) ;
229  addRange(ixmid,ixhi,nbins) ;
230  } else {
231  for (Int_t j=ixlo+1 ; j<ixmid ; j++) {
232  _ay[j] = _ay[ixlo] + (_ay[ixmid]-_ay[ixlo])*(j-ixlo)/(ixmid-ixlo) ;
233  }
234  for (Int_t j=ixmid+1 ; j<ixhi ; j++) {
235  _ay[j] = _ay[ixmid] + (_ay[ixhi]-_ay[ixmid])*(j-ixmid)/(ixhi-ixmid) ;
236  }
237  }
238 
239 }
240 
241 
242 ////////////////////////////////////////////////////////////////////////////////
243 /// Sample function at bin ix
244 
246 {
247  hist()->get(ix) ;
248  _self->x = _xx->getVal() ;
249  _ay[ix] = _self->func.arg().getVal(*_xx) ;
250 
251 }
252 
253 
254 
255 ////////////////////////////////////////////////////////////////////////////////
256 /// Fill the cache object by calling its calculate() method
257 
259 {
260  RICacheElem& riCache = static_cast<RICacheElem&>(cache) ;
261  riCache.calculate(kFALSE) ;
262 }
263 
264 
265 
266 ////////////////////////////////////////////////////////////////////////////////
267 /// Return observable in nset to be cached by RooAbsCachedPdf
268 /// this is always the x observable that is integrated
269 
271 {
272  RooArgSet* ret = new RooArgSet ;
273  ret->add(x.arg()) ;
274  return ret ;
275 }
276 
277 
278 
279 ////////////////////////////////////////////////////////////////////////////////
280 /// Return the parameters of the cache created by RooAbsCachedPdf.
281 /// These are always the input functions parameter, but never the
282 /// integrated variable x.
283 
285 {
287  ret->remove(x.arg(),kTRUE,kTRUE) ;
288  return ret ;
289 }
290 
291 
292 ////////////////////////////////////////////////////////////////////////////////
293 /// Create custom cache element for running integral calculations
294 
296 {
297  return new RICacheElem(*const_cast<RooNumRunningInt*>(this),nset) ;
298 }
299 
300 
301 ////////////////////////////////////////////////////////////////////////////////
302 /// Dummy function that is never called
303 
305 {
306  cout << "RooNumRunningInt::evaluate(" << GetName() << ")" << endl ;
307  return 0 ;
308 }
309 
310 
RooNumRunningInt::fillCacheObject
virtual void fillCacheObject(FuncCacheElem &cacheFunc) const
Fill the cache object by calling its calculate() method.
Definition: RooNumRunningInt.cxx:258
RooNumRunningInt::RICacheElem::calculate
void calculate(Bool_t cdfmode)
Calculate the numeric running integral and store the result in the cache histogram provided by RooAbs...
Definition: RooNumRunningInt.cxx:150
RooNumRunningInt::~RooNumRunningInt
virtual ~RooNumRunningInt()
Destructor.
Definition: RooNumRunningInt.cxx:80
RooNumRunningInt::RICacheElem::addRange
void addRange(Int_t ixlo, Int_t ixhi, Int_t nbins)
Fill all empty histogram bins in the range [ixlo,ixhi] where nbins is the total number of histogram b...
Definition: RooNumRunningInt.cxx:198
kTRUE
const Bool_t kTRUE
Definition: RtypesCore.h:91
e
#define e(i)
Definition: RSha256.hxx:121
RooAbsReal.h
RooNumRunningInt::RICacheElem::_ax
Double_t * _ax
Definition: RooNumRunningInt.h:48
RooMsgService.h
RooNumRunningInt::x
RooRealProxy x
Definition: RooNumRunningInt.h:66
RooNumRunningInt::actualParameters
virtual RooArgSet * actualParameters(const RooArgSet &nset) const
Return the parameters of the cache created by RooAbsCachedPdf.
Definition: RooNumRunningInt.cxx:284
RooTemplateProxy::arg
const T & arg() const
Return reference to object held in proxy.
Definition: RooTemplateProxy.h:271
ClassImp
#define ClassImp(name)
Definition: Rtypes.h:364
RooArgList
Definition: RooArgList.h:21
RooAbsReal::getVal
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:91
RooDataHist::get
virtual const RooArgSet * get() const
Definition: RooDataHist.h:78
RooAbsCollection::remove
virtual Bool_t remove(const RooAbsArg &var, Bool_t silent=kFALSE, Bool_t matchByNameOnly=kFALSE)
Remove the specified argument from our list.
Definition: RooAbsCollection.cxx:584
RooAbsCollection::find
RooAbsArg * find(const char *name) const
Find object with given name in list.
Definition: RooAbsCollection.cxx:813
x
Double_t x[n]
Definition: legend1.C:17
RooArgSet::add
virtual Bool_t add(const RooAbsCollection &col, Bool_t silent=kFALSE)
Add a collection of arguments to this collection by calling add() for each element in the source coll...
Definition: RooArgSet.h:88
RooAbsReal
Definition: RooAbsReal.h:61
RooNumRunningInt::RICacheElem::~RICacheElem
~RICacheElem()
Destructor.
Definition: RooNumRunningInt.cxx:123
bool
RooAbsArg::getParameters
RooArgSet * getParameters(const RooAbsData *data, Bool_t stripDisconnected=kTRUE) const
Create a list of leaf nodes in the arg tree starting with ourself as top node that don't match any of...
Definition: RooAbsArg.cxx:546
RooNumRunningInt::evaluate
virtual Double_t evaluate() const
Dummy function that is never called.
Definition: RooNumRunningInt.cxx:304
hi
float type_of_call hi(const int &, const int &)
RooNumRunningInt::RICacheElem::_xx
RooRealVar * _xx
Definition: RooNumRunningInt.h:50
ROOT::Math::fabs
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
Definition: UnaryOperators.h:131
RooNumRunningInt::actualObservables
virtual RooArgSet * actualObservables(const RooArgSet &nset) const
Return observable in nset to be cached by RooAbsCachedPdf this is always the x observable that is int...
Definition: RooNumRunningInt.cxx:270
RooAbsPdf.h
RooAbsCachedReal::FuncCacheElem::containedArgs
virtual RooArgList containedArgs(Action)
Return list of contained RooAbsArg objects.
Definition: RooAbsCachedReal.cxx:308
RooNumRunningInt::RICacheElem::containedArgs
virtual RooArgList containedArgs(Action)
Return all RooAbsArg components contained in cache element.
Definition: RooNumRunningInt.cxx:134
kFALSE
const Bool_t kFALSE
Definition: RtypesCore.h:92
RooDataHist.h
RooAbsArg::RooArgSet
friend class RooArgSet
Definition: RooAbsArg.h:579
RooNumRunningInt::RICacheElem
Definition: RooNumRunningInt.h:38
RooAbsCachedReal::FuncCacheElem::hist
RooDataHist * hist()
Definition: RooAbsCachedReal.h:69
RooAbsCollection::add
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Add the specified argument to list.
Definition: RooAbsCollection.cxx:437
RooNumRunningInt.h
RooRealVar.h
RooNumRunningInt
Definition: RooNumRunningInt.h:20
RooNumRunningInt::RooNumRunningInt
RooNumRunningInt(const char *name, const char *title, RooAbsReal &_func, RooRealVar &_x, const char *binningName="cache")
Construct running integral of function '_func' over x_print from the lower bound on _x to the present...
Definition: RooNumRunningInt.cxx:52
RooAbsCachedReal::setInterpolationOrder
void setInterpolationOrder(Int_t order)
Set interpolation order of RooHistFunct representing cache histogram.
Definition: RooAbsCachedReal.cxx:291
RooHistPdf.h
Double_t
double Double_t
Definition: RtypesCore.h:59
RooNumRunningInt::RICacheElem
friend class RICacheElem
Definition: RooNumRunningInt.h:54
RooNumRunningInt::RICacheElem::_ay
Double_t * _ay
Definition: RooNumRunningInt.h:49
RooDataHist::numEntries
virtual Int_t numEntries() const
Return the number of bins.
Definition: RooDataHist.cxx:1694
RooAbsCacheElement::Action
Action
Definition: RooAbsCacheElement.h:39
RooNumRunningInt::createCache
virtual FuncCacheElem * createCache(const RooArgSet *nset) const
Create custom cache element for running integral calculations.
Definition: RooNumRunningInt.cxx:295
name
char name[80]
Definition: TGX11.cxx:110
RooNumRunningInt::RICacheElem::RICacheElem
RICacheElem(const RooNumRunningInt &ri, const RooArgSet *nset)
Construct RunningIntegral CacheElement.
Definition: RooNumRunningInt.cxx:102
TNamed::GetName
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:53
RooNumRunningInt::inputBaseName
virtual const char * inputBaseName() const
Return unique name for RooAbsCachedPdf cache components constructed from input function name.
Definition: RooNumRunningInt.cxx:89
RooRealVar
Definition: RooRealVar.h:35
Riostream.h
RooAbsCachedReal::FuncCacheElem
Definition: RooAbsCachedReal.h:59
RooAbsCachedReal
Definition: RooAbsCachedReal.h:24
RooNumRunningInt::RICacheElem::addPoint
void addPoint(Int_t ix)
Sample function at bin ix.
Definition: RooNumRunningInt.cxx:245
RooArgSet
Definition: RooArgSet.h:28
RooNumRunningInt::func
RooRealProxy func
Definition: RooNumRunningInt.h:65
int