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
17Class RooNumRunningInt is an implementation of RooAbsCachedReal that represents a running integral
18\f[ RI(f(x)) = \int_{xlow}^{x} f(x') dx' \f]
19that is calculated internally with a numeric technique: The input function
20is first sampled into a histogram, which is then numerically integrated.
21The output function is an interpolated version of the integrated histogram.
22The sampling density is controlled by the binning named "cache" in the observable x.
23The shape of the p.d.f is always calculated for the entire domain in x and
24cached in a histogram. The cache histogram is automatically recalculated
25when 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
38using 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
52RooNumRunningInt::RooNumRunningInt(const char *name, const char *title, RooAbsReal& _func, RooRealVar& _x, const char* bname) :
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
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
#define e(i)
Definition: RSha256.hxx:103
const Bool_t kFALSE
Definition: RtypesCore.h:90
const Bool_t kTRUE
Definition: RtypesCore.h:89
#define ClassImp(name)
Definition: Rtypes.h:361
char name[80]
Definition: TGX11.cxx:109
float type_of_call hi(const int &, const int &)
friend class RooArgSet
Definition: RooAbsArg.h:572
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:544
virtual RooArgList containedArgs(Action)
Return list of contained RooAbsArg objects.
RooAbsCachedReal is the abstract base class for functions that need or want to cache their evaluate()...
void setInterpolationOrder(Int_t order)
Set interpolation order of RooHistFunct representing cache histogram.
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Add the specified argument to list.
virtual Bool_t remove(const RooAbsArg &var, Bool_t silent=kFALSE, Bool_t matchByNameOnly=kFALSE)
Remove the specified argument from our list.
RooAbsArg * find(const char *name) const
Find object with given name in list.
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:60
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:90
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgList.h:21
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
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
virtual Int_t numEntries() const
Return the number of bins.
virtual const RooArgSet * get() const
Definition: RooDataHist.h:79
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...
virtual RooArgList containedArgs(Action)
Return all RooAbsArg components contained in cache element.
void addPoint(Int_t ix)
Sample function at bin ix.
void calculate(Bool_t cdfmode)
Calculate the numeric running integral and store the result in the cache histogram provided by RooAbs...
RICacheElem(const RooNumRunningInt &ri, const RooArgSet *nset)
Construct RunningIntegral CacheElement.
Class RooNumRunningInt is an implementation of RooAbsCachedReal that represents a running integral.
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 FuncCacheElem * createCache(const RooArgSet *nset) const
Create custom cache element for running integral calculations.
virtual Double_t evaluate() const
Dummy function that is never called.
RooRealProxy func
virtual RooArgSet * actualParameters(const RooArgSet &nset) const
Return the parameters of the cache created by RooAbsCachedPdf.
virtual void fillCacheObject(FuncCacheElem &cacheFunc) const
Fill the cache object by calling its calculate() method.
friend class RICacheElem
virtual ~RooNumRunningInt()
Destructor.
virtual const char * inputBaseName() const
Return unique name for RooAbsCachedPdf cache components constructed from input function name.
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:35
const T & arg() const
Return reference to object held in proxy.
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
Double_t x[n]
Definition: legend1.C:17
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)