Logo ROOT   6.10/09
Reference Guide
RooSimGenContext.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 RooSimGenContext.cxx
19 \class RooSimGenContext
20 \ingroup Roofitcore
21 
22 RooSimGenContext is an efficient implementation of the generator context
23 specific for RooSimultaneous PDFs when generating more than one of the
24 component pdfs.
25 **/
26 
27 #include "RooFit.h"
28 #include "Riostream.h"
29 
30 #include "RooSimGenContext.h"
31 #include "RooSimultaneous.h"
32 #include "RooRealProxy.h"
33 #include "RooDataSet.h"
34 #include "Roo1DTable.h"
35 #include "RooCategory.h"
36 #include "RooMsgService.h"
37 #include "RooRandom.h"
38 #include "RooGlobalFunc.h"
39 
40 using namespace RooFit ;
41 
42 #include <string>
43 
44 using namespace std;
45 
47 ;
48 
49 
50 ////////////////////////////////////////////////////////////////////////////////
51 /// Constructor of specialized generator context for RooSimultaneous p.d.f.s. This
52 /// context creates a dedicated context for each component p.d.f.s and delegates
53 /// generation of events to the appropriate component generator context
54 
56  const RooDataSet *prototype, const RooArgSet* auxProto, Bool_t verbose) :
57  RooAbsGenContext(model,vars,prototype,auxProto,verbose), _pdf(&model), _protoData(0)
58 {
59  // Determine if we are requested to generate the index category
60  RooAbsCategory *idxCat = (RooAbsCategory*) model._indexCat.absArg() ;
61  RooArgSet pdfVars(vars) ;
62 
63  RooArgSet allPdfVars(pdfVars) ;
64  if (prototype) allPdfVars.add(*prototype->get(),kTRUE) ;
65 
66  if (!idxCat->isDerived()) {
67  pdfVars.remove(*idxCat,kTRUE,kTRUE) ;
68  Bool_t doGenIdx = allPdfVars.find(idxCat->GetName())?kTRUE:kFALSE ;
69 
70  if (!doGenIdx) {
71  oocoutE(_pdf,Generation) << "RooSimGenContext::ctor(" << GetName() << ") ERROR: This context must"
72  << " generate the index category" << endl ;
73  _isValid = kFALSE ;
74  _numPdf = 0 ;
76  return ;
77  }
78  } else {
79  TIterator* sIter = idxCat->serverIterator() ;
80  RooAbsArg* server ;
81  Bool_t anyServer(kFALSE), allServers(kTRUE) ;
82  while((server=(RooAbsArg*)sIter->Next())) {
83  if (vars.find(server->GetName())) {
84  anyServer=kTRUE ;
85  pdfVars.remove(*server,kTRUE,kTRUE) ;
86  } else {
87  allServers=kFALSE ;
88  }
89  }
90  delete sIter ;
91 
92  if (anyServer && !allServers) {
93  oocoutE(_pdf,Generation) << "RooSimGenContext::ctor(" << GetName() << ") ERROR: This context must"
94  << " generate all components of a derived index category" << endl ;
95  _isValid = kFALSE ;
96  _numPdf = 0 ;
98  return ;
99  }
100  }
101 
102  // We must either have the prototype or extended likelihood to determined
103  // the relative fractions of the components
104  _haveIdxProto = prototype ? kTRUE : kFALSE ;
105  _idxCatName = idxCat->GetName() ;
106  if (!_haveIdxProto && !model.canBeExtended()) {
107  oocoutE(_pdf,Generation) << "RooSimGenContext::ctor(" << GetName() << ") ERROR: Need either extended mode"
108  << " or prototype data to calculate number of events per category" << endl ;
109  _isValid = kFALSE ;
110  _numPdf = 0 ;
111  return ;
112  }
113 
114  // Initialize fraction threshold array (used only in extended mode)
115  _numPdf = model._pdfProxyList.GetSize() ;
116  _fracThresh = new Double_t[_numPdf+1] ;
117  _fracThresh[0] = 0 ;
118 
119  // Generate index category and all registered PDFS
121  _allVarsPdf.add(allPdfVars) ;
122  RooRealProxy* proxy ;
123  RooAbsPdf* pdf ;
124  Int_t i(1) ;
125  while((proxy=(RooRealProxy*)_proxyIter->Next())) {
126  pdf=(RooAbsPdf*)proxy->absArg() ;
127 
128  // Create generator context for this PDF
129  RooAbsGenContext* cx = pdf->genContext(pdfVars,prototype,auxProto,verbose) ;
130 
131  // Name the context after the associated state and add to list
132  cx->SetName(proxy->name()) ;
133  _gcList.push_back(cx) ;
134  _gcIndex.push_back(idxCat->lookupType(proxy->name())->getVal()) ;
135 
136  // Fill fraction threshold array
137  _fracThresh[i] = _fracThresh[i-1] + (_haveIdxProto?0:pdf->expectedEvents(&allPdfVars)) ;
138  i++ ;
139  }
140 
141  // Normalize fraction threshold array
142  if (!_haveIdxProto) {
143  for(i=0 ; i<_numPdf ; i++)
144  _fracThresh[i] /= _fracThresh[_numPdf] ;
145  }
146 
147 
148  // Clone the index category
149  _idxCatSet = (RooArgSet*) RooArgSet(model._indexCat.arg()).snapshot(kTRUE) ;
150  if (!_idxCatSet) {
151  oocoutE(_pdf,Generation) << "RooSimGenContext::RooSimGenContext(" << GetName() << ") Couldn't deep-clone index category, abort," << endl ;
152  throw std::string("RooSimGenContext::RooSimGenContext() Couldn't deep-clone index category, abort") ;
153  }
154 
156 }
157 
158 
159 
160 ////////////////////////////////////////////////////////////////////////////////
161 /// Destructor. Delete all owned subgenerator contexts
162 
164 {
165  delete[] _fracThresh ;
166  delete _idxCatSet ;
167  for (vector<RooAbsGenContext*>::iterator iter = _gcList.begin() ; iter!=_gcList.end() ; ++iter) {
168  delete (*iter) ;
169  }
170  delete _proxyIter ;
171  if (_protoData) delete _protoData ;
172 }
173 
174 
175 
176 ////////////////////////////////////////////////////////////////////////////////
177 /// Attach the index category clone to the given event buffer
178 
180 {
181  if (_idxCat->isDerived()) {
183  }
184 
185  // Forward initGenerator call to all components
186  for (vector<RooAbsGenContext*>::iterator iter = _gcList.begin() ; iter!=_gcList.end() ; ++iter) {
187  (*iter)->attach(args) ;
188  }
189 
190 }
191 
192 
193 ////////////////////////////////////////////////////////////////////////////////
194 /// Perform one-time initialization of generator context
195 
197 {
198  // Attach the index category clone to the event
199  if (_idxCat->isDerived()) {
201  } else {
202  _idxCat = (RooAbsCategoryLValue*) theEvent.find(_idxCat->GetName()) ;
203  }
204 
205  // Update fractions reflecting possible new parameter values
206  updateFractions() ;
207 
208  // Forward initGenerator call to all components
209  for (vector<RooAbsGenContext*>::iterator iter = _gcList.begin() ; iter!=_gcList.end() ; ++iter) {
210  (*iter)->initGenerator(theEvent) ;
211  }
212 
213 }
214 
215 
216 ////////////////////////////////////////////////////////////////////////////////
217 /// Create an empty dataset to hold the events that will be generated
218 
219 RooDataSet* RooSimGenContext::createDataSet(const char* name, const char* title, const RooArgSet& obs)
220 {
221 
222  // If the observables do not contain the index, make a plain dataset
223  if (!obs.contains(*_idxCat)) {
224  return new RooDataSet(name,title,obs) ;
225  }
226 
227  if (!_protoData) {
228  map<string,RooAbsData*> dmap ;
229  RooCatType* state ;
230  TIterator* iter = _idxCat->typeIterator() ;
231  while((state=(RooCatType*)iter->Next())) {
232  RooAbsPdf* slicePdf = _pdf->getPdf(state->GetName()) ;
233  RooArgSet* sliceObs = slicePdf->getObservables(obs) ;
234  std::string sliceName = Form("%s_slice_%s",name,state->GetName()) ;
235  std::string sliceTitle = Form("%s (index slice %s)",title,state->GetName()) ;
236  RooDataSet* dset = new RooDataSet(sliceName.c_str(),sliceTitle.c_str(),*sliceObs) ;
237  dmap[state->GetName()] = dset ;
238  delete sliceObs ;
239  }
240  delete iter ;
241  _protoData = new RooDataSet(name, title, obs, Index((RooCategory&)*_idxCat), Link(dmap), OwnLinked()) ;
242 
243 // RooDataSet* tmp = _protoData ;
244 // _protoData = 0 ;
245 // return tmp ;
246  }
247 
248  RooDataSet* emptyClone = new RooDataSet(*_protoData,name) ;
249 
250  return emptyClone ;
251 }
252 
253 
254 
255 
256 
257 ////////////////////////////////////////////////////////////////////////////////
258 /// Generate event appropriate for current index state.
259 /// The index state is taken either from the prototype
260 /// or is generated from the fraction threshold table.
261 
263 {
264  if (_haveIdxProto) {
265 
266  // Lookup pdf from selected prototype index state
267  Int_t gidx(0), cidx =_idxCat->getIndex() ;
268  for (Int_t i=0 ; i<(Int_t)_gcIndex.size() ; i++) {
269  if (_gcIndex[i]==cidx) { gidx = i ; break ; }
270  }
271  RooAbsGenContext* cx = _gcList[gidx] ;
272  if (cx) {
273  cx->generateEvent(theEvent,remaining) ;
274  } else {
275  oocoutW(_pdf,Generation) << "RooSimGenContext::generateEvent: WARNING, no PDF to generate event of type " << cidx << endl ;
276  }
277 
278 
279  } else {
280 
281  // Throw a random number and select PDF from fraction threshold table
282  Double_t rand = RooRandom::uniform() ;
283  Int_t i=0 ;
284  for (i=0 ; i<_numPdf ; i++) {
285  if (rand>_fracThresh[i] && rand<_fracThresh[i+1]) {
286  RooAbsGenContext* gen=_gcList[i] ;
287  gen->generateEvent(theEvent,remaining) ;
289  return ;
290  }
291  }
292 
293  }
294 }
295 
296 
297 
298 ////////////////////////////////////////////////////////////////////////////////
299 /// No action needed if we have a proto index
300 
302 {
303  if (_haveIdxProto) return ;
304 
305  // Generate index category and all registered PDFS
306  RooRealProxy* proxy ;
307  RooAbsPdf* pdf ;
308  Int_t i(1) ;
309  _proxyIter->Reset() ;
310  while((proxy=(RooRealProxy*)_proxyIter->Next())) {
311  pdf=(RooAbsPdf*)proxy->absArg() ;
312 
313  // Fill fraction threshold array
315  i++ ;
316  }
317 
318  // Normalize fraction threshold array
319  if (!_haveIdxProto) {
320  for(i=0 ; i<_numPdf ; i++)
321  _fracThresh[i] /= _fracThresh[_numPdf] ;
322  }
323 
324 }
325 
326 
327 
328 ////////////////////////////////////////////////////////////////////////////////
329 /// Set the traversal order of the prototype data to that in the
330 /// given lookup table. This information is passed to all
331 /// component generator contexts
332 
334 {
336 
337  for (vector<RooAbsGenContext*>::iterator iter = _gcList.begin() ; iter!=_gcList.end() ; ++iter) {
338  (*iter)->setProtoDataOrder(lut) ;
339  }
340 }
341 
342 
343 ////////////////////////////////////////////////////////////////////////////////
344 /// Detailed printing interface
345 
346 void RooSimGenContext::printMultiline(ostream &os, Int_t content, Bool_t verbose, TString indent) const
347 {
348  RooAbsGenContext::printMultiline(os,content,verbose,indent) ;
349  os << indent << "--- RooSimGenContext ---" << endl ;
350  os << indent << "Using PDF ";
352  os << indent << "List of component generators" << endl ;
353 
354  TString indent2(indent) ;
355  indent2.Append(" ") ;
356 
357  for (vector<RooAbsGenContext*>::const_iterator iter = _gcList.begin() ; iter!=_gcList.end() ; ++iter) {
358  (*iter)->printMultiline(os,content,verbose,indent2);
359  }
360 }
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
TIterator * _proxyIter
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;.
RooCategoryProxy _indexCat
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
virtual void Reset()=0
virtual ~RooSimGenContext()
Destructor. Delete all owned subgenerator contexts.
RooDataSet * createDataSet(const char *name, const char *title, const RooArgSet &obs)
Create an empty dataset to hold the events that will be generated.
RooArgSet * getObservables(const RooArgSet &set, Bool_t valueOnly=kTRUE) const
Definition: RooAbsArg.h:194
virtual void SetName(const char *name)
Set the name of the TNamed.
Definition: TNamed.cxx:131
RooCmdArg Link(const char *state, RooAbsData &data)
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
virtual void setIndexFast(Int_t index)
virtual Int_t getIndex() const
Return index number of current state.
STL namespace.
RooArgSet _allVarsPdf
Prototype dataset.
Iterator abstract base class.
Definition: TIterator.h:30
void updateFractions()
No action needed if we have a proto index.
virtual void generateEvent(RooArgSet &theEvent, Int_t remaining)
Generate event appropriate for current index state.
RooArgSet * _idxCatSet
#define oocoutE(o, a)
Definition: RooMsgService.h:47
const RooAbsCategory & arg() const
RooSimGenContext(const RooSimultaneous &model, const RooArgSet &vars, const RooDataSet *prototype=0, const RooArgSet *auxProto=0, Bool_t _verbose=kFALSE)
Constructor of specialized generator context for RooSimultaneous p.d.f.s.
RooCatType is an auxilary class for RooAbsCategory and defines a a single category state...
Definition: RooCatType.h:22
virtual void attach(const RooArgSet &params)
Attach the index category clone to the given event buffer.
RooAbsGenContext is the abstract base class for generator contexts of RooAbsPdf objects.
RooAbsCategoryLValue is the common abstract base class for objects that represent a discrete value th...
virtual Double_t expectedEvents(const RooArgSet *nset) const
Return expected number of events from this p.d.f for use in extended likelihood calculations.
Definition: RooAbsPdf.cxx:2930
RooCmdArg OwnLinked()
virtual TIterator * MakeIterator(Bool_t dir=kIterForward) const
Return a list iterator.
Definition: TList.cxx:608
virtual const char * name() const
Definition: RooArgProxy.h:42
std::vector< int > _gcIndex
RooSimGenContext is an efficient implementation of the generator context specific for RooSimultaneous...
TIterator * typeIterator() const
Return iterator over all defined states.
virtual const Text_t * GetName() const
Returns name of object.
Definition: RooCatType.h:44
const RooCatType * lookupType(Int_t index, Bool_t printError=kFALSE) const
Find our type corresponding to the specified index, or return 0 for no match.
bool verbose
char * Form(const char *fmt,...)
RooAbsArg * absArg() const
Definition: RooArgProxy.h:37
virtual void generateEvent(RooArgSet &theEvent, Int_t remaining)=0
virtual void initGenerator(const RooArgSet &theEvent)
Perform one-time initialization of generator context.
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:29
RooCategory represents a fundamental (non-derived) discrete value object.
Definition: RooCategory.h:24
Bool_t canBeExtended() const
Definition: RooAbsPdf.h:216
TIterator * serverIterator() const
Definition: RooAbsArg.h:112
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
Double_t * _fracThresh
RooCmdArg Index(RooCategory &icat)
#define ClassImp(name)
Definition: Rtypes.h:336
static Double_t uniform(TRandom *generator=randomGenerator())
Return a number uniformly distributed from (0,1)
Definition: RooRandom.cxx:84
RooAbsArg * find(const char *name) const
Find object with given name in list.
double Double_t
Definition: RtypesCore.h:55
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.
const RooSimultaneous * _pdf
virtual void printMultiline(std::ostream &os, Int_t contents, Bool_t verbose=kFALSE, TString indent="") const
Interface for multi-line printing.
virtual void printMultiline(std::ostream &os, Int_t content, Bool_t verbose=kFALSE, TString indent="") const
Detailed printing interface.
#define oocoutW(o, a)
Definition: RooMsgService.h:46
std::vector< RooAbsGenContext * > _gcList
virtual Bool_t remove(const RooAbsArg &var, Bool_t silent=kFALSE, Bool_t matchByNameOnly=kFALSE)
Remove the specified argument from our list.
RooAbsPdf * getPdf(const char *catName) const
Return the p.d.f associated with the given index category name.
RooAbsPdf is the abstract interface for all probability density functions The class provides hybrid a...
Definition: RooAbsPdf.h:41
virtual Bool_t isDerived() const
Definition: RooAbsArg.h:81
virtual TObject * Next()=0
RooAbsCategoryLValue * _idxCat
RooAbsCategory is the common abstract base class for objects that represent a discrete value with a f...
RooRealProxy is the concrete proxy for RooAbsReal objects A RooRealProxy is the general mechanism to ...
Definition: RooRealProxy.h:23
Bool_t contains(const RooAbsArg &var) const
RooAbsArg is the common abstract base class for objects that represent a value (of arbitrary type) an...
Definition: RooAbsArg.h:66
virtual Int_t GetSize() const
Definition: TCollection.h:89
Bool_t recursiveRedirectServers(const RooAbsCollection &newServerList, Bool_t mustReplaceAll=kFALSE, Bool_t nameChange=kFALSE, Bool_t recurseInNewSet=kTRUE)
Definition: RooAbsArg.cxx:1088
RooDataSet * _protoData
virtual void setProtoDataOrder(Int_t *lut)
Set the traversal order of the prototype data to that in the given lookup table.
RooSimultaneous facilitates simultaneous fitting of multiple PDFs to subsets of a given dataset...
const Bool_t kTRUE
Definition: RtypesCore.h:91