Logo ROOT  
Reference Guide
RooAcceptReject.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 RooAcceptReject.cxx
19\class RooAcceptReject
20\ingroup Roofitcore
21
22Class RooAcceptReject is a generic toy monte carlo generator implement
23the accept/reject sampling technique on any positively valued function.
24The RooAcceptReject generator is used by the various generator context
25classes to take care of generation of observables for which p.d.fs
26do not define internal methods
27**/
28
29#include "Riostream.h"
30
31#include "RooAcceptReject.h"
32#include "RooAbsReal.h"
33#include "RooCategory.h"
34#include "RooRealVar.h"
35#include "RooDataSet.h"
36#include "RooRandom.h"
37#include "RooErrorHandler.h"
38
39#include "RooMsgService.h"
40#include "TFoam.h"
41#include "RooRealBinding.h"
42#include "RooNumGenFactory.h"
43#include "RooNumGenConfig.h"
44
45#include <assert.h>
46
47using namespace std;
48
50 ;
51
52
53////////////////////////////////////////////////////////////////////////////////
54/// Register RooIntegrator1D, is parameters and capabilities with RooNumIntFactory
55
57{
58 RooRealVar nTrial0D("nTrial0D","Number of trial samples for cat-only generation",100,0,1e9) ;
59 RooRealVar nTrial1D("nTrial1D","Number of trial samples for 1-dim generation",1000,0,1e9) ;
60 RooRealVar nTrial2D("nTrial2D","Number of trial samples for 2-dim generation",100000,0,1e9) ;
61 RooRealVar nTrial3D("nTrial3D","Number of trial samples for N-dim generation",10000000,0,1e9) ;
62
64 fact.storeProtoSampler(proto,RooArgSet(nTrial0D,nTrial1D,nTrial2D,nTrial3D)) ;
65}
66
67
68
69////////////////////////////////////////////////////////////////////////////////
70/// Initialize an accept-reject generator for the specified distribution function,
71/// which must be non-negative but does not need to be normalized over the
72/// variables to be generated, genVars. The function and its dependents are
73/// cloned and so will not be disturbed during the generation process.
74
75RooAcceptReject::RooAcceptReject(const RooAbsReal &func, const RooArgSet &genVars, const RooNumGenConfig& config, bool verbose, const RooAbsReal* maxFuncVal) :
76 RooAbsNumGenerator(func,genVars,verbose,maxFuncVal)
77{
78 _minTrialsArray[0] = static_cast<Int_t>(config.getConfigSection("RooAcceptReject").getRealValue("nTrial0D")) ;
79 _minTrialsArray[1] = static_cast<Int_t>(config.getConfigSection("RooAcceptReject").getRealValue("nTrial1D")) ;
80 _minTrialsArray[2] = static_cast<Int_t>(config.getConfigSection("RooAcceptReject").getRealValue("nTrial2D")) ;
81 _minTrialsArray[3] = static_cast<Int_t>(config.getConfigSection("RooAcceptReject").getRealValue("nTrial3D")) ;
82
84 _catSampleMult = 1 ;
85 for (auto * cat : static_range_cast<RooAbsCategory*>(_catVars)) {
86 _catSampleMult *= cat->numTypes() ;
87 }
88
89
90 // calculate the minimum number of trials needed to estimate our integral and max value
91 if (!_funcMaxVal) {
92
93 if(_realSampleDim > 3) {
95 coutW(Generation) << fName << "::" << ClassName() << ": WARNING: generating " << _realSampleDim
96 << " variables with accept-reject may not be accurate" << endl;
97 }
98 else {
100 }
101 if (_realSampleDim > 1) {
102 coutW(Generation) << "RooAcceptReject::ctor(" << fName
103 << ") WARNING: performing accept/reject sampling on a p.d.f in "
104 << _realSampleDim << " dimensions without prior knowledge on maximum value "
105 << "of p.d.f. Determining maximum value by taking " << _minTrials
106 << " trial samples. If p.d.f contains sharp peaks smaller than average "
107 << "distance between trial sampling points these may be missed and p.d.f. "
108 << "may be sampled incorrectly." << endl ;
109 }
110 } else {
111 // No trials needed if we know the maximum a priori
112 _minTrials=0 ;
113 }
114
115 // Need to fix some things here
116 if (_minTrials>10000) {
117 coutW(Generation) << "RooAcceptReject::ctor(" << fName << "): WARNING: " << _minTrials << " trial samples requested by p.d.f for "
118 << _realSampleDim << "-dimensional accept/reject sampling, this may take some time" << endl ;
119 }
120
121 // print a verbose summary of our configuration, if requested
122 if(_verbose) {
123 coutI(Generation) << fName << "::" << ClassName() << ":" << endl
124 << " Initializing accept-reject generator for" << endl << " ";
126 if (_funcMaxVal) {
127 ccoutI(Generation) << " Function maximum provided, no trial sampling performed" << endl ;
128 } else {
129 ccoutI(Generation) << " Real sampling dimension is " << _realSampleDim << endl;
130 ccoutI(Generation) << " Category sampling multiplier is " << _catSampleMult << endl ;
131 ccoutI(Generation) << " Min sampling trials is " << _minTrials << endl;
132 }
133 if (_catVars.getSize()>0) {
134 ccoutI(Generation) << " Will generate category vars "<< _catVars << endl ;
135 }
136 if (_realVars.getSize()>0) {
137 ccoutI(Generation) << " Will generate real vars " << _realVars << endl ;
138 }
139 }
140
141 // initialize our statistics
142 _maxFuncVal= 0;
143 _funcSum= 0;
144 _totalEvents= 0;
145 _eventsUsed= 0;
146}
147
148
149////////////////////////////////////////////////////////////////////////////////
150/// Return a pointer to a generated event. The caller does not own the event and it
151/// will be overwritten by a subsequent call. The input parameter 'remaining' should
152/// contain your best guess at the total number of subsequent events you will request.
153
154const RooArgSet *RooAcceptReject::generateEvent(UInt_t remaining, double& resampleRatio)
155{
156 // are we actually generating anything? (the cache always contains at least our function value)
157 const RooArgSet *event= _cache->get();
158 if(event->getSize() == 1) return event;
159
160 if (!_funcMaxVal) {
161 // Generation with empirical maximum determination
162
163 // first generate enough events to get reasonable estimates for the integral and
164 // maximum function value
165
166 while(_totalEvents < _minTrials) {
168
169 // Limit cache size to 1M events
170 if (_cache->numEntries()>1000000) {
171 coutI(Generation) << "RooAcceptReject::generateEvent: resetting event cache" << endl ;
172 _cache->reset() ;
173 _eventsUsed = 0 ;
174 }
175 }
176
177 event= 0;
178 double oldMax2(_maxFuncVal);
179 while(0 == event) {
180 // Use any cached events first
181 if (_maxFuncVal>oldMax2) {
182 cxcoutD(Generation) << "RooAcceptReject::generateEvent maxFuncVal has changed, need to resample already accepted events by factor"
183 << oldMax2 << "/" << _maxFuncVal << "=" << oldMax2/_maxFuncVal << endl ;
184 resampleRatio=oldMax2/_maxFuncVal ;
185 }
186 event= nextAcceptedEvent();
187 if(event) break;
188 // When we have used up the cache, start a new cache and add
189 // some more events to it.
190 _cache->reset();
191 _eventsUsed= 0;
192 // Calculate how many more events to generate using our best estimate of our efficiency.
193 // Always generate at least one more event so we don't get stuck.
194 if(_totalEvents*_maxFuncVal <= 0) {
195 coutE(Generation) << "RooAcceptReject::generateEvent: cannot estimate efficiency...giving up" << endl;
196 return 0;
197 }
198
199 double eff= _funcSum/(_totalEvents*_maxFuncVal);
200 Long64_t extra= 1 + (Long64_t)(1.05*remaining/eff);
201 cxcoutD(Generation) << "RooAcceptReject::generateEvent: adding " << extra << " events to the cache, eff = " << eff << endl;
202 double oldMax(_maxFuncVal);
203 while(extra--) {
205 if((_maxFuncVal > oldMax)) {
206 cxcoutD(Generation) << "RooAcceptReject::generateEvent: estimated function maximum increased from "
207 << oldMax << " to " << _maxFuncVal << endl;
208 oldMax = _maxFuncVal ;
209 // Trim cache here
210 }
211 }
212 }
213
214 // Limit cache size to 1M events
215 if (_eventsUsed>1000000) {
216 _cache->reset() ;
217 _eventsUsed = 0 ;
218 }
219
220 } else {
221 // Generation with a priori maximum knowledge
223
224 // Generate enough trials to produce a single accepted event
225 event = 0 ;
226 while(0==event) {
228 event = nextAcceptedEvent() ;
229 }
230
231 }
232 return event;
233}
234
235
236
237////////////////////////////////////////////////////////////////////////////////
238/// Scan through events in the cache which have not been used yet,
239/// looking for the first accepted one which is added to the specified
240/// container. Return a pointer to the accepted event, or else zero
241/// if we use up the cache before we accept an event. The caller does
242/// not own the event and it will be overwritten by a subsequent call.
243
245{
246 const RooArgSet *event = 0;
247 while((event= _cache->get(_eventsUsed))) {
248 _eventsUsed++ ;
249 // accept this cached event?
250 double r= RooRandom::uniform();
252 //cout << " event number " << _eventsUsed << " has been rejected" << endl ;
253 continue;
254 }
255 //cout << " event number " << _eventsUsed << " has been accepted" << endl ;
256 // copy this event into the output container
257 if(_verbose && (_eventsUsed%1000==0)) {
258 cerr << "RooAcceptReject: accepted event (used " << _eventsUsed << " of "
259 << _cache->numEntries() << " so far)" << endl;
260 }
261 break;
262 }
263 //cout << "accepted event " << _eventsUsed << " of " << _cache->numEntries() << endl ;
264 return event;
265}
266
267
268
269////////////////////////////////////////////////////////////////////////////////
270/// Add a trial event to our cache and update our estimates
271/// of the function maximum value and integral.
272
274{
275 // randomize each discrete argument
276 for(auto * cat : static_range_cast<RooCategory*>(_catVars)) cat->randomize();
277
278 // randomize each real argument
279 for(auto * real : static_range_cast<RooRealVar*>(_realVars)) real->randomize();
280
281 // calculate and store our function value at this new point
282 double val= _funcClone->getVal();
283 _funcValPtr->setVal(val);
284
285 // Update the estimated integral and maximum value. Increase our
286 // maximum estimate slightly to give a safety margin with a
287 // corresponding loss of efficiency.
288 if(val > _maxFuncVal) _maxFuncVal= 1.05*val;
289 _funcSum+= val;
290
291 // fill a new entry in our cache dataset for this point
292 _cache->fill();
293 _totalEvents++;
294
295 if (_verbose &&_totalEvents%10000==0) {
296 cerr << "RooAcceptReject: generated " << _totalEvents << " events so far." << endl ;
297 }
298
299}
300
302{
303 // Empirically determine maximum value of function by taking a large number
304 // of samples. The actual number depends on the number of dimensions in which
305 // the sampling occurs
306
307 // Generate the minimum required number of samples for a reliable maximum estimate
308 while(_totalEvents < _minTrials) {
310
311 // Limit cache size to 1M events
312 if (_cache->numEntries()>1000000) {
313 coutI(Generation) << "RooAcceptReject::getFuncMax: resetting event cache" << endl ;
314 _cache->reset() ;
315 _eventsUsed = 0 ;
316 }
317 }
318
319 return _maxFuncVal ;
320}
321
#define coutI(a)
Definition: RooMsgService.h:34
#define cxcoutD(a)
Definition: RooMsgService.h:85
#define coutW(a)
Definition: RooMsgService.h:36
#define coutE(a)
Definition: RooMsgService.h:37
#define ccoutI(a)
Definition: RooMsgService.h:42
long long Long64_t
Definition: RtypesCore.h:80
#define ClassImp(name)
Definition: Rtypes.h:375
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t r
const char * proto
Definition: civetweb.c:17502
double getRealValue(const char *name, double defVal=0.0, bool verbose=false) const
Get value of a RooAbsReal stored in set with given name.
Int_t getSize() const
Return the number of elements in the collection.
virtual void reset()
Definition: RooAbsData.cxx:381
virtual void fill()
Definition: RooAbsData.cxx:367
virtual Int_t numEntries() const
Return number of entries in dataset, i.e., count unweighted entries.
Definition: RooAbsData.cxx:374
Class RooAbsNumGenerator is the abstract base class for MC event generator implementations like RooAc...
RooArgSet _realVars
Sets of discrete and real valued observabeles.
RooDataSet * _cache
Dataset holding generared values of observables.
const RooAbsReal * _funcMaxVal
Container for maximum function value.
RooRealVar * _funcValPtr
RRVs storing function value in context and in output dataset.
RooAbsReal * _funcClone
Pointer to top level node of cloned function.
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:62
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:91
Class RooAcceptReject is a generic toy monte carlo generator implement the accept/reject sampling tec...
UInt_t _minTrialsArray[4]
Minimum number of trials samples for 1,2,3 dimensional problems.
UInt_t _minTrials
Minimum number of max.finding trials, total number of samples.
static void registerSampler(RooNumGenFactory &fact)
Register RooIntegrator1D, is parameters and capabilities with RooNumIntFactory.
double getFuncMax() override
const RooArgSet * nextAcceptedEvent()
Scan through events in the cache which have not been used yet, looking for the first accepted one whi...
void addEventToCache()
Add a trial event to our cache and update our estimates of the function maximum value and integral.
double _funcSum
Maximum function value found, and sum of all samples made.
UInt_t _totalEvents
Total number of function samples.
UInt_t _catSampleMult
Number of real and discrete dimensions to be sampled.
const RooArgSet * generateEvent(UInt_t remaining, double &resampleRatio) override
Return a pointer to a generated event.
UInt_t _eventsUsed
Accepted number of function samples.
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:56
const RooArgSet * get(Int_t index) const override
Return RooArgSet with coordinates of event 'index'.
RooNumGenConfig holds the configuration parameters of the various numeric integrators used by RooReal...
const RooArgSet & getConfigSection(const char *name) const
Retrieve configuration information specific to integrator with given name.
RooNumGenFactory is a factory to instantiate numeric integrators from a given function binding and a ...
bool storeProtoSampler(RooAbsNumGenerator *proto, const RooArgSet &defConfig)
Method accepting registration of a prototype numeric integrator along with a RooArgSet of its default...
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
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:40
void setVal(double value) override
Set value of variable to 'value'.
Definition: RooRealVar.cxx:254
TString fName
Definition: TNamed.h:32
virtual const char * ClassName() const
Returns name of class to which the object belongs.
Definition: TObject.cxx:207
@ Generation
Definition: RooGlobalFunc.h:61