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 std::cout, std::endl, std::string;
39
41
42
43
44////////////////////////////////////////////////////////////////////////////////
45/// Construct running integral of function '_func' over x_print from
46/// the lower bound on _x to the present value of _x using a numeric
47/// sampling technique. The sampling frequency is controlled by the
48/// binning named 'bname' and a default second order interpolation
49/// is applied to smooth the histogram-based c.d.f.
50
51RooNumRunningInt::RooNumRunningInt(const char *name, const char *title, RooAbsReal& _func, RooRealVar& _x, const char* bname) :
53 func("func","func",this,_func),
54 x("x","x",this,_x),
55 _binningName(bname?bname:"cache")
56 {
58 }
59
60
61
62
63////////////////////////////////////////////////////////////////////////////////
64/// Copy constructor
65
68 func("func",this,other.func),
69 x("x",this,other.x),
70 _binningName(other._binningName)
71 {
72 }
73
74
75
76////////////////////////////////////////////////////////////////////////////////
77/// Destructor
78
80{
81}
82
83
84////////////////////////////////////////////////////////////////////////////////
85/// Return unique name for RooAbsCachedPdf cache components
86/// constructed from input function name
87
89{
90 static string ret ;
91 ret = func.arg().GetName() ;
92 ret += "_NUMRUNINT" ;
93 return ret.c_str() ;
94} ;
95
96
97
98////////////////////////////////////////////////////////////////////////////////
99/// Construct RunningIntegral CacheElement
100
102 : FuncCacheElem(self, nset),
103 _self(&const_cast<RooNumRunningInt &>(self)),
104 _xx(static_cast<RooRealVar *>(hist()->get()->find(self.x.arg().GetName())))
105{
106 // Instantiate temp arrays
107 _ax.resize(hist()->numEntries());
108 _ay.resize(hist()->numEntries());
109
110 // Copy X values from histo
111
112 for (int i=0 ; i<hist()->numEntries() ; i++) {
113 hist()->get(i) ;
114 _ax[i] = _xx->getVal() ;
115 _ay[i] = -1 ;
116 }
117
118}
119
120
121////////////////////////////////////////////////////////////////////////////////
122/// Return all RooAbsArg components contained in cache element
123
125{
126 RooArgList ret ;
127 ret.add(FuncCacheElem::containedArgs(action)) ;
128 ret.add(*_self) ;
129 ret.add(*_xx) ;
130 return ret ;
131}
132
133
134
135////////////////////////////////////////////////////////////////////////////////
136/// Calculate the numeric running integral and store
137/// the result in the cache histogram provided
138/// by RooAbsCachedPdf
139
141{
142 // Update contents of histogram
143 Int_t nbins = hist()->numEntries() ;
144
145 double xsave = _self->x ;
146
147 Int_t lastHi=0 ;
148 Int_t nInitRange=32 ;
149 for (int i=1 ; i<=nInitRange ; i++) {
150 Int_t hi = (i*nbins)/nInitRange -1 ;
151 Int_t lo = lastHi ;
152 addRange(lo,hi,nbins) ;
153 lastHi=hi ;
154 }
155
156 // Perform numeric integration
157 for (int i=1 ; i<nbins ; i++) {
158 _ay[i] += _ay[i-1] ;
159 }
160
161 // Normalize and transfer to cache histogram
162 double binv = (_self->x.max()-_self->x.min())/nbins ;
163 for (int i=0 ; i<nbins ; i++) {
164 hist()->get(i) ;
165 if (cdfmode) {
166 hist()->set(i, _ay[i]/_ay[nbins-1], 0.);
167 } else {
168 hist()->set(i, _ay[i]*binv, 0.);
169 }
170 }
171
172 if (cdfmode) {
173 func()->setCdfBoundaries(true) ;
174 }
175 _self->x = xsave ;
176}
177
178
179
180////////////////////////////////////////////////////////////////////////////////
181/// Fill all empty histogram bins in the range [ixlo,ixhi] where nbins is the
182/// total number of histogram bins. This method samples the mid-point of the
183/// range and if the mid-point value is within small tolerance of the interpolated
184/// mid-point value fills all remaining elements through linear interpolation.
185/// If the tolerance is exceeded, the algorithm is recursed on the two subranges
186/// [xlo,xmid] and [xmid,xhi]
187
189{
190 // Add first and last point, if not there already
191 if (_ay[ixlo]<0) {
192 addPoint(ixlo) ;
193 }
194 if (_ay[ixhi]<0) {
195 addPoint(ixhi) ;
196 }
197
198 // Terminate here if there is no gap
199 if (ixhi-ixlo==1) {
200 return ;
201 }
202
203 // If gap size is one, simply fill gap and return
204 if (ixhi-ixlo==2) {
205 addPoint(ixlo+1) ;
206 return ;
207 }
208
209 // Add mid-point
210 Int_t ixmid = (ixlo+ixhi)/2 ;
211 addPoint(ixmid) ;
212
213 // Calculate difference of mid-point w.r.t interpolated value
214 double yInt = _ay[ixlo] + (_ay[ixhi]-_ay[ixlo])*(ixmid-ixlo)/(ixhi-ixlo) ;
215
216 // If relative deviation is greater than tolerance divide and iterate
217 if (std::abs(yInt-_ay[ixmid])*(_ax[nbins-1]-_ax[0])>1e-6) {
218 addRange(ixlo,ixmid,nbins) ;
219 addRange(ixmid,ixhi,nbins) ;
220 } else {
221 for (Int_t j=ixlo+1 ; j<ixmid ; j++) {
222 _ay[j] = _ay[ixlo] + (_ay[ixmid]-_ay[ixlo])*(j-ixlo)/(ixmid-ixlo) ;
223 }
224 for (Int_t j=ixmid+1 ; j<ixhi ; j++) {
225 _ay[j] = _ay[ixmid] + (_ay[ixhi]-_ay[ixmid])*(j-ixmid)/(ixhi-ixmid) ;
226 }
227 }
228
229}
230
231
232////////////////////////////////////////////////////////////////////////////////
233/// Sample function at bin ix
234
236{
237 hist()->get(ix) ;
238 _self->x = _xx->getVal() ;
239 _ay[ix] = _self->func.arg().getVal(*_xx) ;
240
241}
242
243
244
245////////////////////////////////////////////////////////////////////////////////
246/// Fill the cache object by calling its calculate() method
247
249{
250 RICacheElem& riCache = static_cast<RICacheElem&>(cache) ;
251 riCache.calculate(false) ;
252}
253
254
255
256////////////////////////////////////////////////////////////////////////////////
257/// Return observable in nset to be cached by RooAbsCachedPdf
258/// this is always the x observable that is integrated
259
261{
262 RooArgSet* ret = new RooArgSet ;
263 ret->add(x.arg()) ;
265}
266
267
268
269////////////////////////////////////////////////////////////////////////////////
270/// Return the parameters of the cache created by RooAbsCachedPdf.
271/// These are always the input functions parameter, but never the
272/// integrated variable x.
273
275{
276 auto ret = func->getParameters(RooArgSet()) ;
277 ret->remove(x.arg(),true,true) ;
278 return ret;
279}
280
281
282////////////////////////////////////////////////////////////////////////////////
283/// Create custom cache element for running integral calculations
284
286{
287 return new RICacheElem(*const_cast<RooNumRunningInt*>(this),nset) ;
288}
289
290
291////////////////////////////////////////////////////////////////////////////////
292/// Dummy function that is never called
293
295{
296 cout << "RooNumRunningInt::evaluate(" << GetName() << ")" << endl ;
297 return 0 ;
298}
299
300
#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