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 <math.h>
38
39#include "RooGenProdProj.h"
40#include "RooAbsReal.h"
41#include "RooAbsPdf.h"
42#include "RooErrorHandler.h"
43#include "RooProduct.h"
44
45using namespace std;
46
48
49
50
51////////////////////////////////////////////////////////////////////////////////
52/// Default constructor
53
55 _compSetOwnedN(0),
56 _compSetOwnedD(0),
57 _haveD(kFALSE)
58{
59}
60
61
62////////////////////////////////////////////////////////////////////////////////
63/// Constructor for a normalization projection of the product of p.d.f.s _prodSet
64/// integrated over _intSet in range isetRangeName while normalized over _normSet
65
66RooGenProdProj::RooGenProdProj(const char *name, const char *title, const RooArgSet& _prodSet, const RooArgSet& _intSet,
67 const RooArgSet& _normSet, const char* isetRangeName, const char* normRangeName, Bool_t doFactorize) :
68 RooAbsReal(name, title),
69 _compSetOwnedN(0),
70 _compSetOwnedD(0),
71 _compSetN("compSetN","Set of integral components owned by numerator",this,kFALSE),
72 _compSetD("compSetD","Set of integral components owned by denominator",this,kFALSE),
73 _intList("intList","List of integrals",this,kTRUE),
74 _haveD(kFALSE)
75{
76 // Set expensive object cache to that of first item in prodSet
78
79 // Create owners of components created in ctor
82
83 RooAbsReal* numerator = makeIntegral("numerator",_prodSet,_intSet,*_compSetOwnedN,isetRangeName,doFactorize) ;
84 RooAbsReal* denominator = makeIntegral("denominator",_prodSet,_normSet,*_compSetOwnedD,normRangeName,doFactorize) ;
85
86// cout << "RooGenProdPdf::ctor(" << GetName() << ") numerator = " << numerator->GetName() << endl ;
87// numerator->printComponentTree() ;
88// cout << "RooGenProdPdf::ctor(" << GetName() << ") denominator = " << denominator->GetName() << endl ;
89// denominator->printComponentTree() ;
90
91 // Copy all components in (non-owning) set proxy
94
95 _intList.add(*numerator) ;
96 if (denominator) {
97 _intList.add(*denominator) ;
98 _haveD = kTRUE ;
99 }
100}
101
102
103
104////////////////////////////////////////////////////////////////////////////////
105/// Copy constructor
106
108 RooAbsReal(other, name),
109 _compSetOwnedN(0),
110 _compSetOwnedD(0),
111 _compSetN("compSetN","Set of integral components owned by numerator",this),
112 _compSetD("compSetD","Set of integral components owned by denominator",this),
113 _intList("intList","List of integrals",this)
114{
115 // Explicitly remove all server links at this point
116 TIterator* iter = serverIterator() ;
117 RooAbsArg* server ;
118 while((server=(RooAbsArg*)iter->Next())) {
119 removeServer(*server,kTRUE) ;
120 }
121 delete iter ;
122
123 // Copy constructor
126
129
130 RooAbsArg* arg ;
132 while((arg=(RooAbsArg*)nIter->Next())) {
133// cout << "ownedN elem " << arg->GetName() << "(" << arg << ")" << endl ;
134 arg->setOperMode(_operMode) ;
135 }
136 delete nIter ;
138 while((arg=(RooAbsArg*)dIter->Next())) {
139// cout << "ownedD elem " << arg->GetName() << "(" << arg << ")" << endl ;
140 arg->setOperMode(_operMode) ;
141 }
142 delete dIter ;
143
144 // Fill _intList
145 _haveD = other._haveD ;
146 _intList.add(*_compSetN.find(other._intList.at(0)->GetName())) ;
147 if (other._haveD) {
148 _intList.add(*_compSetD.find(other._intList.at(1)->GetName())) ;
149 }
150}
151
152
153
154////////////////////////////////////////////////////////////////////////////////
155/// Destructor
156
158{
159 if (_compSetOwnedN) delete _compSetOwnedN ;
160 if (_compSetOwnedD) delete _compSetOwnedD ;
161}
162
163
164
165////////////////////////////////////////////////////////////////////////////////
166/// Utility function to create integral for product over certain observables.
167/// \param[in] name Name of integral to be created.
168/// \param[in] compSet All components of the product.
169/// \param[in] intSet Observables to be integrated.
170/// \param[in] isetRangeName Integral range.
171/// \param[out] saveSet All component objects needed to represent the product integral are added as owned members to saveSet.
172/// The caller should take ownership of this set.
173/// \return A RooAbsReal object representing the requested integral.
174///
175/// The integration is factorized into components as much as possible and done analytically as far as possible.
176
177RooAbsReal* RooGenProdProj::makeIntegral(const char* name, const RooArgSet& compSet, const RooArgSet& intSet,
178 RooArgSet& saveSet, const char* isetRangeName, Bool_t doFactorize)
179{
180 RooArgSet anaIntSet, numIntSet ;
181
182 // First determine subset of observables in intSet that are factorizable
183 TIterator* compIter = compSet.createIterator() ;
184 TIterator* intIter = intSet.createIterator() ;
185 RooAbsPdf* pdf ;
186 RooAbsArg* arg ;
187 while((arg=(RooAbsArg*)intIter->Next())) {
188 Int_t count(0) ;
189 compIter->Reset() ;
190 while((pdf=(RooAbsPdf*)compIter->Next())) {
191 if (pdf->dependsOn(*arg)) count++ ;
192 }
193
194 if (count==0) {
195 } else if (count==1) {
196 anaIntSet.add(*arg) ;
197 } else {
198 }
199 }
200
201 // Determine which of the factorizable integrals can be done analytically
202 RooArgSet prodSet ;
203 numIntSet.add(intSet) ;
204
205 compIter->Reset() ;
206 while((pdf=(RooAbsPdf*)compIter->Next())) {
207 if (doFactorize && pdf->dependsOn(anaIntSet)) {
208 RooArgSet anaSet ;
209 Int_t code = pdf->getAnalyticalIntegralWN(anaIntSet,anaSet,0,isetRangeName) ;
210 if (code!=0) {
211 // Analytical integral, create integral object
212 RooAbsReal* pai = pdf->createIntegral(anaSet,isetRangeName) ;
213 pai->setOperMode(_operMode) ;
214
215 // Add to integral to product
216 prodSet.add(*pai) ;
217
218 // Remove analytically integratable observables from numeric integration list
219 numIntSet.remove(anaSet) ;
220
221 // Declare ownership of integral
222 saveSet.addOwned(*pai) ;
223 } else {
224 // Analytic integration of factorizable observable not possible, add straight pdf to product
225 prodSet.add(*pdf) ;
226 }
227 } else {
228 // Non-factorizable observables, add straight pdf to product
229 prodSet.add(*pdf) ;
230 }
231 }
232
233 //cout << "RooGenProdProj::makeIntegral(" << GetName() << ") prodset = " << prodSet << endl ;
234
235 // Create product of (partial) analytical integrals
236 TString prodName ;
237 if (isetRangeName) {
238 prodName = Form("%s_%s_Range[%s]",GetName(),name,isetRangeName) ;
239 } else {
240 prodName = Form("%s_%s",GetName(),name) ;
241 }
242 RooProduct* prod = new RooProduct(prodName,"product",prodSet) ;
244 prod->setOperMode(_operMode) ;
245
246 // Declare owndership of product
247 saveSet.addOwned(*prod) ;
248
249 // Create integral performing remaining numeric integration over (partial) analytic product
250 RooAbsReal* ret = prod->createIntegral(numIntSet,isetRangeName) ;
251// cout << "verbose print of integral object" << endl ;
252// ret->Print("v") ;
253 ret->setOperMode(_operMode) ;
254 saveSet.addOwned(*ret) ;
255
256 delete compIter ;
257 delete intIter ;
258
259 // Caller owners returned master integral object
260 return ret ;
261}
262
263
264
265////////////////////////////////////////////////////////////////////////////////
266/// Calculate and return value of normalization projection
267
269{
270 Double_t nom = ((RooAbsReal*)_intList.at(0))->getVal() ;
271
272 if (!_haveD) return nom ;
273
274 Double_t den = ((RooAbsReal*)_intList.at(1))->getVal() ;
275
276 //cout << "RooGenProdProj::eval(" << GetName() << ") nom = " << nom << " den = " << den << endl ;
277
278 return nom / den ;
279}
280
281
282
283////////////////////////////////////////////////////////////////////////////////
284/// Intercept cache mode operation changes and propagate them to the components
285
287{
288 // WVE use cache manager here!
289
290 RooAbsArg* arg ;
292 while((arg=(RooAbsArg*)nIter->Next())) {
293 arg->setOperMode(_operMode) ;
294 }
295 delete nIter ;
296
298 while((arg=(RooAbsArg*)dIter->Next())) {
299 arg->setOperMode(_operMode) ;
300 }
301 delete dIter ;
302
304 if (_haveD) _intList.at(1)->setOperMode(Auto) ; // Denominator always stays in Auto mode (normalization integral)
305}
306
307
308
309
310
311
312
int Int_t
Definition: RtypesCore.h:43
const Bool_t kFALSE
Definition: RtypesCore.h:90
bool Bool_t
Definition: RtypesCore.h:61
double Double_t
Definition: RtypesCore.h:57
const Bool_t kTRUE
Definition: RtypesCore.h:89
#define ClassImp(name)
Definition: Rtypes.h:361
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:73
RooExpensiveObjectCache & expensiveObjectCache() const
Definition: RooAbsArg.cxx:2168
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:730
friend class RooArgSet
Definition: RooAbsArg.h:572
virtual void setExpensiveObjectCache(RooExpensiveObjectCache &cache)
Definition: RooAbsArg.h:521
void setOperMode(OperMode mode, Bool_t recurseADirty=kTRUE)
Change cache operation mode to given mode.
Definition: RooAbsArg.cxx:1718
OperMode _operMode
Mark batches as dirty (only meaningful for RooAbsReal).
Definition: RooAbsArg.h:656
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:405
TIterator * serverIterator() const
Definition: RooAbsArg.h:141
RooAbsArg * first() const
virtual Bool_t remove(const RooAbsArg &var, Bool_t silent=kFALSE, Bool_t matchByNameOnly=kFALSE)
Remove the specified argument from our list.
TIterator * createIterator(Bool_t dir=kIterForward) const
TIterator-style iteration over contained elements.
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:60
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:371
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:90
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:560
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:126
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:30
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