Logo ROOT  
Reference Guide
RooGenProdProj.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 RooGenProdProj.cxx
19\class RooGenProdProj
20\ingroup Roofitcore
21
22
23RooGenProdProj is an auxiliary class for RooProdPdf that calculates
24a general normalised projection of a product of non-factorising PDFs, e.g.
25\f[
26 P_{x,xy} = \frac{\int ( P1 * P2 * \ldots) \mathrm{d}x}{\int ( P1 * P2 * \ldots ) \mathrm{d}x \mathrm{d}y}
27\f]
28
29Partial integrals, which factorise and can be calculated, are calculated
30analytically. Remaining non-factorising observables are integrated numerically.
31**/
32
33
34#include "RooFit.h"
35
36#include "Riostream.h"
37#include "Riostream.h"
38#include <math.h>
39
40#include "RooGenProdProj.h"
41#include "RooAbsReal.h"
42#include "RooAbsPdf.h"
43#include "RooErrorHandler.h"
44#include "RooProduct.h"
45
46using namespace std;
47
49;
50
51
52////////////////////////////////////////////////////////////////////////////////
53/// Default constructor
54
56 _compSetOwnedN(0),
57 _compSetOwnedD(0),
58 _haveD(kFALSE)
59{
60}
61
62
63////////////////////////////////////////////////////////////////////////////////
64/// Constructor for a normalization projection of the product of p.d.f.s _prodSet
65/// integrated over _intSet in range isetRangeName while normalized over _normSet
66
67RooGenProdProj::RooGenProdProj(const char *name, const char *title, const RooArgSet& _prodSet, const RooArgSet& _intSet,
68 const RooArgSet& _normSet, const char* isetRangeName, const char* normRangeName, Bool_t doFactorize) :
69 RooAbsReal(name, title),
70 _compSetOwnedN(0),
71 _compSetOwnedD(0),
72 _compSetN("compSetN","Set of integral components owned by numerator",this,kFALSE),
73 _compSetD("compSetD","Set of integral components owned by denominator",this,kFALSE),
74 _intList("intList","List of integrals",this,kTRUE),
75 _haveD(kFALSE)
76{
77 // Set expensive object cache to that of first item in prodSet
79
80 // Create owners of components created in ctor
83
84 RooAbsReal* numerator = makeIntegral("numerator",_prodSet,_intSet,*_compSetOwnedN,isetRangeName,doFactorize) ;
85 RooAbsReal* denominator = makeIntegral("denominator",_prodSet,_normSet,*_compSetOwnedD,normRangeName,doFactorize) ;
86
87// cout << "RooGenProdPdf::ctor(" << GetName() << ") numerator = " << numerator->GetName() << endl ;
88// numerator->printComponentTree() ;
89// cout << "RooGenProdPdf::ctor(" << GetName() << ") denominator = " << denominator->GetName() << endl ;
90// denominator->printComponentTree() ;
91
92 // Copy all components in (non-owning) set proxy
95
96 _intList.add(*numerator) ;
97 if (denominator) {
98 _intList.add(*denominator) ;
99 _haveD = kTRUE ;
100 }
101}
102
103
104
105////////////////////////////////////////////////////////////////////////////////
106/// Copy constructor
107
109 RooAbsReal(other, name),
110 _compSetOwnedN(0),
111 _compSetOwnedD(0),
112 _compSetN("compSetN","Set of integral components owned by numerator",this),
113 _compSetD("compSetD","Set of integral components owned by denominator",this),
114 _intList("intList","List of integrals",this)
115{
116 // Explicitly remove all server links at this point
117 TIterator* iter = serverIterator() ;
118 RooAbsArg* server ;
119 while((server=(RooAbsArg*)iter->Next())) {
120 removeServer(*server,kTRUE) ;
121 }
122 delete iter ;
123
124 // Copy constructor
127
130
131 RooAbsArg* arg ;
133 while((arg=(RooAbsArg*)nIter->Next())) {
134// cout << "ownedN elem " << arg->GetName() << "(" << arg << ")" << endl ;
135 arg->setOperMode(_operMode) ;
136 }
137 delete nIter ;
139 while((arg=(RooAbsArg*)dIter->Next())) {
140// cout << "ownedD elem " << arg->GetName() << "(" << arg << ")" << endl ;
141 arg->setOperMode(_operMode) ;
142 }
143 delete dIter ;
144
145 // Fill _intList
146 _haveD = other._haveD ;
147 _intList.add(*_compSetN.find(other._intList.at(0)->GetName())) ;
148 if (other._haveD) {
149 _intList.add(*_compSetD.find(other._intList.at(1)->GetName())) ;
150 }
151}
152
153
154
155////////////////////////////////////////////////////////////////////////////////
156/// Destructor
157
159{
160 if (_compSetOwnedN) delete _compSetOwnedN ;
161 if (_compSetOwnedD) delete _compSetOwnedD ;
162}
163
164
165
166////////////////////////////////////////////////////////////////////////////////
167/// Utility function to create integral for product over certain observables.
168/// \param[in] name Name of integral to be created.
169/// \param[in] compSet All components of the product.
170/// \param[in] intSet Observables to be integrated.
171/// \param[in] isetRangeName Integral range.
172/// \param[out] saveSet All component objects needed to represent the product integral are added as owned members to saveSet.
173/// The caller should take ownership of this set.
174/// \return A RooAbsReal object representing the requested integral.
175///
176/// The integration is factorized into components as much as possible and done analytically as far as possible.
177
178RooAbsReal* RooGenProdProj::makeIntegral(const char* name, const RooArgSet& compSet, const RooArgSet& intSet,
179 RooArgSet& saveSet, const char* isetRangeName, Bool_t doFactorize)
180{
181 RooArgSet anaIntSet, numIntSet ;
182
183 // First determine subset of observables in intSet that are factorizable
184 TIterator* compIter = compSet.createIterator() ;
185 TIterator* intIter = intSet.createIterator() ;
186 RooAbsPdf* pdf ;
187 RooAbsArg* arg ;
188 while((arg=(RooAbsArg*)intIter->Next())) {
189 Int_t count(0) ;
190 compIter->Reset() ;
191 while((pdf=(RooAbsPdf*)compIter->Next())) {
192 if (pdf->dependsOn(*arg)) count++ ;
193 }
194
195 if (count==0) {
196 } else if (count==1) {
197 anaIntSet.add(*arg) ;
198 } else {
199 }
200 }
201
202 // Determine which of the factorizable integrals can be done analytically
203 RooArgSet prodSet ;
204 numIntSet.add(intSet) ;
205
206 compIter->Reset() ;
207 while((pdf=(RooAbsPdf*)compIter->Next())) {
208 if (doFactorize && pdf->dependsOn(anaIntSet)) {
209 RooArgSet anaSet ;
210 Int_t code = pdf->getAnalyticalIntegralWN(anaIntSet,anaSet,0,isetRangeName) ;
211 if (code!=0) {
212 // Analytical integral, create integral object
213 RooAbsReal* pai = pdf->createIntegral(anaSet,isetRangeName) ;
214 pai->setOperMode(_operMode) ;
215
216 // Add to integral to product
217 prodSet.add(*pai) ;
218
219 // Remove analytically integratable observables from numeric integration list
220 numIntSet.remove(anaSet) ;
221
222 // Declare ownership of integral
223 saveSet.addOwned(*pai) ;
224 } else {
225 // Analytic integration of factorizable observable not possible, add straight pdf to product
226 prodSet.add(*pdf) ;
227 }
228 } else {
229 // Non-factorizable observables, add straight pdf to product
230 prodSet.add(*pdf) ;
231 }
232 }
233
234 //cout << "RooGenProdProj::makeIntegral(" << GetName() << ") prodset = " << prodSet << endl ;
235
236 // Create product of (partial) analytical integrals
237 TString prodName ;
238 if (isetRangeName) {
239 prodName = Form("%s_%s_Range[%s]",GetName(),name,isetRangeName) ;
240 } else {
241 prodName = Form("%s_%s",GetName(),name) ;
242 }
243 RooProduct* prod = new RooProduct(prodName,"product",prodSet) ;
245 prod->setOperMode(_operMode) ;
246
247 // Declare owndership of product
248 saveSet.addOwned(*prod) ;
249
250 // Create integral performing remaining numeric integration over (partial) analytic product
251 RooAbsReal* ret = prod->createIntegral(numIntSet,isetRangeName) ;
252// cout << "verbose print of integral object" << endl ;
253// ret->Print("v") ;
254 ret->setOperMode(_operMode) ;
255 saveSet.addOwned(*ret) ;
256
257 delete compIter ;
258 delete intIter ;
259
260 // Caller owners returned master integral object
261 return ret ;
262}
263
264
265
266////////////////////////////////////////////////////////////////////////////////
267/// Calculate and return value of normalization projection
268
270{
271 Double_t nom = ((RooAbsReal*)_intList.at(0))->getVal() ;
272
273 if (!_haveD) return nom ;
274
275 Double_t den = ((RooAbsReal*)_intList.at(1))->getVal() ;
276
277 //cout << "RooGenProdProj::eval(" << GetName() << ") nom = " << nom << " den = " << den << endl ;
278
279 return nom / den ;
280}
281
282
283
284////////////////////////////////////////////////////////////////////////////////
285/// Intercept cache mode operation changes and propagate them to the components
286
288{
289 // WVE use cache manager here!
290
291 RooAbsArg* arg ;
293 while((arg=(RooAbsArg*)nIter->Next())) {
294 arg->setOperMode(_operMode) ;
295 }
296 delete nIter ;
297
299 while((arg=(RooAbsArg*)dIter->Next())) {
300 arg->setOperMode(_operMode) ;
301 }
302 delete dIter ;
303
305 if (_haveD) _intList.at(1)->setOperMode(Auto) ; // Denominator always stays in Auto mode (normalization integral)
306}
307
308
309
310
311
312
313
int Int_t
Definition: RtypesCore.h:41
const Bool_t kFALSE
Definition: RtypesCore.h:88
bool Bool_t
Definition: RtypesCore.h:59
double Double_t
Definition: RtypesCore.h:55
const Bool_t kTRUE
Definition: RtypesCore.h:87
#define ClassImp(name)
Definition: Rtypes.h:365
char name[80]
Definition: TGX11.cxx:109
char * Form(const char *fmt,...)
RooAbsArg is the common abstract base class for objects that represent a value (of arbitrary type) an...
Definition: RooAbsArg.h:71
RooExpensiveObjectCache & expensiveObjectCache() const
Definition: RooAbsArg.cxx:2176
Bool_t dependsOn(const RooAbsCollection &serverList, const RooAbsArg *ignoreArg=0, Bool_t valueOnly=kFALSE) const
Test whether we depend on (ie, are served by) any object in the specified collection.
Definition: RooAbsArg.cxx:734
friend class RooArgSet
Definition: RooAbsArg.h:551
virtual void setExpensiveObjectCache(RooExpensiveObjectCache &cache)
Definition: RooAbsArg.h:500
void setOperMode(OperMode mode, Bool_t recurseADirty=kTRUE)
Change cache operation mode to given mode.
Definition: RooAbsArg.cxx:1726
OperMode _operMode
Mark batches as dirty (only meaningful for RooAbsReal).
Definition: RooAbsArg.h:635
void removeServer(RooAbsArg &server, Bool_t force=kFALSE)
Unregister another RooAbsArg as a server to us, ie, declare that we no longer depend on its value and...
Definition: RooAbsArg.cxx:409
TIterator * serverIterator() const R__SUGGEST_ALTERNATIVE("Use servers() and begin()
RooAbsArg * first() const
TIterator * createIterator(Bool_t dir=kIterForward) const R__SUGGEST_ALTERNATIVE("begin()
TIterator-style iteration over contained elements.
virtual Bool_t remove(const RooAbsArg &var, Bool_t silent=kFALSE, Bool_t matchByNameOnly=kFALSE)
Remove the specified argument from our list.
RooAbsArg * find(const char *name) const
Find object with given name in list.
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:59
virtual Int_t getAnalyticalIntegralWN(RooArgSet &allVars, RooArgSet &analVars, const RooArgSet *normSet, const char *rangeName=0) const
Variant of getAnalyticalIntegral that is also passed the normalization set that should be applied to ...
Definition: RooAbsReal.cxx:373
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:87
RooAbsReal * createIntegral(const RooArgSet &iset, const RooCmdArg &arg1, const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none(), const RooCmdArg &arg6=RooCmdArg::none(), const RooCmdArg &arg7=RooCmdArg::none(), const RooCmdArg &arg8=RooCmdArg::none()) const
Create an object that represents the integral of the function over one or more observables listed in ...
Definition: RooAbsReal.cxx:562
RooAbsArg * at(Int_t idx) const
Return object at given index, or nullptr if index is out of range.
Definition: RooArgList.h:74
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
RooArgSet * snapshot(bool deepCopy=true) const
Use RooAbsCollection::snapshot(), but return as RooArgSet.
Definition: RooArgSet.h:134
virtual Bool_t addOwned(const RooAbsCollection &col, Bool_t silent=kFALSE)
Add a collection of arguments to this collection by calling addOwned() for each element in the source...
Definition: RooArgSet.h:92
virtual Bool_t add(const RooAbsCollection &col, Bool_t silent=kFALSE)
Add a collection of arguments to this collection by calling add() for each element in the source coll...
Definition: RooArgSet.h:88
RooGenProdProj is an auxiliary class for RooProdPdf that calculates a general normalised projection o...
virtual ~RooGenProdProj()
Destructor.
RooGenProdProj()
Default constructor.
RooSetProxy _compSetD
virtual void operModeHook()
Intercept cache mode operation changes and propagate them to the components.
Double_t evaluate() const
Calculate and return value of normalization projection.
RooListProxy _intList
RooArgSet * _compSetOwnedD
RooAbsReal * makeIntegral(const char *name, const RooArgSet &compSet, const RooArgSet &intSet, RooArgSet &saveSet, const char *isetRangeName, Bool_t doFactorize)
Utility function to create integral for product over certain observables.
RooArgSet * _compSetOwnedN
RooSetProxy _compSetN
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Reimplementation of standard RooArgList::add()
A RooProduct represents the product of a given set of RooAbsReal objects.
Definition: RooProduct.h:32
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Overloaded RooArgSet::add() method inserts 'var' into set and registers 'var' as server to owner with...
Iterator abstract base class.
Definition: TIterator.h:30
virtual void Reset()=0
virtual TObject * Next()=0
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
Basic string class.
Definition: TString.h:131