ROOT  6.07/01
Reference Guide
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
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 
void setInterpolationOrder(Int_t order)
Set interpolation order of RooHistFunct representing cache histogram.
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...
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...
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Add the specified argument to list.
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...
Class RooNumRunningInt is an implementation of RooAbsCachedReal that represents a running integral t...
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
const Bool_t kFALSE
Definition: Rtypes.h:92
int nbins[3]
ClassImp(RooNumRunningInt)
const RooAbsReal & arg() const
Definition: RooRealProxy.h:43
Double_t x[n]
Definition: legend1.C:17
virtual FuncCacheElem * createCache(const RooArgSet *nset) const
Create custom cache element for running integral calculations.
friend class RICacheElem
friend class RooArgSet
Definition: RooAbsArg.h:469
Double_t getVal(const RooArgSet *set=0) const
Definition: RooAbsReal.h:64
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:560
RooRealVar represents a fundamental (non-derived) real valued object.
Definition: RooRealVar.h:37
virtual ~RooNumRunningInt()
Destructor.
virtual void fillCacheObject(FuncCacheElem &cacheFunc) const
Fill the cache object by calling its calculate() method.
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
RooAbsCachedReal is the abstract base class for functions that need or want to cache their evaluate()...
RooRealProxy func
return
Definition: TBase64.cxx:62
TPaveLabel title(3, 27.1, 15, 28.7,"ROOT Environment and Tools")
void addPoint(Int_t ix)
Sample function at bin ix.
virtual RooArgList containedArgs(Action)
Return list of contained RooAbsArg objects.
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:51
double Double_t
Definition: RtypesCore.h:55
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:53
void calculate(Bool_t cdfmode)
Calculate the numeric running integral and store the result in the cache histogram provided by RooAbs...
double func(double *x, double *p)
Definition: stressTF1.cxx:213
virtual const char * inputBaseName() const
Return unique name for RooAbsCachedPdf cache components constructed from input function name...
#define name(a, b)
Definition: linkTestLib0.cpp:5
virtual Bool_t remove(const RooAbsArg &var, Bool_t silent=kFALSE, Bool_t matchByNameOnly=kFALSE)
Remove the specified argument from our list.
virtual RooArgSet * actualParameters(const RooArgSet &nset) const
Return the parameters of the cache created by RooAbsCachedPdf.
RICacheElem(const RooNumRunningInt &ri, const RooArgSet *nset)
Construct RunningIntegral CacheElement.
float type_of_call hi(const int &, const int &)
const Bool_t kTRUE
Definition: Rtypes.h:91
virtual RooArgList containedArgs(Action)
Return all RooAbsArg components contained in cache element.
virtual Double_t evaluate() const
Dummy function that is never called.
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Add element to non-owning set.
Definition: RooArgSet.cxx:448