Logo 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 *
10 * and Stanford University. All rights reserved. *
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 "Riostream.h"
60#include "TH2F.h"
61#include "RooNumConvolution.h"
62#include "RooArgList.h"
63#include "RooRealVar.h"
64#include "RooCustomizer.h"
66#include "RooNumIntFactory.h"
67#include "RooGenContext.h"
68#include "RooConvGenContext.h"
69#include "RooMsgService.h"
70
71
72using namespace std;
73
75
76
77
78////////////////////////////////////////////////////////////////////////////////
79
81 _init(false),
82 _integrand(0),
83 _integrator(0),
84 _cloneVar(0),
85 _clonePdf(0),
86 _cloneModel(0),
87 _useWindow(false),
88 _windowScale(1),
89 _verboseThresh(2000),
90 _doProf(false),
91 _callHist(0)
92{
93}
94
95
96
97////////////////////////////////////////////////////////////////////////////////
98/// Constructor of convolution operator PDF
99///
100/// convVar : convolution variable (on which both pdf and resmodel should depend)
101/// pdf : input 'physics' pdf
102/// resmodel : input 'resultion' pdf
103///
104/// output is pdf(x) (X) resmodel(x) = Int [ pdf(x') resmodel (x-x') ] dx'
105///
106
107RooNumConvolution::RooNumConvolution(const char *name, const char *title, RooRealVar& convVar, RooAbsReal& inPdf, RooAbsReal& resmodel, const RooNumConvolution* proto) :
108 RooAbsReal(name,title),
109 _init(false),
110 _convIntConfig(RooNumIntConfig::defaultConfig()),
111 _integrand(0),
112 _integrator(0),
113 _origVar("origVar","Original Convolution variable",this,convVar),
114 _origPdf("origPdf","Original Input PDF",this,inPdf),
115 _origModel("origModel","Original Resolution model",this,resmodel),
116 _ownedClonedPdfSet("ownedClonePdfSet"),
117 _ownedClonedModelSet("ownedCloneModelSet"),
118 _cloneVar(0),
119 _clonePdf(0),
120 _cloneModel(0),
121 _useWindow(false),
122 _windowScale(1),
123 _windowParam("windowParam","Convolution window parameter",this,false),
124 _verboseThresh(2000),
125 _doProf(false),
126 _callHist(0)
127{
128 // Use Adaptive Gauss-Kronrod integration by default for the convolution integral
129 _convIntConfig.method1D().setLabel("RooAdaptiveGaussKronrodIntegrator1D") ;
130 _convIntConfig.method1DOpen().setLabel("RooAdaptiveGaussKronrodIntegrator1D") ;
131
132 if (proto) {
133 convIntConfig() = proto->convIntConfig() ;
134 if (proto->_useWindow) {
135 setConvolutionWindow((RooAbsReal&)*proto->_windowParam.at(0),(RooAbsReal&)*proto->_windowParam.at(1),proto->_windowScale) ;
136 }
137 }
138}
139
140
141
142////////////////////////////////////////////////////////////////////////////////
143/// Copy constructor
144
146 RooAbsReal(other,name),
147 _init(false),
148 _convIntConfig(other._convIntConfig),
149 _integrand(0),
150 _integrator(0),
151 _origVar("origVar",this,other._origVar),
152 _origPdf("origPdf",this,other._origPdf),
153 _origModel("origModel",this,other._origModel),
154 _ownedClonedPdfSet("ownedClonePdfSet"),
155 _ownedClonedModelSet("ownedCloneModelSet"),
156 _cloneVar(0),
157 _clonePdf(0),
158 _cloneModel(0),
159 _useWindow(other._useWindow),
160 _windowScale(other._windowScale),
161 _windowParam("windowParam",this,other._windowParam),
162 _verboseThresh(other._verboseThresh),
163 _doProf(other._doProf),
164 _callHist(other._callHist)
165{
166}
167
168
169
170////////////////////////////////////////////////////////////////////////////////
171/// One-time initialization of object
172
174{
175 // Initialization function -- create clone of convVar (x') and deep-copy clones of pdf and
176 // model that are connected to x' rather than x (convVar)
177
178 // Start out clean
181
182 if (_cloneVar) delete _cloneVar ;
183
184 // Customize a copy of origPdf that is connected to x' rather than x
185 // store all cloned components in _clonePdfSet as well as x' itself
186 _cloneVar = new RooRealVar(Form("%s_prime",_origVar.arg().GetName()),"Convolution Variable",0) ;
187
188 RooCustomizer mgr1(pdf(),"NumConv_PdfClone") ;
190 mgr1.replaceArg(var(),*_cloneVar) ;
191 _clonePdf = (RooAbsReal*) mgr1.build() ;
192
193 RooCustomizer mgr2(model(),"NumConv_ModelClone") ;
195 mgr2.replaceArg(var(),*_cloneVar) ;
196 _cloneModel = (RooAbsReal*) mgr2.build() ;
197
198 // Change name back to original name
200
201 // Create Convolution integrand
203
204 // Instantiate integrator for convolution integrand
207
208 _init = true ;
209}
210
211
212
213
214////////////////////////////////////////////////////////////////////////////////
215/// Destructor
216
218{
219}
220
221
222
223////////////////////////////////////////////////////////////////////////////////
224/// Calculate convolution integral
225
227{
228 // Check if deferred initialization has occurred
229 if (!_init) initialize() ;
230
231 // Retrieve current value of convolution variable
232 double x = _origVar ;
233
234 // Propagate current normalization set to integrand
236
237 // Adjust convolution integration window
238 if (_useWindow) {
239 double center = ((RooAbsReal*)_windowParam.at(0))->getVal() ;
240 double width = _windowScale * ((RooAbsReal*)_windowParam.at(1))->getVal() ;
241 _integrator->setLimits(x-center-width,x-center+width) ;
242 } else {
244 }
245
246 // Calculate convolution for present x
248 double ret = _integrator->integral(&x) ;
249 if (_doProf) {
252 coutW(Integration) << "RooNumConvolution::eveluate(" << GetName() << ") WARNING convolution integral at x=" << x
253 << " required " << _integrand->numCall() << " function evaluations" << endl ;
254 }
255 }
256
257 return ret ;
258}
259
260
261
262////////////////////////////////////////////////////////////////////////////////
263/// Intercept server redirects. Throw away cache, as figuring out redirections on the cache is an unsolvable problem.
264
265bool RooNumConvolution::redirectServersHook(const RooAbsCollection& newServerList, bool mustReplaceAll,
266 bool nameChange, bool isRecursive)
267{
268 _init = false ;
269 return RooAbsReal::redirectServersHook(newServerList, mustReplaceAll, nameChange, isRecursive);
270}
271
272
273
274////////////////////////////////////////////////////////////////////////////////
275/// Removes previously defined convolution window, reverting to convolution from -inf to +inf
276
278{
279 _useWindow = false ;
281}
282
283
284
285////////////////////////////////////////////////////////////////////////////////
286/// Restrict convolution integral to finite range [ x - C - S*W, x - C + S*W ]
287/// where x is current value of convolution variablem, C = centerParam, W=widthParam and S = widthScaleFactor
288/// Inputs centerParam and withParam can be function expressions (RooAbsReal, RooFormulaVar) etc.
289
290void RooNumConvolution::setConvolutionWindow(RooAbsReal& centerParam, RooAbsReal& widthParam, double widthScaleFactor)
291{
292 _useWindow = true ;
294 _windowParam.add(centerParam) ;
295 _windowParam.add(widthParam) ;
296 _windowScale = widthScaleFactor ;
297}
298
299
300
301////////////////////////////////////////////////////////////////////////////////
302/// Activate warning messages if number of function calls needed for evaluation of convolution integral
303/// exceeds given threshold
304
306{
307 if (threshold<0) {
308 coutE(InputArguments) << "RooNumConvolution::setCallWarning(" << GetName() << ") ERROR: threshold must be positive, value unchanged" << endl ;
309 return ;
310 }
311 _verboseThresh = threshold ;
312}
313
314
315////////////////////////////////////////////////////////////////////////////////
316/// Activate call profile if flag is set to true. A 2-D histogram is kept that stores the required number
317/// of function calls versus the value of x, the convolution variable
318///
319/// All clones of RooNumConvolution objects will keep logging to the histogram of the original class
320/// so that performance of temporary object clones, such as used in e.g. fitting, plotting and generating
321/// are all logged in a single place.
322///
323/// Function caller should take ownership of profiling histogram as it is not deleted at the RooNumConvolution dtor
324///
325/// Calling this function with flag set to false will deactivate call profiling and delete the profiling histogram
326
327void RooNumConvolution::setCallProfiling(bool flag, Int_t nbinX, Int_t nbinCall, Int_t nCallHigh)
328{
329 if (flag) {
330 if (_doProf) {
331 delete _callHist ;
332 }
333 _callHist = new TH2F(Form("callHist_%s",GetName()),Form("Call Profiling of RooNumConvolution %s",GetTitle()),
334 nbinX,_origVar.min(),_origVar.max(),
335 nbinCall,0,nCallHigh) ;
336 _doProf=true ;
337
338 } else if (_doProf) {
339
340 delete _callHist ;
341 _callHist = 0 ;
342 _doProf = false ;
343 }
344
345}
346
347
348
349////////////////////////////////////////////////////////////////////////////////
350/// Hook function to intercept printCompactTree() calls so that it can print out
351/// the content of its private cache in the print sequence
352
354{
355 os << indent << "RooNumConvolution begin cache" << endl ;
356
357 if (_init) {
358 _cloneVar->printCompactTree(os,Form("%s[Var]",indent)) ;
359 _clonePdf->printCompactTree(os,Form("%s[Pdf]",indent)) ;
360 _cloneModel->printCompactTree(os,Form("%s[Mod]",indent)) ;
361 }
362
363 os << indent << "RooNumConvolution end cache" << endl ;
364}
365
366
#define coutW(a)
Definition: RooMsgService.h:36
#define coutE(a)
Definition: RooMsgService.h:37
#define ClassImp(name)
Definition: Rtypes.h:375
static void indent(ostringstream &buf, int indent_level)
Option_t Option_t width
char name[80]
Definition: TGX11.cxx:110
char * Form(const char *fmt,...)
Formats a string in a circular formatting buffer.
Definition: TString.cxx:2452
const char * proto
Definition: civetweb.c:17502
void SetName(const char *name) override
Set the name of the TNamed.
Definition: RooAbsArg.cxx:2314
void printCompactTree(const char *indent="", const char *fileName=nullptr, const char *namePat=nullptr, RooAbsArg *client=nullptr)
Print tree structure of expression tree on stdout, or to file if filename is specified.
Definition: RooAbsArg.cxx:1893
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
Reset function call counter.
Definition: RooAbsFunc.h:52
Int_t numCall() const
Return number of function calls since last reset.
Definition: RooAbsFunc.h:47
virtual bool setUseIntegrandLimits(bool flag)
Interface function that allows to defer limit definition to integrand definition.
virtual double integral(const double *yvec=nullptr)=0
virtual bool setLimits(double *, double *)
const RooArgSet * nset() const
Definition: RooAbsProxy.h:48
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
bool redirectServersHook(const RooAbsCollection &newServerList, bool mustReplaceAll, bool nameChange, bool isRecursiveStep) override
A buffer for reading values from trees.
RooAbsArg * at(Int_t idx) const
Return object at given index, or nullptr if index is out of range.
Definition: RooArgList.h:110
bool setLabel(const char *label, bool printError=true) override
Set value by specifying the name of the desired state.
void removeAll() override
Remove all argument inset using remove(const RooAbsArg&).
bool add(const RooAbsArg &var, bool valueServer, bool shapeServer, bool silent)
Overloaded RooCollection_t::add() method insert object into set and registers object as server to own...
Implementation of RooAbsFunc that represent the integrand of a generic (numeric) convolution A (x) B ...
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:35
void setCloneBranchSet(RooArgSet &cloneBranchSet)
Releases ownership of list of cloned branch nodes.
void replaceArg(const RooAbsArg &orig, const RooAbsArg &subst)
Replace any occurence of arg 'orig' with arg 'subst'.
RooAbsArg * build(const char *masterCatState, bool verbose=false)
Build a clone of the prototype executing all registered 'replace' rules and 'split' rules for the mas...
Numeric 1-dimensional convolution operator PDF.
void clearConvolutionWindow()
Removes previously defined convolution window, reverting to convolution from -inf to +inf.
RooArgSet _ownedClonedModelSet
Owning set of cloned model components.
RooAbsReal * _cloneModel
Pointer to cloned model.
void setCallProfiling(bool flag, Int_t nbinX=40, Int_t nbinCall=40, Int_t nCallHigh=1000)
Activate call profile if flag is set to true.
~RooNumConvolution() override
Destructor.
bool redirectServersHook(const RooAbsCollection &newServerList, bool mustReplaceAll, bool nameChange, bool isRecursive) override
Intercept server redirects. Throw away cache, as figuring out redirections on the cache is an unsolva...
void initialize() const
One-time initialization of object.
RooRealProxy _origVar
Original convolution variable.
Int_t _verboseThresh
Call count threshold for verbose printing.
RooConvIntegrandBinding * _integrand
! Binding of Convolution Integrand function
RooAbsIntegrator * _integrator
! Numeric integrator of convolution integrand
void printCompactTreeHook(std::ostream &os, const char *indent="") override
Hook function to intercept printCompactTree() calls so that it can print out the content of its priva...
double evaluate() const override
Calculate convolution integral.
RooListProxy _windowParam
Holder for optional convolution integration window scaling parameter.
double _windowScale
Scale factor for window parameter.
RooArgSet _ownedClonedPdfSet
Owning set of cloned PDF components.
RooRealVar & var() const
bool _doProf
Switch to activate profiling option.
RooAbsReal * _cloneVar
Pointer to cloned convolution variable.
RooAbsReal * _clonePdf
Pointer to cloned PDF.
RooNumIntConfig & convIntConfig()
void setCallWarning(Int_t threshold=2000)
Activate warning messages if number of function calls needed for evaluation of convolution integral e...
RooNumIntConfig _convIntConfig
Configuration of numeric convolution integral ;.
RooAbsReal & model() const
RooAbsReal & pdf() const
TH2 * _callHist
! Histogram recording number of calls per convolution integral calculation
bool _useWindow
Switch to activate window convolution.
void setConvolutionWindow(RooAbsReal &centerParam, RooAbsReal &widthParam, double widthScaleFactor=1)
Restrict convolution integral to finite range [ x - C - S*W, x - C + S*W ] where x is current value o...
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 isBinned=false) 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 infinity()
Return internal infinity representation.
Definition: RooNumber.cxx:48
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:40
double max(const char *rname=nullptr) 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.
double min(const char *rname=nullptr) const
Query lower limit of range. This requires the payload to be RooAbsRealLValue or derived.
2-D histogram with a float per channel (see TH1 documentation)}
Definition: TH2.h:257
Int_t Fill(Double_t) override
Invalid Fill method.
Definition: TH2.cxx:347
const char * GetName() const override
Returns name of object.
Definition: TNamed.h:47
const char * GetTitle() const override
Returns title of object.
Definition: TNamed.h:48
Double_t x[n]
Definition: legend1.C:17
@ InputArguments
Definition: RooGlobalFunc.h:62
@ Integration
Definition: RooGlobalFunc.h:61