Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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
17Implementation 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),
104 _self(&const_cast<RooNumRunningInt &>(self)),
105 _xx(static_cast<RooRealVar *>(hist()->get()->find(self.x.arg().GetName())))
106{
107 // Instantiate temp arrays
108 _ax.resize(hist()->numEntries());
109 _ay.resize(hist()->numEntries());
110
111 // Copy X values from histo
112
113 for (int i=0 ; i<hist()->numEntries() ; i++) {
114 hist()->get(i) ;
115 _ax[i] = _xx->getVal() ;
116 _ay[i] = -1 ;
117 }
118
119}
120
121
122////////////////////////////////////////////////////////////////////////////////
123/// Return all RooAbsArg components contained in cache element
124
126{
127 RooArgList ret ;
128 ret.add(FuncCacheElem::containedArgs(action)) ;
129 ret.add(*_self) ;
130 ret.add(*_xx) ;
131 return ret ;
132}
133
134
135
136////////////////////////////////////////////////////////////////////////////////
137/// Calculate the numeric running integral and store
138/// the result in the cache histogram provided
139/// by RooAbsCachedPdf
140
142{
143 // Update contents of histogram
144 Int_t nbins = hist()->numEntries() ;
145
146 double xsave = _self->x ;
147
148 Int_t lastHi=0 ;
149 Int_t nInitRange=32 ;
150 for (int i=1 ; i<=nInitRange ; i++) {
151 Int_t hi = (i*nbins)/nInitRange -1 ;
152 Int_t lo = lastHi ;
153 addRange(lo,hi,nbins) ;
154 lastHi=hi ;
155 }
156
157 // Perform numeric integration
158 for (int i=1 ; i<nbins ; i++) {
159 _ay[i] += _ay[i-1] ;
160 }
161
162 // Normalize and transfer to cache histogram
163 double binv = (_self->x.max()-_self->x.min())/nbins ;
164 for (int i=0 ; i<nbins ; i++) {
165 hist()->get(i) ;
166 if (cdfmode) {
167 hist()->set(i, _ay[i]/_ay[nbins-1], 0.);
168 } else {
169 hist()->set(i, _ay[i]*binv, 0.);
170 }
171 }
172
173 if (cdfmode) {
174 func()->setCdfBoundaries(true) ;
175 }
176 _self->x = xsave ;
177}
178
179
180
181////////////////////////////////////////////////////////////////////////////////
182/// Fill all empty histogram bins in the range [ixlo,ixhi] where nbins is the
183/// total number of histogram bins. This method samples the mid-point of the
184/// range and if the mid-point value is within small tolerance of the interpolated
185/// mid-point value fills all remaining elements through linear interpolation.
186/// If the tolerance is exceeded, the algorithm is recursed on the two subranges
187/// [xlo,xmid] and [xmid,xhi]
188
190{
191 // Add first and last point, if not there already
192 if (_ay[ixlo]<0) {
193 addPoint(ixlo) ;
194 }
195 if (_ay[ixhi]<0) {
196 addPoint(ixhi) ;
197 }
198
199 // Terminate here if there is no gap
200 if (ixhi-ixlo==1) {
201 return ;
202 }
203
204 // If gap size is one, simply fill gap and return
205 if (ixhi-ixlo==2) {
206 addPoint(ixlo+1) ;
207 return ;
208 }
209
210 // Add mid-point
211 Int_t ixmid = (ixlo+ixhi)/2 ;
212 addPoint(ixmid) ;
213
214 // Calculate difference of mid-point w.r.t interpolated value
215 double yInt = _ay[ixlo] + (_ay[ixhi]-_ay[ixlo])*(ixmid-ixlo)/(ixhi-ixlo) ;
216
217 // If relative deviation is greater than tolerance divide and iterate
218 if (std::abs(yInt-_ay[ixmid])*(_ax[nbins-1]-_ax[0])>1e-6) {
219 addRange(ixlo,ixmid,nbins) ;
220 addRange(ixmid,ixhi,nbins) ;
221 } else {
222 for (Int_t j=ixlo+1 ; j<ixmid ; j++) {
223 _ay[j] = _ay[ixlo] + (_ay[ixmid]-_ay[ixlo])*(j-ixlo)/(ixmid-ixlo) ;
224 }
225 for (Int_t j=ixmid+1 ; j<ixhi ; j++) {
226 _ay[j] = _ay[ixmid] + (_ay[ixhi]-_ay[ixmid])*(j-ixmid)/(ixhi-ixmid) ;
227 }
228 }
229
230}
231
232
233////////////////////////////////////////////////////////////////////////////////
234/// Sample function at bin ix
235
237{
238 hist()->get(ix) ;
239 _self->x = _xx->getVal() ;
240 _ay[ix] = _self->func.arg().getVal(*_xx) ;
241
242}
243
244
245
246////////////////////////////////////////////////////////////////////////////////
247/// Fill the cache object by calling its calculate() method
248
250{
251 RICacheElem& riCache = static_cast<RICacheElem&>(cache) ;
252 riCache.calculate(false) ;
253}
254
255
256
257////////////////////////////////////////////////////////////////////////////////
258/// Return observable in nset to be cached by RooAbsCachedPdf
259/// this is always the x observable that is integrated
260
262{
263 RooArgSet* ret = new RooArgSet ;
264 ret->add(x.arg()) ;
266}
267
268
269
270////////////////////////////////////////////////////////////////////////////////
271/// Return the parameters of the cache created by RooAbsCachedPdf.
272/// These are always the input functions parameter, but never the
273/// integrated variable x.
274
276{
277 auto ret = func->getParameters(RooArgSet()) ;
278 ret->remove(x.arg(),true,true) ;
279 return ret;
280}
281
282
283////////////////////////////////////////////////////////////////////////////////
284/// Create custom cache element for running integral calculations
285
287{
288 return new RICacheElem(*const_cast<RooNumRunningInt*>(this),nset) ;
289}
290
291
292////////////////////////////////////////////////////////////////////////////////
293/// Dummy function that is never called
294
296{
297 cout << "RooNumRunningInt::evaluate(" << GetName() << ")" << endl ;
298 return 0 ;
299}
300
301
#define e(i)
Definition RSha256.hxx:103
RooAbsReal * _func
Pointer to original input function.
#define ClassImp(name)
Definition Rtypes.h:377
char name[80]
Definition TGX11.cxx:110
#define hi
RooFit::OwningPtr< 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...
Abstract base class for functions that need or want to cache their evaluate() output in a RooHistFunc...
void setInterpolationOrder(Int_t order)
Set interpolation order of RooHistFunct representing cache histogram.
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to list.
virtual Int_t numEntries() const
Return number of entries in dataset, i.e., count unweighted entries.
Abstract base class for objects that represent a real value and implements functionality common to al...
Definition RooAbsReal.h:59
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition RooAbsReal.h:103
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:55
const RooArgSet * get() const override
Get bin centre of current bin.
Definition RooDataHist.h:76
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.
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.
RooFit::OwningPtr< RooArgSet > actualParameters(const RooArgSet &nset) const override
Return the parameters of the cache created by RooAbsCachedPdf.
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.
RooFit::OwningPtr< 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...
~RooNumRunningInt() override
Destructor.
RooRealProxy x
Integrated observable.
void fillCacheObject(FuncCacheElem &cacheFunc) const override
Fill the cache object by calling its calculate() method.
Variable that can be changed from the outside.
Definition RooRealVar.h:37
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
T * OwningPtr
An alias for raw pointers for indicating that the return type of a RooFit function is an owning point...
Definition Config.h:35