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.resize(hist()->numEntries());
107 _ay.resize(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/// Return all RooAbsArg components contained in cache element
122
124{
125 RooArgList ret ;
126 ret.add(FuncCacheElem::containedArgs(action)) ;
127 ret.add(*_self) ;
128 ret.add(*_xx) ;
129 return ret ;
130}
131
132
133
134////////////////////////////////////////////////////////////////////////////////
135/// Calculate the numeric running integral and store
136/// the result in the cache histogram provided
137/// by RooAbsCachedPdf
138
140{
141 // Update contents of histogram
142 Int_t nbins = hist()->numEntries() ;
143
144 double xsave = _self->x ;
145
146 Int_t lastHi=0 ;
147 Int_t nInitRange=32 ;
148 for (int i=1 ; i<=nInitRange ; i++) {
149 Int_t hi = (i*nbins)/nInitRange -1 ;
150 Int_t lo = lastHi ;
151 addRange(lo,hi,nbins) ;
152 lastHi=hi ;
153 }
154
155 // Perform numeric integration
156 for (int i=1 ; i<nbins ; i++) {
157 _ay[i] += _ay[i-1] ;
158 }
159
160 // Normalize and transfer to cache histogram
161 double binv = (_self->x.max()-_self->x.min())/nbins ;
162 for (int i=0 ; i<nbins ; i++) {
163 hist()->get(i) ;
164 if (cdfmode) {
165 hist()->set(i, _ay[i]/_ay[nbins-1], 0.);
166 } else {
167 hist()->set(i, _ay[i]*binv, 0.);
168 }
169 }
170
171 if (cdfmode) {
172 func()->setCdfBoundaries(true) ;
173 }
174 _self->x = xsave ;
175}
176
177
178
179////////////////////////////////////////////////////////////////////////////////
180/// Fill all empty histogram bins in the range [ixlo,ixhi] where nbins is the
181/// total number of histogram bins. This method samples the mid-point of the
182/// range and if the mid-point value is within small tolerance of the interpolated
183/// mid-point value fills all remaining elements through linear interpolation.
184/// If the tolerance is exceeded, the algorithm is recursed on the two subranges
185/// [xlo,xmid] and [xmid,xhi]
186
188{
189 // Add first and last point, if not there already
190 if (_ay[ixlo]<0) {
191 addPoint(ixlo) ;
192 }
193 if (_ay[ixhi]<0) {
194 addPoint(ixhi) ;
195 }
196
197 // Terminate here if there is no gap
198 if (ixhi-ixlo==1) {
199 return ;
200 }
201
202 // If gap size is one, simply fill gap and return
203 if (ixhi-ixlo==2) {
204 addPoint(ixlo+1) ;
205 return ;
206 }
207
208 // Add mid-point
209 Int_t ixmid = (ixlo+ixhi)/2 ;
210 addPoint(ixmid) ;
211
212 // Calculate difference of mid-point w.r.t interpolated value
213 double yInt = _ay[ixlo] + (_ay[ixhi]-_ay[ixlo])*(ixmid-ixlo)/(ixhi-ixlo) ;
214
215 // If relative deviation is greater than tolerance divide and iterate
216 if (fabs(yInt-_ay[ixmid])*(_ax[nbins-1]-_ax[0])>1e-6) {
217 addRange(ixlo,ixmid,nbins) ;
218 addRange(ixmid,ixhi,nbins) ;
219 } else {
220 for (Int_t j=ixlo+1 ; j<ixmid ; j++) {
221 _ay[j] = _ay[ixlo] + (_ay[ixmid]-_ay[ixlo])*(j-ixlo)/(ixmid-ixlo) ;
222 }
223 for (Int_t j=ixmid+1 ; j<ixhi ; j++) {
224 _ay[j] = _ay[ixmid] + (_ay[ixhi]-_ay[ixmid])*(j-ixmid)/(ixhi-ixmid) ;
225 }
226 }
227
228}
229
230
231////////////////////////////////////////////////////////////////////////////////
232/// Sample function at bin ix
233
235{
236 hist()->get(ix) ;
237 _self->x = _xx->getVal() ;
238 _ay[ix] = _self->func.arg().getVal(*_xx) ;
239
240}
241
242
243
244////////////////////////////////////////////////////////////////////////////////
245/// Fill the cache object by calling its calculate() method
246
248{
249 RICacheElem& riCache = static_cast<RICacheElem&>(cache) ;
250 riCache.calculate(false) ;
251}
252
253
254
255////////////////////////////////////////////////////////////////////////////////
256/// Return observable in nset to be cached by RooAbsCachedPdf
257/// this is always the x observable that is integrated
258
260{
261 RooArgSet* ret = new RooArgSet ;
262 ret->add(x.arg()) ;
263 return ret ;
264}
265
266
267
268////////////////////////////////////////////////////////////////////////////////
269/// Return the parameters of the cache created by RooAbsCachedPdf.
270/// These are always the input functions parameter, but never the
271/// integrated variable x.
272
274{
276 ret->remove(x.arg(),true,true) ;
277 return ret ;
278}
279
280
281////////////////////////////////////////////////////////////////////////////////
282/// Create custom cache element for running integral calculations
283
285{
286 return new RICacheElem(*const_cast<RooNumRunningInt*>(this),nset) ;
287}
288
289
290////////////////////////////////////////////////////////////////////////////////
291/// Dummy function that is never called
292
294{
295 cout << "RooNumRunningInt::evaluate(" << GetName() << ")" << endl ;
296 return 0 ;
297}
298
299
#define e(i)
Definition: RSha256.hxx:103
#define ClassImp(name)
Definition: Rtypes.h:375
char name[80]
Definition: TGX11.cxx:110
#define hi
Definition: THbookFile.cxx:128
RooArgSet * getParameters(const RooAbsData *data, bool stripDisconnected=true) 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:541
RooArgList containedArgs(Action) override
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 remove(const RooAbsArg &var, bool silent=false, bool matchByNameOnly=false)
Remove the specified argument from our list.
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to 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:62
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:91
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgList.h:22
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:56
Int_t numEntries() const override
Return the number of bins.
const RooArgSet * get() const override
Get bin centre of current bin.
Definition: RooDataHist.h:82
std::vector< double > _ay
std::vector< double > _ax
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...
RooArgList containedArgs(Action) override
Return all RooAbsArg components contained in cache element.
void addPoint(Int_t ix)
Sample function at bin ix.
void calculate(bool 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...
double evaluate() const override
Dummy function that is never called.
const char * inputBaseName() const override
Return unique name for RooAbsCachedPdf cache components constructed from input function name.
RooRealProxy func
Proxy to functions whose running integral is calculated.
FuncCacheElem * createCache(const RooArgSet *nset) const override
Create custom cache element for running integral calculations.
~RooNumRunningInt() override
Destructor.
RooArgSet * actualParameters(const RooArgSet &nset) const override
Return the parameters of the cache created by RooAbsCachedPdf.
friend class RICacheElem
RooArgSet * actualObservables(const RooArgSet &nset) const override
Return observable in nset to be cached by RooAbsCachedPdf this is always the x observable that is int...
RooRealProxy x
Intergrated observable.
void fillCacheObject(FuncCacheElem &cacheFunc) const override
Fill the cache object by calling its calculate() method.
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:40
const T & arg() const
Return reference to object held in proxy.
const char * GetName() const override
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)