Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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 std::endl, std::ostream;
73
74
75
76
77////////////////////////////////////////////////////////////////////////////////
78
80 _init(false),
81 _integrand(nullptr),
82 _cloneVar(nullptr),
83 _clonePdf(nullptr),
84 _cloneModel(nullptr),
85 _useWindow(false),
86 _windowScale(1),
87 _verboseThresh(2000),
88 _doProf(false),
89 _callHist(nullptr)
90{
91}
92
93
94
95////////////////////////////////////////////////////////////////////////////////
96/// Constructor of convolution operator PDF
97///
98/// convVar : convolution variable (on which both pdf and resmodel should depend)
99/// pdf : input 'physics' pdf
100/// resmodel : input 'resolution' pdf
101///
102/// output is pdf(x) (X) resmodel(x) = Int [ pdf(x') resmodel (x-x') ] dx'
103///
104
106 RooAbsReal(name,title),
107 _init(false),
108 _convIntConfig(RooNumIntConfig::defaultConfig()),
109 _integrand(nullptr),
110 _origVar("origVar","Original Convolution variable",this,convVar),
111 _origPdf("origPdf","Original Input PDF",this,inPdf),
112 _origModel("origModel","Original Resolution model",this,resmodel),
113 _ownedClonedPdfSet("ownedClonePdfSet"),
114 _ownedClonedModelSet("ownedCloneModelSet"),
115 _cloneVar(nullptr),
116 _clonePdf(nullptr),
117 _cloneModel(nullptr),
118 _useWindow(false),
119 _windowScale(1),
120 _windowParam("windowParam","Convolution window parameter",this,false),
121 _verboseThresh(2000),
122 _doProf(false),
123 _callHist(nullptr)
124{
125 // Use Adaptive Gauss-Kronrod integration by default for the convolution integral
126 _convIntConfig.method1D().setLabel("RooAdaptiveGaussKronrodIntegrator1D") ;
127 _convIntConfig.method1DOpen().setLabel("RooAdaptiveGaussKronrodIntegrator1D") ;
128
129 if (proto) {
130 convIntConfig() = proto->convIntConfig() ;
131 if (proto->_useWindow) {
132 setConvolutionWindow(static_cast<RooAbsReal&>(*proto->_windowParam.at(0)),static_cast<RooAbsReal&>(*proto->_windowParam.at(1)),proto->_windowScale) ;
133 }
134 }
135}
136
137
138
139////////////////////////////////////////////////////////////////////////////////
140/// Copy constructor
141
144 _init(false),
145 _convIntConfig(other._convIntConfig),
146 _integrand(nullptr),
147 _origVar("origVar",this,other._origVar),
148 _origPdf("origPdf",this,other._origPdf),
149 _origModel("origModel",this,other._origModel),
150 _ownedClonedPdfSet("ownedClonePdfSet"),
151 _ownedClonedModelSet("ownedCloneModelSet"),
152 _cloneVar(nullptr),
153 _clonePdf(nullptr),
154 _cloneModel(nullptr),
155 _useWindow(other._useWindow),
156 _windowScale(other._windowScale),
157 _windowParam("windowParam",this,other._windowParam),
158 _verboseThresh(other._verboseThresh),
159 _doProf(other._doProf),
160 _callHist(other._callHist)
161{
162}
163
164
165
166////////////////////////////////////////////////////////////////////////////////
167/// One-time initialization of object
168
170{
171 // Initialization function -- create clone of convVar (x') and deep-copy clones of pdf and
172 // model that are connected to x' rather than x (convVar)
173
174 // Start out clean
177
178 if (_cloneVar) delete _cloneVar ;
179
180 // Customize a copy of origPdf that is connected to x' rather than x
181 // store all cloned components in _clonePdfSet as well as x' itself
182 _cloneVar = new RooRealVar(Form("%s_prime",_origVar.arg().GetName()),"Convolution Variable",0) ;
183
184 RooCustomizer mgr1(pdf(),"NumConv_PdfClone") ;
185 mgr1.setCloneBranchSet(_ownedClonedPdfSet) ;
186 mgr1.replaceArg(var(),*_cloneVar) ;
187 _clonePdf = static_cast<RooAbsReal*>(mgr1.build()) ;
188
189 RooCustomizer mgr2(model(),"NumConv_ModelClone") ;
190 mgr2.setCloneBranchSet(_ownedClonedModelSet) ;
191 mgr2.replaceArg(var(),*_cloneVar) ;
192 _cloneModel = static_cast<RooAbsReal*>(mgr2.build()) ;
193
194 // Change name back to original name
196
197 // Create Convolution integrand
199
200 // Instantiate integrator for convolution integrand
202 _integrator->setUseIntegrandLimits(false) ;
203
204 _init = true ;
205}
206
207
208
209
210////////////////////////////////////////////////////////////////////////////////
211/// Destructor
212
216
217
218
219////////////////////////////////////////////////////////////////////////////////
220/// Calculate convolution integral
221
223{
224 // Check if deferred initialization has occurred
225 if (!_init) initialize() ;
226
227 // Retrieve current value of convolution variable
228 double x = _origVar ;
229
230 // Propagate current normalization set to integrand
231 _integrand->setNormalizationSet(_origVar.nset()) ;
232
233 // Adjust convolution integration window
234 if (_useWindow) {
235 double center = (static_cast<RooAbsReal*>(_windowParam.at(0)))->getVal() ;
236 double width = _windowScale * (static_cast<RooAbsReal*>(_windowParam.at(1)))->getVal() ;
237 _integrator->setLimits(x-center-width,x-center+width) ;
238 } else {
240 }
241
242 // Calculate convolution for present x
243 if (_doProf) _integrand->resetNumCall() ;
244 double ret = _integrator->integral(&x) ;
245 if (_doProf) {
246 _callHist->Fill(x,_integrand->numCall()) ;
247 if (_integrand->numCall()>_verboseThresh) {
248 coutW(Integration) << "RooNumConvolution::evaluate(" << GetName() << ") WARNING convolution integral at x=" << x
249 << " required " << _integrand->numCall() << " function evaluations" << std::endl ;
250 }
251 }
252
253 return ret ;
254}
255
256
257
258////////////////////////////////////////////////////////////////////////////////
259/// Intercept server redirects. Throw away cache, as figuring out redirections on the cache is an unsolvable problem.
260
267
268
269
270////////////////////////////////////////////////////////////////////////////////
271/// Removes previously defined convolution window, reverting to convolution from -inf to +inf
272
278
279
280
281////////////////////////////////////////////////////////////////////////////////
282/// Restrict convolution integral to finite range [ x - C - S*W, x - C + S*W ]
283/// where x is current value of convolution variablem, C = centerParam, W=widthParam and S = widthScaleFactor
284/// Inputs centerParam and withParam can be function expressions (RooAbsReal, RooFormulaVar) etc.
285
294
295
296
297////////////////////////////////////////////////////////////////////////////////
298/// Activate warning messages if number of function calls needed for evaluation of convolution integral
299/// exceeds given threshold
300
302{
303 if (threshold<0) {
304 coutE(InputArguments) << "RooNumConvolution::setCallWarning(" << GetName() << ") ERROR: threshold must be positive, value unchanged" << std::endl ;
305 return ;
306 }
308}
309
310
311////////////////////////////////////////////////////////////////////////////////
312/// Activate call profile if flag is set to true. A 2-D histogram is kept that stores the required number
313/// of function calls versus the value of x, the convolution variable
314///
315/// All clones of RooNumConvolution objects will keep logging to the histogram of the original class
316/// so that performance of temporary object clones, such as used in e.g. fitting, plotting and generating
317/// are all logged in a single place.
318///
319/// Function caller should take ownership of profiling histogram as it is not deleted at the RooNumConvolution dtor
320///
321/// Calling this function with flag set to false will deactivate call profiling and delete the profiling histogram
322
324{
325 if (flag) {
326 if (_doProf) {
327 delete _callHist ;
328 }
329 _callHist = new TH2F(Form("callHist_%s",GetName()),Form("Call Profiling of RooNumConvolution %s",GetTitle()),
332 _doProf=true ;
333
334 } else if (_doProf) {
335
336 delete _callHist ;
337 _callHist = nullptr ;
338 _doProf = false ;
339 }
340
341}
342
343
344
345////////////////////////////////////////////////////////////////////////////////
346/// Hook function to intercept printCompactTree() calls so that it can print out
347/// the content of its private cache in the print sequence
348
350{
351 os << indent << "RooNumConvolution begin cache" << std::endl ;
352
353 if (_init) {
354 _cloneVar->printCompactTree(os,Form("%s[Var]",indent)) ;
355 _clonePdf->printCompactTree(os,Form("%s[Pdf]",indent)) ;
356 _cloneModel->printCompactTree(os,Form("%s[Mod]",indent)) ;
357 }
358
359 os << indent << "RooNumConvolution end cache" << std::endl ;
360}
361
362
#define coutW(a)
#define coutE(a)
static void indent(ostringstream &buf, int indent_level)
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
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:2489
const char * proto
Definition civetweb.c:17535
void SetName(const char *name) override
Set the name of the TNamed.
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.
Abstract container object that can hold multiple RooAbsArg objects.
virtual void removeAll()
Remove all arguments from our set, deleting them if we own them.
const RooArgSet * nset() const
Definition RooAbsProxy.h:52
Abstract base class for objects that represent a real value and implements functionality common to al...
Definition RooAbsReal.h:59
bool redirectServersHook(const RooAbsCollection &newServerList, bool mustReplaceAll, bool nameChange, bool isRecursiveStep) override
Function that is called at the end of redirectServers().
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...
RooCustomizer is a factory class to produce clones of a prototype composite PDF object with the same ...
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.
std::unique_ptr< RooAbsIntegrator > _integrator
! Numeric integrator of convolution integrand
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
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...
Holds the configuration parameters of the various numeric integrators used by RooRealIntegral.
RooCategory & method1D()
RooCategory & method1DOpen()
static RooNumIntFactory & instance()
Static method returning reference to singleton instance of factory.
static constexpr double infinity()
Return internal infinity representation.
Definition RooNumber.h:25
Variable that can be changed from the outside.
Definition RooRealVar.h:37
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:307
Int_t Fill(Double_t) override
Invalid Fill method.
Definition TH2.cxx:357
const char * GetName() const override
Returns name of object.
Definition TNamed.h:49
const char * GetTitle() const override
Returns title of object.
Definition TNamed.h:50
Double_t x[n]
Definition legend1.C:17