Logo ROOT   6.14/05
Reference Guide
RooNumConvPdf.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 RooNumConvPdf.cxx
19 \class RooNumConvPdf
20 \ingroup Roofitcore
21 
22 Numeric 1-dimensional convolution operator PDF. This class can convolve any PDF
23 with any other PDF using a straightforward numeric calculation of the
24 convolution integral
25 This class should be used as last resort as numeric convolution calculated
26 this way is computationally intensive and prone to stability fitting problems.
27 <b>The preferred way to compute numeric convolutions is RooFFTConvPdf</b>,
28 which calculates convolutions using Fourier Transforms (requires external free
29 FFTW3 package)
30 RooNumConvPdf implements reasonable defaults that should convolve most
31 functions reasonably well, but results strongly depend on the shape of your
32 input PDFS so always check your result.
33 The default integration engine for the numeric convolution is the
34 adaptive Gauss-Kronrod method, which empirically seems the most robust
35 for this task. You can override the convolution integration settings via
36 the RooNumIntConfig object reference returned by the convIntConfig() member
37 function
38 By default the numeric convolution is integrated from -infinity to
39 +infinity through a <pre>x -> 1/x</pre> coordinate transformation of the
40 tails. For convolution with a very small bandwidth it may be
41 advantageous (for both CPU consumption and stability) if the
42 integration domain is limited to a finite range. The function
43 setConvolutionWindow(mean,width,scale) allows to set a sliding
44 window around the x value to be calculated taking a RooAbsReal
45 expression for an offset and a width to be taken around the x
46 value. These input expression can be RooFormulaVars or other
47 function objects although the 3d 'scale' argument 'scale'
48 multiplies the width RooAbsReal expression given in the 2nd
49 argument, allowing for an appropriate window definition for most
50 cases without need for a RooFormulaVar object: e.g. a Gaussian
51 resolution PDF do setConvolutionWindow(gaussMean,gaussSigma,5)
52 Note that for a 'wide' Gaussian the -inf to +inf integration
53 may converge more quickly than that over a finite range!
54 The default numeric precision is 1e-7, i.e. the global default for
55 numeric integration but you should experiment with this value to
56 see if it is sufficient for example by studying the number of function
57 calls that MINUIT needs to fit your function as function of the
58 convolution precision.
59 **/
60 
61 #include "RooFit.h"
62 
63 #include "Riostream.h"
64 #include "Riostream.h"
65 #include "TH2F.h"
66 #include "RooNumConvPdf.h"
67 #include "RooArgList.h"
68 #include "RooRealVar.h"
69 #include "RooFormulaVar.h"
70 #include "RooCustomizer.h"
72 #include "RooNumIntFactory.h"
73 #include "RooGenContext.h"
74 #include "RooConvGenContext.h"
75 
76 
77 
78 using namespace std;
79 
81 ;
82 
83 
84 
85 ////////////////////////////////////////////////////////////////////////////////
86 
88  _init(kFALSE),
89  _conv(0)
90 {
91 }
92 
93 
94 
95 
96 //_____________________________________________________________________________R
97 RooNumConvPdf::RooNumConvPdf(const char *name, const char *title, RooRealVar& convVar, RooAbsPdf& inPdf, RooAbsPdf& resmodel) :
98  RooAbsPdf(name,title),
99  _init(kFALSE),
100  _conv(0),
101  _origVar("!origVar","Original Convolution variable",this,convVar),
102  _origPdf("!origPdf","Original Input PDF",this,inPdf),
103  _origModel("!origModel","Original Resolution model",this,resmodel)
104 {
105  // Constructor of convolution operator PDF
106  //
107  // convVar : convolution variable (on which both pdf and resmodel should depend)
108  // pdf : input 'physics' pdf
109  // resmodel : input 'resultion' pdf
110  //
111  // output is pdf(x) (X) resmodel(x) = Int [ pdf(x') resmodel (x-x') ] dx'
112  //
113 }
114 
115 
116 
117 ////////////////////////////////////////////////////////////////////////////////
118 /// Copy constructor
119 
120 RooNumConvPdf::RooNumConvPdf(const RooNumConvPdf& other, const char* name) :
121  RooAbsPdf(other,name),
122  _init(kFALSE),
123  _origVar("!origVar",this,other._origVar),
124  _origPdf("!origPdf",this,other._origPdf),
125  _origModel("!origModel",this,other._origModel)
126 {
127  // Make temporary clone of original convolution to preserve configuration information
128  // This information will be propagated to a newly create convolution in a subsequent
129  // call to initialize()
130  if (other._conv) {
131  _conv = new RooNumConvolution(*other._conv,Form("%s_CONV",name?name:GetName())) ;
132  } else {
133  _conv = 0 ;
134  }
135 }
136 
137 
138 
139 ////////////////////////////////////////////////////////////////////////////////
140 /// Destructor
141 
143 {
144  if (_init) {
145  delete _conv ;
146  }
147 }
148 
149 
150 
151 ////////////////////////////////////////////////////////////////////////////////
152 /// Calculate and return value of p.d.f
153 
155 {
156  if (!_init) initialize() ;
157 
158  return _conv->evaluate() ;
159 }
160 
161 
162 
163 ////////////////////////////////////////////////////////////////////////////////
164 /// One-time initialization of object
165 
167 {
168  // Save pointer to any prototype convolution object (only present if this object is made through
169  // a copy constructor)
170  RooNumConvolution* protoConv = _conv ;
171 
172  // Optionally pass along configuration data from prototype object
173  _conv = new RooNumConvolution(Form("%s_CONV",GetName()),GetTitle(),var(),pdf(),model(),protoConv) ;
174 
175  // Delete prototype object now
176  if (protoConv) {
177  delete protoConv ;
178  }
179 
180  _init = kTRUE ;
181 }
182 
183 
184 
185 ////////////////////////////////////////////////////////////////////////////////
186 /// Return appropriate generator context for this convolved p.d.f. If both pdf and resolution
187 /// model support internal generation return and optimization convolution generation context
188 /// that uses a smearing algorithm. Otherwise return a standard accept/reject sampling
189 /// context on the convoluted shape.
190 
192  const RooArgSet* auxProto, Bool_t verbose) const
193 {
194  if (!_init) initialize() ;
195 
196  // Check if physics PDF and resolution model can both directly generate the convolution variable
197  RooArgSet* modelDep = _conv->model().getObservables(&vars) ;
198  modelDep->remove(_conv->var(),kTRUE,kTRUE) ;
199  Int_t numAddDep = modelDep->getSize() ;
200  delete modelDep ;
201 
202  RooArgSet dummy ;
203  Bool_t pdfCanDir = (((RooAbsPdf&)_conv->pdf()).getGenerator(_conv->var(),dummy) != 0 && \
205  Bool_t resCanDir = (((RooAbsPdf&)_conv->model()).getGenerator(_conv->var(),dummy) !=0 &&
207 
208  if (numAddDep>0 || !pdfCanDir || !resCanDir) {
209  // Any resolution model with more dependents than the convolution variable
210  // or pdf or resmodel do not support direct generation
211  return new RooGenContext(*this,vars,prototype,auxProto,verbose) ;
212  }
213 
214  // Any other resolution model: use specialized generator context
215  return new RooConvGenContext(*this,vars,prototype,auxProto,verbose) ;
216 }
217 
218 
219 
220 ////////////////////////////////////////////////////////////////////////////////
221 /// Customized printing of arguments of a RooNumConvPdf to more intuitively reflect the contents of the
222 /// product operator construction
223 
224 void RooNumConvPdf::printMetaArgs(ostream& os) const
225 {
226  os << _origPdf.arg().GetName() << "(" << _origVar.arg().GetName() << ") (*) " << _origModel.arg().GetName() << "(" << _origVar.arg().GetName() << ") " ;
227 }
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
virtual ~RooNumConvPdf()
Destructor.
RooArgSet * getObservables(const RooArgSet &set, Bool_t valueOnly=kTRUE) const
Definition: RooAbsArg.h:194
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
STL namespace.
Double_t evaluate() const
Calculate convolution integral.
friend class RooConvGenContext
Definition: RooNumConvPdf.h:76
RooAbsReal & model() const
Definition: RooNumConvPdf.h:54
Numeric 1-dimensional convolution operator PDF.
void initialize() const
do not persist
RooAbsGenContext is the abstract base class for generator contexts of RooAbsPdf objects.
void printMetaArgs(std::ostream &os) const
Customized printing of arguments of a RooNumConvPdf to more intuitively reflect the contents of the p...
RooRealVar represents a fundamental (non-derived) real valued object.
Definition: RooRealVar.h:36
virtual Double_t evaluate() const
Calculate and return value of p.d.f.
Int_t getSize() const
RooAbsReal & pdf() const
char * Form(const char *fmt,...)
RooAbsReal & pdf() const
Definition: RooNumConvPdf.h:53
virtual RooAbsGenContext * genContext(const RooArgSet &vars, const RooDataSet *prototype=0, const RooArgSet *auxProto=0, Bool_t verbose=kFALSE) const
Return appropriate generator context for this convolved p.d.f.
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:29
Class RooGenContext implement a universal generator context for all RooAbsPdf classes that do not hav...
Definition: RooGenContext.h:30
virtual Bool_t isDirectGenSafe(const RooAbsArg &arg) const
Check if given observable can be safely generated using the pdfs internal generator mechanism (if tha...
Definition: RooAbsPdf.cxx:2159
const Bool_t kFALSE
Definition: RtypesCore.h:88
RooAbsReal & model() const
#define ClassImp(name)
Definition: Rtypes.h:359
double Double_t
Definition: RtypesCore.h:55
virtual Int_t getGenerator(const RooArgSet &directVars, RooArgSet &generateVars, Bool_t staticInitOK=kTRUE) const
Load generatedVars with the subset of directVars that we can generate events for, and return a code t...
Definition: RooAbsPdf.cxx:2124
static RooMathCoreReg dummy
RooRealProxy _origPdf
Definition: RooNumConvPdf.h:70
RooRealProxy _origVar
Actual convolution calculation.
Definition: RooNumConvPdf.h:69
virtual Bool_t remove(const RooAbsArg &var, Bool_t silent=kFALSE, Bool_t matchByNameOnly=kFALSE)
Remove the specified argument from our list.
RooRealProxy _origModel
Definition: RooNumConvPdf.h:71
RooAbsPdf is the abstract interface for all probability density functions The class provides hybrid a...
Definition: RooAbsPdf.h:41
RooRealVar & var() const
Definition: RooNumConvPdf.h:52
RooRealVar & var() const
const RooAbsReal & arg() const
Definition: RooRealProxy.h:43
Numeric 1-dimensional convolution operator PDF.
Definition: RooNumConvPdf.h:26
RooNumConvolution * _conv
Definition: RooNumConvPdf.h:67
const Bool_t kTRUE
Definition: RtypesCore.h:87
char name[80]
Definition: TGX11.cxx:109
virtual const char * GetTitle() const
Returns title of object.
Definition: TNamed.h:48