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 
23 RooGenProdProj is an auxiliary class for RooProdPdf that calculates
24 a 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 
29 Partial integrals, which factorise and can be calculated, are calculated
30 analytically. 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 
45 using 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 
66 RooGenProdProj::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 
177 RooAbsReal* 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 
RooAbsArg::serverIterator
TIterator * serverIterator() const
Definition: RooAbsArg.h:141
RooSetProxy::add
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...
Definition: RooSetProxy.cxx:165
RooGenProdProj::_compSetOwnedN
RooArgSet * _compSetOwnedN
Definition: RooGenProdProj.h:45
kTRUE
const Bool_t kTRUE
Definition: RtypesCore.h:91
RooAbsReal.h
RooAbsCollection::first
RooAbsArg * first() const
Definition: RooAbsCollection.h:176
RooFit.h
RooAbsArg::Auto
@ Auto
Definition: RooAbsArg.h:491
RooAbsReal::createIntegral
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:561
ClassImp
#define ClassImp(name)
Definition: Rtypes.h:364
Form
char * Form(const char *fmt,...)
RooAbsReal::getVal
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:91
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
RooAbsCollection::find
RooAbsArg * find(const char *name) const
Find object with given name in list.
Definition: RooAbsCollection.cxx:813
RooAbsArg::setOperMode
void setOperMode(OperMode mode, Bool_t recurseADirty=kTRUE)
Change cache operation mode to given mode.
Definition: RooAbsArg.cxx:1720
RooArgSet::add
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
RooAbsReal
Definition: RooAbsReal.h:61
RooArgList::at
RooAbsArg * at(Int_t idx) const
Return object at given index, or nullptr if index is out of range.
Definition: RooArgList.h:88
RooAbsArg::setExpensiveObjectCache
virtual void setExpensiveObjectCache(RooExpensiveObjectCache &cache)
Definition: RooAbsArg.h:528
TString
Definition: TString.h:136
RooAbsArg::expensiveObjectCache
RooExpensiveObjectCache & expensiveObjectCache() const
Definition: RooAbsArg.cxx:2170
RooGenProdProj
Definition: RooGenProdProj.h:26
bool
TIterator
Definition: TIterator.h:30
RooArgSet::addOwned
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
RooGenProdProj.h
RooArgSet::snapshot
RooArgSet * snapshot(bool deepCopy=true) const
Use RooAbsCollection::snapshot(), but return as RooArgSet.
Definition: RooArgSet.h:134
RooProduct
Definition: RooProduct.h:30
RooAbsPdf.h
RooAbsArg::removeServer
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:407
RooAbsArg::_operMode
OperMode _operMode
Mark batches as dirty (only meaningful for RooAbsReal).
Definition: RooAbsArg.h:663
RooGenProdProj::_compSetOwnedD
RooArgSet * _compSetOwnedD
Definition: RooGenProdProj.h:46
kFALSE
const Bool_t kFALSE
Definition: RtypesCore.h:92
RooAbsCollection::createIterator
TIterator * createIterator(Bool_t dir=kIterForward) const
TIterator-style iteration over contained elements.
Definition: RooAbsCollection.h:118
RooAbsArg::RooArgSet
friend class RooArgSet
Definition: RooAbsArg.h:579
RooListProxy::add
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Reimplementation of standard RooArgList::add()
Definition: RooListProxy.cxx:104
RooGenProdProj::~RooGenProdProj
virtual ~RooGenProdProj()
Destructor.
Definition: RooGenProdProj.cxx:157
RooProduct.h
RooGenProdProj::makeIntegral
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.
Definition: RooGenProdProj.cxx:177
TIterator::Next
virtual TObject * Next()=0
TIterator::Reset
virtual void Reset()=0
RooGenProdProj::RooGenProdProj
RooGenProdProj()
Default constructor.
Definition: RooGenProdProj.cxx:54
RooAbsReal::getAnalyticalIntegralWN
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:372
Double_t
double Double_t
Definition: RtypesCore.h:59
RooGenProdProj::_compSetN
RooSetProxy _compSetN
Definition: RooGenProdProj.h:47
RooGenProdProj::_intList
RooListProxy _intList
Definition: RooGenProdProj.h:49
name
char name[80]
Definition: TGX11.cxx:110
RooAbsArg::dependsOn
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:732
RooErrorHandler.h
RooAbsArg
Definition: RooAbsArg.h:73
RooAbsPdf
Definition: RooAbsPdf.h:40
RooGenProdProj::_haveD
Bool_t _haveD
Definition: RooGenProdProj.h:50
TNamed::GetName
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:53
RooGenProdProj::_compSetD
RooSetProxy _compSetD
Definition: RooGenProdProj.h:48
RooGenProdProj::evaluate
Double_t evaluate() const
Calculate and return value of normalization projection.
Definition: RooGenProdProj.cxx:268
RooGenProdProj::operModeHook
virtual void operModeHook()
Intercept cache mode operation changes and propagate them to the components.
Definition: RooGenProdProj.cxx:286
Riostream.h
RooArgSet
Definition: RooArgSet.h:28
int