Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
RooGenContext.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 RooGenContext.cxx
19\class RooGenContext
20\ingroup Roofitcore
21
22Implements a universal generator context for all
23RooAbsPdf classes that do not have or need a specialized generator
24context. This generator context queries the input p.d.f which observables
25it can generate internally and delegates generation of those observables
26to the p.d.f if it deems that safe. The other observables are generated
27use a RooAcceptReject sampling technique.
28**/
29
30#include "RooMsgService.h"
31#include "Riostream.h"
32
33#include "RooGenContext.h"
34#include "RooAbsPdf.h"
35#include "RooDataSet.h"
36#include "RooRealIntegral.h"
37#include "RooAcceptReject.h"
38#include "RooRealVar.h"
39#include "RooDataHist.h"
40#include "RooErrorHandler.h"
41#include "RooNumGenConfig.h"
42#include "RooNumGenFactory.h"
43
44#include "TString.h"
45
46
47using std::endl, std::string, std::ostream;
48
50
51
52////////////////////////////////////////////////////////////////////////////////
53/// Initialize a new context for generating events with the specified
54/// variables, using the specified PDF model. A prototype dataset (if provided)
55/// is not cloned and still belongs to the caller. The contents and shape
56/// of this dataset can be changed between calls to generate() as long as the
57/// expected columns to be copied to the generated dataset are present.
58/// Any argument supplied in the forceDirect RooArgSet are always offered
59/// for internal generation to the p.d.f., even if this is deemed unsafe by
60/// the logic of RooGenContext.
61
63 const RooDataSet *prototype, const RooArgSet* auxProto,
64 bool verbose, const RooArgSet* forceDirect) :
65 RooAbsGenContext(model,vars,prototype,auxProto,verbose),
66 _updateFMaxPerEvent(0)
67{
68 cxcoutI(Generation) << "RooGenContext::ctor() setting up event generator context for p.d.f. " << model.GetName()
69 << " for generation of observable(s) " << vars ;
70 if (prototype) ccxcoutI(Generation) << " with prototype data for " << *prototype->get() ;
71 if (auxProto && !auxProto->empty()) ccxcoutI(Generation) << " with auxiliary prototypes " << *auxProto ;
72 if (forceDirect && !forceDirect->empty()) ccxcoutI(Generation) << " with internal generation forced for observables " << *forceDirect ;
73 ccxcoutI(Generation) << endl ;
74
75
76 // Clone the model and all nodes that it depends on so that this context
77 // is independent of any existing objects.
78 RooArgSet nodes(model,model.GetName());
79 if (nodes.snapshot(_cloneSet, true)) {
80 coutE(Generation) << "RooGenContext::RooGenContext(" << GetName() << ") Couldn't deep-clone PDF, abort," << endl ;
82 }
83
84 // Find the clone in the snapshot list
85 _pdfClone = static_cast<RooAbsPdf*>(_cloneSet.find(model.GetName()));
87
88 // Optionally fix RooAddPdf normalizations
89 if (prototype&&_pdfClone->dependsOn(*prototype->get())) {
90 RooArgSet fullNormSet(vars) ;
91 fullNormSet.add(*prototype->get()) ;
93 }
94
95 // Analyze the list of variables to generate...
96 _isValid= true;
97 const RooAbsArg *arg = nullptr;
98 for(RooAbsArg * tmp : vars) {
99 if(!_isValid) break;
100
101 // is this argument derived?
102 if(!tmp->isFundamental()) {
103 coutE(Generation) << "RooGenContext::ctor(): cannot generate values for derived \"" << tmp->GetName() << "\"" << endl;
104 _isValid= false;
105 continue;
106 }
107 // lookup this argument in the cloned set of PDF dependents
108 arg= static_cast<RooAbsArg const*>(_cloneSet.find(tmp->GetName()));
109 if(nullptr == arg) {
110 coutI(Generation) << "RooGenContext::ctor() WARNING model does not depend on \"" << tmp->GetName()
111 << "\" which will have uniform distribution" << endl;
112 _uniformVars.add(*tmp);
113 }
114 else {
115
116 // does the model depend on this variable directly, ie, like "x" in
117 // f(x) or f(x,g(x,y)) or even f(x,x) ?
118 const RooAbsArg *direct= arg ;
119 if (forceDirect==nullptr || !forceDirect->find(direct->GetName())) {
120 if (!_pdfClone->isDirectGenSafe(*arg)) {
121 cxcoutD(Generation) << "RooGenContext::ctor() observable " << arg->GetName() << " has been determined to be unsafe for internal generation" << endl;
122 direct=nullptr ;
123 }
124 }
125
126 // does the model depend indirectly on this variable through an lvalue chain?
127 // otherwise, this variable will have to be generated with accept/reject
128 if(direct) {
129 _directVars.add(*arg);
130 } else {
131 _otherVars.add(*arg);
132 }
133 }
134 }
135
136 if(!isValid()) {
137 coutE(Generation) << "RooGenContext::ctor() constructor failed with errors" << endl;
138 return;
139 }
140
141 // At this point directVars are all variables that are safe to be generated directly
142 // otherVars are all variables that are _not_ safe to be generated directly
143 if (!_directVars.empty()) {
144 cxcoutD(Generation) << "RooGenContext::ctor() observables " << _directVars << " are safe for internal generator (if supported by p.d.f)" << endl ;
145 }
146 if (!_otherVars.empty()) {
147 cxcoutD(Generation) << "RooGenContext::ctor() observables " << _otherVars << " are NOT safe for internal generator (if supported by p.d.f)" << endl ;
148 }
149
150 // If PDF depends on prototype data, direct generator cannot use static initialization
151 // in initGenerator()
152 bool staticInitOK = !_pdfClone->dependsOn(_protoVars) ;
153 if (!staticInitOK) {
154 cxcoutD(Generation) << "RooGenContext::ctor() Model depends on supplied protodata observables, static initialization of internal generator will not be allowed" << endl ;
155 }
156
157 // Can the model generate any of the direct variables itself?
158 RooArgSet generatedVars;
159 if (!_directVars.empty()) {
160 _code= _pdfClone->getGenerator(_directVars,generatedVars,staticInitOK);
161
162 cxcoutD(Generation) << "RooGenContext::ctor() Model indicates that it can internally generate observables "
163 << generatedVars << " with configuration identifier " << _code << endl ;
164 } else {
165 _code = 0 ;
166 }
167
168 // Move variables which cannot be generated into the list to be generated with accept/reject
169 _directVars.remove(generatedVars);
171
172 // Update _directVars to only include variables that will actually be directly generated
174 _directVars.add(generatedVars);
175
176 cxcoutI(Generation) << "RooGenContext::ctor() Context will" ;
177 if (!_directVars.empty()) ccxcoutI(Generation) << " let model internally generate observables " << _directVars ;
178 if (!_directVars.empty() && !_otherVars.empty()) ccxcoutI(Generation) << " and" ;
179 if (!_otherVars.empty()) ccxcoutI(Generation) << " generate variables " << _otherVars << " with accept/reject sampling" ;
180 if (!_uniformVars.empty()) ccxcoutI(Generation) << ", non-observable variables " << _uniformVars << " will be generated with flat distribution" ;
181 ccxcoutI(Generation) << endl ;
182
183 // initialize the accept-reject generator
184 RooArgSet depList;
186 depList.remove(_otherVars);
187
188 TString nname(_pdfClone->GetName()) ;
189 nname.Append("_AccRej") ;
190 TString ntitle(_pdfClone->GetTitle()) ;
191 ntitle.Append(" (Accept/Reject)") ;
192
193
194 RooArgSet protoDeps;
195 model.getObservables(&_protoVars, protoDeps);
196
197
198 if (_protoVars.empty()) {
199
200 // No prototype variables
201
202 if(depList.empty()) {
203 // All variable are generated with accept-reject
204
205 // Check if PDF supports maximum finding
206 Int_t maxFindCode = _pdfClone->getMaxVal(_otherVars) ;
207 if (maxFindCode != 0) {
208 coutI(Generation) << "RooGenContext::ctor() no prototype data provided, all observables are generated with numerically and "
209 << "model supports analytical maximum findin:, can provide analytical pdf maximum to numeric generator" << endl ;
210 double maxVal = _pdfClone->maxVal(maxFindCode) / _pdfClone->getNorm(&_theEvent) ;
211 _maxVar = std::make_unique<RooRealVar>("funcMax","function maximum",maxVal) ;
212 cxcoutD(Generation) << "RooGenContext::ctor() maximum value returned by RooAbsPdf::maxVal() is " << maxVal << endl ;
213 }
214 }
215
216 if (!_otherVars.empty()) {
217 _pdfClone->getVal(&vars) ; // WVE debug
218 _acceptRejectFunc= std::make_unique<RooRealIntegral>(nname,ntitle,*_pdfClone,depList,&vars);
219 cxcoutI(Generation) << "RooGenContext::ctor() accept/reject sampling function is " << _acceptRejectFunc->GetName() << endl ;
220 } else {
221 _acceptRejectFunc = nullptr;
222 }
223
224 } else {
225
226 // Generation _with_ prototype variable
227 depList.remove(_protoVars,true,true);
228 _acceptRejectFunc= std::unique_ptr<RooAbsReal>{_pdfClone->createIntegral(depList,vars)};
229 cxcoutI(Generation) << "RooGenContext::ctor() accept/reject sampling function is " << _acceptRejectFunc->GetName() << endl ;
230
231 // Check if PDF supports maximum finding for the entire phase space
232 RooArgSet allVars(_otherVars);
233 allVars.add(_directVars);
234 Int_t maxFindCode = _pdfClone->getMaxVal(allVars) ;
235 if (maxFindCode != 0) {
236 // Special case: PDF supports max-finding in otherVars, no need to scan other+proto space for maximum
237 coutI(Generation) << "RooGenContext::ctor() prototype data provided, and "
238 << "model supports analytical maximum finding in the full phase space: "
239 << "can provide analytical pdf maximum to numeric generator" << endl ;
240 double maxVal = _pdfClone->maxVal(maxFindCode) / _pdfClone->getNorm(&allVars);
241 _maxVar = std::make_unique<RooRealVar>("funcMax", "function maximum", maxVal);
242 } else {
243 maxFindCode = _pdfClone->getMaxVal(_otherVars) ;
244 if (maxFindCode != 0) {
245 _updateFMaxPerEvent = maxFindCode ;
246 coutI(Generation) << "RooGenContext::ctor() prototype data provided, and "
247 << "model supports analytical maximum finding in the variables which are not"
248 << " internally generated. Can provide analytical pdf maximum to numeric "
249 << "generator" << endl;
250 cxcoutD(Generation) << "RooGenContext::ctor() maximum value must be reevaluated for each "
251 << "event with configuration code " << maxFindCode << endl ;
252 _maxVar = std::make_unique<RooRealVar>("funcMax","function maximum",1) ;
253 }
254 }
255
256 if (!_maxVar) {
257
258 // Regular case: First find maximum in other+proto space
259 RooArgSet otherAndProto(_otherVars) ;
260
261 otherAndProto.add(protoDeps) ;
262
263 if (!_otherVars.empty()) {
264
265 cxcoutD(Generation) << "RooGenContext::ctor() prototype data provided, observables are generated numerically no "
266 << "analytical estimate of maximum function value provided by model, must determine maximum value through initial sampling space "
267 << "of accept/reject observables plus prototype observables: " << otherAndProto << endl ;
268
269 // Calculate maximum in other+proto space if there are any accept/reject generated observables
270 std::unique_ptr<RooAbsNumGenerator> maxFinder{RooNumGenFactory::instance().createSampler(*_acceptRejectFunc,otherAndProto,RooArgSet(_protoVars),
271 *model.getGeneratorConfig(),_verbose)};
272// RooAcceptReject maxFinder(*_acceptRejectFunc,otherAndProto,RooNumGenConfig::defaultConfig(),_verbose) ;
273 double max = maxFinder->getFuncMax() ;
274 _maxVar = std::make_unique<RooRealVar>("funcMax","function maximum",max) ;
275
276 if (max==0) {
277 oocoutE(nullptr, Generation) << "RooGenContext::ctor(" << model.GetName()
278 << ") ERROR: generating conditional p.d.f. which requires prior knowledge of function maximum, "
279 << "but chosen numeric generator (" << maxFinder->generatorName() << ") does not support maximum finding" << endl ;
280 throw string("RooGenContext::ctor()") ;
281 }
282
283 cxcoutD(Generation) << "RooGenContext::ctor() maximum function value found through initial sampling is " << max << endl ;
284 }
285 }
286
287 }
288
291 cxcoutD(Generation) << "RooGenContext::ctor() creating MC sampling generator " << _generator->generatorName() << " from function for observables " << _otherVars << endl ;
292 //_generator= new RooAcceptReject(*_acceptRejectFunc,_otherVars,RooNumGenConfig::defaultConfig(),_verbose,_maxVar);
293 } else {
294 _generator = nullptr ;
295 }
296
298}
299
301
302////////////////////////////////////////////////////////////////////////////////
303/// Attach the cloned model to the event buffer we will be filling.
304
306{
308 if (_acceptRejectFunc) {
309 _acceptRejectFunc->recursiveRedirectServers(args,false) ; // WVE DEBUG
310 }
311
312 // Attach the RooAcceptReject generator the event buffer
313 if (_generator) {
314 _generator->attachParameters(args) ;
315 }
316
317}
318
319
320
321////////////////////////////////////////////////////////////////////////////////
322/// Perform one-time initialization of the generator context
323
325{
326 for (auto* arg : theEvent) {
327 arg->setOperMode(RooAbsArg::ADirty) ;
328 }
329
330 attach(theEvent) ;
331
332 // Reset the cloned model's error counters.
334
335 // Initialize the PDFs internal generator
336 if (!_directVars.empty()) {
337 cxcoutD(Generation) << "RooGenContext::initGenerator() initializing internal generator of model with code " << _code << endl ;
339 }
340}
341
342
343////////////////////////////////////////////////////////////////////////////////
344/// Generate one event. The 'remaining' integer is not used other than
345/// for printing messages
346
348{
349 if(!_otherVars.empty()) {
350 // call the accept-reject generator to generate its variables
351
352 if (_updateFMaxPerEvent!=0) {
354 cxcoutD(Generation) << "RooGenContext::initGenerator() reevaluation of maximum function value is required for each event, new value is " << max << endl ;
355 _maxVar->setVal(max) ;
356 }
357
358 if (_generator) {
359 double resampleRatio(1) ;
360 const RooArgSet *subEvent= _generator->generateEvent(remaining,resampleRatio);
361 if (resampleRatio<1) {
362 coutI(Generation) << "RooGenContext::generateEvent INFO: accept/reject generator requests resampling of previously produced events by factor "
363 << resampleRatio << " due to increased maximum weight" << endl ;
364 resampleData(resampleRatio) ;
365 }
366 if(nullptr == subEvent) {
367 coutE(Generation) << "RooGenContext::generateEvent ERROR accept/reject generator failed" << endl ;
368 return;
369 }
370 theEvent.assignValueOnly(*subEvent) ;
371 //theEvent= *subEvent;
372
373 }
374 }
375
376 // Use the model's optimized generator, if one is available.
377 // The generator writes directly into our local 'event' since we attached it above.
378 if(!_directVars.empty()) {
380 }
381
382 // Generate uniform variables (non-dependents)
383 for(auto * uniVar : _uniformVars) {
384 auto * arglv = dynamic_cast<RooAbsLValue*>(uniVar);
385 if (!arglv) {
386 coutE(Generation) << "RooGenContext::generateEvent(" << GetName() << ") ERROR: uniform variable " << uniVar->GetName() << " is not an lvalue" << endl ;
388 }
389 arglv->randomize() ;
390 }
391 theEvent.assign(_uniformVars) ;
392
393}
394
395
396////////////////////////////////////////////////////////////////////////////////
397/// Printing interface
398
399void RooGenContext::printMultiline(ostream &os, Int_t content, bool verbose, TString indent) const
400{
401 RooAbsGenContext::printMultiline(os,content,verbose,indent);
402 os << indent << " --- RooGenContext --- " << endl ;
403 os << indent << "Using PDF ";
405
406 if(verbose) {
407 os << indent << "Use PDF generator for " << _directVars << endl ;
408 os << indent << "Use MC sampling generator " << (_generator ? _generator->generatorName() : "<none>") << " for " << _otherVars << endl ;
409 if (!_protoVars.empty()) {
410 os << indent << "Prototype observables are " << _protoVars << endl ;
411 }
412 }
413}
#define coutI(a)
#define cxcoutI(a)
#define cxcoutD(a)
#define oocoutE(o, a)
#define coutE(a)
#define ccxcoutI(a)
#define ClassImp(name)
Definition Rtypes.h:382
static void indent(ostringstream &buf, int indent_level)
winID h direct
Common abstract base class for objects that represent a value and a "shape" in RooFit.
Definition RooAbsArg.h:79
bool dependsOn(const RooAbsCollection &serverList, const RooAbsArg *ignoreArg=nullptr, bool valueOnly=false) const
Test whether we depend on (ie, are served by) any object in the specified collection.
bool recursiveRedirectServers(const RooAbsCollection &newServerList, bool mustReplaceAll=false, bool nameChange=false, bool recurseInNewSet=true)
Recursively replace all servers with the new servers in newSet.
void setOperMode(OperMode mode, bool recurseADirty=true)
Set the operation mode of this node.
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 void removeAll()
Remove all arguments from our set, deleting them if we own them.
virtual bool remove(const RooAbsArg &var, bool silent=false, bool matchByNameOnly=false)
Remove the specified argument from our list.
RooAbsCollection & assignValueOnly(const RooAbsCollection &other, bool forceIfSizeOne=false)
Sets the value of any argument in our set that also appears in the other set.
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to list.
void assign(const RooAbsCollection &other) const
Sets the value, cache and constant attribute of any argument in our set that also appears in the othe...
RooAbsArg * find(const char *name) const
Find object with given name in list.
Abstract base class for generator contexts of RooAbsPdf objects.
RooArgSet _theEvent
Pointer to observable event being generated.
void printMultiline(std::ostream &os, Int_t contents, bool verbose=false, TString indent="") const override
Interface for multi-line printing.
RooArgSet _protoVars
Prototype observables.
bool _verbose
Verbose messaging?
bool _isValid
Is context in valid state?
bool isValid() const
void resampleData(double &ratio)
Rescale existing output buffer with given ratio.
Abstract base class for objects that are lvalues, i.e.
virtual std::string const & generatorName() const =0
Return unique name of generator implementation.
Abstract interface for all probability density functions.
Definition RooAbsPdf.h:40
double getNorm(const RooArgSet &nset) const
Get normalisation term needed to normalise the raw values returned by getVal().
Definition RooAbsPdf.h:195
virtual void generateEvent(Int_t code)
Interface for generation of an event using the algorithm corresponding to the specified code.
virtual void resetErrorCounters(Int_t resetValue=10)
Reset error counter to given value, limiting the number of future error messages for this pdf to 'res...
virtual bool isDirectGenSafe(const RooAbsArg &arg) const
Check if given observable can be safely generated using the pdfs internal generator mechanism (if tha...
const RooNumGenConfig * getGeneratorConfig() const
Return the numeric MC generator configuration used for this object.
virtual void initGenerator(Int_t code)
Interface for one-time initialization to setup the generator for the specified code.
virtual Int_t getGenerator(const RooArgSet &directVars, RooArgSet &generateVars, bool staticInitOK=true) const
Load generatedVars with the subset of directVars that we can generate events for, and return a code t...
RooFit::OwningPtr< RooAbsReal > createIntegral(const RooArgSet &iset, const RooCmdArg &arg1, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={}, const RooCmdArg &arg6={}, const RooCmdArg &arg7={}, const RooCmdArg &arg8={}) const
Create an object that represents the integral of the function over one or more observables listed in ...
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition RooAbsReal.h:103
virtual void fixAddCoefNormalization(const RooArgSet &addNormSet=RooArgSet(), bool force=true)
Fix the interpretation of the coefficient of any RooAddPdf component in the expression tree headed by...
virtual double maxVal(Int_t code) const
Return maximum value for set of observables identified by code assigned in getMaxVal.
virtual Int_t getMaxVal(const RooArgSet &vars) const
Advertise capability to determine maximum value of function for given set of observables.
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:24
RooArgSet * snapshot(bool deepCopy=true) const
Use RooAbsCollection::snapshot(), but return as RooArgSet.
Definition RooArgSet.h:154
Container class to hold unbinned data.
Definition RooDataSet.h:33
const RooArgSet * get(Int_t index) const override
Return RooArgSet with coordinates of event 'index'.
static void softAbort()
Soft abort function that interrupts macro execution but doesn't kill ROOT.
Implements a universal generator context for all RooAbsPdf classes that do not have or need a special...
void printMultiline(std::ostream &os, Int_t content, bool verbose=false, TString indent="") const override
Printing interface.
void generateEvent(RooArgSet &theEvent, Int_t remaining) override
Generate one event.
RooArgSet _otherVars
List of observables generated internally, randomly, and by accept/reject sampling.
Int_t _code
Internal generation code.
RooArgSet _uniformVars
void initGenerator(const RooArgSet &theEvent) override
Perform one-time initialization of the generator context.
~RooGenContext() override
std::unique_ptr< RooAbsNumGenerator > _generator
MC sampling generation engine.
RooArgSet _directVars
RooAbsPdf * _pdfClone
Clone of input p.d.f.
RooArgSet _cloneSet
Clone of all nodes of input p.d.f.
Int_t _updateFMaxPerEvent
If true, maximum p.d.f value needs to be recalculated for each event.
std::unique_ptr< RooRealVar > _maxVar
Variable holding maximum value of p.d.f.
std::unique_ptr< RooAbsReal > _acceptRejectFunc
Projection function to be passed to accept/reject sampler.
void attach(const RooArgSet &params) override
Attach the cloned model to the event buffer we will be filling.
RooGenContext(const RooAbsPdf &model, const RooArgSet &vars, const RooDataSet *prototype=nullptr, const RooArgSet *auxProto=nullptr, bool verbose=false, const RooArgSet *forceDirect=nullptr)
Initialize a new context for generating events with the specified variables, using the specified PDF ...
RooAbsNumGenerator * createSampler(RooAbsReal &func, const RooArgSet &genVars, const RooArgSet &condVars, const RooNumGenConfig &config, bool verbose=false, RooAbsReal *maxFuncVal=nullptr)
Construct a numeric integrator instance that operates on function 'func' and is configured with 'conf...
static RooNumGenFactory & instance()
Static method returning reference to singleton instance of factory.
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,...
const char * GetName() const override
Returns name of object.
Definition TNamed.h:47
const char * GetTitle() const override
Returns title of object.
Definition TNamed.h:48
Basic string class.
Definition TString.h:139
TString & Append(const char *cs)
Definition TString.h:572