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[out] saveSet All component objects needed to represent the product integral are added as owned members to saveSet.
171 /// \note The set owns new components that are created for the integral.
172 /// \param[in] isetRangeName Integral range.
173 /// \param[in] doFactorize
174 ///
175 /// \return A RooAbsReal object representing the requested integral. The object is owned by `saveSet`.
176 ///
177 /// The integration is factorized into components as much as possible and done analytically as far as possible.
178 RooAbsReal* 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  for (const auto arg : intSet) {
185  auto count = std::count_if(compSet.begin(), compSet.end(), [arg](const RooAbsArg* pdfAsArg){
186  auto pdf = static_cast<const RooAbsPdf*>(pdfAsArg);
187  return (pdf->dependsOn(*arg));
188  });
189 
190  if (count==1) {
191  anaIntSet.add(*arg) ;
192  }
193  }
194 
195  // Determine which of the factorizable integrals can be done analytically
196  RooArgSet prodSet ;
197  numIntSet.add(intSet) ;
198 
199  for (const auto pdfAsArg : compSet) {
200  auto pdf = static_cast<const RooAbsPdf*>(pdfAsArg);
201 
202  if (doFactorize && pdf->dependsOn(anaIntSet)) {
203  RooArgSet anaSet ;
204  Int_t code = pdf->getAnalyticalIntegralWN(anaIntSet,anaSet,0,isetRangeName) ;
205  if (code!=0) {
206  // Analytical integral, create integral object
207  RooAbsReal* pai = pdf->createIntegral(anaSet,isetRangeName) ;
208  pai->setOperMode(_operMode) ;
209 
210  // Add to integral to product
211  prodSet.add(*pai) ;
212 
213  // Remove analytically integratable observables from numeric integration list
214  numIntSet.remove(anaSet) ;
215 
216  // Declare ownership of integral
217  saveSet.addOwned(*pai) ;
218  } else {
219  // Analytic integration of factorizable observable not possible, add straight pdf to product
220  prodSet.add(*pdf) ;
221  }
222  } else {
223  // Non-factorizable observables, add straight pdf to product
224  prodSet.add(*pdf) ;
225  }
226  }
227 
228  // Create product of (partial) analytical integrals
229  TString prodName ;
230  if (isetRangeName) {
231  prodName = Form("%s_%s_Range[%s]",GetName(),name,isetRangeName) ;
232  } else {
233  prodName = Form("%s_%s",GetName(),name) ;
234  }
235 
236  // Create clones of the elements in prodSet. These need to be cloned
237  // because when caching optimisation lvl 2 is activated, pre-computed
238  // values are side-loaded into the elements.
239  // Those pre-cached values already contain normalisation constants, so
240  // the integral comes out wrongly. Therefore, we create here nodes that
241  // don't participate in any caching, which are used to compute integrals.
242  RooArgSet prodSetClone;
243  prodSet.snapshot(prodSetClone, false);
244  saveSet.addOwned(prodSetClone);
245  prodSetClone.releaseOwnership();
246 
247  RooProduct* prod = new RooProduct(prodName, "product", prodSetClone);
249  prod->setOperMode(_operMode) ;
250 
251  // Declare ownership of product
252  saveSet.addOwned(*prod);
253 
254  // Create integral performing remaining numeric integration over (partial) analytic product
255  RooAbsReal* ret = prod->createIntegral(numIntSet,isetRangeName) ;
256  ret->setOperMode(_operMode) ;
257  saveSet.addOwned(*ret) ;
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:140
RooGenProdProj::_compSetOwnedN
RooArgSet * _compSetOwnedN
Definition: RooGenProdProj.h:45
kTRUE
const Bool_t kTRUE
Definition: RtypesCore.h:100
RooAbsReal.h
RooAbsCollection::first
RooAbsArg * first() const
Definition: RooAbsCollection.h:236
RooFit.h
RooAbsArg::Auto
@ Auto
Definition: RooAbsArg.h:390
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:548
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:585
RooAbsCollection::find
RooAbsArg * find(const char *name) const
Find object with given name in list.
Definition: RooAbsCollection.cxx:810
RooAbsArg::setOperMode
void setOperMode(OperMode mode, Bool_t recurseADirty=kTRUE)
Set the operation mode of this node.
Definition: RooAbsArg.cxx:1817
RooAbsReal
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
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:72
RooAbsArg::setExpensiveObjectCache
virtual void setExpensiveObjectCache(RooExpensiveObjectCache &cache)
Definition: RooAbsArg.h:503
TString
Basic string class.
Definition: TString.h:136
RooAbsArg::expensiveObjectCache
RooExpensiveObjectCache & expensiveObjectCache() const
Definition: RooAbsArg.cxx:2277
RooGenProdProj
RooGenProdProj is an auxiliary class for RooProdPdf that calculates a general normalised projection o...
Definition: RooGenProdProj.h:26
bool
TIterator
Iterator abstract base class.
Definition: TIterator.h:30
RooGenProdProj.h
RooArgSet::snapshot
RooArgSet * snapshot(bool deepCopy=true) const
Use RooAbsCollection::snapshot(), but return as RooArgSet.
Definition: RooArgSet.h:115
RooProduct
A RooProduct represents the product of a given set of RooAbsReal objects.
Definition: RooProduct.h:29
RooAbsPdf.h
Double_t
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:439
RooAbsArg::_operMode
OperMode _operMode
Mark batches as dirty (only meaningful for RooAbsReal).
Definition: RooAbsArg.h:679
RooGenProdProj::_compSetOwnedD
RooArgSet * _compSetOwnedD
Definition: RooGenProdProj.h:46
kFALSE
const Bool_t kFALSE
Definition: RtypesCore.h:101
RooAbsCollection::createIterator
TIterator * createIterator(Bool_t dir=kIterForward) const
TIterator-style iteration over contained elements.
Definition: RooAbsCollection.h:178
RooAbsArg::RooArgSet
friend class RooArgSet
Definition: RooAbsArg.h:599
RooAbsCollection::add
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Add the specified argument to list.
Definition: RooAbsCollection.cxx:455
RooGenProdProj::~RooGenProdProj
virtual ~RooGenProdProj()
Destructor.
Definition: RooGenProdProj.cxx:157
RooProduct.h
RooAbsCollection::end
const_iterator end() const
Definition: RooAbsCollection.h:202
RooAbsCollection::addOwned
virtual Bool_t addOwned(RooAbsArg &var, Bool_t silent=kFALSE)
Add the specified argument to list.
Definition: RooAbsCollection.cxx:403
RooAbsCollection::releaseOwnership
void releaseOwnership()
Definition: RooAbsCollection.h:299
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:178
TIterator::Next
virtual TObject * Next()=0
RooGenProdProj::RooGenProdProj
RooGenProdProj()
Default constructor.
Definition: RooGenProdProj.cxx:54
RooAbsCollection::begin
const_iterator begin() const
Definition: RooAbsCollection.h:198
RooGenProdProj::_compSetN
RooSetProxy _compSetN
Definition: RooGenProdProj.h:47
RooGenProdProj::_intList
RooListProxy _intList
Definition: RooGenProdProj.h:49
name
char name[80]
Definition: TGX11.cxx:110
RooErrorHandler.h
RooAbsArg
RooAbsArg is the common abstract base class for objects that represent a value and a "shape" in RooFi...
Definition: RooAbsArg.h:72
RooListProxy::add
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE) override
Reimplementation of standard RooArgList::add()
Definition: RooListProxy.cxx:104
RooAbsPdf
Definition: RooAbsPdf.h:41
RooGenProdProj::_haveD
Bool_t _haveD
Definition: RooGenProdProj.h:50
TNamed::GetName
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
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
RooSetProxy::add
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE) override
Overloaded RooArgSet::add() method inserts 'var' into set and registers 'var' as server to owner with...
Definition: RooSetProxy.cxx:165
RooArgSet
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:33
int