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(false)
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 doFactorize) :
66 RooAbsReal(name, title),
67 _compSetOwnedN(0),
68 _compSetOwnedD(0),
69 _compSetN("compSetN","Set of integral components owned by numerator",this,false),
70 _compSetD("compSetD","Set of integral components owned by denominator",this,false),
71 _intList("intList","List of integrals",this,true),
72 _haveD(false)
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 = true ;
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 for(RooAbsArg * server : servers()) {
115 removeServer(*server,true) ;
116 }
117
118 // Copy constructor
121
124
125 for (RooAbsArg * arg : *_compSetOwnedN) {
126 arg->setOperMode(_operMode) ;
127 }
128 for (RooAbsArg * arg : *_compSetOwnedD) {
129 arg->setOperMode(_operMode) ;
130 }
131
132 // Fill _intList
133 _haveD = other._haveD ;
134 _intList.add(*_compSetN.find(other._intList.at(0)->GetName())) ;
135 if (other._haveD) {
136 _intList.add(*_compSetD.find(other._intList.at(1)->GetName())) ;
137 }
138}
139
140
141
142////////////////////////////////////////////////////////////////////////////////
143/// Destructor
144
146{
147 if (_compSetOwnedN) delete _compSetOwnedN ;
148 if (_compSetOwnedD) delete _compSetOwnedD ;
149}
150
151
152
153////////////////////////////////////////////////////////////////////////////////
154/// Utility function to create integral for product over certain observables.
155/// \param[in] name Name of integral to be created.
156/// \param[in] compSet All components of the product.
157/// \param[in] intSet Observables to be integrated.
158/// \param[out] saveSet All component objects needed to represent the product integral are added as owned members to saveSet.
159/// \note The set owns new components that are created for the integral.
160/// \param[in] isetRangeName Integral range.
161/// \param[in] doFactorize
162///
163/// \return A RooAbsReal object representing the requested integral. The object is owned by saveSet.
164///
165/// The integration is factorized into components as much as possible and done analytically as far as possible.
166RooAbsReal* RooGenProdProj::makeIntegral(const char* name, const RooArgSet& compSet, const RooArgSet& intSet,
167 RooArgSet& saveSet, const char* isetRangeName, bool doFactorize)
168{
169 RooArgSet anaIntSet, numIntSet ;
170
171 // First determine subset of observables in intSet that are factorizable
172 for (const auto arg : intSet) {
173 auto count = std::count_if(compSet.begin(), compSet.end(), [arg](const RooAbsArg* pdfAsArg){
174 auto pdf = static_cast<const RooAbsPdf*>(pdfAsArg);
175 return (pdf->dependsOn(*arg));
176 });
177
178 if (count==1) {
179 anaIntSet.add(*arg) ;
180 }
181 }
182
183 // Determine which of the factorizable integrals can be done analytically
184 RooArgSet prodSet ;
185 numIntSet.add(intSet) ;
186
187 for (const auto pdfAsArg : compSet) {
188 auto pdf = static_cast<const RooAbsPdf*>(pdfAsArg);
189
190 if (doFactorize && pdf->dependsOn(anaIntSet)) {
191 RooArgSet anaSet ;
192 Int_t code = pdf->getAnalyticalIntegralWN(anaIntSet,anaSet,0,isetRangeName) ;
193 if (code!=0) {
194 // Analytical integral, create integral object
195 RooAbsReal* pai = pdf->createIntegral(anaSet,isetRangeName) ;
196 pai->setOperMode(_operMode) ;
197
198 // Add to integral to product
199 prodSet.add(*pai) ;
200
201 // Remove analytically integratable observables from numeric integration list
202 numIntSet.remove(anaSet) ;
203
204 // Declare ownership of integral
205 saveSet.addOwned(*pai) ;
206 } else {
207 // Analytic integration of factorizable observable not possible, add straight pdf to product
208 prodSet.add(*pdf) ;
209 }
210 } else {
211 // Non-factorizable observables, add straight pdf to product
212 prodSet.add(*pdf) ;
213 }
214 }
215
216 // Create product of (partial) analytical integrals
217 TString prodName ;
218 if (isetRangeName) {
219 prodName = Form("%s_%s_Range[%s]",GetName(),name,isetRangeName) ;
220 } else {
221 prodName = Form("%s_%s",GetName(),name) ;
222 }
223
224 // Create clones of the elements in prodSet. These need to be cloned
225 // because when caching optimisation lvl 2 is activated, pre-computed
226 // values are side-loaded into the elements.
227 // Those pre-cached values already contain normalisation constants, so
228 // the integral comes out wrongly. Therefore, we create here nodes that
229 // don't participate in any caching, which are used to compute integrals.
230 RooArgSet prodSetClone;
231 prodSet.snapshot(prodSetClone, false);
232
233 auto prod = std::make_unique<RooProduct>(prodName, "product", prodSetClone);
234 prod->setExpensiveObjectCache(expensiveObjectCache()) ;
235 prod->setOperMode(_operMode) ;
236
237 // Create integral performing remaining numeric integration over (partial) analytic product
238 std::unique_ptr<RooAbsReal> integral{prod->createIntegral(numIntSet,isetRangeName)};
239 integral->setOperMode(_operMode) ;
240 auto ret = integral.get();
241
242 // Declare ownership of prodSet, product, and integral
243 saveSet.addOwned(std::move(prodSetClone));
244 saveSet.addOwned(std::move(prod));
245 saveSet.addOwned(std::move(integral)) ;
246
247
248 // Caller owners returned master integral object
249 return ret ;
250}
251
252
253
254////////////////////////////////////////////////////////////////////////////////
255/// Calculate and return value of normalization projection
256
258{
259 RooArgSet const* nset = _intList.nset();
260
261 double nom = static_cast<RooAbsReal*>(_intList.at(0))->getVal(nset);
262
263 if (!_haveD) return nom ;
264
265 double den = static_cast<RooAbsReal*>(_intList.at(1))->getVal(nset);
266
267 //cout << "RooGenProdProj::eval(" << GetName() << ") nom = " << nom << " den = " << den << endl ;
268
269 return nom / den ;
270}
271
272
273
274////////////////////////////////////////////////////////////////////////////////
275/// Intercept cache mode operation changes and propagate them to the components
276
278{
279 // WVE use cache manager here!
280
281 for(RooAbsArg * arg : *_compSetOwnedN) {
282 arg->setOperMode(_operMode) ;
283 }
284
285 for(RooAbsArg * arg : *_compSetOwnedD) {
286 arg->setOperMode(_operMode) ;
287 }
288
290 if (_haveD) _intList.at(1)->setOperMode(Auto) ; // Denominator always stays in Auto mode (normalization integral)
291}
292
293
294
295
296
297
298
#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:2456
RooAbsArg is the common abstract base class for objects that represent a value and a "shape" in RooFi...
Definition: RooAbsArg.h:71
RooExpensiveObjectCache & expensiveObjectCache() const
Definition: RooAbsArg.cxx:2271
void setOperMode(OperMode mode, bool recurseADirty=true)
Set the operation mode of this node.
Definition: RooAbsArg.cxx:1866
void removeServer(RooAbsArg &server, bool force=false)
Unregister another RooAbsArg as a server to us, ie, declare that we no longer depend on its value and...
Definition: RooAbsArg.cxx:402
virtual void setExpensiveObjectCache(RooExpensiveObjectCache &cache)
Definition: RooAbsArg.h:503
const RefCountList_t & servers() const
List of all servers of this object.
Definition: RooAbsArg.h:198
OperMode _operMode
Definition: RooAbsArg.h:712
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.
const_iterator end() const
RooAbsArg * first() const
virtual bool addOwned(RooAbsArg &var, bool silent=false)
Add an argument and transfer the ownership to the collection.
const_iterator begin() const
RooAbsArg * find(const char *name) const
Find object with given name in list.
const RooArgSet * nset() const
Definition: RooAbsProxy.h:48
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:62
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:91
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 std::liste...
Definition: RooAbsReal.cxx:522
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:56
RooArgSet * snapshot(bool deepCopy=true) const
Use RooAbsCollection::snapshot(), but return as RooArgSet.
Definition: RooArgSet.h:179
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.
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.
~RooGenProdProj() override
Destructor.
RooSetProxy _compSetD
Set proxy for denominator components.
RooListProxy _intList
Master integrals representing numerator and denominator.
RooArgSet * _compSetOwnedD
Owner of denominator components.
bool _haveD
Do we have a denominator term?
RooArgSet * _compSetOwnedN
Owner of numerator components.
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:136