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