Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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 <cmath>
36
37#include "RooGenProdProj.h"
38#include "RooAbsReal.h"
39#include "RooAbsPdf.h"
40#include "RooErrorHandler.h"
41#include "RooProduct.h"
42
43
44////////////////////////////////////////////////////////////////////////////////
45/// Constructor for a normalization projection of the product of p.d.f.s _prodSet
46/// integrated over _intSet in range isetRangeName while normalized over _normSet
47
48RooGenProdProj::RooGenProdProj(const char *name, const char *title, const RooArgSet& _prodSet, const RooArgSet& _intSet,
49 const RooArgSet& _normSet, const char* isetRangeName, const char* normRangeName, bool doFactorize) :
50 RooAbsReal(name, title),
51 _compSetN("compSetN","Set of integral components owned by numerator",this,false),
52 _compSetD("compSetD","Set of integral components owned by denominator",this,false),
53 _intList("intList","List of integrals",this,true)
54{
55 // Set expensive object cache to that of first item in prodSet
57
58 // Create owners of components created in constructor
59 _compSetOwnedN = std::make_unique<RooArgSet>();
60 _compSetOwnedD = std::make_unique<RooArgSet>();
61
62 RooAbsReal* numerator = makeIntegral("numerator",_prodSet,_intSet,*_compSetOwnedN,isetRangeName,doFactorize) ;
63 RooAbsReal* denominator = makeIntegral("denominator",_prodSet,_normSet,*_compSetOwnedD,normRangeName,doFactorize) ;
64
65// cout << "RooGenProdPdf::ctor(" << GetName() << ") numerator = " << numerator->GetName() << endl ;
66// numerator->printComponentTree() ;
67// cout << "RooGenProdPdf::ctor(" << GetName() << ") denominator = " << denominator->GetName() << endl ;
68// denominator->printComponentTree() ;
69
70 // Copy all components in (non-owning) set proxy
73
74 _intList.add(*numerator) ;
75 if (denominator) {
76 _intList.add(*denominator) ;
77 _haveD = true ;
78 }
79}
80
81
82
83////////////////////////////////////////////////////////////////////////////////
84/// Copy constructor
85
87 : RooAbsReal(other, name),
88 _compSetN("compSetN", "Set of integral components owned by numerator", this),
89 _compSetD("compSetD", "Set of integral components owned by denominator", this),
90 _intList("intList", "List of integrals", this),
91 _haveD(other._haveD)
92{
93 // Copy constructor
94 _compSetOwnedN = std::make_unique<RooArgSet>();
97
98 _compSetOwnedD = std::make_unique<RooArgSet>();
101
102 for (RooAbsArg * arg : *_compSetOwnedN) {
103 arg->setOperMode(_operMode) ;
104 }
105 for (RooAbsArg * arg : *_compSetOwnedD) {
106 arg->setOperMode(_operMode) ;
107 }
108
109 // Fill _intList
110
111 _intList.add(*_compSetN.find(other._intList.at(0)->GetName())) ;
112 if (other._haveD) {
113 _intList.add(*_compSetD.find(other._intList.at(1)->GetName())) ;
114 }
115}
116
117
118////////////////////////////////////////////////////////////////////////////////
119/// Utility function to create integral for product over certain observables.
120/// \param[in] name Name of integral to be created.
121/// \param[in] compSet All components of the product.
122/// \param[in] intSet Observables to be integrated.
123/// \param[out] saveSet All component objects needed to represent the product integral are added as owned members to saveSet.
124/// \note The set owns new components that are created for the integral.
125/// \param[in] isetRangeName Integral range.
126/// \param[in] doFactorize
127///
128/// \return A RooAbsReal object representing the requested integral. The object is owned by `saveSet`.
129///
130/// The integration is factorized into components as much as possible and done analytically as far as possible.
131RooAbsReal* RooGenProdProj::makeIntegral(const char* name, const RooArgSet& compSet, const RooArgSet& intSet,
132 RooArgSet& saveSet, const char* isetRangeName, bool doFactorize)
133{
134 RooArgSet anaIntSet;
135 RooArgSet numIntSet;
136
137 // First determine subset of observables in intSet that are factorizable
138 for (const auto arg : intSet) {
139 auto count = std::count_if(compSet.begin(), compSet.end(), [arg](const RooAbsArg* pdfAsArg){
140 auto pdf = static_cast<const RooAbsPdf*>(pdfAsArg);
141 return (pdf->dependsOn(*arg));
142 });
143
144 if (count==1) {
145 anaIntSet.add(*arg) ;
146 }
147 }
148
149 // Determine which of the factorizable integrals can be done analytically
150 RooArgSet prodSet ;
151 numIntSet.add(intSet) ;
152
153 // The idea of the RooGenProdProj is that we divide two integral objects each
154 // created with this makeIntegral() function to get the normalized integral of
155 // a product. Therefore, we don't need to normalize the numerater and
156 // denominator integrals themselves. Doing the normalization would be
157 // expensive and it would cancel out anyway. However, if we don't specify an
158 // explicit normalization integral in createIntegral(), the last-used
159 // normalization set might be used to normalize the pdf, resulting in
160 // redundant computations.
161 //
162 // For this reason, the normalization set of the integrated pdfs is fixed to
163 // an empty set in this case. Note that in RooFit, a nullptr normalization
164 // set and an empty normalization set is not equivalent. The former implies
165 // taking the last-used normalization set, and the latter means explicitly no
166 // normalization.
167 RooArgSet emptyNormSet{};
168
169 RooArgSet keepAlive;
170
171 for (const auto pdfAsArg : compSet) {
172 auto pdf = static_cast<const RooAbsPdf*>(pdfAsArg);
173
174 if (doFactorize && pdf->dependsOn(anaIntSet)) {
175 RooArgSet anaSet ;
176 Int_t code = pdf->getAnalyticalIntegralWN(anaIntSet,anaSet,nullptr,isetRangeName) ;
177 if (code!=0) {
178 // Analytical integral, create integral object
179 std::unique_ptr<RooAbsReal> pai{pdf->createIntegral(anaSet,emptyNormSet,isetRangeName)};
180 pai->setOperMode(_operMode) ;
181
182 // Add to integral to product
183 prodSet.add(*pai) ;
184
185 // Remove analytically integratable observables from numeric integration list
186 numIntSet.remove(anaSet) ;
187
188 // Keep integral alive until the prodSet is cloned later
189 keepAlive.addOwned(std::move(pai));
190 } else {
191 // Analytic integration of factorizable observable not possible, add straight pdf to product
192 prodSet.add(*pdf) ;
193 }
194 } else {
195 // Non-factorizable observables, add straight pdf to product
196 prodSet.add(*pdf) ;
197 }
198 }
199
200 // Create product of (partial) analytical integrals
201 TString prodName ;
202 if (isetRangeName) {
203 prodName = Form("%s_%s_Range[%s]",GetName(),name,isetRangeName) ;
204 } else {
205 prodName = Form("%s_%s",GetName(),name) ;
206 }
207
208 // Create clones of the elements in prodSet. These need to be cloned
209 // because when caching optimisation lvl 2 is activated, pre-computed
210 // values are side-loaded into the elements.
211 // Those pre-cached values already contain normalisation constants, so
212 // the integral comes out wrongly. Therefore, we create here nodes that
213 // don't participate in any caching, which are used to compute integrals.
214 RooArgSet prodSetClone;
215 prodSet.snapshot(prodSetClone, false);
216
217 auto prod = std::make_unique<RooProduct>(prodName, "product", prodSetClone);
218 prod->setExpensiveObjectCache(expensiveObjectCache()) ;
219 prod->setOperMode(_operMode) ;
220
221 // Create integral performing remaining numeric integration over (partial) analytic product
222 std::unique_ptr<RooAbsReal> integral{prod->createIntegral(numIntSet,emptyNormSet,isetRangeName)};
223 integral->setOperMode(_operMode) ;
224 auto ret = integral.get();
225
226 // Declare ownership of prodSet, product, and integral
227 saveSet.addOwned(std::move(prodSetClone));
228 saveSet.addOwned(std::move(prod));
229 saveSet.addOwned(std::move(integral)) ;
230
231
232 // Caller owners returned master integral object
233 return ret ;
234}
235
236
237
238////////////////////////////////////////////////////////////////////////////////
239/// Calculate and return value of normalization projection
240
242{
243 RooArgSet const* nset = _intList.nset();
244
245 double nom = static_cast<RooAbsReal*>(_intList.at(0))->getVal(nset);
246
247 if (!_haveD) return nom ;
248
249 double den = static_cast<RooAbsReal*>(_intList.at(1))->getVal(nset);
250
251 //cout << "RooGenProdProj::eval(" << GetName() << ") nom = " << nom << " den = " << den << endl ;
252
253 return nom / den ;
254}
255
256
257
258////////////////////////////////////////////////////////////////////////////////
259/// Intercept cache mode operation changes and propagate them to the components
260
262{
263 // WVE use cache manager here!
264
265 for(RooAbsArg * arg : *_compSetOwnedN) {
266 arg->setOperMode(_operMode) ;
267 }
268
269 for(RooAbsArg * arg : *_compSetOwnedD) {
270 arg->setOperMode(_operMode) ;
271 }
272
274 if (_haveD) _intList.at(1)->setOperMode(Auto) ; // Denominator always stays in Auto mode (normalization integral)
275}
RooArgSet * _normSet
Pointer to set with observables used for normalization.
char name[80]
Definition TGX11.cxx:110
char * Form(const char *fmt,...)
Formats a string in a circular formatting buffer.
Definition TString.cxx:2489
Common abstract base class for objects that represent a value and a "shape" in RooFit.
Definition RooAbsArg.h:77
RooExpensiveObjectCache & expensiveObjectCache() const
void setOperMode(OperMode mode, bool recurseADirty=true)
Set the operation mode of this node.
virtual void setExpensiveObjectCache(RooExpensiveObjectCache &cache)
Definition RooAbsArg.h:501
OperMode _operMode
Definition RooAbsArg.h:714
virtual bool remove(const RooAbsArg &var, bool silent=false, bool matchByNameOnly=false)
Remove the specified argument from our list.
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to list.
TIterator begin()
RooAbsArg * first() const
virtual bool addOwned(RooAbsArg &var, bool silent=false)
Add an argument and transfer the ownership to the collection.
TIterator end() and range-based for loops.")
RooAbsArg * find(const char *name) const
Find object with given name in list.
Abstract interface for all probability density functions.
Definition RooAbsPdf.h:40
const RooArgSet * nset() const
Definition RooAbsProxy.h:52
Abstract base class for objects that represent a real value and implements functionality common to al...
Definition RooAbsReal.h:59
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition RooAbsReal.h:103
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:55
RooArgSet * snapshot(bool deepCopy=true) const
Use RooAbsCollection::snapshot(), but return as RooArgSet.
Definition RooArgSet.h:178
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...
General form of projected integral of product of PDFs, utility class for RooProdPdf.
void operModeHook() override
Intercept cache mode operation changes and propagate them to the components.
double 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 doFactorize)
Utility function to create integral for product over certain observables.
RooSetProxy _compSetD
Set proxy for denominator components.
RooListProxy _intList
Master integrals representing numerator and denominator.
std::unique_ptr< RooArgSet > _compSetOwnedN
Owner of numerator components.
std::unique_ptr< RooArgSet > _compSetOwnedD
Owner of denominator components.
bool _haveD
Do we have a denominator term?
RooGenProdProj(const char *name, const char *title, const RooArgSet &_prodSet, const RooArgSet &_intSet, const RooArgSet &_normSet, const char *isetRangeName, const char *normRangeName=nullptr, bool doFactorize=true)
Constructor for a normalization projection of the product of p.d.f.s _prodSet integrated over _intSet...
RooSetProxy _compSetN
Set proxy for numerator components.
const char * GetName() const override
Returns name of object.
Definition TNamed.h:47
Basic string class.
Definition TString.h:139