Logo ROOT  
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
22RooSimGenContext is an efficient implementation of the generator context
23specific for RooSimultaneous PDFs when generating more than one of the
24component pdfs.
25It runs in two modes:
26- Proto data with category entries are given: An event from the same category as
27in the proto data is created.
28- No proto data: A category is chosen randomly.
29\note This requires that the PDFs are extended, to determine the relative probabilities
30that an event originates from a certain category.
31**/
32
33#include "RooFit.h"
34
35#include "RooSimGenContext.h"
36#include "RooSimultaneous.h"
37#include "RooRealProxy.h"
38#include "RooDataSet.h"
39#include "Roo1DTable.h"
40#include "RooCategory.h"
41#include "RooMsgService.h"
42#include "RooRandom.h"
43#include "RooGlobalFunc.h"
44
45using namespace RooFit;
46
47#include <iostream>
48#include <string>
49
50using namespace std;
51
53;
54
55
56////////////////////////////////////////////////////////////////////////////////
57/// Constructor of specialized generator context for RooSimultaneous p.d.f.s. This
58/// context creates a dedicated context for each component p.d.f.s and delegates
59/// generation of events to the appropriate component generator context
60
62 const RooDataSet *prototype, const RooArgSet* auxProto, Bool_t verbose) :
63 RooAbsGenContext(model,vars,prototype,auxProto,verbose), _pdf(&model), _protoData(0)
64{
65 // Determine if we are requested to generate the index category
66 RooAbsCategory *idxCat = (RooAbsCategory*) model._indexCat.absArg() ;
67 RooArgSet pdfVars(vars) ;
68
69 RooArgSet allPdfVars(pdfVars) ;
70 if (prototype) allPdfVars.add(*prototype->get(),kTRUE) ;
71
72 if (!idxCat->isDerived()) {
73 pdfVars.remove(*idxCat,kTRUE,kTRUE) ;
74 Bool_t doGenIdx = allPdfVars.find(idxCat->GetName())?kTRUE:kFALSE ;
75
76 if (!doGenIdx) {
77 oocoutE(_pdf,Generation) << "RooSimGenContext::ctor(" << GetName() << ") ERROR: This context must"
78 << " generate the index category" << endl ;
80 _numPdf = 0 ;
82 return ;
83 }
84 } else {
85 TIterator* sIter = idxCat->serverIterator() ;
86 RooAbsArg* server ;
87 Bool_t anyServer(kFALSE), allServers(kTRUE) ;
88 while((server=(RooAbsArg*)sIter->Next())) {
89 if (vars.find(server->GetName())) {
90 anyServer=kTRUE ;
91 pdfVars.remove(*server,kTRUE,kTRUE) ;
92 } else {
93 allServers=kFALSE ;
94 }
95 }
96 delete sIter ;
97
98 if (anyServer && !allServers) {
99 oocoutE(_pdf,Generation) << "RooSimGenContext::ctor(" << GetName() << ") ERROR: This context must"
100 << " generate all components of a derived index category" << endl ;
101 _isValid = kFALSE ;
102 _numPdf = 0 ;
104 return ;
105 }
106 }
107
108 // We must either have the prototype or extended likelihood to determined
109 // the relative fractions of the components
110 _haveIdxProto = prototype ? kTRUE : kFALSE ;
111 _idxCatName = idxCat->GetName() ;
112 if (!_haveIdxProto && !model.canBeExtended()) {
113 oocoutE(_pdf,Generation) << "RooSimGenContext::ctor(" << GetName() << ") ERROR: Need either extended mode"
114 << " or prototype data to calculate number of events per category" << endl ;
115 _isValid = kFALSE ;
116 _numPdf = 0 ;
117 return ;
118 }
119
120 // Initialize fraction threshold array (used only in extended mode)
121 _numPdf = model._pdfProxyList.GetSize() ;
122 _fracThresh = new Double_t[_numPdf+1] ;
123 _fracThresh[0] = 0 ;
124
125 // Generate index category and all registered PDFS
127 _allVarsPdf.add(allPdfVars) ;
128 RooRealProxy* proxy ;
129 RooAbsPdf* pdf ;
130 Int_t i(1) ;
131 while((proxy=(RooRealProxy*)_proxyIter->Next())) {
132 pdf=(RooAbsPdf*)proxy->absArg() ;
133
134 // Create generator context for this PDF
135 RooAbsGenContext* cx = pdf->genContext(pdfVars,prototype,auxProto,verbose) ;
136
137 // Name the context after the associated state and add to list
138 cx->SetName(proxy->name()) ;
139 _gcList.push_back(cx) ;
140 _gcIndex.push_back(idxCat->lookupIndex(proxy->name()));
141
142 // Fill fraction threshold array
143 _fracThresh[i] = _fracThresh[i-1] + (_haveIdxProto?0:pdf->expectedEvents(&allPdfVars)) ;
144 i++ ;
145 }
146
147 // Normalize fraction threshold array
148 if (!_haveIdxProto) {
149 for(i=0 ; i<_numPdf ; i++)
151 }
152
153
154 // Clone the index category
155 _idxCatSet = (RooArgSet*) RooArgSet(model._indexCat.arg()).snapshot(kTRUE) ;
156 if (!_idxCatSet) {
157 oocoutE(_pdf,Generation) << "RooSimGenContext::RooSimGenContext(" << GetName() << ") Couldn't deep-clone index category, abort," << endl ;
158 throw std::string("RooSimGenContext::RooSimGenContext() Couldn't deep-clone index category, abort") ;
159 }
160
162}
163
164
165
166////////////////////////////////////////////////////////////////////////////////
167/// Destructor. Delete all owned subgenerator contexts
168
170{
171 delete[] _fracThresh ;
172 delete _idxCatSet ;
173 for (vector<RooAbsGenContext*>::iterator iter = _gcList.begin() ; iter!=_gcList.end() ; ++iter) {
174 delete (*iter) ;
175 }
176 delete _proxyIter ;
177 if (_protoData) delete _protoData ;
178}
179
180
181
182////////////////////////////////////////////////////////////////////////////////
183/// Attach the index category clone to the given event buffer
184
186{
187 if (_idxCat->isDerived()) {
189 }
190
191 // Forward initGenerator call to all components
192 for (vector<RooAbsGenContext*>::iterator iter = _gcList.begin() ; iter!=_gcList.end() ; ++iter) {
193 (*iter)->attach(args) ;
194 }
195
196}
197
198
199////////////////////////////////////////////////////////////////////////////////
200/// Perform one-time initialization of generator context
201
203{
204 // Attach the index category clone to the event
205 if (_idxCat->isDerived()) {
207 } else {
209 }
210
211 // Update fractions reflecting possible new parameter values
213
214 // Forward initGenerator call to all components
215 for (vector<RooAbsGenContext*>::iterator iter = _gcList.begin() ; iter!=_gcList.end() ; ++iter) {
216 (*iter)->initGenerator(theEvent) ;
217 }
218
219}
220
221
222////////////////////////////////////////////////////////////////////////////////
223/// Create an empty dataset to hold the events that will be generated
224
225RooDataSet* RooSimGenContext::createDataSet(const char* name, const char* title, const RooArgSet& obs)
226{
227
228 // If the observables do not contain the index, make a plain dataset
229 if (!obs.contains(*_idxCat)) {
230 return new RooDataSet(name,title,obs) ;
231 }
232
233 if (!_protoData) {
234 map<string,RooAbsData*> dmap ;
235 for (const auto& nameIdx : *_idxCat) {
236 RooAbsPdf* slicePdf = _pdf->getPdf(nameIdx.first.c_str());
237 RooArgSet* sliceObs = slicePdf->getObservables(obs) ;
238 std::string sliceName = Form("%s_slice_%s", name, nameIdx.first.c_str());
239 std::string sliceTitle = Form("%s (index slice %s)", title, nameIdx.first.c_str());
240 RooDataSet* dset = new RooDataSet(sliceName.c_str(),sliceTitle.c_str(),*sliceObs) ;
241 dmap[nameIdx.first] = dset ;
242 delete sliceObs ;
243 }
244 _protoData = new RooDataSet(name, title, obs, Index((RooCategory&)*_idxCat), Link(dmap), OwnLinked()) ;
245 }
246
247 RooDataSet* emptyClone = new RooDataSet(*_protoData,name) ;
248
249 return emptyClone ;
250}
251
252
253
254
255
256////////////////////////////////////////////////////////////////////////////////
257/// Generate event appropriate for current index state.
258/// The index state is taken either from the prototype
259/// or is generated from the fraction threshold table.
260
262{
263 if (_haveIdxProto) {
264
265 // Lookup pdf from selected prototype index state
266 Int_t gidx(0), cidx =_idxCat->getCurrentIndex() ;
267 for (Int_t i=0 ; i<(Int_t)_gcIndex.size() ; i++) {
268 if (_gcIndex[i]==cidx) { gidx = i ; break ; }
269 }
270 RooAbsGenContext* cx = _gcList[gidx] ;
271 if (cx) {
272 cx->generateEvent(theEvent,remaining) ;
273 } else {
274 oocoutW(_pdf,Generation) << "RooSimGenContext::generateEvent: WARNING, no PDF to generate event of type " << cidx << endl ;
275 }
276
277
278 } else {
279
280 // Throw a random number and select PDF from fraction threshold table
282 Int_t i=0 ;
283 for (i=0 ; i<_numPdf ; i++) {
284 if (rand>_fracThresh[i] && rand<_fracThresh[i+1]) {
285 RooAbsGenContext* gen=_gcList[i] ;
286 gen->generateEvent(theEvent,remaining) ;
287 //Write through to sub-categories because they might be written to dataset:
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++)
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
347{
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}
#define oocoutW(o, a)
Definition: RooMsgService.h:47
#define oocoutE(o, a)
Definition: RooMsgService.h:48
int Int_t
Definition: RtypesCore.h:43
const Bool_t kFALSE
Definition: RtypesCore.h:90
bool Bool_t
Definition: RtypesCore.h:61
double Double_t
Definition: RtypesCore.h:57
const Bool_t kTRUE
Definition: RtypesCore.h:89
#define ClassImp(name)
Definition: Rtypes.h:361
static void indent(ostringstream &buf, int indent_level)
char name[80]
Definition: TGX11.cxx:109
char * Form(const char *fmt,...)
RooAbsArg is the common abstract base class for objects that represent a value (of arbitrary type) an...
Definition: RooAbsArg.h:73
RooArgSet * getObservables(const RooArgSet &set, Bool_t valueOnly=kTRUE) const
Return the observables of this pdf given a set of observables.
Definition: RooAbsArg.h:276
virtual Bool_t isDerived() const
Does value or shape of this arg depend on any other arg?
Definition: RooAbsArg.h:93
TIterator * serverIterator() const
Definition: RooAbsArg.h:141
Bool_t recursiveRedirectServers(const RooAbsCollection &newServerList, Bool_t mustReplaceAll=kFALSE, Bool_t nameChange=kFALSE, Bool_t recurseInNewSet=kTRUE)
Definition: RooAbsArg.cxx:1064
RooAbsCategoryLValue is the common abstract base class for objects that represent a discrete value th...
virtual bool setIndex(value_type index, bool printError=true)=0
Change category state by specifying the index code of the desired state.
RooAbsCategory is the base class for objects that represent a discrete value with a finite number of ...
virtual value_type getCurrentIndex() const
Return index number of current state.
value_type lookupIndex(const std::string &stateName) const
Find the index number corresponding to the state name.
Bool_t contains(const RooAbsArg &var) const
Check if collection contains an argument with the same name as var.
virtual Bool_t remove(const RooAbsArg &var, Bool_t silent=kFALSE, Bool_t matchByNameOnly=kFALSE)
Remove the specified argument from our list.
RooAbsArg * find(const char *name) const
Find object with given name in list.
RooAbsGenContext is the abstract base class for generator contexts of RooAbsPdf objects.
virtual void printMultiline(std::ostream &os, Int_t contents, Bool_t verbose=kFALSE, TString indent="") const
Interface for multi-line printing.
virtual void generateEvent(RooArgSet &theEvent, Int_t remaining)=0
virtual void setProtoDataOrder(Int_t *lut)
Set the traversal order of prototype data to that in the lookup tables passed as argument.
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:1998
Bool_t canBeExtended() const
Definition: RooAbsPdf.h:236
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:3283
RooAbsArg * absArg() const
Definition: RooArgProxy.h:37
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
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:88
RooCategory is an object to represent discrete states.
Definition: RooCategory.h:23
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:33
virtual const RooArgSet * get(Int_t index) const override
Return RooArgSet with coordinates of event 'index'.
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,...
static Double_t uniform(TRandom *generator=randomGenerator())
Return a number uniformly distributed from (0,1)
Definition: RooRandom.cxx:83
RooSimGenContext is an efficient implementation of the generator context specific for RooSimultaneous...
RooAbsCategoryLValue * _idxCat
RooArgSet _allVarsPdf
Prototype dataset.
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.
virtual void attach(const RooArgSet &params)
Attach the index category clone to the given event buffer.
RooArgSet * _idxCatSet
virtual void setProtoDataOrder(Int_t *lut)
Set the traversal order of the prototype data to that in the given lookup table.
void updateFractions()
No action needed if we have a proto index.
RooDataSet * _protoData
virtual void generateEvent(RooArgSet &theEvent, Int_t remaining)
Generate event appropriate for current index state.
RooDataSet * createDataSet(const char *name, const char *title, const RooArgSet &obs)
Create an empty dataset to hold the events that will be generated.
const RooSimultaneous * _pdf
std::vector< RooAbsGenContext * > _gcList
virtual ~RooSimGenContext()
Destructor. Delete all owned subgenerator contexts.
Double_t * _fracThresh
std::vector< int > _gcIndex
virtual void initGenerator(const RooArgSet &theEvent)
Perform one-time initialization of generator context.
virtual void printMultiline(std::ostream &os, Int_t content, Bool_t verbose=kFALSE, TString indent="") const
Detailed printing interface.
TIterator * _proxyIter
RooSimultaneous facilitates simultaneous fitting of multiple PDFs to subsets of a given dataset.
RooCategoryProxy _indexCat
RooAbsPdf * getPdf(const char *catName) const
Return the p.d.f associated with the given index category name.
const T & arg() const
Return reference to object held in proxy.
virtual Int_t GetSize() const
Return the capacity of the collection, i.e.
Definition: TCollection.h:182
Iterator abstract base class.
Definition: TIterator.h:30
virtual void Reset()=0
virtual TObject * Next()=0
virtual TIterator * MakeIterator(Bool_t dir=kIterForward) const
Return a list iterator.
Definition: TList.cxx:721
virtual void SetName(const char *name)
Set the name of the TNamed.
Definition: TNamed.cxx:140
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
Basic string class.
Definition: TString.h:131
TString & Append(const char *cs)
Definition: TString.h:559
RooCmdArg OwnLinked()
RooCmdArg Index(RooCategory &icat)
RooCmdArg Link(const char *state, RooAbsData &data)
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
@ Generation
Definition: RooGlobalFunc.h:67