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 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(),true) ;
69
70 if (!idxCat->isDerived()) {
71 pdfVars.remove(*idxCat,true,true) ;
72 bool doGenIdx = allPdfVars.find(idxCat->GetName())?true:false ;
73
74 if (!doGenIdx) {
75 oocoutE(_pdf,Generation) << "RooSimGenContext::ctor(" << GetName() << ") ERROR: This context must"
76 << " generate the index category" << endl ;
77 _isValid = false ;
78 _numPdf = 0 ;
79 _haveIdxProto = false ;
80 return ;
81 }
82 } else {
83 bool anyServer(false), allServers(true) ;
84 for(RooAbsArg* server : idxCat->servers()) {
85 if (vars.find(server->GetName())) {
86 anyServer=true ;
87 pdfVars.remove(*server,true,true) ;
88 } else {
89 allServers=false ;
90 }
91 }
92
93 if (anyServer && !allServers) {
94 oocoutE(_pdf,Generation) << "RooSimGenContext::ctor(" << GetName() << ") ERROR: This context must"
95 << " generate all components of a derived index category" << endl ;
96 _isValid = false ;
97 _numPdf = 0 ;
98 _haveIdxProto = false ;
99 return ;
100 }
101 }
102
103 // We must either have the prototype or extended likelihood to determined
104 // the relative fractions of the components
105 _haveIdxProto = prototype ? true : false ;
106 _idxCatName = idxCat->GetName() ;
107 if (!_haveIdxProto && !model.canBeExtended()) {
108 oocoutE(_pdf,Generation) << "RooSimGenContext::ctor(" << GetName() << ") ERROR: Need either extended mode"
109 << " or prototype data to calculate number of events per category" << endl ;
110 _isValid = false ;
111 _numPdf = 0 ;
112 return ;
113 }
114
115 // Initialize fraction threshold array (used only in extended mode)
116 _numPdf = model._pdfProxyList.GetSize() ;
117 _fracThresh = new double[_numPdf+1] ;
118 _fracThresh[0] = 0 ;
119
120 // Generate index category and all registered PDFS
121 _allVarsPdf.add(allPdfVars) ;
122 Int_t i(1) ;
123 for(auto * proxy : static_range_cast<RooRealProxy*>(model._pdfProxyList)) {
124 auto* pdf = static_cast<RooAbsPdf*>(proxy->absArg());
125
126 // Create generator context for this PDF
127 RooAbsGenContext* cx = pdf->genContext(pdfVars,prototype,auxProto,verbose) ;
128
129 // Name the context after the associated state and add to list
130 cx->SetName(proxy->name()) ;
131 _gcList.push_back(cx) ;
132 _gcIndex.push_back(idxCat->lookupIndex(proxy->name()));
133
134 // Fill fraction threshold array
135 _fracThresh[i] = _fracThresh[i-1] + (_haveIdxProto?0:pdf->expectedEvents(&allPdfVars)) ;
136 i++ ;
137 }
138
139 // Normalize fraction threshold array
140 if (!_haveIdxProto) {
141 for(i=0 ; i<_numPdf ; i++)
143 }
144
145
146 // Clone the index category
147 _idxCatSet = (RooArgSet*) RooArgSet(model._indexCat.arg()).snapshot(true) ;
148 if (!_idxCatSet) {
149 oocoutE(_pdf,Generation) << "RooSimGenContext::RooSimGenContext(" << GetName() << ") Couldn't deep-clone index category, abort," << endl ;
150 throw std::string("RooSimGenContext::RooSimGenContext() Couldn't deep-clone index category, abort") ;
151 }
152
154}
155
156
157
158////////////////////////////////////////////////////////////////////////////////
159/// Destructor. Delete all owned subgenerator contexts
160
162{
163 delete[] _fracThresh ;
164 delete _idxCatSet ;
165 for (vector<RooAbsGenContext*>::iterator iter = _gcList.begin() ; iter!=_gcList.end() ; ++iter) {
166 delete (*iter) ;
167 }
168 if (_protoData) delete _protoData ;
169}
170
171
172
173////////////////////////////////////////////////////////////////////////////////
174/// Attach the index category clone to the given event buffer
175
177{
178 if (_idxCat->isDerived()) {
180 }
181
182 // Forward initGenerator call to all components
183 for (vector<RooAbsGenContext*>::iterator iter = _gcList.begin() ; iter!=_gcList.end() ; ++iter) {
184 (*iter)->attach(args) ;
185 }
186
187}
188
189
190////////////////////////////////////////////////////////////////////////////////
191/// Perform one-time initialization of generator context
192
194{
195 // Attach the index category clone to the event
196 if (_idxCat->isDerived()) {
197 _idxCat->recursiveRedirectServers(theEvent,true) ;
198 } else {
200 }
201
202 // Update fractions reflecting possible new parameter values
204
205 // Forward initGenerator call to all components
206 for (vector<RooAbsGenContext*>::iterator iter = _gcList.begin() ; iter!=_gcList.end() ; ++iter) {
207 (*iter)->initGenerator(theEvent) ;
208 }
209
210}
211
212
213////////////////////////////////////////////////////////////////////////////////
214/// Create an empty dataset to hold the events that will be generated
215
216RooDataSet* RooSimGenContext::createDataSet(const char* name, const char* title, const RooArgSet& obs)
217{
218
219 // If the observables do not contain the index, make a plain dataset
220 if (!obs.contains(*_idxCat)) {
221 return new RooDataSet(name,title,obs) ;
222 }
223
224 if (!_protoData) {
225 map<string,RooAbsData*> dmap ;
226 for (const auto& nameIdx : *_idxCat) {
227 RooAbsPdf* slicePdf = _pdf->getPdf(nameIdx.first.c_str());
228 RooArgSet* sliceObs = slicePdf->getObservables(obs) ;
229 std::string sliceName = Form("%s_slice_%s", name, nameIdx.first.c_str());
230 std::string sliceTitle = Form("%s (index slice %s)", title, nameIdx.first.c_str());
231 RooDataSet* dset = new RooDataSet(sliceName.c_str(),sliceTitle.c_str(),*sliceObs) ;
232 dmap[nameIdx.first] = dset ;
233 delete sliceObs ;
234 }
235 _protoData = new RooDataSet(name, title, obs, Index((RooCategory&)*_idxCat), Link(dmap), OwnLinked()) ;
236 }
237
238 RooDataSet* emptyClone = new RooDataSet(*_protoData,name) ;
239
240 return emptyClone ;
241}
242
243
244
245
246
247////////////////////////////////////////////////////////////////////////////////
248/// Generate event appropriate for current index state.
249/// The index state is taken either from the prototype
250/// or is generated from the fraction threshold table.
251
253{
254 if (_haveIdxProto) {
255
256 // Lookup pdf from selected prototype index state
257 Int_t gidx(0), cidx =_idxCat->getCurrentIndex() ;
258 for (Int_t i=0 ; i<(Int_t)_gcIndex.size() ; i++) {
259 if (_gcIndex[i]==cidx) { gidx = i ; break ; }
260 }
261 RooAbsGenContext* cx = _gcList[gidx] ;
262 if (cx) {
263 cx->generateEvent(theEvent,remaining) ;
264 } else {
265 oocoutW(_pdf,Generation) << "RooSimGenContext::generateEvent: WARNING, no PDF to generate event of type " << cidx << endl ;
266 }
267
268
269 } else {
270
271 // Throw a random number and select PDF from fraction threshold table
272 double rand = RooRandom::uniform() ;
273 Int_t i=0 ;
274 for (i=0 ; i<_numPdf ; i++) {
275 if (rand>_fracThresh[i] && rand<_fracThresh[i+1]) {
276 RooAbsGenContext* gen=_gcList[i] ;
277 gen->generateEvent(theEvent,remaining) ;
278 //Write through to sub-categories because they might be written to dataset:
280 return ;
281 }
282 }
283
284 }
285}
286
287
288
289////////////////////////////////////////////////////////////////////////////////
290/// No action needed if we have a proto index
291
293{
294 if (_haveIdxProto) return ;
295
296 // Generate index category and all registered PDFS
297 Int_t i(1) ;
298 for(auto * proxy : static_range_cast<RooRealProxy*>(_pdf->_pdfProxyList)) {
299 auto* pdf = static_cast<RooAbsPdf*>(proxy->absArg());
300
301 // Fill fraction threshold array
302 _fracThresh[i] = _fracThresh[i-1] + (_haveIdxProto?0:pdf->expectedEvents(&_allVarsPdf)) ;
303 i++ ;
304 }
305
306 // Normalize fraction threshold array
307 if (!_haveIdxProto) {
308 for(i=0 ; i<_numPdf ; i++)
310 }
311
312}
313
314
315
316////////////////////////////////////////////////////////////////////////////////
317/// Set the traversal order of the prototype data to that in the
318/// given lookup table. This information is passed to all
319/// component generator contexts
320
322{
324
325 for (vector<RooAbsGenContext*>::iterator iter = _gcList.begin() ; iter!=_gcList.end() ; ++iter) {
326 (*iter)->setProtoDataOrder(lut) ;
327 }
328}
329
330
331////////////////////////////////////////////////////////////////////////////////
332/// Detailed printing interface
333
334void RooSimGenContext::printMultiline(ostream &os, Int_t content, bool verbose, TString indent) const
335{
337 os << indent << "--- RooSimGenContext ---" << endl ;
338 os << indent << "Using PDF ";
340 os << indent << "List of component generators" << endl ;
341
342 TString indent2(indent) ;
343 indent2.Append(" ") ;
344
345 for (vector<RooAbsGenContext*>::const_iterator iter = _gcList.begin() ; iter!=_gcList.end() ; ++iter) {
346 (*iter)->printMultiline(os,content,verbose,indent2);
347 }
348}
#define oocoutW(o, a)
Definition: RooMsgService.h:51
#define oocoutE(o, a)
Definition: RooMsgService.h:52
int Int_t
Definition: RtypesCore.h:45
#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:2452
RooAbsArg is the common abstract base class for objects that represent a value and a "shape" in RooFi...
Definition: RooAbsArg.h:71
bool recursiveRedirectServers(const RooAbsCollection &newServerList, bool mustReplaceAll=false, bool nameChange=false, bool recurseInNewSet=true)
Recursively replace all servers with the new servers in newSet.
Definition: RooAbsArg.cxx:1165
RooArgSet * getObservables(const RooArgSet &set, bool valueOnly=true) const
Given a set of possible observables, return the observables that this PDF depends on.
Definition: RooAbsArg.h:293
const RefCountList_t & servers() const
List of all servers of this object.
Definition: RooAbsArg.h:198
virtual bool isDerived() const
Does value or shape of this arg depend on any other arg?
Definition: RooAbsArg.h:91
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 remove(const RooAbsArg &var, bool silent=false, bool matchByNameOnly=false)
Remove the specified argument from our list.
bool contains(const RooAbsArg &var) const
Check if collection contains an argument with the same name as var.
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to 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 verbose=false, TString indent="") const override
Interface for multi-line printing.
virtual void generateEvent(RooArgSet &theEvent, Int_t remaining)=0
bool _isValid
Is context in valid state?
virtual void setProtoDataOrder(Int_t *lut)
Set the traversal order of prototype data to that in the lookup tables passed as argument.
bool canBeExtended() const
If true, PDF can provide extended likelihood term.
Definition: RooAbsPdf.h:264
RooAbsArg * absArg() const
Return pointer to contained argument.
Definition: RooArgProxy.h:47
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:56
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 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.
void printMultiline(std::ostream &os, Int_t content, bool verbose=false, TString indent="") const override
Detailed printing interface.
RooArgSet _allVarsPdf
All pdf variables.
double * _fracThresh
[_numPdf] Fraction threshold array
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
const RooSimultaneous * _pdf
Original PDF.
RooSimGenContext(const RooSimultaneous &model, const RooArgSet &vars, const RooDataSet *prototype=nullptr, const RooArgSet *auxProto=nullptr, bool _verbose=false)
Constructor of specialized generator context for RooSimultaneous p.d.f.s.
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.
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.
bool _haveIdxProto
Flag set if generation of index is requested.
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
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:61