ROOT   Reference Guide
RooNumConvolution.cxx
Go to the documentation of this file.
1/*****************************************************************************
2 * Project: RooFit *
3 * Package: RooFitCore *
4 * @(#)root/roofitcore:$Id$
5 * Authors: *
6 * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu *
7 * DK, David Kirkby, UC Irvine, dkirkby@uci.edu *
8 * *
9 * Copyright (c) 2000-2005, Regents of the University of California *
11 * *
12 * Redistribution and use in source and binary forms, *
13 * with or without modification, are permitted according to the terms *
14 * listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
15 *****************************************************************************/
16
17/**
18\file RooNumConvolution.cxx
19\class RooNumConvolution
20\ingroup Roofitcore
21
22Numeric 1-dimensional convolution operator PDF. This class can convolve any PDF
23with any other PDF
24This class should not be used blindly as numeric convolution is computing
25intensive and prone to stability fitting problems. If an analytic convolution
26can be calculated, you should use that or implement it if not available.
27RooNumConvolution implements reasonable defaults that should convolve most
28functions reasonably well, but results strongly depend on the shape of your
29input PDFS so always check your result.
30
31The default integration engine for the numeric convolution is the
32adaptive Gauss-Kronrod method, which empirically seems the most robust
33for this task. You can override the convolution integration settings via
34the RooNumIntConfig object reference returned by the convIntConfig() member
35function
36By default the numeric convolution is integrated from -infinity to
37+infinity through a <pre>x -> 1/x</pre> coordinate transformation of the
38tails. For convolution with a very small bandwidth it may be
39advantageous (for both CPU consumption and stability) if the
40integration domain is limited to a finite range. The function
41setConvolutionWindow(mean,width,scale) allows to set a sliding
42window around the x value to be calculated taking a RooAbsReal
43expression for an offset and a width to be taken around the x
44value. These input expression can be RooFormulaVars or other
45function objects although the 3d 'scale' argument 'scale'
46multiplies the width RooAbsReal expression given in the 2nd
47argument, allowing for an appropriate window definition for most
48cases without need for a RooFormulaVar object: e.g. a Gaussian
49resolution PDF do setConvolutionWindow(gaussMean,gaussSigma,5)
50Note that for a 'wide' Gaussian the -inf to +inf integration
51may converge more quickly than that over a finite range!
52The default numeric precision is 1e-7, i.e. the global default for
53numeric integration but you should experiment with this value to
54see if it is sufficient for example by studying the number of function
55calls that MINUIT needs to fit your function as function of the
56convolution precision.
57**/
58
59#include "RooFit.h"
60
61#include "Riostream.h"
62#include "TH2F.h"
63#include "RooNumConvolution.h"
64#include "RooArgList.h"
65#include "RooRealVar.h"
66#include "RooFormulaVar.h"
67#include "RooCustomizer.h"
69#include "RooNumIntFactory.h"
70#include "RooGenContext.h"
71#include "RooConvGenContext.h"
72#include "RooMsgService.h"
73
74
75using namespace std;
76
78
79
80
81////////////////////////////////////////////////////////////////////////////////
82
84 _init(kFALSE),
85 _integrand(0),
86 _integrator(0),
87 _cloneVar(0),
88 _clonePdf(0),
89 _cloneModel(0),
90 _useWindow(kFALSE),
91 _windowScale(1),
92 _verboseThresh(2000),
93 _doProf(kFALSE),
94 _callHist(0)
95{
96}
97
98
99
100////////////////////////////////////////////////////////////////////////////////
101/// Constructor of convolution operator PDF
102///
103/// convVar : convolution variable (on which both pdf and resmodel should depend)
104/// pdf : input 'physics' pdf
105/// resmodel : input 'resultion' pdf
106///
107/// output is pdf(x) (X) resmodel(x) = Int [ pdf(x') resmodel (x-x') ] dx'
108///
109
110RooNumConvolution::RooNumConvolution(const char *name, const char *title, RooRealVar& convVar, RooAbsReal& inPdf, RooAbsReal& resmodel, const RooNumConvolution* proto) :
111 RooAbsReal(name,title),
112 _init(kFALSE),
113 _convIntConfig(RooNumIntConfig::defaultConfig()),
114 _integrand(0),
115 _integrator(0),
116 _origVar("origVar","Original Convolution variable",this,convVar),
117 _origPdf("origPdf","Original Input PDF",this,inPdf),
118 _origModel("origModel","Original Resolution model",this,resmodel),
119 _ownedClonedPdfSet("ownedClonePdfSet"),
120 _ownedClonedModelSet("ownedCloneModelSet"),
121 _cloneVar(0),
122 _clonePdf(0),
123 _cloneModel(0),
124 _useWindow(kFALSE),
125 _windowScale(1),
126 _windowParam("windowParam","Convolution window parameter",this,kFALSE),
127 _verboseThresh(2000),
128 _doProf(kFALSE),
129 _callHist(0)
130{
131 // Use Adaptive Gauss-Kronrod integration by default for the convolution integral
134
135 if (proto) {
136 convIntConfig() = proto->convIntConfig() ;
137 if (proto->_useWindow) {
138 setConvolutionWindow((RooAbsReal&)*proto->_windowParam.at(0),(RooAbsReal&)*proto->_windowParam.at(1),proto->_windowScale) ;
139 }
140 }
141}
142
143
144
145////////////////////////////////////////////////////////////////////////////////
146/// Copy constructor
147
149 RooAbsReal(other,name),
150 _init(kFALSE),
151 _convIntConfig(other._convIntConfig),
152 _integrand(0),
153 _integrator(0),
154 _origVar("origVar",this,other._origVar),
155 _origPdf("origPdf",this,other._origPdf),
156 _origModel("origModel",this,other._origModel),
157 _ownedClonedPdfSet("ownedClonePdfSet"),
158 _ownedClonedModelSet("ownedCloneModelSet"),
159 _cloneVar(0),
160 _clonePdf(0),
161 _cloneModel(0),
162 _useWindow(other._useWindow),
163 _windowScale(other._windowScale),
164 _windowParam("windowParam",this,other._windowParam),
165 _verboseThresh(other._verboseThresh),
166 _doProf(other._doProf),
167 _callHist(other._callHist)
168{
169}
170
171
172
173////////////////////////////////////////////////////////////////////////////////
174/// One-time initialization of object
175
177{
178 // Initialization function -- create clone of convVar (x') and deep-copy clones of pdf and
179 // model that are connected to x' rather than x (convVar)
180
181 // Start out clean
184
185 if (_cloneVar) delete _cloneVar ;
186
187 // Customize a copy of origPdf that is connected to x' rather than x
188 // store all cloned components in _clonePdfSet as well as x' itself
189 _cloneVar = new RooRealVar(Form("%s_prime",_origVar.arg().GetName()),"Convolution Variable",0) ;
190
191 RooCustomizer mgr1(pdf(),"NumConv_PdfClone") ;
193 mgr1.replaceArg(var(),*_cloneVar) ;
194 _clonePdf = (RooAbsReal*) mgr1.build() ;
195
196 RooCustomizer mgr2(model(),"NumConv_ModelClone") ;
198 mgr2.replaceArg(var(),*_cloneVar) ;
199 _cloneModel = (RooAbsReal*) mgr2.build() ;
200
201 // Change name back to original name
203
204 // Create Convolution integrand
206
207 // Instantiate integrator for convolution integrand
210
211 _init = kTRUE ;
212}
213
214
215
216
217////////////////////////////////////////////////////////////////////////////////
218/// Destructor
219
221{
222}
223
224
225
226////////////////////////////////////////////////////////////////////////////////
227/// Calculate convolution integral
228
230{
231 // Check if deferred initialization has occurred
232 if (!_init) initialize() ;
233
234 // Retrieve current value of convolution variable
236
237 // Propagate current normalization set to integrand
239
240 // Adjust convolution integration window
241 if (_useWindow) {
242 Double_t center = ((RooAbsReal*)_windowParam.at(0))->getVal() ;
244 _integrator->setLimits(x-center-width,x-center+width) ;
245 } else {
247 }
248
249 // Calculate convolution for present x
251 Double_t ret = _integrator->integral(&x) ;
252 if (_doProf) {
255 coutW(Integration) << "RooNumConvolution::eveluate(" << GetName() << ") WARNING convolution integral at x=" << x
256 << " required " << _integrand->numCall() << " function evaluations" << endl ;
257 }
258 }
259
260 return ret ;
261}
262
263
264
265////////////////////////////////////////////////////////////////////////////////
266/// Intercept server redirects. Throw away cache, as figuring out redirections on the cache is an unsolvable problem.
267
268Bool_t RooNumConvolution::redirectServersHook(const RooAbsCollection& /*newServerList*/, Bool_t /*mustReplaceAll*/,
269 Bool_t /*nameChange*/, Bool_t /*isRecursive*/)
270{
271 _init = kFALSE ;
272 return kFALSE ;
273}
274
275
276
277////////////////////////////////////////////////////////////////////////////////
278/// Removes previously defined convolution window, reverting to convolution from -inf to +inf
279
281{
284}
285
286
287
288////////////////////////////////////////////////////////////////////////////////
289/// Restrict convolution integral to finite range [ x - C - S*W, x - C + S*W ]
290/// where x is current value of convolution variablem, C = centerParam, W=widthParam and S = widthScaleFactor
291/// Inputs centerParam and withParam can be function expressions (RooAbsReal, RooFormulaVar) etc.
292
293void RooNumConvolution::setConvolutionWindow(RooAbsReal& centerParam, RooAbsReal& widthParam, Double_t widthScaleFactor)
294{
295 _useWindow = kTRUE ;
299 _windowScale = widthScaleFactor ;
300}
301
302
303
304////////////////////////////////////////////////////////////////////////////////
305/// Activate warning messages if number of function calls needed for evaluation of convolution integral
306/// exceeds given threshold
307
309{
310 if (threshold<0) {
311 coutE(InputArguments) << "RooNumConvolution::setCallWarning(" << GetName() << ") ERROR: threshold must be positive, value unchanged" << endl ;
312 return ;
313 }
314 _verboseThresh = threshold ;
315}
316
317
318////////////////////////////////////////////////////////////////////////////////
319/// Activate call profile if flag is set to true. A 2-D histogram is kept that stores the required number
320/// of function calls versus the value of x, the convolution variable
321///
322/// All clones of RooNumConvolution objects will keep logging to the histogram of the original class
323/// so that performance of temporary object clones, such as used in e.g. fitting, plotting and generating
324/// are all logged in a single place.
325///
326/// Function caller should take ownership of profiling histogram as it is not deleted at the RooNumConvolution dtor
327///
328/// Calling this function with flag set to false will deactivate call profiling and delete the profiling histogram
329
330void RooNumConvolution::setCallProfiling(Bool_t flag, Int_t nbinX, Int_t nbinCall, Int_t nCallHigh)
331{
332 if (flag) {
333 if (_doProf) {
334 delete _callHist ;
335 }
336 _callHist = new TH2F(Form("callHist_%s",GetName()),Form("Call Profiling of RooNumConvolution %s",GetTitle()),
337 nbinX,_origVar.min(),_origVar.max(),
338 nbinCall,0,nCallHigh) ;
339 _doProf=kTRUE ;
340
341 } else if (_doProf) {
342
343 delete _callHist ;
344 _callHist = 0 ;
345 _doProf = kFALSE ;
346 }
347
348}
349
350
351
352////////////////////////////////////////////////////////////////////////////////
353/// Hook function to intercept printCompactTree() calls so that it can print out
354/// the content of its private cache in the print sequence
355
356void RooNumConvolution::printCompactTreeHook(ostream& os, const char* indent)
357{
358 os << indent << "RooNumConvolution begin cache" << endl ;
359
360 if (_init) {
361 _cloneVar->printCompactTree(os,Form("%s[Var]",indent)) ;
362 _clonePdf->printCompactTree(os,Form("%s[Pdf]",indent)) ;
363 _cloneModel->printCompactTree(os,Form("%s[Mod]",indent)) ;
364 }
365
366 os << indent << "RooNumConvolution end cache" << endl ;
367}
368
369
#define coutW(a)
Definition: RooMsgService.h:32
#define coutE(a)
Definition: RooMsgService.h:33
const Bool_t kFALSE
Definition: RtypesCore.h:90
const Bool_t kTRUE
Definition: RtypesCore.h:89
#define ClassImp(name)
Definition: Rtypes.h:361
static void indent(ostringstream &buf, int indent_level)
include TDocParser_001 C image html pict1_TDocParser_001 png width
Definition: TDocParser.cxx:121
char name[80]
Definition: TGX11.cxx:109
char * Form(const char *fmt,...)
const char * proto
Definition: civetweb.c:16604
void printCompactTree(const char *indent="", const char *fileName=0, const char *namePat=0, RooAbsArg *client=0)
Print tree structure of expression tree on stdout, or to file if filename is specified.
Definition: RooAbsArg.cxx:1745
void SetName(const char *name)
Set the name of the TNamed.
Definition: RooAbsArg.cxx:2216
RooAbsCollection is an abstract container object that can hold multiple RooAbsArg objects.
virtual void removeAll()
Remove all arguments from our set, deleting them if we own them.
void resetNumCall() const
Definition: RooAbsFunc.h:46
Int_t numCall() const
Definition: RooAbsFunc.h:42
virtual Double_t integral(const Double_t *yvec=0)=0
virtual Bool_t setUseIntegrandLimits(Bool_t flag)
Interface function that allows to defer limit definition to integrand definition.
virtual Bool_t setLimits(Double_t *, Double_t *)
const RooArgSet * nset() const
Definition: RooAbsProxy.h:46
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
RooAbsArg * at(Int_t idx) const
Return object at given index, or nullptr if index is out of range.
Definition: RooArgList.h:74
virtual Bool_t setLabel(const char *label, bool printError=true) override
Set value by specifying the name of the desired state.
Implementation of RooAbsFunc that represent the the integrand of a generic (numeric) convolution A (x...
void setNormalizationSet(const RooArgSet *nset)
RooCustomizer is a factory class to produce clones of a prototype composite PDF object with the same ...
Definition: RooCustomizer.h:32
void setCloneBranchSet(RooArgSet &cloneBranchSet)
Install the input RooArgSet as container in which all cloned branches will be stored.
void replaceArg(const RooAbsArg &orig, const RooAbsArg &subst)
Replace any occurence of arg 'orig' with arg 'subst'.
RooAbsArg * build(const char *masterCatState, Bool_t verbose=kFALSE)
Build a clone of the prototype executing all registered 'replace' rules and 'split' rules for the mas...
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Reimplementation of standard RooArgList::add()
virtual void removeAll()
Reimplementation of standard RooArgList::removeAll()
Numeric 1-dimensional convolution operator PDF.
void clearConvolutionWindow()
Removes previously defined convolution window, reverting to convolution from -inf to +inf.
RooArgSet _ownedClonedModelSet
RooAbsReal * _cloneModel
Bool_t redirectServersHook(const RooAbsCollection &newServerList, Bool_t mustReplaceAll, Bool_t nameChange, Bool_t isRecursive)
Intercept server redirects. Throw away cache, as figuring out redirections on the cache is an unsolva...
void setConvolutionWindow(RooAbsReal &centerParam, RooAbsReal &widthParam, Double_t widthScaleFactor=1)
Restrict convolution integral to finite range [ x - C - S*W, x - C + S*W ] where x is current value o...
Double_t evaluate() const
Calculate convolution integral.
void initialize() const
One-time initialization of object.
RooRealProxy _origVar
Numeric integrator of convolution integrand.
virtual void printCompactTreeHook(std::ostream &os, const char *indent="")
Hook function to intercept printCompactTree() calls so that it can print out the content of its priva...
RooConvIntegrandBinding * _integrand
RooAbsIntegrator * _integrator
Binding of Convolution Integrand function.
void setCallProfiling(Bool_t flag, Int_t nbinX=40, Int_t nbinCall=40, Int_t nCallHigh=1000)
Activate call profile if flag is set to true.
RooListProxy _windowParam
RooArgSet _ownedClonedPdfSet
RooRealVar & var() const
RooAbsReal * _cloneVar
RooAbsReal * _clonePdf
RooNumIntConfig & convIntConfig()
virtual ~RooNumConvolution()
Destructor.
void setCallWarning(Int_t threshold=2000)
Activate warning messages if number of function calls needed for evaluation of convolution integral e...
RooNumIntConfig _convIntConfig
RooAbsReal & model() const
RooAbsReal & pdf() const
RooNumIntConfig holds the configuration parameters of the various numeric integrators used by RooReal...
RooCategory & method1D()
RooCategory & method1DOpen()
RooAbsIntegrator * createIntegrator(RooAbsFunc &func, const RooNumIntConfig &config, Int_t ndim=0, Bool_t isBinned=kFALSE) const
Construct a numeric integrator instance that operates on function 'func' and is configured with 'conf...
static RooNumIntFactory & instance()
Static method returning reference to singleton instance of factory.
static Double_t infinity()
Return internal infinity representation.
Definition: RooNumber.cxx:49
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:35
double min(const char *rname=0) const
Query lower limit of range. This requires the payload to be RooAbsRealLValue or derived.
double max(const char *rname=0) const
Query upper limit of range. This requires the payload to be RooAbsRealLValue or derived.
const T & arg() const
Return reference to object held in proxy.
2-D histogram with a float per channel (see TH1 documentation)}
Definition: TH2.h:251
Int_t Fill(Double_t)
Invalid Fill method.
Definition: TH2.cxx:294
virtual const char * GetTitle() const
Returns title of object.
Definition: TNamed.h:48
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
Double_t x[n]
Definition: legend1.C:17
@ InputArguments
Definition: RooGlobalFunc.h:68
@ Integration
Definition: RooGlobalFunc.h:67