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 "RooSimGenContext.h"
34#include "RooSimultaneous.h"
35#include "RooRealProxy.h"
36#include "RooDataSet.h"
37#include "Roo1DTable.h"
38#include "RooCategory.h"
39#include "RooMsgService.h"
40#include "RooRandom.h"
41#include "RooGlobalFunc.h"
42
43using namespace RooFit;
44
45#include <iostream>
46#include <string>
47
48using namespace std;
49
51;
52
53
54////////////////////////////////////////////////////////////////////////////////
55/// Constructor of specialized generator context for RooSimultaneous p.d.f.s. This
56/// context creates a dedicated context for each component p.d.f.s and delegates
57/// generation of events to the appropriate component generator context
58
60 const RooDataSet *prototype, const RooArgSet* auxProto, Bool_t verbose) :
61 RooAbsGenContext(model,vars,prototype,auxProto,verbose), _pdf(&model), _protoData(0)
62{
63 // Determine if we are requested to generate the index category
64 RooAbsCategory *idxCat = (RooAbsCategory*) model._indexCat.absArg() ;
65 RooArgSet pdfVars(vars) ;
66
67 RooArgSet allPdfVars(pdfVars) ;
68 if (prototype) allPdfVars.add(*prototype->get(),kTRUE) ;
69
70 if (!idxCat->isDerived()) {
71 pdfVars.remove(*idxCat,kTRUE,kTRUE) ;
72 Bool_t doGenIdx = allPdfVars.find(idxCat->GetName())?kTRUE:kFALSE ;
73
74 if (!doGenIdx) {
75 oocoutE(_pdf,Generation) << "RooSimGenContext::ctor(" << GetName() << ") ERROR: This context must"
76 << " generate the index category" << endl ;
78 _numPdf = 0 ;
80 return ;
81 }
82 } else {
83 TIterator* sIter = idxCat->serverIterator() ;
84 RooAbsArg* server ;
85 Bool_t anyServer(kFALSE), allServers(kTRUE) ;
86 while((server=(RooAbsArg*)sIter->Next())) {
87 if (vars.find(server->GetName())) {
88 anyServer=kTRUE ;
89 pdfVars.remove(*server,kTRUE,kTRUE) ;
90 } else {
91 allServers=kFALSE ;
92 }
93 }
94 delete sIter ;
95
96 if (anyServer && !allServers) {
97 oocoutE(_pdf,Generation) << "RooSimGenContext::ctor(" << GetName() << ") ERROR: This context must"
98 << " generate all components of a derived index category" << endl ;
100 _numPdf = 0 ;
102 return ;
103 }
104 }
105
106 // We must either have the prototype or extended likelihood to determined
107 // the relative fractions of the components
108 _haveIdxProto = prototype ? kTRUE : kFALSE ;
109 _idxCatName = idxCat->GetName() ;
110 if (!_haveIdxProto && !model.canBeExtended()) {
111 oocoutE(_pdf,Generation) << "RooSimGenContext::ctor(" << GetName() << ") ERROR: Need either extended mode"
112 << " or prototype data to calculate number of events per category" << endl ;
113 _isValid = kFALSE ;
114 _numPdf = 0 ;
115 return ;
116 }
117
118 // Initialize fraction threshold array (used only in extended mode)
119 _numPdf = model._pdfProxyList.GetSize() ;
120 _fracThresh = new Double_t[_numPdf+1] ;
121 _fracThresh[0] = 0 ;
122
123 // Generate index category and all registered PDFS
125 _allVarsPdf.add(allPdfVars) ;
126 RooRealProxy* proxy ;
127 RooAbsPdf* pdf ;
128 Int_t i(1) ;
129 while((proxy=(RooRealProxy*)_proxyIter->Next())) {
130 pdf=(RooAbsPdf*)proxy->absArg() ;
131
132 // Create generator context for this PDF
133 RooAbsGenContext* cx = pdf->genContext(pdfVars,prototype,auxProto,verbose) ;
134
135 // Name the context after the associated state and add to list
136 cx->SetName(proxy->name()) ;
137 _gcList.push_back(cx) ;
138 _gcIndex.push_back(idxCat->lookupIndex(proxy->name()));
139
140 // Fill fraction threshold array
141 _fracThresh[i] = _fracThresh[i-1] + (_haveIdxProto?0:pdf->expectedEvents(&allPdfVars)) ;
142 i++ ;
143 }
144
145 // Normalize fraction threshold array
146 if (!_haveIdxProto) {
147 for(i=0 ; i<_numPdf ; i++)
149 }
150
151
152 // Clone the index category
153 _idxCatSet = (RooArgSet*) RooArgSet(model._indexCat.arg()).snapshot(kTRUE) ;
154 if (!_idxCatSet) {
155 oocoutE(_pdf,Generation) << "RooSimGenContext::RooSimGenContext(" << GetName() << ") Couldn't deep-clone index category, abort," << endl ;
156 throw std::string("RooSimGenContext::RooSimGenContext() Couldn't deep-clone index category, abort") ;
157 }
158
160}
161
162
163
164////////////////////////////////////////////////////////////////////////////////
165/// Destructor. Delete all owned subgenerator contexts
166
168{
169 delete[] _fracThresh ;
170 delete _idxCatSet ;
171 for (vector<RooAbsGenContext*>::iterator iter = _gcList.begin() ; iter!=_gcList.end() ; ++iter) {
172 delete (*iter) ;
173 }
174 delete _proxyIter ;
175 if (_protoData) delete _protoData ;
176}
177
178
179
180////////////////////////////////////////////////////////////////////////////////
181/// Attach the index category clone to the given event buffer
182
184{
185 if (_idxCat->isDerived()) {
187 }
188
189 // Forward initGenerator call to all components
190 for (vector<RooAbsGenContext*>::iterator iter = _gcList.begin() ; iter!=_gcList.end() ; ++iter) {
191 (*iter)->attach(args) ;
192 }
193
194}
195
196
197////////////////////////////////////////////////////////////////////////////////
198/// Perform one-time initialization of generator context
199
201{
202 // Attach the index category clone to the event
203 if (_idxCat->isDerived()) {
205 } else {
207 }
208
209 // Update fractions reflecting possible new parameter values
211
212 // Forward initGenerator call to all components
213 for (vector<RooAbsGenContext*>::iterator iter = _gcList.begin() ; iter!=_gcList.end() ; ++iter) {
214 (*iter)->initGenerator(theEvent) ;
215 }
216
217}
218
219
220////////////////////////////////////////////////////////////////////////////////
221/// Create an empty dataset to hold the events that will be generated
222
223RooDataSet* RooSimGenContext::createDataSet(const char* name, const char* title, const RooArgSet& obs)
224{
225
226 // If the observables do not contain the index, make a plain dataset
227 if (!obs.contains(*_idxCat)) {
228 return new RooDataSet(name,title,obs) ;
229 }
230
231 if (!_protoData) {
232 map<string,RooAbsData*> dmap ;
233 for (const auto& nameIdx : *_idxCat) {
234 RooAbsPdf* slicePdf = _pdf->getPdf(nameIdx.first.c_str());
235 RooArgSet* sliceObs = slicePdf->getObservables(obs) ;
236 std::string sliceName = Form("%s_slice_%s", name, nameIdx.first.c_str());
237 std::string sliceTitle = Form("%s (index slice %s)", title, nameIdx.first.c_str());
238 RooDataSet* dset = new RooDataSet(sliceName.c_str(),sliceTitle.c_str(),*sliceObs) ;
239 dmap[nameIdx.first] = dset ;
240 delete sliceObs ;
241 }
242 _protoData = new RooDataSet(name, title, obs, Index((RooCategory&)*_idxCat), Link(dmap), OwnLinked()) ;
243 }
244
245 RooDataSet* emptyClone = new RooDataSet(*_protoData,name) ;
246
247 return emptyClone ;
248}
249
250
251
252
253
254////////////////////////////////////////////////////////////////////////////////
255/// Generate event appropriate for current index state.
256/// The index state is taken either from the prototype
257/// or is generated from the fraction threshold table.
258
260{
261 if (_haveIdxProto) {
262
263 // Lookup pdf from selected prototype index state
264 Int_t gidx(0), cidx =_idxCat->getCurrentIndex() ;
265 for (Int_t i=0 ; i<(Int_t)_gcIndex.size() ; i++) {
266 if (_gcIndex[i]==cidx) { gidx = i ; break ; }
267 }
268 RooAbsGenContext* cx = _gcList[gidx] ;
269 if (cx) {
270 cx->generateEvent(theEvent,remaining) ;
271 } else {
272 oocoutW(_pdf,Generation) << "RooSimGenContext::generateEvent: WARNING, no PDF to generate event of type " << cidx << endl ;
273 }
274
275
276 } else {
277
278 // Throw a random number and select PDF from fraction threshold table
280 Int_t i=0 ;
281 for (i=0 ; i<_numPdf ; i++) {
282 if (rand>_fracThresh[i] && rand<_fracThresh[i+1]) {
283 RooAbsGenContext* gen=_gcList[i] ;
284 gen->generateEvent(theEvent,remaining) ;
285 //Write through to sub-categories because they might be written to dataset:
287 return ;
288 }
289 }
290
291 }
292}
293
294
295
296////////////////////////////////////////////////////////////////////////////////
297/// No action needed if we have a proto index
298
300{
301 if (_haveIdxProto) return ;
302
303 // Generate index category and all registered PDFS
304 RooRealProxy* proxy ;
305 RooAbsPdf* pdf ;
306 Int_t i(1) ;
307 _proxyIter->Reset() ;
308 while((proxy=(RooRealProxy*)_proxyIter->Next())) {
309 pdf=(RooAbsPdf*)proxy->absArg() ;
310
311 // Fill fraction threshold array
313 i++ ;
314 }
315
316 // Normalize fraction threshold array
317 if (!_haveIdxProto) {
318 for(i=0 ; i<_numPdf ; i++)
320 }
321
322}
323
324
325
326////////////////////////////////////////////////////////////////////////////////
327/// Set the traversal order of the prototype data to that in the
328/// given lookup table. This information is passed to all
329/// component generator contexts
330
332{
334
335 for (vector<RooAbsGenContext*>::iterator iter = _gcList.begin() ; iter!=_gcList.end() ; ++iter) {
336 (*iter)->setProtoDataOrder(lut) ;
337 }
338}
339
340
341////////////////////////////////////////////////////////////////////////////////
342/// Detailed printing interface
343
345{
347 os << indent << "--- RooSimGenContext ---" << endl ;
348 os << indent << "Using PDF ";
350 os << indent << "List of component generators" << endl ;
351
352 TString indent2(indent) ;
353 indent2.Append(" ") ;
354
355 for (vector<RooAbsGenContext*>::const_iterator iter = _gcList.begin() ; iter!=_gcList.end() ; ++iter) {
356 (*iter)->printMultiline(os,content,verbose,indent2);
357 }
358}
#define oocoutW(o, a)
Definition: RooMsgService.h:47
#define oocoutE(o, a)
Definition: RooMsgService.h:48
int Int_t
Definition: RtypesCore.h:45
const Bool_t kFALSE
Definition: RtypesCore.h:101
const Bool_t kTRUE
Definition: RtypesCore.h:100
#define ClassImp(name)
Definition: Rtypes.h:375
static void indent(ostringstream &buf, int indent_level)
char name[80]
Definition: TGX11.cxx:110
char * Form(const char *fmt,...)
Formats a string in a circular formatting buffer.
Definition: TString.cxx:2447
RooAbsArg is the common abstract base class for objects that represent a value and a "shape" in RooFi...
Definition: RooAbsArg.h:78
RooArgSet * getObservables(const RooArgSet &set, Bool_t valueOnly=kTRUE) const
Given a set of possible observables, return the observables that this PDF depends on.
Definition: RooAbsArg.h:318
virtual Bool_t isDerived() const
Does value or shape of this arg depend on any other arg?
Definition: RooAbsArg.h:98
TIterator * serverIterator() const
Definition: RooAbsArg.h:146
Bool_t recursiveRedirectServers(const RooAbsCollection &newServerList, Bool_t mustReplaceAll=kFALSE, Bool_t nameChange=kFALSE, Bool_t recurseInNewSet=kTRUE)
Recursively replace all servers with the new servers in newSet.
Definition: RooAbsArg.cxx:1197
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.
A space to attach TBranches.
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.
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Add the specified argument to list.
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.
void printMultiline(std::ostream &os, Int_t contents, Bool_t verbose=kFALSE, TString indent="") const override
Interface for multi-line printing.
Bool_t _isValid
Is context in valid state?
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:1903
Bool_t canBeExtended() const
If true, PDF can provide extended likelihood term.
Definition: RooAbsPdf.h:263
virtual Double_t expectedEvents(const RooArgSet *nset) const
Return expected number of events to be used in calculation of extended likelihood.
Definition: RooAbsPdf.cxx:3231
RooAbsArg * absArg() const
Return pointer to contained argument.
Definition: RooArgProxy.h:40
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:57
RooCategory is an object to represent discrete states.
Definition: RooCategory.h:28
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:55
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:81
RooSimGenContext is an efficient implementation of the generator context specific for RooSimultaneous...
RooAbsCategoryLValue * _idxCat
Clone of index category.
void initGenerator(const RooArgSet &theEvent) override
Perform one-time initialization of generator context.
RooArgSet _allVarsPdf
All pdf variables.
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.
void printMultiline(std::ostream &os, Int_t content, Bool_t verbose=kFALSE, TString indent="") const override
Detailed printing interface.
RooArgSet * _idxCatSet
Owner of index category components.
void setProtoDataOrder(Int_t *lut) override
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
! Prototype dataset
Bool_t _haveIdxProto
Flag set if generation of index is requested.
const RooSimultaneous * _pdf
Original PDF.
RooDataSet * createDataSet(const char *name, const char *title, const RooArgSet &obs) override
Create an empty dataset to hold the events that will be generated.
TString _idxCatName
Name of index category.
std::vector< RooAbsGenContext * > _gcList
List of component generator contexts.
void generateEvent(RooArgSet &theEvent, Int_t remaining) override
Generate event appropriate for current index state.
Int_t _numPdf
Number of generated PDFs.
Double_t * _fracThresh
[_numPdf] Fraction threshold array
std::vector< int > _gcIndex
Index value corresponding to component.
void attach(const RooArgSet &params) override
Attach the index category clone to the given event buffer.
~RooSimGenContext() override
Destructor. Delete all owned subgenerator contexts.
TIterator * _proxyIter
Iterator over pdf proxies.
RooSimultaneous facilitates simultaneous fitting of multiple PDFs to subsets of a given dataset.
TList _pdfProxyList
List of PDF proxies (named after applicable category state)
RooCategoryProxy _indexCat
Index category.
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:184
Iterator abstract base class.
Definition: TIterator.h:30
virtual void Reset()=0
virtual TObject * Next()=0
TIterator * MakeIterator(Bool_t dir=kIterForward) const override
Return a list iterator.
Definition: TList.cxx:722
const char * GetName() const override
Returns name of object.
Definition: TNamed.h:47
virtual void SetName(const char *name)
Set the name of the TNamed.
Definition: TNamed.cxx:140
Basic string class.
Definition: TString.h:136
TString & Append(const char *cs)
Definition: TString.h:564
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...
Definition: Common.h:18
@ Generation
Definition: RooGlobalFunc.h:63