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 "Riostream.h"
35#include <math.h>
36
37#include "RooGenProdProj.h"
38#include "RooAbsReal.h"
39#include "RooAbsPdf.h"
40#include "RooErrorHandler.h"
41#include "RooProduct.h"
42
43using namespace std;
44
46
47
48
49////////////////////////////////////////////////////////////////////////////////
50/// Default constructor
51
53 _compSetOwnedN(0),
54 _compSetOwnedD(0),
55 _haveD(kFALSE)
56{
57}
58
59
60////////////////////////////////////////////////////////////////////////////////
61/// Constructor for a normalization projection of the product of p.d.f.s _prodSet
62/// integrated over _intSet in range isetRangeName while normalized over _normSet
63
64RooGenProdProj::RooGenProdProj(const char *name, const char *title, const RooArgSet& _prodSet, const RooArgSet& _intSet,
65 const RooArgSet& _normSet, const char* isetRangeName, const char* normRangeName, Bool_t doFactorize) :
66 RooAbsReal(name, title),
67 _compSetOwnedN(0),
68 _compSetOwnedD(0),
69 _compSetN("compSetN","Set of integral components owned by numerator",this,kFALSE),
70 _compSetD("compSetD","Set of integral components owned by denominator",this,kFALSE),
71 _intList("intList","List of integrals",this,kTRUE),
72 _haveD(kFALSE)
73{
74 // Set expensive object cache to that of first item in prodSet
76
77 // Create owners of components created in ctor
80
81 RooAbsReal* numerator = makeIntegral("numerator",_prodSet,_intSet,*_compSetOwnedN,isetRangeName,doFactorize) ;
82 RooAbsReal* denominator = makeIntegral("denominator",_prodSet,_normSet,*_compSetOwnedD,normRangeName,doFactorize) ;
83
84// cout << "RooGenProdPdf::ctor(" << GetName() << ") numerator = " << numerator->GetName() << endl ;
85// numerator->printComponentTree() ;
86// cout << "RooGenProdPdf::ctor(" << GetName() << ") denominator = " << denominator->GetName() << endl ;
87// denominator->printComponentTree() ;
88
89 // Copy all components in (non-owning) set proxy
92
93 _intList.add(*numerator) ;
94 if (denominator) {
95 _intList.add(*denominator) ;
96 _haveD = kTRUE ;
97 }
98}
99
100
101
102////////////////////////////////////////////////////////////////////////////////
103/// Copy constructor
104
106 RooAbsReal(other, name),
107 _compSetOwnedN(0),
108 _compSetOwnedD(0),
109 _compSetN("compSetN","Set of integral components owned by numerator",this),
110 _compSetD("compSetD","Set of integral components owned by denominator",this),
111 _intList("intList","List of integrals",this)
112{
113 // Explicitly remove all server links at this point
114 TIterator* iter = serverIterator() ;
115 RooAbsArg* server ;
116 while((server=(RooAbsArg*)iter->Next())) {
117 removeServer(*server,kTRUE) ;
118 }
119 delete iter ;
120
121 // Copy constructor
124
127
128 RooAbsArg* arg ;
130 while((arg=(RooAbsArg*)nIter->Next())) {
131// cout << "ownedN elem " << arg->GetName() << "(" << arg << ")" << endl ;
132 arg->setOperMode(_operMode) ;
133 }
134 delete nIter ;
136 while((arg=(RooAbsArg*)dIter->Next())) {
137// cout << "ownedD elem " << arg->GetName() << "(" << arg << ")" << endl ;
138 arg->setOperMode(_operMode) ;
139 }
140 delete dIter ;
141
142 // Fill _intList
143 _haveD = other._haveD ;
144 _intList.add(*_compSetN.find(other._intList.at(0)->GetName())) ;
145 if (other._haveD) {
146 _intList.add(*_compSetD.find(other._intList.at(1)->GetName())) ;
147 }
148}
149
150
151
152////////////////////////////////////////////////////////////////////////////////
153/// Destructor
154
156{
157 if (_compSetOwnedN) delete _compSetOwnedN ;
158 if (_compSetOwnedD) delete _compSetOwnedD ;
159}
160
161
162
163////////////////////////////////////////////////////////////////////////////////
164/// Utility function to create integral for product over certain observables.
165/// \param[in] name Name of integral to be created.
166/// \param[in] compSet All components of the product.
167/// \param[in] intSet Observables to be integrated.
168/// \param[out] saveSet All component objects needed to represent the product integral are added as owned members to saveSet.
169/// \note The set owns new components that are created for the integral.
170/// \param[in] isetRangeName Integral range.
171/// \param[in] doFactorize
172///
173/// \return A RooAbsReal object representing the requested integral. The object is owned by `saveSet`.
174///
175/// The integration is factorized into components as much as possible and done analytically as far as possible.
176RooAbsReal* RooGenProdProj::makeIntegral(const char* name, const RooArgSet& compSet, const RooArgSet& intSet,
177 RooArgSet& saveSet, const char* isetRangeName, Bool_t doFactorize)
178{
179 RooArgSet anaIntSet, numIntSet ;
180
181 // First determine subset of observables in intSet that are factorizable
182 for (const auto arg : intSet) {
183 auto count = std::count_if(compSet.begin(), compSet.end(), [arg](const RooAbsArg* pdfAsArg){
184 auto pdf = static_cast<const RooAbsPdf*>(pdfAsArg);
185 return (pdf->dependsOn(*arg));
186 });
187
188 if (count==1) {
189 anaIntSet.add(*arg) ;
190 }
191 }
192
193 // Determine which of the factorizable integrals can be done analytically
194 RooArgSet prodSet ;
195 numIntSet.add(intSet) ;
196
197 for (const auto pdfAsArg : compSet) {
198 auto pdf = static_cast<const RooAbsPdf*>(pdfAsArg);
199
200 if (doFactorize && pdf->dependsOn(anaIntSet)) {
201 RooArgSet anaSet ;
202 Int_t code = pdf->getAnalyticalIntegralWN(anaIntSet,anaSet,0,isetRangeName) ;
203 if (code!=0) {
204 // Analytical integral, create integral object
205 RooAbsReal* pai = pdf->createIntegral(anaSet,isetRangeName) ;
206 pai->setOperMode(_operMode) ;
207
208 // Add to integral to product
209 prodSet.add(*pai) ;
210
211 // Remove analytically integratable observables from numeric integration list
212 numIntSet.remove(anaSet) ;
213
214 // Declare ownership of integral
215 saveSet.addOwned(*pai) ;
216 } else {
217 // Analytic integration of factorizable observable not possible, add straight pdf to product
218 prodSet.add(*pdf) ;
219 }
220 } else {
221 // Non-factorizable observables, add straight pdf to product
222 prodSet.add(*pdf) ;
223 }
224 }
225
226 // Create product of (partial) analytical integrals
227 TString prodName ;
228 if (isetRangeName) {
229 prodName = Form("%s_%s_Range[%s]",GetName(),name,isetRangeName) ;
230 } else {
231 prodName = Form("%s_%s",GetName(),name) ;
232 }
233
234 // Create clones of the elements in prodSet. These need to be cloned
235 // because when caching optimisation lvl 2 is activated, pre-computed
236 // values are side-loaded into the elements.
237 // Those pre-cached values already contain normalisation constants, so
238 // the integral comes out wrongly. Therefore, we create here nodes that
239 // don't participate in any caching, which are used to compute integrals.
240 RooArgSet prodSetClone;
241 prodSet.snapshot(prodSetClone, false);
242
243 auto prod = std::make_unique<RooProduct>(prodName, "product", prodSetClone);
244 prod->setExpensiveObjectCache(expensiveObjectCache()) ;
245 prod->setOperMode(_operMode) ;
246
247 // Create integral performing remaining numeric integration over (partial) analytic product
248 std::unique_ptr<RooAbsReal> integral{prod->createIntegral(numIntSet,isetRangeName)};
249 integral->setOperMode(_operMode) ;
250 auto ret = integral.get();
251
252 // Declare ownership of prodSet, product, and integral
253 saveSet.addOwned(std::move(prodSetClone));
254 saveSet.addOwned(std::move(prod));
255 saveSet.addOwned(std::move(integral)) ;
256
257
258 // Caller owners returned master integral object
259 return ret ;
260}
261
262
263
264////////////////////////////////////////////////////////////////////////////////
265/// Calculate and return value of normalization projection
266
268{
269 Double_t nom = ((RooAbsReal*)_intList.at(0))->getVal() ;
270
271 if (!_haveD) return nom ;
272
273 Double_t den = ((RooAbsReal*)_intList.at(1))->getVal() ;
274
275 //cout << "RooGenProdProj::eval(" << GetName() << ") nom = " << nom << " den = " << den << endl ;
276
277 return nom / den ;
278}
279
280
281
282////////////////////////////////////////////////////////////////////////////////
283/// Intercept cache mode operation changes and propagate them to the components
284
286{
287 // WVE use cache manager here!
288
289 RooAbsArg* arg ;
291 while((arg=(RooAbsArg*)nIter->Next())) {
292 arg->setOperMode(_operMode) ;
293 }
294 delete nIter ;
295
297 while((arg=(RooAbsArg*)dIter->Next())) {
298 arg->setOperMode(_operMode) ;
299 }
300 delete dIter ;
301
303 if (_haveD) _intList.at(1)->setOperMode(Auto) ; // Denominator always stays in Auto mode (normalization integral)
304}
305
306
307
308
309
310
311
const Bool_t kFALSE
Definition: RtypesCore.h:101
const Bool_t kTRUE
Definition: RtypesCore.h:100
#define ClassImp(name)
Definition: Rtypes.h:375
char name[80]
Definition: TGX11.cxx:110
char * Form(const char *fmt,...)
Formats a string in a circular formatting buffer.
Definition: TString.cxx:2447
RooAbsArg is the common abstract base class for objects that represent a value and a "shape" in RooFi...
Definition: RooAbsArg.h:78
RooExpensiveObjectCache & expensiveObjectCache() const
Definition: RooAbsArg.cxx:2354
friend class RooArgSet
Definition: RooAbsArg.h:651
virtual void setExpensiveObjectCache(RooExpensiveObjectCache &cache)
Definition: RooAbsArg.h:527
void setOperMode(OperMode mode, Bool_t recurseADirty=kTRUE)
Set the operation mode of this node.
Definition: RooAbsArg.cxx:1876
OperMode _operMode
Mark batches as dirty (only meaningful for RooAbsReal).
Definition: RooAbsArg.h:738
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:430
TIterator * serverIterator() const
Definition: RooAbsArg.h:146
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Add the specified argument to list.
const_iterator end() const
virtual Bool_t addOwned(RooAbsArg &var, Bool_t silent=kFALSE)
Add an argument and transfer the ownership to the collection.
RooAbsArg * first() const
const_iterator begin() 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:63
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:553
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:93
RooAbsArg * at(Int_t idx) const
Return object at given index, or nullptr if index is out of range.
Definition: RooArgList.h:110
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:57
RooArgSet * snapshot(bool deepCopy=true) const
Use RooAbsCollection::snapshot(), but return as RooArgSet.
Definition: RooArgSet.h:180
bool add(const RooAbsArg &var, bool valueServer, bool shapeServer, bool silent)
Overloaded RooCollection_t::add() method insert object into set and registers object as server to own...
RooGenProdProj is an auxiliary class for RooProdPdf that calculates a general normalised projection o...
void operModeHook() override
Intercept cache mode operation changes and propagate them to the components.
RooGenProdProj()
Default constructor.
~RooGenProdProj() override
Destructor.
RooSetProxy _compSetD
Set proxy for denominator components.
RooListProxy _intList
Master integrals representing numerator and denominator.
Bool_t _haveD
Do we have a denominator term?
RooArgSet * _compSetOwnedD
Owner of denominator components.
Double_t evaluate() const override
Calculate and return value of normalization projection.
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
Owner of numerator components.
RooSetProxy _compSetN
Set proxy for numerator components.
Iterator abstract base class.
Definition: TIterator.h:30
virtual TObject * Next()=0
const char * GetName() const override
Returns name of object.
Definition: TNamed.h:47
Basic string class.
Definition: TString.h:136