Logo ROOT  
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 "TH2F.h"
65 #include "RooNumConvPdf.h"
66 #include "RooArgList.h"
67 #include "RooRealVar.h"
68 #include "RooFormulaVar.h"
69 #include "RooCustomizer.h"
71 #include "RooNumIntFactory.h"
72 #include "RooGenContext.h"
73 #include "RooConvGenContext.h"
74 
75 
76 
77 using namespace std;
78 
80 
81 
82 
83 
84 ////////////////////////////////////////////////////////////////////////////////
85 
87  _init(kFALSE),
88  _conv(0)
89 {
90 }
91 
92 
93 
94 
95 //_____________________________________________________________________________R
96 RooNumConvPdf::RooNumConvPdf(const char *name, const char *title, RooRealVar& convVar, RooAbsPdf& inPdf, RooAbsPdf& resmodel) :
97  RooAbsPdf(name,title),
98  _init(kFALSE),
99  _conv(0),
100  _origVar("!origVar","Original Convolution variable",this,convVar),
101  _origPdf("!origPdf","Original Input PDF",this,inPdf),
102  _origModel("!origModel","Original Resolution model",this,resmodel)
103 {
104  // Constructor of convolution operator PDF
105  //
106  // convVar : convolution variable (on which both pdf and resmodel should depend)
107  // pdf : input 'physics' pdf
108  // resmodel : input 'resultion' pdf
109  //
110  // output is pdf(x) (X) resmodel(x) = Int [ pdf(x') resmodel (x-x') ] dx'
111  //
112 }
113 
114 
115 
116 ////////////////////////////////////////////////////////////////////////////////
117 /// Copy constructor
118 
119 RooNumConvPdf::RooNumConvPdf(const RooNumConvPdf& other, const char* name) :
120  RooAbsPdf(other,name),
121  _init(kFALSE),
122  _origVar("!origVar",this,other._origVar),
123  _origPdf("!origPdf",this,other._origPdf),
124  _origModel("!origModel",this,other._origModel)
125 {
126  // Make temporary clone of original convolution to preserve configuration information
127  // This information will be propagated to a newly create convolution in a subsequent
128  // call to initialize()
129  if (other._conv) {
130  _conv = new RooNumConvolution(*other._conv,Form("%s_CONV",name?name:GetName())) ;
131  } else {
132  _conv = 0 ;
133  }
134 }
135 
136 
137 
138 ////////////////////////////////////////////////////////////////////////////////
139 /// Destructor
140 
142 {
143  if (_init) {
144  delete _conv ;
145  }
146 }
147 
148 
149 
150 ////////////////////////////////////////////////////////////////////////////////
151 /// Calculate and return value of p.d.f
152 
154 {
155  if (!_init) initialize() ;
156 
157  return _conv->evaluate() ;
158 }
159 
160 
161 
162 ////////////////////////////////////////////////////////////////////////////////
163 /// One-time initialization of object
164 
166 {
167  // Save pointer to any prototype convolution object (only present if this object is made through
168  // a copy constructor)
169  RooNumConvolution* protoConv = _conv ;
170 
171  // Optionally pass along configuration data from prototype object
172  _conv = new RooNumConvolution(Form("%s_CONV",GetName()),GetTitle(),var(),pdf(),model(),protoConv) ;
173 
174  // Delete prototype object now
175  if (protoConv) {
176  delete protoConv ;
177  }
178 
179  _init = kTRUE ;
180 }
181 
182 
183 
184 ////////////////////////////////////////////////////////////////////////////////
185 /// Return appropriate generator context for this convolved p.d.f. If both pdf and resolution
186 /// model support internal generation return and optimization convolution generation context
187 /// that uses a smearing algorithm. Otherwise return a standard accept/reject sampling
188 /// context on the convoluted shape.
189 
191  const RooArgSet* auxProto, Bool_t verbose) const
192 {
193  if (!_init) initialize() ;
194 
195  // Check if physics PDF and resolution model can both directly generate the convolution variable
196  RooArgSet* modelDep = _conv->model().getObservables(&vars) ;
197  modelDep->remove(_conv->var(),kTRUE,kTRUE) ;
198  Int_t numAddDep = modelDep->getSize() ;
199  delete modelDep ;
200 
201  RooArgSet dummy ;
202  Bool_t pdfCanDir = (((RooAbsPdf&)_conv->pdf()).getGenerator(_conv->var(),dummy) != 0 && \
204  Bool_t resCanDir = (((RooAbsPdf&)_conv->model()).getGenerator(_conv->var(),dummy) !=0 &&
206 
207  if (numAddDep>0 || !pdfCanDir || !resCanDir) {
208  // Any resolution model with more dependents than the convolution variable
209  // or pdf or resmodel do not support direct generation
210  return new RooGenContext(*this,vars,prototype,auxProto,verbose) ;
211  }
212 
213  // Any other resolution model: use specialized generator context
214  return new RooConvGenContext(*this,vars,prototype,auxProto,verbose) ;
215 }
216 
217 
218 
219 ////////////////////////////////////////////////////////////////////////////////
220 /// Customized printing of arguments of a RooNumConvPdf to more intuitively reflect the contents of the
221 /// product operator construction
222 
223 void RooNumConvPdf::printMetaArgs(ostream& os) const
224 {
225  os << _origPdf.arg().GetName() << "(" << _origVar.arg().GetName() << ") (*) " << _origModel.arg().GetName() << "(" << _origVar.arg().GetName() << ") " ;
226 }
RooNumConvPdf::RooConvGenContext
friend class RooConvGenContext
Definition: RooNumConvPdf.h:76
RooFormulaVar.h
kTRUE
const Bool_t kTRUE
Definition: RtypesCore.h:91
RooNumConvPdf::_init
Bool_t _init
Definition: RooNumConvPdf.h:65
RooFit.h
RooTemplateProxy::arg
const T & arg() const
Return reference to object held in proxy.
Definition: RooTemplateProxy.h:271
RooNumConvPdf::printMetaArgs
void printMetaArgs(std::ostream &os) const
Customized printing of arguments of a RooNumConvPdf to more intuitively reflect the contents of the p...
Definition: RooNumConvPdf.cxx:223
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
RooGenContext
Definition: RooGenContext.h:30
RooAbsCollection::remove
virtual Bool_t remove(const RooAbsArg &var, Bool_t silent=kFALSE, Bool_t matchByNameOnly=kFALSE)
Remove the specified argument from our list.
Definition: RooAbsCollection.cxx:584
RooAbsPdf::getGenerator
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:2413
RooNumConvPdf::_origPdf
RooRealProxy _origPdf
Definition: RooNumConvPdf.h:70
RooNumConvPdf::_origVar
RooRealProxy _origVar
Actual convolution calculation.
Definition: RooNumConvPdf.h:69
RooNumConvPdf::_origModel
RooRealProxy _origModel
Definition: RooNumConvPdf.h:71
RooNumConvPdf::var
RooRealVar & var() const
Definition: RooNumConvPdf.h:52
RooNumConvPdf::RooNumConvPdf
RooNumConvPdf()
Definition: RooNumConvPdf.cxx:86
RooAbsPdf::isDirectGenSafe
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:2448
bool
RooNumConvPdf.h
RooNumConvolution::evaluate
Double_t evaluate() const
Calculate convolution integral.
Definition: RooNumConvolution.cxx:229
RooNumConvolution
Definition: RooNumConvolution.h:29
RooCustomizer.h
RooNumConvolution::pdf
RooAbsReal & pdf() const
Definition: RooNumConvolution.h:56
kFALSE
const Bool_t kFALSE
Definition: RtypesCore.h:92
RooNumConvPdf::~RooNumConvPdf
virtual ~RooNumConvPdf()
Destructor.
Definition: RooNumConvPdf.cxx:141
RooAbsGenContext
Definition: RooAbsGenContext.h:26
RooGenContext.h
RooNumConvPdf::_conv
RooNumConvolution * _conv
Definition: RooNumConvPdf.h:67
RooNumConvolution::model
RooAbsReal & model() const
Definition: RooNumConvolution.h:57
RooRealVar.h
TH2F.h
RooNumConvolution::var
RooRealVar & var() const
Definition: RooNumConvolution.h:55
Double_t
double Double_t
Definition: RtypesCore.h:59
RooAbsArg::getObservables
RooArgSet * getObservables(const RooArgSet &set, Bool_t valueOnly=kTRUE) const
Return the observables of this pdf given a set of observables.
Definition: RooAbsArg.h:294
dummy
static RooMathCoreReg dummy
Definition: RooMathCoreReg.cxx:27
RooNumConvPdf::evaluate
virtual Double_t evaluate() const
Calculate and return value of p.d.f.
Definition: RooNumConvPdf.cxx:153
RooNumConvPdf
Definition: RooNumConvPdf.h:26
name
char name[80]
Definition: TGX11.cxx:110
RooNumConvPdf::model
RooAbsReal & model() const
Definition: RooNumConvPdf.h:54
genreflex::verbose
bool verbose
Definition: rootcling_impl.cxx:133
RooConvIntegrandBinding.h
RooDataSet
Definition: RooDataSet.h:33
RooAbsPdf
Definition: RooAbsPdf.h:40
TNamed::GetName
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:53
RooNumConvPdf::initialize
void initialize() const
do not persist
Definition: RooNumConvPdf.cxx:165
RooRealVar
Definition: RooRealVar.h:35
Riostream.h
RooArgList.h
RooNumConvPdf::pdf
RooAbsReal & pdf() const
Definition: RooNumConvPdf.h:53
RooConvGenContext.h
RooNumConvPdf::genContext
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.
Definition: RooNumConvPdf.cxx:190
RooNumIntFactory.h
RooAbsCollection::getSize
Int_t getSize() const
Definition: RooAbsCollection.h:171
RooArgSet
Definition: RooArgSet.h:28
int