Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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
22Efficient 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
43#include <iostream>
44#include <string>
45
47
48////////////////////////////////////////////////////////////////////////////////
49/// Constructor of specialized generator context for RooSimultaneous p.d.f.s. This
50/// context creates a dedicated context for each component p.d.f.s and delegates
51/// generation of events to the appropriate component generator context
52
54 const RooDataSet *prototype, const RooArgSet* auxProto, bool verbose) :
55 RooAbsGenContext(model,vars,prototype,auxProto,verbose), _pdf(&model)
56{
57 // Determine if we are requested to generate the index category
58 RooAbsCategoryLValue const& idxCat = model.indexCat();
59 RooArgSet pdfVars(vars) ;
60
61 RooArgSet allPdfVars(pdfVars) ;
62 if (prototype) allPdfVars.add(*prototype->get(),true) ;
63
64 RooArgSet catsAmongAllVars;
65 allPdfVars.selectCommon(model.flattenedCatList(), catsAmongAllVars);
66
67 if(catsAmongAllVars.size() != model.flattenedCatList().size()) {
68 oocoutE(_pdf,Generation) << "RooSimGenContext::ctor(" << GetName() << ") ERROR: This context must"
69 << " generate all components of the index category" << std::endl ;
70 _isValid = false ;
71 _numPdf = 0 ;
72 _haveIdxProto = false ;
73 return ;
74 }
75
76 // We must either have the prototype or extended likelihood to determined
77 // the relative fractions of the components
78 _haveIdxProto = prototype ? true : false ;
79 _idxCatName = idxCat.GetName() ;
80 if (!_haveIdxProto && !model.canBeExtended()) {
81 oocoutE(_pdf,Generation) << "RooSimGenContext::ctor(" << GetName() << ") ERROR: Need either extended mode"
82 << " or prototype data to calculate number of events per category" << std::endl ;
83 _isValid = false ;
84 _numPdf = 0 ;
85 return ;
86 }
87
88 // Initialize fraction threshold array (used only in extended mode)
89 _numPdf = model._pdfProxyList.GetSize() ;
90 _fracThresh = new double[_numPdf+1] ;
91 _fracThresh[0] = 0 ;
92
93 // Generate index category and all registered PDFS
94 _allVarsPdf.add(allPdfVars) ;
95 Int_t i(1) ;
96 for(auto * proxy : static_range_cast<RooRealProxy*>(model._pdfProxyList)) {
97 auto* pdf = static_cast<RooAbsPdf*>(proxy->absArg());
98
99 // Create generator context for this PDF
100 RooAbsGenContext* cx = pdf->genContext(pdfVars,prototype,auxProto,verbose) ;
101
102 // Name the context after the associated state and add to list
103 cx->SetName(proxy->name()) ;
104 _gcList.push_back(cx) ;
105 _gcIndex.push_back(idxCat.lookupIndex(proxy->name()));
106
107 // Fill fraction threshold array
108 _fracThresh[i] = _fracThresh[i-1] + (_haveIdxProto?0:pdf->expectedEvents(&allPdfVars)) ;
109 i++ ;
110 }
111
112 // Normalize fraction threshold array
113 if (!_haveIdxProto) {
114 for(i=0 ; i<_numPdf ; i++)
116 }
117
118
119 // Clone the index category
120 _idxCatSet = new RooArgSet;
121 RooArgSet(model.indexCat()).snapshot(*_idxCatSet, true);
122 if (!_idxCatSet) {
123 oocoutE(_pdf,Generation) << "RooSimGenContext::RooSimGenContext(" << GetName() << ") Couldn't deep-clone index category, abort," << std::endl ;
124 throw std::string("RooSimGenContext::RooSimGenContext() Couldn't deep-clone index category, abort") ;
125 }
126
127 _idxCat = static_cast<RooAbsCategoryLValue*>(_idxCatSet->find(model.indexCat().GetName()));
128}
129
130
131
132////////////////////////////////////////////////////////////////////////////////
133/// Destructor. Delete all owned subgenerator contexts
134
136{
137 delete[] _fracThresh ;
138 delete _idxCatSet ;
139 for(RooAbsGenContext *item : _gcList) {
140 delete item;
141 }
142 if (_protoData) delete _protoData ;
143}
144
145
146
147////////////////////////////////////////////////////////////////////////////////
148/// Attach the index category clone to the given event buffer
149
151{
152 if (_idxCat->isDerived()) {
154 }
155
156 // Forward initGenerator call to all components
157 for(RooAbsGenContext *item : _gcList) {
158 item->attach(args) ;
159 }
160
161}
162
163
164////////////////////////////////////////////////////////////////////////////////
165/// Perform one-time initialization of generator context
166
168{
169 // Attach the index category clone to the event
170 if (_idxCat->isDerived()) {
172 } else {
173 _idxCat = static_cast<RooAbsCategoryLValue*>(theEvent.find(_idxCat->GetName())) ;
174 }
175
176 // Update fractions reflecting possible new parameter values
178
179 // Forward initGenerator call to all components
180 for(RooAbsGenContext *item : _gcList) {
181 item->initGenerator(theEvent) ;
182 }
183
184}
185
186
187////////////////////////////////////////////////////////////////////////////////
188/// Create an empty dataset to hold the events that will be generated
189
190RooDataSet* RooSimGenContext::createDataSet(const char* name, const char* title, const RooArgSet& obs)
191{
192
193 // If the observables do not contain the index, make a plain dataset
194 if (!obs.contains(*_idxCat)) {
195 return new RooDataSet(name,title,obs) ;
196 }
197
198 if (!_protoData) {
199 std::map<std::string,RooAbsData*> dmap ;
200 for (const auto& nameIdx : *_idxCat) {
201 RooAbsPdf* slicePdf = _pdf->getPdf(nameIdx.first.c_str());
202 std::unique_ptr<RooArgSet> sliceObs{slicePdf->getObservables(obs)};
203 std::string sliceName = Form("%s_slice_%s", name, nameIdx.first.c_str());
204 std::string sliceTitle = Form("%s (index slice %s)", title, nameIdx.first.c_str());
205 dmap[nameIdx.first] = new RooDataSet(sliceName,sliceTitle,*sliceObs);
206 }
207 using namespace RooFit;
208 _protoData = new RooDataSet(name, title, obs, Index(static_cast<RooCategory&>(*_idxCat)), Link(dmap), OwnLinked()) ;
209 }
210
211 RooDataSet* emptyClone = new RooDataSet(*_protoData,name) ;
212
213 return emptyClone ;
214}
215
216
217
218
219
220////////////////////////////////////////////////////////////////////////////////
221/// Generate event appropriate for current index state.
222/// The index state is taken either from the prototype
223/// or is generated from the fraction threshold table.
224
226{
227 if (_haveIdxProto) {
228
229 // Lookup pdf from selected prototype index state
230 Int_t gidx(0);
231 Int_t cidx = _idxCat->getCurrentIndex();
232 for (Int_t i=0 ; i<(Int_t)_gcIndex.size() ; i++) {
233 if (_gcIndex[i]==cidx) { gidx = i ; break ; }
234 }
235 RooAbsGenContext* cx = _gcList[gidx] ;
236 if (cx) {
237 cx->generateEvent(theEvent,remaining) ;
238 } else {
239 oocoutW(_pdf,Generation) << "RooSimGenContext::generateEvent: WARNING, no PDF to generate event of type " << cidx << std::endl ;
240 }
241
242
243 } else {
244
245 // Throw a random number and select PDF from fraction threshold table
246 double rand = RooRandom::uniform() ;
247 Int_t i=0 ;
248 for (i=0 ; i<_numPdf ; i++) {
249 if (rand>_fracThresh[i] && rand<_fracThresh[i+1]) {
250 RooAbsGenContext* gen=_gcList[i] ;
251 gen->generateEvent(theEvent,remaining) ;
252 //Write through to sub-categories because they might be written to dataset:
254 return ;
255 }
256 }
257
258 }
259}
260
261
262
263////////////////////////////////////////////////////////////////////////////////
264/// No action needed if we have a proto index
265
267{
268 if (_haveIdxProto) return ;
269
270 // Generate index category and all registered PDFS
271 Int_t i(1) ;
272 for(auto * proxy : static_range_cast<RooRealProxy*>(_pdf->_pdfProxyList)) {
273 auto* pdf = static_cast<RooAbsPdf*>(proxy->absArg());
274
275 // Fill fraction threshold array
276 _fracThresh[i] = _fracThresh[i-1] + (_haveIdxProto?0:pdf->expectedEvents(&_allVarsPdf)) ;
277 i++ ;
278 }
279
280 // Normalize fraction threshold array
281 if (!_haveIdxProto) {
282 for(i=0 ; i<_numPdf ; i++)
284 }
285
286}
287
288
289
290////////////////////////////////////////////////////////////////////////////////
291/// Set the traversal order of the prototype data to that in the
292/// given lookup table. This information is passed to all
293/// component generator contexts
294
296{
298
299 for (RooAbsGenContext *item : _gcList) {
300 item->setProtoDataOrder(lut) ;
301 }
302}
303
304
305////////////////////////////////////////////////////////////////////////////////
306/// Detailed printing interface
307
308void RooSimGenContext::printMultiline(std::ostream &os, Int_t content, bool verbose, TString indent) const
309{
310 RooAbsGenContext::printMultiline(os,content,verbose,indent) ;
311 os << indent << "--- RooSimGenContext ---" << std::endl ;
312 os << indent << "Using PDF ";
314 os << indent << "List of component generators" << std::endl ;
315
316 TString indent2(indent) ;
317 indent2.Append(" ") ;
318
319 for (RooAbsGenContext *item : _gcList) {
320 item->printMultiline(os,content,verbose,indent2);
321 }
322}
#define oocoutW(o, a)
#define oocoutE(o, a)
int Int_t
Definition RtypesCore.h:45
#define ClassImp(name)
Definition Rtypes.h:377
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:2489
bool recursiveRedirectServers(const RooAbsCollection &newServerList, bool mustReplaceAll=false, bool nameChange=false, bool recurseInNewSet=true)
Recursively replace all servers with the new servers in newSet.
RooFit::OwningPtr< RooArgSet > getObservables(const RooArgSet &set, bool valueOnly=true) const
Given a set of possible observables, return the observables that this PDF depends on.
virtual bool isDerived() const
Does value or shape of this arg depend on any other arg?
Definition RooAbsArg.h:97
Abstract base class for objects that represent a discrete value that can be set from the outside,...
virtual bool setIndex(value_type index, bool printError=true)=0
Change category state by specifying the index code of the desired state.
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 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.
Storage_t::size_type size() const
bool selectCommon(const RooAbsCollection &refColl, RooAbsCollection &outColl) const
Create a subset of the current collection, consisting only of those elements that are contained as we...
RooAbsArg * find(const char *name) const
Find object with given name in list.
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.
Abstract interface for all probability density functions.
Definition RooAbsPdf.h:40
bool canBeExtended() const
If true, PDF can provide extended likelihood term.
Definition RooAbsPdf.h:219
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:55
RooArgSet * snapshot(bool deepCopy=true) const
Use RooAbsCollection::snapshot(), but return as RooArgSet.
Definition RooArgSet.h:178
Object to represent discrete states.
Definition RooCategory.h:28
Container class to hold unbinned data.
Definition RooDataSet.h:57
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:78
Efficient implementation of the generator context specific for RooSimultaneous PDFs when generating m...
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.
Facilitates simultaneous fitting of multiple PDFs to subsets of a given dataset.
TList _pdfProxyList
List of PDF proxies (named after applicable category state)
RooArgSet const & flattenedCatList() const
Internal utility function to get a list of all category components for this RooSimultaneous.
RooAbsPdf * getPdf(RooStringView catName) const
Return the p.d.f associated with the given index category name.
const RooAbsCategoryLValue & indexCat() const
virtual Int_t GetSize() const
Return the capacity of the collection, i.e.
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:139
TString & Append(const char *cs)
Definition TString.h:572
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Definition JSONIO.h:26