Logo ROOT   6.14/05
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 
23 RooGenProdProj is an auxiliary class for RooProdPdf that calculates
24 a general normalized projection of a product of non-factorizing PDFs, e.g.
25 
26  Int ( P1 * P2 * ....) dx
27  P_x_xy = -------------------------------
28  Int (P1 * P2 * ... ) dx dy
29 
30 Partial integrals that factorize that can be calculated are calculated
31 analytically. Remaining non-factorizing observables are integrated numerically.
32 **/
33 
34 
35 #include "RooFit.h"
36 
37 #include "Riostream.h"
38 #include "Riostream.h"
39 #include <math.h>
40 
41 #include "RooGenProdProj.h"
42 #include "RooAbsReal.h"
43 #include "RooAbsPdf.h"
44 #include "RooErrorHandler.h"
45 #include "RooProduct.h"
46 
47 using namespace std;
48 
50 ;
51 
52 
53 ////////////////////////////////////////////////////////////////////////////////
54 /// Default constructor
55 
57  _compSetOwnedN(0),
58  _compSetOwnedD(0),
59  _haveD(kFALSE)
60 {
61 }
62 
63 
64 ////////////////////////////////////////////////////////////////////////////////
65 /// Constructor for a normalization projection of the product of p.d.f.s _prodSet
66 /// integrated over _intSet in range isetRangeName while normalized over _normSet
67 
68 RooGenProdProj::RooGenProdProj(const char *name, const char *title, const RooArgSet& _prodSet, const RooArgSet& _intSet,
69  const RooArgSet& _normSet, const char* isetRangeName, const char* normRangeName, Bool_t doFactorize) :
70  RooAbsReal(name, title),
71  _compSetOwnedN(0),
72  _compSetOwnedD(0),
73  _compSetN("compSetN","Set of integral components owned by numerator",this,kFALSE),
74  _compSetD("compSetD","Set of integral components owned by denominator",this,kFALSE),
75  _intList("intList","List of integrals",this,kTRUE),
76  _haveD(kFALSE)
77 {
78  // Set expensive object cache to that of first item in prodSet
80 
81  // Create owners of components created in ctor
84 
85  RooAbsReal* numerator = makeIntegral("numerator",_prodSet,_intSet,*_compSetOwnedN,isetRangeName,doFactorize) ;
86  RooAbsReal* denominator = makeIntegral("denominator",_prodSet,_normSet,*_compSetOwnedD,normRangeName,doFactorize) ;
87 
88 // cout << "RooGenProdPdf::ctor(" << GetName() << ") numerator = " << numerator->GetName() << endl ;
89 // numerator->printComponentTree() ;
90 // cout << "RooGenProdPdf::ctor(" << GetName() << ") denominator = " << denominator->GetName() << endl ;
91 // denominator->printComponentTree() ;
92 
93  // Copy all components in (non-owning) set proxy
94  _compSetN.add(*_compSetOwnedN) ;
96 
97  _intList.add(*numerator) ;
98  if (denominator) {
99  _intList.add(*denominator) ;
100  _haveD = kTRUE ;
101  }
102 }
103 
104 
105 
106 ////////////////////////////////////////////////////////////////////////////////
107 /// Copy constructor
108 
110  RooAbsReal(other, name),
111  _compSetOwnedN(0),
112  _compSetOwnedD(0),
113  _compSetN("compSetN","Set of integral components owned by numerator",this),
114  _compSetD("compSetD","Set of integral components owned by denominator",this),
115  _intList("intList","List of integrals",this)
116 {
117  // Explicitly remove all server links at this point
118  TIterator* iter = serverIterator() ;
119  RooAbsArg* server ;
120  while((server=(RooAbsArg*)iter->Next())) {
121  removeServer(*server,kTRUE) ;
122  }
123  delete iter ;
124 
125  // Copy constructor
128 
131 
132  RooAbsArg* arg ;
134  while((arg=(RooAbsArg*)nIter->Next())) {
135 // cout << "ownedN elem " << arg->GetName() << "(" << arg << ")" << endl ;
136  arg->setOperMode(_operMode) ;
137  }
138  delete nIter ;
140  while((arg=(RooAbsArg*)dIter->Next())) {
141 // cout << "ownedD elem " << arg->GetName() << "(" << arg << ")" << endl ;
142  arg->setOperMode(_operMode) ;
143  }
144  delete dIter ;
145 
146  // Fill _intList
147  _haveD = other._haveD ;
148  _intList.add(*_compSetN.find(other._intList.at(0)->GetName())) ;
149  if (other._haveD) {
150  _intList.add(*_compSetD.find(other._intList.at(1)->GetName())) ;
151  }
152 }
153 
154 
155 
156 ////////////////////////////////////////////////////////////////////////////////
157 /// Destructor
158 
160 {
161  if (_compSetOwnedN) delete _compSetOwnedN ;
162  if (_compSetOwnedD) delete _compSetOwnedD ;
163 }
164 
165 
166 
167 ////////////////////////////////////////////////////////////////////////////////
168 /// Utility function to create integral over observables intSet in range isetRangeName over product of p.d.fs in compSet.
169 /// The integration is factorized into components as much as possible and done analytically as far as possible.
170 /// All component object needed to represent product integral are added as owned members to saveSet.
171 /// The return value is a RooAbsReal object representing the requested integral
172 
173 RooAbsReal* RooGenProdProj::makeIntegral(const char* name, const RooArgSet& compSet, const RooArgSet& intSet,
174  RooArgSet& saveSet, const char* isetRangeName, Bool_t doFactorize)
175 {
176  RooArgSet anaIntSet, numIntSet ;
177 
178  // First determine subset of observables in intSet that are factorizable
179  TIterator* compIter = compSet.createIterator() ;
180  TIterator* intIter = intSet.createIterator() ;
181  RooAbsPdf* pdf ;
182  RooAbsArg* arg ;
183  while((arg=(RooAbsArg*)intIter->Next())) {
184  Int_t count(0) ;
185  compIter->Reset() ;
186  while((pdf=(RooAbsPdf*)compIter->Next())) {
187  if (pdf->dependsOn(*arg)) count++ ;
188  }
189 
190  if (count==0) {
191  } else if (count==1) {
192  anaIntSet.add(*arg) ;
193  } else {
194  }
195  }
196 
197  // Determine which of the factorizable integrals can be done analytically
198  RooArgSet prodSet ;
199  numIntSet.add(intSet) ;
200 
201  compIter->Reset() ;
202  while((pdf=(RooAbsPdf*)compIter->Next())) {
203  if (doFactorize && pdf->dependsOn(anaIntSet)) {
204  RooArgSet anaSet ;
205  Int_t code = pdf->getAnalyticalIntegralWN(anaIntSet,anaSet,0,isetRangeName) ;
206  if (code!=0) {
207  // Analytical integral, create integral object
208  RooAbsReal* pai = pdf->createIntegral(anaSet,isetRangeName) ;
209  pai->setOperMode(_operMode) ;
210 
211  // Add to integral to product
212  prodSet.add(*pai) ;
213 
214  // Remove analytically integratable observables from numeric integration list
215  numIntSet.remove(anaSet) ;
216 
217  // Declare ownership of integral
218  saveSet.addOwned(*pai) ;
219  } else {
220  // Analytic integration of factorizable observable not possible, add straight pdf to product
221  prodSet.add(*pdf) ;
222  }
223  } else {
224  // Non-factorizable observables, add straight pdf to product
225  prodSet.add(*pdf) ;
226  }
227  }
228 
229  //cout << "RooGenProdProj::makeIntegral(" << GetName() << ") prodset = " << prodSet << endl ;
230 
231  // Create product of (partial) analytical integrals
232  TString prodName ;
233  if (isetRangeName) {
234  prodName = Form("%s_%s_Range[%s]",GetName(),name,isetRangeName) ;
235  } else {
236  prodName = Form("%s_%s",GetName(),name) ;
237  }
238  RooProduct* prod = new RooProduct(prodName,"product",prodSet) ;
240  prod->setOperMode(_operMode) ;
241 
242  // Declare owndership of product
243  saveSet.addOwned(*prod) ;
244 
245  // Create integral performing remaining numeric integration over (partial) analytic product
246  RooAbsReal* ret = prod->createIntegral(numIntSet,isetRangeName) ;
247 // cout << "verbose print of integral object" << endl ;
248 // ret->Print("v") ;
249  ret->setOperMode(_operMode) ;
250  saveSet.addOwned(*ret) ;
251 
252  delete compIter ;
253  delete intIter ;
254 
255  // Caller owners returned master integral object
256  return ret ;
257 }
258 
259 
260 
261 ////////////////////////////////////////////////////////////////////////////////
262 /// Calculate and return value of normalization projection
263 
265 {
266  Double_t nom = ((RooAbsReal*)_intList.at(0))->getVal() ;
267 
268  if (!_haveD) return nom ;
269 
270  Double_t den = ((RooAbsReal*)_intList.at(1))->getVal() ;
271 
272  //cout << "RooGenProdProj::eval(" << GetName() << ") nom = " << nom << " den = " << den << endl ;
273 
274  return nom / den ;
275 }
276 
277 
278 
279 ////////////////////////////////////////////////////////////////////////////////
280 /// Intercept cache mode operation changes and propagate them to the components
281 
283 {
284  // WVE use cache manager here!
285 
286  RooAbsArg* arg ;
288  while((arg=(RooAbsArg*)nIter->Next())) {
289  arg->setOperMode(_operMode) ;
290  }
291  delete nIter ;
292 
294  while((arg=(RooAbsArg*)dIter->Next())) {
295  arg->setOperMode(_operMode) ;
296  }
297  delete dIter ;
298 
300  if (_haveD) _intList.at(1)->setOperMode(Auto) ; // Denominator always stays in Auto mode (normalization integral)
301 }
302 
303 
304 
305 
306 
307 
308 
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
TIterator * createIterator(Bool_t dir=kIterForward) const
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:86
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:735
virtual void Reset()=0
Double_t getVal(const RooArgSet *set=0) const
Definition: RooAbsReal.h:64
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
STL namespace.
virtual void operModeHook()
Intercept cache mode operation changes and propagate them to the components.
Iterator abstract base class.
Definition: TIterator.h:30
virtual ~RooGenProdProj()
Destructor.
friend class RooArgSet
Definition: RooAbsArg.h:471
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Reimplementation of standard RooArgList::add()
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:501
RooAbsCollection * snapshot(Bool_t deepCopy=kTRUE) const
Take a snap shot of current collection contents: An owning collection is returned containing clones o...
RooExpensiveObjectCache & expensiveObjectCache() const
Definition: RooAbsArg.cxx:2333
RooAbsArg * at(Int_t idx) const
Definition: RooArgList.h:84
RooAbsArg * first() const
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:318
RooListProxy _intList
RooGenProdProj()
Default constructor.
char * Form(const char *fmt,...)
RooArgSet * _compSetOwnedN
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:90
OperMode _operMode
Definition: RooAbsArg.h:575
RooSetProxy _compSetD
TIterator * serverIterator() const
Definition: RooAbsArg.h:112
RooProduct a RooAbsReal implementation that represent the product of a given set of other RooAbsReal ...
Definition: RooProduct.h:32
const Bool_t kFALSE
Definition: RtypesCore.h:88
virtual void setExpensiveObjectCache(RooExpensiveObjectCache &cache)
Definition: RooAbsArg.h:499
RooArgSet * _compSetOwnedD
#define ClassImp(name)
Definition: Rtypes.h:359
RooAbsArg * find(const char *name) const
Find object with given name in list.
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
RooGenProdProj is an auxiliary class for RooProdPdf that calculates a general normalized projection o...
RooSetProxy _compSetN
virtual Bool_t remove(const RooAbsArg &var, Bool_t silent=kFALSE, Bool_t matchByNameOnly=kFALSE)
Remove the specified argument from our list.
RooAbsPdf is the abstract interface for all probability density functions The class provides hybrid a...
Definition: RooAbsPdf.h:41
RooAbsReal * makeIntegral(const char *name, const RooArgSet &compSet, const RooArgSet &intSet, RooArgSet &saveSet, const char *isetRangeName, Bool_t doFactorize)
Utility function to create integral over observables intSet in range isetRangeName over product of p...
virtual TObject * Next()=0
Double_t evaluate() const
Calculate and return value of normalization projection.
void setOperMode(OperMode mode, Bool_t recurseADirty=kTRUE)
Change cache operation mode to given mode.
Definition: RooAbsArg.cxx:1745
RooAbsArg is the common abstract base class for objects that represent a value (of arbitrary type) an...
Definition: RooAbsArg.h:66
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:385
const Bool_t kTRUE
Definition: RtypesCore.h:87
char name[80]
Definition: TGX11.cxx:109
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Overloaded RooArgSet::add() method inserts &#39;var&#39; into set and registers &#39;var&#39; as server to owner with...