Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
RooProdGenContext.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 RooProdGenContext.cxx
19\class RooProdGenContext
20\ingroup Roofitcore
21
22Efficient implementation of the generator context
23specific for RooProdPdf PDFs. The sim-context owns a list of
24component generator contexts that are used to generate the dependents
25for each component PDF sequentially.
26**/
27
28#include "Riostream.h"
29#include "RooMsgService.h"
30
31#include "RooProdGenContext.h"
32#include "RooProdPdf.h"
33#include "RooDataSet.h"
34#include "RooRealVar.h"
35#include "RooGlobalFunc.h"
36
37
38
39using namespace std;
40
42
43
44////////////////////////////////////////////////////////////////////////////////
45
47 const RooDataSet *prototype, const RooArgSet* auxProto, bool verbose) :
48 RooAbsGenContext(model,vars,prototype,auxProto,verbose), _pdf(&model)
49{
50 // Constructor of optimization generator context for RooProdPdf objects
51
52 //Build an array of generator contexts for each product component PDF
53 cxcoutI(Generation) << "RooProdGenContext::ctor() setting up event special generator context for product p.d.f. " << model.GetName()
54 << " for generation of observable(s) " << vars ;
55 if (prototype) ccxcoutI(Generation) << " with prototype data for " << *prototype->get() ;
56 if (auxProto && !auxProto->empty()) ccxcoutI(Generation) << " with auxiliary prototypes " << *auxProto ;
57 ccxcoutI(Generation) << endl ;
58
59 // Make full list of dependents (generated & proto)
60 RooArgSet deps(vars) ;
61 if (prototype) {
62 deps.remove(*std::unique_ptr<RooArgSet>{model.getObservables(*prototype->get())},true,true) ;
63 }
64
65 // Factorize product in irreducible terms
66 RooLinkedList termList;
67 RooLinkedList depsList;
68 RooLinkedList impDepList;
69 RooLinkedList crossDepList;
70 RooLinkedList intList;
71 model.factorizeProduct(deps,RooArgSet(),termList,depsList,impDepList,crossDepList,intList) ;
72
73 if (dologD(Generation)) {
74 cxcoutD(Generation) << "RooProdGenContext::ctor() factorizing product expression in irriducible terms " ;
75 for(auto * t : static_range_cast<RooArgSet*>(termList)) {
76 ccxcoutD(Generation) << *t ;
77 }
78 ccxcoutD(Generation) << std::endl;
79 }
80
81 RooArgSet genDeps ;
82 // First add terms that do not import observables
83
84 bool anyAction = true ;
85 bool go=true ;
86 while(go) {
87
88 auto termIter = termList.begin();
89 auto impIter = impDepList.begin();
90 auto normIter = depsList.begin();
91
92 bool anyPrevAction=anyAction ;
93 anyAction=false ;
94
95 if (termList.empty()) {
96 break ;
97 }
98
99 while(termIter != termList.end()) {
100
101 auto * term = static_cast<RooArgSet*>(*termIter);
102 auto * impDeps = static_cast<RooArgSet*>(*impIter);
103 auto * termDeps = static_cast<RooArgSet*>(*normIter);
104 if (impDeps==nullptr || termDeps==nullptr) {
105 break ;
106 }
107
108 cxcoutD(Generation) << "RooProdGenContext::ctor() analyzing product term " << *term << " with observable(s) " << *termDeps ;
109 if (!impDeps->empty()) {
110 ccxcoutD(Generation) << " which has dependence of external observable(s) " << *impDeps << " that to be generated first by other terms" ;
111 }
112 ccxcoutD(Generation) << std::endl;
113
114 // Add this term if we have no imported dependents, or imported dependents are already generated
115 RooArgSet neededDeps(*impDeps) ;
116 neededDeps.remove(genDeps,true,true) ;
117
118 if (!neededDeps.empty()) {
119 if (!anyPrevAction) {
120 cxcoutD(Generation) << "RooProdGenContext::ctor() no convergence in single term analysis loop, terminating loop and process remainder of terms as single unit " << endl ;
121 go=false ;
122 break ;
123 }
124 cxcoutD(Generation) << "RooProdGenContext::ctor() skipping this term for now because it needs imported dependents that are not generated yet" << endl ;
125 ++termIter;
126 ++impIter;
127 ++normIter;
128 continue ;
129 }
130
131 // Check if this component has any dependents that need to be generated
132 // e.g. it can happen that there are none if all dependents of this component are prototyped
133 if (termDeps->empty()) {
134 cxcoutD(Generation) << "RooProdGenContext::ctor() term has no observables requested to be generated, removing it" << endl ;
135
136 // Increment the iterators first, because Removing the corresponding element
137 // would invalidate them otherwise.
138 ++termIter;
139 ++normIter;
140 ++impIter;
141 termList.Remove(term);
142 depsList.Remove(termDeps);
143 impDepList.Remove(impDeps);
144
145 delete term ;
146 delete termDeps ;
147 delete impDeps ;
148 anyAction=true ;
149 continue ;
150 }
151
152 if (term->size()==1) {
153 // Simple term
154
155 auto pdf = static_cast<RooAbsPdf*>((*term)[0]);
156 std::unique_ptr<RooArgSet> pdfDep{pdf->getObservables(termDeps)};
157 if (!pdfDep->empty()) {
158 coutI(Generation) << "RooProdGenContext::ctor() creating subcontext for generation of observables " << *pdfDep << " from model " << pdf->GetName() << endl ;
159 std::unique_ptr<RooArgSet> auxProto2{pdf->getObservables(impDeps)};
160 _gcList.emplace_back(pdf->genContext(*pdfDep,prototype,auxProto2.get(),verbose)) ;
161 }
162
163// cout << "adding following dependents to list of generated observables: " ; pdfDep->Print("1") ;
164 genDeps.add(*pdfDep) ;
165
166 } else {
167
168 // Composite term
169 if (!termDeps->empty()) {
170 const std::string name = model.makeRGPPName("PRODGEN_",*term,RooArgSet(),RooArgSet(),nullptr) ;
171
172 // Construct auxiliary PDF expressing product of composite terms,
173 // following Conditional component specification of input model
174 RooLinkedList cmdList ;
175 RooLinkedList pdfSetList ;
176 RooArgSet fullPdfSet ;
177 for(auto * pdf : static_range_cast<RooAbsPdf*>(*term)) {
178
179 RooArgSet* pdfnset = model.findPdfNSet(*pdf) ;
180 RooArgSet* pdfSet = new RooArgSet(*pdf) ;
181 pdfSetList.Add(pdfSet) ;
182
183 if (pdfnset && !pdfnset->empty()) {
184 // This PDF requires a Conditional() construction
185 cmdList.Add(RooFit::Conditional(*pdfSet,*pdfnset).Clone()) ;
186// cout << "Conditional " << pdf->GetName() << " " ; pdfnset->Print("1") ;
187 } else {
188 fullPdfSet.add(*pdfSet) ;
189 }
190
191 }
192 auto multiPdf = std::make_unique<RooProdPdf>(name.c_str(),name.c_str(),fullPdfSet,cmdList) ;
193 cmdList.Delete() ;
194 pdfSetList.Delete() ;
195
196 multiPdf->setOperMode(RooAbsArg::ADirty,true) ;
197 multiPdf->useDefaultGen(true) ;
198
199 coutI(Generation) << "RooProdGenContext()::ctor creating subcontext for generation of observables " << *termDeps
200 << "for irriducuble composite term using sub-product object " << multiPdf->GetName() ;
201 _gcList.emplace_back(multiPdf->genContext(*termDeps,prototype,auxProto,verbose));
202 _ownedMultiProds.addOwned(std::move(multiPdf));
203
204 genDeps.add(*termDeps) ;
205
206 }
207 }
208
209 // Increment the iterators first, because Removing the corresponding
210 // element would invalidate them otherwise.
211 ++termIter;
212 ++normIter;
213 ++impIter;
214 termList.Remove(term);
215 depsList.Remove(termDeps);
216 impDepList.Remove(impDeps);
217 delete term ;
218 delete termDeps ;
219 delete impDeps ;
220 anyAction=true ;
221 }
222 }
223
224 // Check if there are any left over terms that cannot be generated
225 // separately due to cross dependency of observables
226 if (!termList.empty()) {
227
228 cxcoutD(Generation) << "RooProdGenContext::ctor() there are left-over terms that need to be generated separately" << endl ;
229
230 // Concatenate remaining terms
231 auto normIter = depsList.begin();
232 RooArgSet trailerTerm ;
233 RooArgSet trailerTermDeps ;
234 for(auto * term : static_range_cast<RooArgSet*>(termList)) {
235 auto* termDeps = static_cast<RooArgSet*>(*normIter);
236 trailerTerm.add(*term) ;
237 trailerTermDeps.add(*termDeps) ;
238 ++normIter;
239 }
240
241 const std::string name = model.makeRGPPName("PRODGEN_",trailerTerm,RooArgSet(),RooArgSet(),nullptr) ;
242
243 // Construct auxiliary PDF expressing product of composite terms,
244 // following Partial/Full component specification of input model
245 RooLinkedList cmdList ;
246 RooLinkedList pdfSetList ;
247 RooArgSet fullPdfSet ;
248
249 for(auto * pdf : static_range_cast<RooAbsPdf*>(trailerTerm)) {
250
251 RooArgSet* pdfnset = model.findPdfNSet(*pdf) ;
252 RooArgSet* pdfSet = new RooArgSet(*pdf) ;
253 pdfSetList.Add(pdfSet) ;
254
255 if (pdfnset && !pdfnset->empty()) {
256 // This PDF requires a Conditional() construction
257 cmdList.Add(RooFit::Conditional(*pdfSet,*pdfnset).Clone()) ;
258 } else {
259 fullPdfSet.add(*pdfSet) ;
260 }
261
262 }
263// cmdList.Print("v") ;
264 auto multiPdf = std::make_unique<RooProdPdf>(name.c_str(),name.c_str(),fullPdfSet,cmdList);
265 cmdList.Delete() ;
266 pdfSetList.Delete() ;
267
268 multiPdf->setOperMode(RooAbsArg::ADirty,true) ;
269 multiPdf->useDefaultGen(true) ;
270
271 cxcoutD(Generation) << "RooProdGenContext(" << model.GetName() << "): creating context for irreducible composite trailer term "
272 << multiPdf->GetName() << " that generates observables " << trailerTermDeps << endl ;
273 _gcList.emplace_back(multiPdf->genContext(trailerTermDeps,prototype,auxProto,verbose));
274
275 _ownedMultiProds.addOwned(std::move(multiPdf));
276 }
277
278 // Now check if the are observables in vars that are not generated by any of the above p.d.f.s
279 // If not, generate uniform distributions for these using a special context
280 _uniObs.add(vars) ;
281 _uniObs.remove(genDeps,true,true) ;
282 if (!_uniObs.empty()) {
283 coutI(Generation) << "RooProdGenContext(" << model.GetName() << "): generating uniform distribution for non-dependent observable(s) " << _uniObs << std::endl;
284 }
285
286
287 // We own contents of lists filled by factorizeProduct()
288 termList.Delete() ;
289 depsList.Delete() ;
290 impDepList.Delete() ;
291 crossDepList.Delete() ;
292 intList.Delete() ;
293
294}
295
296
297
298////////////////////////////////////////////////////////////////////////////////
299/// Destructor. Delete all owned subgenerator contexts
300
302
303
304////////////////////////////////////////////////////////////////////////////////
305/// Attach generator to given event buffer
306
308{
309 //Forward initGenerator call to all components
310 for (auto const& elem : _gcList) {
311 elem->attach(args) ;
312 }
313}
314
315
316////////////////////////////////////////////////////////////////////////////////
317/// One-time initialization of generator context, forward to component generators
318
320{
321 // Forward initGenerator call to all components
322 for (auto const& elem : _gcList) {
323 elem->initGenerator(theEvent) ;
324 }
325}
326
327
328
329////////////////////////////////////////////////////////////////////////////////
330/// Generate a single event of the product by generating the components
331/// of the products sequentially. The subcontext have been order such
332/// that all conditional dependencies are correctly taken into account
333/// when processed in sequential order
334
336{
337 // Loop over the component generators
338
339 for (auto const& elem : _gcList) {
340 elem->generateEvent(theEvent,remaining) ;
341 }
342
343 // Generate uniform variables (non-dependents)
344 if (!_uniObs.empty()) {
345 for(auto * arglv : dynamic_range_cast<RooAbsLValue*>(_uniObs)) {
346 if (arglv) {
347 arglv->randomize() ;
348 }
349 }
350 theEvent.assign(_uniObs) ;
351 }
352}
353
354
355
356////////////////////////////////////////////////////////////////////////////////
357/// Set the traversal order of the prototype dataset by the
358/// given lookup table
359
361{
362 // Forward call to component generators
364
365 for (auto const& elem : _gcList) {
366 elem->setProtoDataOrder(lut) ;
367 }
368}
369
370
371
372////////////////////////////////////////////////////////////////////////////////
373/// Detailed printing interface
374
375void RooProdGenContext::printMultiline(ostream &os, Int_t content, bool verbose, TString indent) const
376{
377 RooAbsGenContext::printMultiline(os,content,verbose,indent) ;
378 os << indent << "--- RooProdGenContext ---" << endl ;
379 os << indent << "Using PDF ";
381 os << indent << "List of component generators" << endl ;
382
383 TString indent2(indent) ;
384 indent2.Append(" ") ;
385
386 for (auto const& elem : _gcList) {
387 elem->printMultiline(os,content,verbose,indent2) ;
388 }
389}
#define coutI(a)
#define cxcoutI(a)
#define cxcoutD(a)
#define dologD(a)
#define ccxcoutD(a)
#define ccxcoutI(a)
#define ClassImp(name)
Definition Rtypes.h:377
static void indent(ostringstream &buf, int indent_level)
char name[80]
Definition TGX11.cxx:110
RooFit::OwningPtr< RooArgSet > getObservables(const RooArgSet &set, bool valueOnly=true) const
Given a set of possible observables, return the observables that this PDF depends on.
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.
void assign(const RooAbsCollection &other) const
Sets the value, cache and constant attribute of any argument in our set that also appears in the othe...
virtual bool addOwned(RooAbsArg &var, bool silent=false)
Add an argument and transfer the ownership to the collection.
Abstract base class for generator contexts of RooAbsPdf objects.
void printMultiline(std::ostream &os, Int_t contents, bool verbose=false, TString indent="") const override
Interface for multi-line printing.
virtual void setProtoDataOrder(Int_t *lut)
Set the traversal order of prototype data to that in the lookup tables passed as argument.
Abstract interface for all probability density functions.
Definition RooAbsPdf.h:40
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:55
Container class to hold unbinned data.
Definition RooDataSet.h:57
const RooArgSet * get(Int_t index) const override
Return RooArgSet with coordinates of event 'index'.
Collection class for internal use, storing a collection of RooAbsArg pointers in a doubly linked list...
RooLinkedListIterImpl end() const
bool empty() const
void Delete(Option_t *o=nullptr) override
Remove all elements in collection and delete all elements NB: Collection does not own elements,...
virtual void Add(TObject *arg)
RooLinkedListIterImpl begin() const
virtual bool Remove(TObject *arg)
Remove object from collection.
virtual void printStream(std::ostream &os, Int_t contents, StyleOption style, TString indent="") const
Print description of object on ostream, printing contents set by contents integer,...
Efficient implementation of the generator context specific for RooProdPdf PDFs.
std::list< std::unique_ptr< RooAbsGenContext > > _gcList
List of component generator contexts.
void setProtoDataOrder(Int_t *lut) override
Set the traversal order of the prototype dataset by the given lookup table.
void attach(const RooArgSet &params) override
Attach generator to given event buffer.
RooArgSet _ownedMultiProds
Owned auxiliary multi-term product PDFs.
const RooProdPdf * _pdf
Original PDF.
void initGenerator(const RooArgSet &theEvent) override
One-time initialization of generator context, forward to component generators.
RooArgSet _uniObs
Observable to be generated with flat distribution.
~RooProdGenContext() override
Destructor. Delete all owned subgenerator contexts.
void generateEvent(RooArgSet &theEvent, Int_t remaining) override
Generate a single event of the product by generating the components of the products sequentially.
RooProdGenContext(const RooProdPdf &model, const RooArgSet &vars, const RooDataSet *prototype=nullptr, const RooArgSet *auxProto=nullptr, bool _verbose=false)
void printMultiline(std::ostream &os, Int_t content, bool verbose=false, TString indent="") const override
Detailed printing interface.
Efficient implementation of a product of PDFs of the form.
Definition RooProdPdf.h:33
void factorizeProduct(const RooArgSet &normSet, const RooArgSet &intSet, RooLinkedList &termList, RooLinkedList &normList, RooLinkedList &impDepList, RooLinkedList &crossDepList, RooLinkedList &intList) const
Factorize product in irreducible terms for given choice of integration/normalization.
std::string makeRGPPName(const char *pfx, const RooArgSet &term, const RooArgSet &iset, const RooArgSet &nset, const char *isetRangeName) const
Make an appropriate automatic name for a RooGenProdProj object in getPartIntList()
RooArgSet * findPdfNSet(RooAbsPdf const &pdf) const
Look up user specified normalization set for given input PDF component.
TObject * Clone(const char *newname="") const override
Make a clone of an object using the Streamer facility.
Definition TNamed.cxx:74
const char * GetName() const override
Returns name of object.
Definition TNamed.h:47
Basic string class.
Definition TString.h:139
TString & Append(const char *cs)
Definition TString.h:574
RooCmdArg Conditional(const RooArgSet &pdfSet, const RooArgSet &depSet, bool depsAreCond=false)