Logo ROOT   6.08/07
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 "Riostream.h"
63 #include "TH2F.h"
64 #include "RooNumConvolution.h"
65 #include "RooArgList.h"
66 #include "RooRealVar.h"
67 #include "RooFormulaVar.h"
68 #include "RooCustomizer.h"
70 #include "RooNumIntFactory.h"
71 #include "RooGenContext.h"
72 #include "RooConvGenContext.h"
73 #include "RooMsgService.h"
74 
75 
76 using namespace std;
77 
79 ;
80 
81 
82 ////////////////////////////////////////////////////////////////////////////////
83 
85  _init(kFALSE),
86  _integrand(0),
87  _integrator(0),
88  _cloneVar(0),
89  _clonePdf(0),
90  _cloneModel(0),
91  _useWindow(kFALSE),
92  _windowScale(1),
93  _verboseThresh(2000),
94  _doProf(kFALSE),
95  _callHist(0)
96 {
97 }
98 
99 
100 
101 ////////////////////////////////////////////////////////////////////////////////
102 /// Constructor of convolution operator PDF
103 ///
104 /// convVar : convolution variable (on which both pdf and resmodel should depend)
105 /// pdf : input 'physics' pdf
106 /// resmodel : input 'resultion' pdf
107 ///
108 /// output is pdf(x) (X) resmodel(x) = Int [ pdf(x') resmodel (x-x') ] dx'
109 ///
110 
111 RooNumConvolution::RooNumConvolution(const char *name, const char *title, RooRealVar& convVar, RooAbsReal& inPdf, RooAbsReal& resmodel, const RooNumConvolution* proto) :
112  RooAbsReal(name,title),
113  _init(kFALSE),
114  _convIntConfig(RooNumIntConfig::defaultConfig()),
115  _integrand(0),
116  _integrator(0),
117  _origVar("origVar","Original Convolution variable",this,convVar),
118  _origPdf("origPdf","Original Input PDF",this,inPdf),
119  _origModel("origModel","Original Resolution model",this,resmodel),
120  _ownedClonedPdfSet("ownedClonePdfSet"),
121  _ownedClonedModelSet("ownedCloneModelSet"),
122  _cloneVar(0),
123  _clonePdf(0),
124  _cloneModel(0),
126  _windowScale(1),
127  _windowParam("windowParam","Convolution window parameter",this,kFALSE),
128  _verboseThresh(2000),
129  _doProf(kFALSE),
130  _callHist(0)
131 {
132  // Use Adaptive Gauss-Kronrod integration by default for the convolution integral
133  _convIntConfig.method1D().setLabel("RooAdaptiveGaussKronrodIntegrator1D") ;
134  _convIntConfig.method1DOpen().setLabel("RooAdaptiveGaussKronrodIntegrator1D") ;
135 
136  if (proto) {
137  convIntConfig() = proto->convIntConfig() ;
138  if (proto->_useWindow) {
140  }
141  }
142 }
143 
144 
145 
146 ////////////////////////////////////////////////////////////////////////////////
147 /// Copy constructor
148 
150  RooAbsReal(other,name),
151  _init(kFALSE),
153  _integrand(0),
154  _integrator(0),
155  _origVar("origVar",this,other._origVar),
156  _origPdf("origPdf",this,other._origPdf),
157  _origModel("origModel",this,other._origModel),
158  _ownedClonedPdfSet("ownedClonePdfSet"),
159  _ownedClonedModelSet("ownedCloneModelSet"),
160  _cloneVar(0),
161  _clonePdf(0),
162  _cloneModel(0),
163  _useWindow(other._useWindow),
164  _windowScale(other._windowScale),
165  _windowParam("windowParam",this,other._windowParam),
167  _doProf(other._doProf),
168  _callHist(other._callHist)
169 {
170 }
171 
172 
173 
174 ////////////////////////////////////////////////////////////////////////////////
175 /// One-time initialization of object
176 
178 {
179  // Initialization function -- create clone of convVar (x') and deep-copy clones of pdf and
180  // model that are connected to x' rather than x (convVar)
181 
182  // Start out clean
185 
186  if (_cloneVar) delete _cloneVar ;
187 
188  // Customize a copy of origPdf that is connected to x' rather than x
189  // store all cloned components in _clonePdfSet as well as x' itself
190  _cloneVar = new RooRealVar(Form("%s_prime",_origVar.arg().GetName()),"Convolution Variable",0) ;
191 
192  RooCustomizer mgr1(pdf(),"NumConv_PdfClone") ;
194  mgr1.replaceArg(var(),*_cloneVar) ;
195  _clonePdf = (RooAbsReal*) mgr1.build() ;
196 
197  RooCustomizer mgr2(model(),"NumConv_ModelClone") ;
199  mgr2.replaceArg(var(),*_cloneVar) ;
200  _cloneModel = (RooAbsReal*) mgr2.build() ;
201 
202  // Change name back to original name
203  _cloneVar->SetName(var().GetName()) ;
204 
205  // Create Convolution integrand
207 
208  // Instantiate integrator for convolution integrand
211 
212  _init = kTRUE ;
213 }
214 
215 
216 
217 
218 ////////////////////////////////////////////////////////////////////////////////
219 /// Destructor
220 
222 {
223 }
224 
225 
226 
227 ////////////////////////////////////////////////////////////////////////////////
228 /// Calculate convolution integral
229 
231 {
232  // Check if deferred initialization has occurred
233  if (!_init) initialize() ;
234 
235  // Retrieve current value of convolution variable
236  Double_t x = _origVar ;
237 
238  // Propagate current normalization set to integrand
240 
241  // Adjust convolution integration window
242  if (_useWindow) {
243  Double_t center = ((RooAbsReal*)_windowParam.at(0))->getVal() ;
244  Double_t width = _windowScale * ((RooAbsReal*)_windowParam.at(1))->getVal() ;
245  _integrator->setLimits(x-center-width,x-center+width) ;
246  } else {
248  }
249 
250  // Calculate convolution for present x
252  Double_t ret = _integrator->integral(&x) ;
253  if (_doProf) {
256  coutW(Integration) << "RooNumConvolution::eveluate(" << GetName() << ") WARNING convolution integral at x=" << x
257  << " required " << _integrand->numCall() << " function evaluations" << endl ;
258  }
259  }
260 
261  return ret ;
262 }
263 
264 
265 
266 ////////////////////////////////////////////////////////////////////////////////
267 /// Intercept server redirects. Throw away cache, as figuring out redirections on the cache is an unsolvable problem.
268 
269 Bool_t RooNumConvolution::redirectServersHook(const RooAbsCollection& /*newServerList*/, Bool_t /*mustReplaceAll*/,
270  Bool_t /*nameChange*/, Bool_t /*isRecursive*/)
271 {
272  _init = kFALSE ;
273  return kFALSE ;
274 }
275 
276 
277 
278 ////////////////////////////////////////////////////////////////////////////////
279 /// Removes previously defined convolution window, reverting to convolution from -inf to +inf
280 
282 {
283  _useWindow = kFALSE ;
285 }
286 
287 
288 
289 ////////////////////////////////////////////////////////////////////////////////
290 /// Restrict convolution integral to finite range [ x - C - S*W, x - C + S*W ]
291 /// where x is current value of convolution variablem, C = centerParam, W=widthParam and S = widthScaleFactor
292 /// Inputs centerParam and withParam can be function expressions (RooAbsReal, RooFormulaVar) etc.
293 
294 void RooNumConvolution::setConvolutionWindow(RooAbsReal& centerParam, RooAbsReal& widthParam, Double_t widthScaleFactor)
295 {
296  _useWindow = kTRUE ;
298  _windowParam.add(centerParam) ;
299  _windowParam.add(widthParam) ;
300  _windowScale = widthScaleFactor ;
301 }
302 
303 
304 
305 ////////////////////////////////////////////////////////////////////////////////
306 /// Activate warning messages if number of function calls needed for evaluation of convolution integral
307 /// exceeds given threshold
308 
310 {
311  if (threshold<0) {
312  coutE(InputArguments) << "RooNumConvolution::setCallWarning(" << GetName() << ") ERROR: threshold must be positive, value unchanged" << endl ;
313  return ;
314  }
315  _verboseThresh = threshold ;
316 }
317 
318 
319 ////////////////////////////////////////////////////////////////////////////////
320 /// Activate call profile if flag is set to true. A 2-D histogram is kept that stores the required number
321 /// of function calls versus the value of x, the convolution variable
322 ///
323 /// All clones of RooNumConvolution objects will keep logging to the histogram of the original class
324 /// so that performance of temporary object clones, such as used in e.g. fitting, plotting and generating
325 /// are all logged in a single place.
326 ///
327 /// Function caller should take ownership of profiling histogram as it is not deleted at the RooNumConvolution dtor
328 ///
329 /// Calling this function with flag set to false will deactivate call profiling and delete the profiling histogram
330 
331 void RooNumConvolution::setCallProfiling(Bool_t flag, Int_t nbinX, Int_t nbinCall, Int_t nCallHigh)
332 {
333  if (flag) {
334  if (_doProf) {
335  delete _callHist ;
336  }
337  _callHist = new TH2F(Form("callHist_%s",GetName()),Form("Call Profiling of RooNumConvolution %s",GetTitle()),
338  nbinX,_origVar.min(),_origVar.max(),
339  nbinCall,0,nCallHigh) ;
340  _doProf=kTRUE ;
341 
342  } else if (_doProf) {
343 
344  delete _callHist ;
345  _callHist = 0 ;
346  _doProf = kFALSE ;
347  }
348 
349 }
350 
351 
352 
353 ////////////////////////////////////////////////////////////////////////////////
354 /// Hook function to intercept printCompactTree() calls so that it can print out
355 /// the content of its private cache in the print sequence
356 
357 void RooNumConvolution::printCompactTreeHook(ostream& os, const char* indent)
358 {
359  os << indent << "RooNumConvolution begin cache" << endl ;
360 
361  if (_init) {
362  _cloneVar->printCompactTree(os,Form("%s[Var]",indent)) ;
363  _clonePdf->printCompactTree(os,Form("%s[Pdf]",indent)) ;
364  _cloneModel->printCompactTree(os,Form("%s[Mod]",indent)) ;
365  }
366 
367  os << indent << "RooNumConvolution end cache" << endl ;
368 }
369 
370 
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:51
RooListProxy _windowParam
#define coutE(a)
Definition: RooMsgService.h:35
RooNumIntConfig & convIntConfig()
RooAbsReal * _cloneModel
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.
RooNumIntConfig holds the configuration parameters of the various numeric integrators used by RooReal...
void initialize() const
One-time initialization of object.
virtual ~RooNumConvolution()
Destructor.
static RooNumIntFactory & instance()
Static method returning reference to singleton instance of factory.
RooRealProxy _origModel
void resetNumCall() const
Definition: RooAbsFunc.h:46
RooAbsReal * _cloneVar
Double_t getVal(const RooArgSet *set=0) const
Definition: RooAbsReal.h:64
RooRealProxy _origVar
Numeric integrator of convolution integrand.
RooNumIntConfig _convIntConfig
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
const Bool_t kFALSE
Definition: Rtypes.h:92
virtual Bool_t setLimits(Double_t *, Double_t *)
STL namespace.
#define coutW(a)
Definition: RooMsgService.h:34
Double_t evaluate() const
Calculate convolution integral.
Int_t numCall() const
Definition: RooAbsFunc.h:42
virtual Bool_t setLabel(const char *label, Bool_t printError=kTRUE)
Set value by specifying the name of the desired state If printError is set, a message will be printed...
void setNormalizationSet(const RooArgSet *nset)
RooCategory & method1D()
Double_t x[n]
Definition: legend1.C:17
RooAbsIntegrator * createIntegrator(RooAbsFunc &func, const RooNumIntConfig &config, Int_t ndim=0, Bool_t isBinned=kFALSE)
Construct a numeric integrator instance that operates on function &#39;func&#39; and is configured with &#39;conf...
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:1783
virtual Double_t integral(const Double_t *yvec=0)=0
RooAbsReal * _clonePdf
virtual void removeAll()
Remove all arguments from our set, deleting them if we own them.
Numeric 1-dimensional convolution operator PDF.
void setCallWarning(Int_t threshold=2000)
Activate warning messages if number of function calls needed for evaluation of convolution integral e...
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Reimplementation of standard RooArgList::add()
const RooArgSet * nset() const
Definition: RooAbsProxy.h:47
RooRealVar represents a fundamental (non-derived) real valued object.
Definition: RooRealVar.h:37
RooConvIntegrandBinding * _integrand
RooAbsArg * at(Int_t idx) const
Definition: RooArgList.h:84
RooAbsReal & pdf() const
void SetName(const char *name)
Set the name of the TNamed.
Definition: RooAbsArg.cxx:2390
char * Form(const char *fmt,...)
static Double_t infinity()
Return internal infinity representation.
Definition: RooNumber.cxx:49
RooArgSet _ownedClonedModelSet
static void indent(ostringstream &buf, int indent_level)
RooAbsReal & model() const
void replaceArg(const RooAbsArg &orig, const RooAbsArg &subst)
Replace any occurence of arg &#39;orig&#39; with arg &#39;subst&#39;.
#define ClassImp(name)
Definition: Rtypes.h:279
Double_t min(const char *rname=0) const
Definition: RooRealProxy.h:56
double Double_t
Definition: RtypesCore.h:55
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:53
virtual Bool_t setUseIntegrandLimits(Bool_t flag)
Interface function that allows to defer limit definition to integrand definition. ...
RooCategory & method1DOpen()
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 setCloneBranchSet(RooArgSet &cloneBranchSet)
Install the input RooArgSet as container in which all cloned branches will be stored.
RooAbsCollection is an abstract container object that can hold multiple RooAbsArg objects...
void clearConvolutionWindow()
Removes previously defined convolution window, reverting to convolution from -inf to +inf...
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 max(const char *rname=0) const
Definition: RooRealProxy.h:57
virtual void removeAll()
Reimplementation of standard RooArgList::removeAll()
RooAbsArg * build(const char *masterCatState, Bool_t verbose=kFALSE)
Build a clone of the prototype executing all registered &#39;replace&#39; rules and &#39;split&#39; rules for the mas...
RooRealProxy _origPdf
RooRealVar & var() const
Implementation of RooAbsFunc that represent the the integrand of a generic (numeric) convolution A (x...
const char * proto
Definition: civetweb.c:11652
const RooAbsReal & arg() const
Definition: RooRealProxy.h:43
RooAbsIntegrator * _integrator
Binding of Convolution Integrand function.
const Bool_t kTRUE
Definition: Rtypes.h:91
Int_t Fill(Double_t)
Invalid Fill method.
Definition: TH2.cxx:292
THist< 2, float, THistStatContent, THistStatUncertainty > TH2F
Definition: THist.hxx:308
RooArgSet _ownedClonedPdfSet
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...
return
Definition: HLFactory.cxx:514
char name[80]
Definition: TGX11.cxx:109
virtual const char * GetTitle() const
Returns title of object.
Definition: TNamed.h:52