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 "RooFormulaVar.h"
65#include "RooCustomizer.h"
67#include "RooNumIntFactory.h"
68#include "RooGenContext.h"
69#include "RooConvGenContext.h"
70#include "RooMsgService.h"
71
72
73using namespace std;
74
76
77
78
79////////////////////////////////////////////////////////////////////////////////
80
82 _init(kFALSE),
83 _integrand(0),
84 _integrator(0),
85 _cloneVar(0),
86 _clonePdf(0),
87 _cloneModel(0),
88 _useWindow(kFALSE),
89 _windowScale(1),
90 _verboseThresh(2000),
91 _doProf(kFALSE),
92 _callHist(0)
93{
94}
95
96
97
98////////////////////////////////////////////////////////////////////////////////
99/// Constructor of convolution operator PDF
100///
101/// convVar : convolution variable (on which both pdf and resmodel should depend)
102/// pdf : input 'physics' pdf
103/// resmodel : input 'resultion' pdf
104///
105/// output is pdf(x) (X) resmodel(x) = Int [ pdf(x') resmodel (x-x') ] dx'
106///
107
108RooNumConvolution::RooNumConvolution(const char *name, const char *title, RooRealVar& convVar, RooAbsReal& inPdf, RooAbsReal& resmodel, const RooNumConvolution* proto) :
109 RooAbsReal(name,title),
110 _init(kFALSE),
111 _convIntConfig(RooNumIntConfig::defaultConfig()),
112 _integrand(0),
113 _integrator(0),
114 _origVar("origVar","Original Convolution variable",this,convVar),
115 _origPdf("origPdf","Original Input PDF",this,inPdf),
116 _origModel("origModel","Original Resolution model",this,resmodel),
117 _ownedClonedPdfSet("ownedClonePdfSet"),
118 _ownedClonedModelSet("ownedCloneModelSet"),
119 _cloneVar(0),
120 _clonePdf(0),
121 _cloneModel(0),
122 _useWindow(kFALSE),
123 _windowScale(1),
124 _windowParam("windowParam","Convolution window parameter",this,kFALSE),
125 _verboseThresh(2000),
126 _doProf(kFALSE),
127 _callHist(0)
128{
129 // Use Adaptive Gauss-Kronrod integration by default for the convolution integral
130 _convIntConfig.method1D().setLabel("RooAdaptiveGaussKronrodIntegrator1D") ;
131 _convIntConfig.method1DOpen().setLabel("RooAdaptiveGaussKronrodIntegrator1D") ;
132
133 if (proto) {
134 convIntConfig() = proto->convIntConfig() ;
135 if (proto->_useWindow) {
136 setConvolutionWindow((RooAbsReal&)*proto->_windowParam.at(0),(RooAbsReal&)*proto->_windowParam.at(1),proto->_windowScale) ;
137 }
138 }
139}
140
141
142
143////////////////////////////////////////////////////////////////////////////////
144/// Copy constructor
145
147 RooAbsReal(other,name),
148 _init(kFALSE),
149 _convIntConfig(other._convIntConfig),
150 _integrand(0),
151 _integrator(0),
152 _origVar("origVar",this,other._origVar),
153 _origPdf("origPdf",this,other._origPdf),
154 _origModel("origModel",this,other._origModel),
155 _ownedClonedPdfSet("ownedClonePdfSet"),
156 _ownedClonedModelSet("ownedCloneModelSet"),
157 _cloneVar(0),
158 _clonePdf(0),
159 _cloneModel(0),
160 _useWindow(other._useWindow),
161 _windowScale(other._windowScale),
162 _windowParam("windowParam",this,other._windowParam),
163 _verboseThresh(other._verboseThresh),
164 _doProf(other._doProf),
165 _callHist(other._callHist)
166{
167}
168
169
170
171////////////////////////////////////////////////////////////////////////////////
172/// One-time initialization of object
173
175{
176 // Initialization function -- create clone of convVar (x') and deep-copy clones of pdf and
177 // model that are connected to x' rather than x (convVar)
178
179 // Start out clean
182
183 if (_cloneVar) delete _cloneVar ;
184
185 // Customize a copy of origPdf that is connected to x' rather than x
186 // store all cloned components in _clonePdfSet as well as x' itself
187 _cloneVar = new RooRealVar(Form("%s_prime",_origVar.arg().GetName()),"Convolution Variable",0) ;
188
189 RooCustomizer mgr1(pdf(),"NumConv_PdfClone") ;
191 mgr1.replaceArg(var(),*_cloneVar) ;
192 _clonePdf = (RooAbsReal*) mgr1.build() ;
193
194 RooCustomizer mgr2(model(),"NumConv_ModelClone") ;
196 mgr2.replaceArg(var(),*_cloneVar) ;
197 _cloneModel = (RooAbsReal*) mgr2.build() ;
198
199 // Change name back to original name
201
202 // Create Convolution integrand
204
205 // Instantiate integrator for convolution integrand
208
209 _init = kTRUE ;
210}
211
212
213
214
215////////////////////////////////////////////////////////////////////////////////
216/// Destructor
217
219{
220}
221
222
223
224////////////////////////////////////////////////////////////////////////////////
225/// Calculate convolution integral
226
228{
229 // Check if deferred initialization has occurred
230 if (!_init) initialize() ;
231
232 // Retrieve current value of convolution variable
234
235 // Propagate current normalization set to integrand
237
238 // Adjust convolution integration window
239 if (_useWindow) {
240 Double_t center = ((RooAbsReal*)_windowParam.at(0))->getVal() ;
242 _integrator->setLimits(x-center-width,x-center+width) ;
243 } else {
245 }
246
247 // Calculate convolution for present x
249 Double_t ret = _integrator->integral(&x) ;
250 if (_doProf) {
253 coutW(Integration) << "RooNumConvolution::eveluate(" << GetName() << ") WARNING convolution integral at x=" << x
254 << " required " << _integrand->numCall() << " function evaluations" << endl ;
255 }
256 }
257
258 return ret ;
259}
260
261
262
263////////////////////////////////////////////////////////////////////////////////
264/// Intercept server redirects. Throw away cache, as figuring out redirections on the cache is an unsolvable problem.
265
266Bool_t RooNumConvolution::redirectServersHook(const RooAbsCollection& /*newServerList*/, Bool_t /*mustReplaceAll*/,
267 Bool_t /*nameChange*/, Bool_t /*isRecursive*/)
268{
269 _init = kFALSE ;
270 return kFALSE ;
271}
272
273
274
275////////////////////////////////////////////////////////////////////////////////
276/// Removes previously defined convolution window, reverting to convolution from -inf to +inf
277
279{
282}
283
284
285
286////////////////////////////////////////////////////////////////////////////////
287/// Restrict convolution integral to finite range [ x - C - S*W, x - C + S*W ]
288/// where x is current value of convolution variablem, C = centerParam, W=widthParam and S = widthScaleFactor
289/// Inputs centerParam and withParam can be function expressions (RooAbsReal, RooFormulaVar) etc.
290
291void RooNumConvolution::setConvolutionWindow(RooAbsReal& centerParam, RooAbsReal& widthParam, Double_t widthScaleFactor)
292{
293 _useWindow = kTRUE ;
295 _windowParam.add(centerParam) ;
296 _windowParam.add(widthParam) ;
297 _windowScale = widthScaleFactor ;
298}
299
300
301
302////////////////////////////////////////////////////////////////////////////////
303/// Activate warning messages if number of function calls needed for evaluation of convolution integral
304/// exceeds given threshold
305
307{
308 if (threshold<0) {
309 coutE(InputArguments) << "RooNumConvolution::setCallWarning(" << GetName() << ") ERROR: threshold must be positive, value unchanged" << endl ;
310 return ;
311 }
312 _verboseThresh = threshold ;
313}
314
315
316////////////////////////////////////////////////////////////////////////////////
317/// Activate call profile if flag is set to true. A 2-D histogram is kept that stores the required number
318/// of function calls versus the value of x, the convolution variable
319///
320/// All clones of RooNumConvolution objects will keep logging to the histogram of the original class
321/// so that performance of temporary object clones, such as used in e.g. fitting, plotting and generating
322/// are all logged in a single place.
323///
324/// Function caller should take ownership of profiling histogram as it is not deleted at the RooNumConvolution dtor
325///
326/// Calling this function with flag set to false will deactivate call profiling and delete the profiling histogram
327
328void RooNumConvolution::setCallProfiling(Bool_t flag, Int_t nbinX, Int_t nbinCall, Int_t nCallHigh)
329{
330 if (flag) {
331 if (_doProf) {
332 delete _callHist ;
333 }
334 _callHist = new TH2F(Form("callHist_%s",GetName()),Form("Call Profiling of RooNumConvolution %s",GetTitle()),
335 nbinX,_origVar.min(),_origVar.max(),
336 nbinCall,0,nCallHigh) ;
337 _doProf=kTRUE ;
338
339 } else if (_doProf) {
340
341 delete _callHist ;
342 _callHist = 0 ;
343 _doProf = kFALSE ;
344 }
345
346}
347
348
349
350////////////////////////////////////////////////////////////////////////////////
351/// Hook function to intercept printCompactTree() calls so that it can print out
352/// the content of its private cache in the print sequence
353
355{
356 os << indent << "RooNumConvolution begin cache" << endl ;
357
358 if (_init) {
359 _cloneVar->printCompactTree(os,Form("%s[Var]",indent)) ;
360 _clonePdf->printCompactTree(os,Form("%s[Pdf]",indent)) ;
361 _cloneModel->printCompactTree(os,Form("%s[Mod]",indent)) ;
362 }
363
364 os << indent << "RooNumConvolution end cache" << endl ;
365}
366
367
#define coutW(a)
Definition: RooMsgService.h:32
#define coutE(a)
Definition: RooMsgService.h:33
const Bool_t kFALSE
Definition: RtypesCore.h:101
const Bool_t kTRUE
Definition: RtypesCore.h:100
#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:2447
const char * proto
Definition: civetweb.c:17493
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:1903
void SetName(const char *name) override
Set the name of the TNamed.
Definition: RooAbsArg.cxx:2399
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 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:45
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:63
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:93
RooAbsArg * at(Int_t idx) const
Return object at given index, or nullptr if index is out of range.
Definition: RooArgList.h:110
Bool_t 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 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: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_t verbose=kFALSE)
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.
~RooNumConvolution() override
Destructor.
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...
void initialize() const
One-time initialization of object.
Double_t _windowScale
Scale factor for window parameter.
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...
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
Holder for optional convolution integration window scaling parameter.
RooArgSet _ownedClonedPdfSet
Owning set of cloned PDF components.
RooRealVar & var() const
Bool_t redirectServersHook(const RooAbsCollection &newServerList, Bool_t mustReplaceAll, Bool_t nameChange, Bool_t isRecursive) override
Intercept server redirects. Throw away cache, as figuring out redirections on the cache is an unsolva...
RooAbsReal * _cloneVar
Pointer to cloned convolution variable.
Bool_t _useWindow
Switch to activate window convolution.
Bool_t _doProf
Switch to activate profiling option.
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
Double_t evaluate() const override
Calculate convolution integral.
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:48
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:40
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:257
Int_t Fill(Double_t) override
Invalid Fill method.
Definition: TH2.cxx:358
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:64
@ Integration
Definition: RooGlobalFunc.h:63