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 
22 Numeric 1-dimensional convolution operator PDF. This class can convolve any PDF
23 with any other PDF
24 This class should not be used blindly as numeric convolution is computing
25 intensive and prone to stability fitting problems. If an analytic convolution
26 can be calculated, you should use that or implement it if not available.
27 RooNumConvolution implements reasonable defaults that should convolve most
28 functions reasonably well, but results strongly depend on the shape of your
29 input PDFS so always check your result.
30 
31 The default integration engine for the numeric convolution is the
32 adaptive Gauss-Kronrod method, which empirically seems the most robust
33 for this task. You can override the convolution integration settings via
34 the RooNumIntConfig object reference returned by the convIntConfig() member
35 function
36 By default the numeric convolution is integrated from -infinity to
37 +infinity through a <pre>x -> 1/x</pre> coordinate transformation of the
38 tails. For convolution with a very small bandwidth it may be
39 advantageous (for both CPU consumption and stability) if the
40 integration domain is limited to a finite range. The function
41 setConvolutionWindow(mean,width,scale) allows to set a sliding
42 window around the x value to be calculated taking a RooAbsReal
43 expression for an offset and a width to be taken around the x
44 value. These input expression can be RooFormulaVars or other
45 function objects although the 3d 'scale' argument 'scale'
46 multiplies the width RooAbsReal expression given in the 2nd
47 argument, allowing for an appropriate window definition for most
48 cases without need for a RooFormulaVar object: e.g. a Gaussian
49 resolution PDF do setConvolutionWindow(gaussMean,gaussSigma,5)
50 Note that for a 'wide' Gaussian the -inf to +inf integration
51 may converge more quickly than that over a finite range!
52 The default numeric precision is 1e-7, i.e. the global default for
53 numeric integration but you should experiment with this value to
54 see if it is sufficient for example by studying the number of function
55 calls that MINUIT needs to fit your function as function of the
56 convolution 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 
75 using 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 
110 RooNumConvolution::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
132  _convIntConfig.method1D().setLabel("RooAdaptiveGaussKronrodIntegrator1D") ;
133  _convIntConfig.method1DOpen().setLabel("RooAdaptiveGaussKronrodIntegrator1D") ;
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
202  _cloneVar->SetName(var().GetName()) ;
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
235  Double_t x = _origVar ;
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 
268 Bool_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 {
282  _useWindow = kFALSE ;
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 
293 void RooNumConvolution::setConvolutionWindow(RooAbsReal& centerParam, RooAbsReal& widthParam, Double_t widthScaleFactor)
294 {
295  _useWindow = kTRUE ;
297  _windowParam.add(centerParam) ;
298  _windowParam.add(widthParam) ;
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 
330 void 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 
356 void 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 
RooFormulaVar.h
RooNumConvolution::setCallProfiling
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.
Definition: RooNumConvolution.cxx:330
RooNumConvolution::_useWindow
Bool_t _useWindow
Definition: RooNumConvolution.h:89
RooNumConvolution::_convIntConfig
RooNumIntConfig _convIntConfig
Definition: RooNumConvolution.h:69
kTRUE
const Bool_t kTRUE
Definition: RtypesCore.h:91
RooMsgService.h
RooCustomizer::setCloneBranchSet
void setCloneBranchSet(RooArgSet &cloneBranchSet)
Install the input RooArgSet as container in which all cloned branches will be stored.
Definition: RooCustomizer.cxx:678
RooFit.h
RooTemplateProxy::arg
const T & arg() const
Return reference to object held in proxy.
Definition: RooTemplateProxy.h:271
RooFit::InputArguments
@ InputArguments
Definition: RooGlobalFunc.h:68
TH2F
2-D histogram with a float per channel (see TH1 documentation)}
Definition: TH2.h:251
ClassImp
#define ClassImp(name)
Definition: Rtypes.h:364
Form
char * Form(const char *fmt,...)
TNamed::GetTitle
virtual const char * GetTitle() const
Returns title of object.
Definition: TNamed.h:54
RooCustomizer
Definition: RooCustomizer.h:35
RooNumIntConfig::method1D
RooCategory & method1D()
Definition: RooNumIntConfig.h:34
coutE
#define coutE(a)
Definition: RooMsgService.h:33
coutW
#define coutW(a)
Definition: RooMsgService.h:32
RooNumConvolution::convIntConfig
RooNumIntConfig & convIntConfig()
Definition: RooNumConvolution.h:44
RooAbsReal::getVal
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:91
RooNumConvolution::_cloneModel
RooAbsReal * _cloneModel
Definition: RooNumConvolution.h:82
RooAbsArg::SetName
void SetName(const char *name)
Set the name of the TNamed.
Definition: RooAbsArg.cxx:2218
RooNumConvolution::initialize
void initialize() const
One-time initialization of object.
Definition: RooNumConvolution.cxx:176
width
include TDocParser_001 C image html pict1_TDocParser_001 png width
Definition: TDocParser.cxx:121
RooNumConvolution::~RooNumConvolution
virtual ~RooNumConvolution()
Destructor.
Definition: RooNumConvolution.cxx:220
x
Double_t x[n]
Definition: legend1.C:17
indent
static void indent(ostringstream &buf, int indent_level)
Definition: TClingCallFunc.cxx:87
RooNumIntFactory::createIntegrator
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...
Definition: RooNumIntFactory.cxx:172
RooAbsArg::printCompactTree
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:1747
RooAbsReal
Definition: RooAbsReal.h:61
RooNumConvolution::_callHist
TH2 * _callHist
Definition: RooNumConvolution.h:95
RooArgList::at
RooAbsArg * at(Int_t idx) const
Return object at given index, or nullptr if index is out of range.
Definition: RooArgList.h:88
RooCustomizer::build
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...
Definition: RooCustomizer.cxx:395
bool
RooCustomizer::replaceArg
void replaceArg(const RooAbsArg &orig, const RooAbsArg &subst)
Replace any occurence of arg 'orig' with arg 'subst'.
Definition: RooCustomizer.cxx:339
RooConvIntegrandBinding::setNormalizationSet
void setNormalizationSet(const RooArgSet *nset)
Definition: RooConvIntegrandBinding.h:35
RooNumConvolution::evaluate
Double_t evaluate() const
Calculate convolution integral.
Definition: RooNumConvolution.cxx:229
RooCategory::setLabel
virtual Bool_t setLabel(const char *label, bool printError=true) override
Set value by specifying the name of the desired state.
Definition: RooCategory.cxx:185
RooNumConvolution
Definition: RooNumConvolution.h:29
RooCustomizer.h
RooAbsIntegrator::setLimits
virtual Bool_t setLimits(Double_t *, Double_t *)
Definition: RooAbsIntegrator.h:75
RooNumConvolution::_clonePdf
RooAbsReal * _clonePdf
Definition: RooNumConvolution.h:81
RooNumConvolution::pdf
RooAbsReal & pdf() const
Definition: RooNumConvolution.h:56
RooTemplateProxy::min
double min(const char *rname=0) const
Query lower limit of range. This requires the payload to be RooAbsRealLValue or derived.
Definition: RooTemplateProxy.h:287
kFALSE
const Bool_t kFALSE
Definition: RtypesCore.h:92
RooNumConvolution::setCallWarning
void setCallWarning(Int_t threshold=2000)
Activate warning messages if number of function calls needed for evaluation of convolution integral e...
Definition: RooNumConvolution.cxx:308
RooAbsFunc::resetNumCall
void resetNumCall() const
Definition: RooAbsFunc.h:46
RooNumConvolution::_integrand
RooConvIntegrandBinding * _integrand
Definition: RooNumConvolution.h:70
RooAbsCollection
Definition: RooAbsCollection.h:30
RooGenContext.h
RooListProxy::add
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Reimplementation of standard RooArgList::add()
Definition: RooListProxy.cxx:104
RooListProxy::removeAll
virtual void removeAll()
Reimplementation of standard RooArgList::removeAll()
Definition: RooListProxy.cxx:157
RooAbsProxy::nset
const RooArgSet * nset() const
Definition: RooAbsProxy.h:59
RooNumConvolution::model
RooAbsReal & model() const
Definition: RooNumConvolution.h:57
RooConvIntegrandBinding
Definition: RooConvIntegrandBinding.h:25
RooRealVar.h
TH2F.h
TH2::Fill
Int_t Fill(Double_t)
Invalid Fill method.
Definition: TH2.cxx:294
RooNumConvolution::var
RooRealVar & var() const
Definition: RooNumConvolution.h:55
proto
const char * proto
Definition: civetweb.c:16604
RooNumConvolution::redirectServersHook
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...
Definition: RooNumConvolution.cxx:268
RooNumber::infinity
static Double_t infinity()
Return internal infinity representation.
Definition: RooNumber.cxx:49
Double_t
double Double_t
Definition: RtypesCore.h:59
RooNumConvolution::_windowScale
Double_t _windowScale
Definition: RooNumConvolution.h:90
RooNumConvolution.h
RooNumConvolution::_ownedClonedModelSet
RooArgSet _ownedClonedModelSet
Definition: RooNumConvolution.h:78
RooNumConvolution::_init
Bool_t _init
Definition: RooNumConvolution.h:63
RooAbsIntegrator::integral
virtual Double_t integral(const Double_t *yvec=0)=0
name
char name[80]
Definition: TGX11.cxx:110
RooAbsFunc::numCall
Int_t numCall() const
Definition: RooAbsFunc.h:42
RooConvIntegrandBinding.h
RooNumConvolution::setConvolutionWindow
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...
Definition: RooNumConvolution.cxx:293
RooNumConvolution::clearConvolutionWindow
void clearConvolutionWindow()
Removes previously defined convolution window, reverting to convolution from -inf to +inf.
Definition: RooNumConvolution.cxx:280
RooTemplateProxy::max
double max(const char *rname=0) const
Query upper limit of range. This requires the payload to be RooAbsRealLValue or derived.
Definition: RooTemplateProxy.h:289
RooNumIntConfig
Definition: RooNumIntConfig.h:25
RooNumConvolution::RooNumConvolution
RooNumConvolution()
Definition: RooNumConvolution.cxx:83
RooNumConvolution::_verboseThresh
Int_t _verboseThresh
Definition: RooNumConvolution.h:93
RooNumConvolution::_ownedClonedPdfSet
RooArgSet _ownedClonedPdfSet
Definition: RooNumConvolution.h:77
RooAbsIntegrator::setUseIntegrandLimits
virtual Bool_t setUseIntegrandLimits(Bool_t flag)
Interface function that allows to defer limit definition to integrand definition.
Definition: RooAbsIntegrator.cxx:90
TNamed::GetName
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:53
RooNumIntFactory::instance
static RooNumIntFactory & instance()
Static method returning reference to singleton instance of factory.
Definition: RooNumIntFactory.cxx:96
RooNumConvolution::printCompactTreeHook
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...
Definition: RooNumConvolution.cxx:356
RooNumIntConfig::method1DOpen
RooCategory & method1DOpen()
Definition: RooNumIntConfig.h:42
RooFit::Integration
@ Integration
Definition: RooGlobalFunc.h:67
RooRealVar
Definition: RooRealVar.h:35
RooAbsCollection::removeAll
virtual void removeAll()
Remove all arguments from our set, deleting them if we own them.
Definition: RooAbsCollection.cxx:643
Riostream.h
RooArgList.h
RooConvGenContext.h
RooNumIntFactory.h
RooNumConvolution::_windowParam
RooListProxy _windowParam
Definition: RooNumConvolution.h:91
RooNumConvolution::_cloneVar
RooAbsReal * _cloneVar
Definition: RooNumConvolution.h:80
RooNumConvolution::_integrator
RooAbsIntegrator * _integrator
Binding of Convolution Integrand function.
Definition: RooNumConvolution.h:71
RooNumConvolution::_doProf
Bool_t _doProf
Definition: RooNumConvolution.h:94
RooNumConvolution::_origVar
RooRealProxy _origVar
Numeric integrator of convolution integrand.
Definition: RooNumConvolution.h:73
int