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
30#include "RooFit.h"
31#include "Riostream.h"
32
33#include "RooAcceptReject.h"
34#include "RooAbsReal.h"
35#include "RooCategory.h"
36#include "RooRealVar.h"
37#include "RooDataSet.h"
38#include "RooRandom.h"
39#include "RooErrorHandler.h"
40
41#include "TString.h"
42#include "TIterator.h"
43#include "RooMsgService.h"
44#include "TClass.h"
45#include "TFoam.h"
46#include "RooRealBinding.h"
47#include "RooNumGenFactory.h"
48#include "RooNumGenConfig.h"
49
50#include <assert.h>
51
52using namespace std;
53
55 ;
56
57
58////////////////////////////////////////////////////////////////////////////////
59/// Register RooIntegrator1D, is parameters and capabilities with RooNumIntFactory
60
62{
63 RooRealVar nTrial0D("nTrial0D","Number of trial samples for cat-only generation",100,0,1e9) ;
64 RooRealVar nTrial1D("nTrial1D","Number of trial samples for 1-dim generation",1000,0,1e9) ;
65 RooRealVar nTrial2D("nTrial2D","Number of trial samples for 2-dim generation",100000,0,1e9) ;
66 RooRealVar nTrial3D("nTrial3D","Number of trial samples for N-dim generation",10000000,0,1e9) ;
67
69 fact.storeProtoSampler(proto,RooArgSet(nTrial0D,nTrial1D,nTrial2D,nTrial3D)) ;
70}
71
72
73
74////////////////////////////////////////////////////////////////////////////////
75/// Initialize an accept-reject generator for the specified distribution function,
76/// which must be non-negative but does not need to be normalized over the
77/// variables to be generated, genVars. The function and its dependents are
78/// cloned and so will not be disturbed during the generation process.
79
80RooAcceptReject::RooAcceptReject(const RooAbsReal &func, const RooArgSet &genVars, const RooNumGenConfig& config, Bool_t verbose, const RooAbsReal* maxFuncVal) :
81 RooAbsNumGenerator(func,genVars,verbose,maxFuncVal), _nextCatVar(0), _nextRealVar(0)
82{
83 _minTrialsArray[0] = static_cast<Int_t>(config.getConfigSection("RooAcceptReject").getRealValue("nTrial0D")) ;
84 _minTrialsArray[1] = static_cast<Int_t>(config.getConfigSection("RooAcceptReject").getRealValue("nTrial1D")) ;
85 _minTrialsArray[2] = static_cast<Int_t>(config.getConfigSection("RooAcceptReject").getRealValue("nTrial2D")) ;
86 _minTrialsArray[3] = static_cast<Int_t>(config.getConfigSection("RooAcceptReject").getRealValue("nTrial3D")) ;
87
90 RooAbsCategory* cat ;
91 _catSampleMult = 1 ;
92 while((cat=(RooAbsCategory*)iter->Next())) {
93 _catSampleMult *= cat->numTypes() ;
94 }
95 delete iter ;
96
97
98 // calculate the minimum number of trials needed to estimate our integral and max value
99 if (!_funcMaxVal) {
100
101 if(_realSampleDim > 3) {
103 coutW(Generation) << fName << "::" << ClassName() << ": WARNING: generating " << _realSampleDim
104 << " variables with accept-reject may not be accurate" << endl;
105 }
106 else {
108 }
109 if (_realSampleDim > 1) {
110 coutW(Generation) << "RooAcceptReject::ctor(" << fName
111 << ") WARNING: performing accept/reject sampling on a p.d.f in "
112 << _realSampleDim << " dimensions without prior knowledge on maximum value "
113 << "of p.d.f. Determining maximum value by taking " << _minTrials
114 << " trial samples. If p.d.f contains sharp peaks smaller than average "
115 << "distance between trial sampling points these may be missed and p.d.f. "
116 << "may be sampled incorrectly." << endl ;
117 }
118 } else {
119 // No trials needed if we know the maximum a priori
120 _minTrials=0 ;
121 }
122
123 // Need to fix some things here
124 if (_minTrials>10000) {
125 coutW(Generation) << "RooAcceptReject::ctor(" << fName << "): WARNING: " << _minTrials << " trial samples requested by p.d.f for "
126 << _realSampleDim << "-dimensional accept/reject sampling, this may take some time" << endl ;
127 }
128
129 // print a verbose summary of our configuration, if requested
130 if(_verbose) {
131 coutI(Generation) << fName << "::" << ClassName() << ":" << endl
132 << " Initializing accept-reject generator for" << endl << " ";
134 if (_funcMaxVal) {
135 ccoutI(Generation) << " Function maximum provided, no trial sampling performed" << endl ;
136 } else {
137 ccoutI(Generation) << " Real sampling dimension is " << _realSampleDim << endl;
138 ccoutI(Generation) << " Category sampling multiplier is " << _catSampleMult << endl ;
139 ccoutI(Generation) << " Min sampling trials is " << _minTrials << endl;
140 }
141 if (_catVars.getSize()>0) {
142 ccoutI(Generation) << " Will generate category vars "<< _catVars << endl ;
143 }
144 if (_realVars.getSize()>0) {
145 ccoutI(Generation) << " Will generate real vars " << _realVars << endl ;
146 }
147 }
148 // create iterators for the new sets
151 assert(0 != _nextCatVar && 0 != _nextRealVar);
152
153 // initialize our statistics
154 _maxFuncVal= 0;
155 _funcSum= 0;
156 _totalEvents= 0;
157 _eventsUsed= 0;
158}
159
160
161
162////////////////////////////////////////////////////////////////////////////////
163/// Destructor
164
166{
167 delete _nextCatVar;
168 delete _nextRealVar;
169}
170
171
172
173////////////////////////////////////////////////////////////////////////////////
174/// Return a pointer to a generated event. The caller does not own the event and it
175/// will be overwritten by a subsequent call. The input parameter 'remaining' should
176/// contain your best guess at the total number of subsequent events you will request.
177
178const RooArgSet *RooAcceptReject::generateEvent(UInt_t remaining, Double_t& resampleRatio)
179{
180 // are we actually generating anything? (the cache always contains at least our function value)
181 const RooArgSet *event= _cache->get();
182 if(event->getSize() == 1) return event;
183
184 if (!_funcMaxVal) {
185 // Generation with empirical maximum determination
186
187 // first generate enough events to get reasonable estimates for the integral and
188 // maximum function value
189
190 while(_totalEvents < _minTrials) {
192
193 // Limit cache size to 1M events
194 if (_cache->numEntries()>1000000) {
195 coutI(Generation) << "RooAcceptReject::generateEvent: resetting event cache" << endl ;
196 _cache->reset() ;
197 _eventsUsed = 0 ;
198 }
199 }
200
201 event= 0;
202 Double_t oldMax2(_maxFuncVal);
203 while(0 == event) {
204 // Use any cached events first
205 if (_maxFuncVal>oldMax2) {
206 cxcoutD(Generation) << "RooAcceptReject::generateEvent maxFuncVal has changed, need to resample already accepted events by factor"
207 << oldMax2 << "/" << _maxFuncVal << "=" << oldMax2/_maxFuncVal << endl ;
208 resampleRatio=oldMax2/_maxFuncVal ;
209 }
210 event= nextAcceptedEvent();
211 if(event) break;
212 // When we have used up the cache, start a new cache and add
213 // some more events to it.
214 _cache->reset();
215 _eventsUsed= 0;
216 // Calculate how many more events to generate using our best estimate of our efficiency.
217 // Always generate at least one more event so we don't get stuck.
218 if(_totalEvents*_maxFuncVal <= 0) {
219 coutE(Generation) << "RooAcceptReject::generateEvent: cannot estimate efficiency...giving up" << endl;
220 return 0;
221 }
222
224 Long64_t extra= 1 + (Long64_t)(1.05*remaining/eff);
225 cxcoutD(Generation) << "RooAcceptReject::generateEvent: adding " << extra << " events to the cache, eff = " << eff << endl;
226 Double_t oldMax(_maxFuncVal);
227 while(extra--) {
229 if((_maxFuncVal > oldMax)) {
230 cxcoutD(Generation) << "RooAcceptReject::generateEvent: estimated function maximum increased from "
231 << oldMax << " to " << _maxFuncVal << endl;
232 oldMax = _maxFuncVal ;
233 // Trim cache here
234 }
235 }
236 }
237
238 // Limit cache size to 1M events
239 if (_eventsUsed>1000000) {
240 _cache->reset() ;
241 _eventsUsed = 0 ;
242 }
243
244 } else {
245 // Generation with a priori maximum knowledge
247
248 // Generate enough trials to produce a single accepted event
249 event = 0 ;
250 while(0==event) {
252 event = nextAcceptedEvent() ;
253 }
254
255 }
256 return event;
257}
258
259
260
261////////////////////////////////////////////////////////////////////////////////
262/// Scan through events in the cache which have not been used yet,
263/// looking for the first accepted one which is added to the specified
264/// container. Return a pointer to the accepted event, or else zero
265/// if we use up the cache before we accept an event. The caller does
266/// not own the event and it will be overwritten by a subsequent call.
267
269{
270 const RooArgSet *event = 0;
271 while((event= _cache->get(_eventsUsed))) {
272 _eventsUsed++ ;
273 // accept this cached event?
276 //cout << " event number " << _eventsUsed << " has been rejected" << endl ;
277 continue;
278 }
279 //cout << " event number " << _eventsUsed << " has been accepted" << endl ;
280 // copy this event into the output container
281 if(_verbose && (_eventsUsed%1000==0)) {
282 cerr << "RooAcceptReject: accepted event (used " << _eventsUsed << " of "
283 << _cache->numEntries() << " so far)" << endl;
284 }
285 break;
286 }
287 //cout << "accepted event " << _eventsUsed << " of " << _cache->numEntries() << endl ;
288 return event;
289}
290
291
292
293////////////////////////////////////////////////////////////////////////////////
294/// Add a trial event to our cache and update our estimates
295/// of the function maximum value and integral.
296
298{
299 // randomize each discrete argument
301 RooCategory *cat = 0;
302 while((cat= (RooCategory*)_nextCatVar->Next())) cat->randomize();
303
304 // randomize each real argument
306 RooRealVar *real = 0;
307 while((real= (RooRealVar*)_nextRealVar->Next())) real->randomize();
308
309 // calculate and store our function value at this new point
310 Double_t val= _funcClone->getVal();
311 _funcValPtr->setVal(val);
312
313 // Update the estimated integral and maximum value. Increase our
314 // maximum estimate slightly to give a safety margin with a
315 // corresponding loss of efficiency.
316 if(val > _maxFuncVal) _maxFuncVal= 1.05*val;
317 _funcSum+= val;
318
319 // fill a new entry in our cache dataset for this point
320 _cache->fill();
321 _totalEvents++;
322
323 if (_verbose &&_totalEvents%10000==0) {
324 cerr << "RooAcceptReject: generated " << _totalEvents << " events so far." << endl ;
325 }
326
327}
328
330{
331 // Empirically determine maximum value of function by taking a large number
332 // of samples. The actual number depends on the number of dimensions in which
333 // the sampling occurs
334
335 // Generate the minimum required number of samples for a reliable maximum estimate
336 while(_totalEvents < _minTrials) {
338
339 // Limit cache size to 1M events
340 if (_cache->numEntries()>1000000) {
341 coutI(Generation) << "RooAcceptReject::getFuncMax: resetting event cache" << endl ;
342 _cache->reset() ;
343 _eventsUsed = 0 ;
344 }
345 }
346
347 return _maxFuncVal ;
348}
349
ROOT::R::TRInterface & r
Definition: Object.C:4
#define coutI(a)
Definition: RooMsgService.h:30
#define cxcoutD(a)
Definition: RooMsgService.h:81
#define coutW(a)
Definition: RooMsgService.h:32
#define coutE(a)
Definition: RooMsgService.h:33
#define ccoutI(a)
Definition: RooMsgService.h:38
int Int_t
Definition: RtypesCore.h:43
unsigned int UInt_t
Definition: RtypesCore.h:44
bool Bool_t
Definition: RtypesCore.h:61
double Double_t
Definition: RtypesCore.h:57
long long Long64_t
Definition: RtypesCore.h:71
#define ClassImp(name)
Definition: Rtypes.h:361
const char * proto
Definition: civetweb.c:16604
virtual void randomize(const char *rangeName=0)
Randomize current value.
RooAbsCategory is the base class for objects that represent a discrete value with a finite number of ...
Int_t numTypes(const char *=0) const
Return number of types defined (in range named rangeName if rangeName!=0)
Int_t getSize() const
Double_t getRealValue(const char *name, Double_t defVal=0, Bool_t verbose=kFALSE) const
Get value of a RooAbsReal stored in set with given name.
TIterator * createIterator(Bool_t dir=kIterForward) const
TIterator-style iteration over contained elements.
virtual void reset()
Definition: RooAbsData.cxx:313
virtual void fill()
Definition: RooAbsData.cxx:299
virtual Int_t numEntries() const
Definition: RooAbsData.cxx:306
Class RooAbsNumGenerator is the abstract base class for MC event generator implementations like RooAc...
const RooAbsReal * _funcMaxVal
RooRealVar * _funcValPtr
virtual void randomize(const char *rangeName=0)
Set a new value sampled from a uniform distribution over the fit range.
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:60
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:90
Class RooAcceptReject is a generic toy monte carlo generator implement the accept/reject sampling tec...
UInt_t _minTrialsArray[4]
virtual ~RooAcceptReject()
Destructor.
TIterator * _nextCatVar
static void registerSampler(RooNumGenFactory &fact)
Register RooIntegrator1D, is parameters and capabilities with RooNumIntFactory.
const RooArgSet * generateEvent(UInt_t remaining, Double_t &resampleRatio)
Return a pointer to a generated event.
Double_t _maxFuncVal
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.
TIterator * _nextRealVar
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
RooCategory is an object to represent discrete states.
Definition: RooCategory.h:23
virtual 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_t 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_t uniform(TRandom *generator=randomGenerator())
Return a number uniformly distributed from (0,1)
Definition: RooRandom.cxx:83
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:35
virtual void setVal(Double_t value)
Set value of variable to 'value'.
Definition: RooRealVar.cxx:225
Iterator abstract base class.
Definition: TIterator.h:30
virtual void Reset()=0
virtual TObject * Next()=0
TString fName
Definition: TNamed.h:32
virtual const char * ClassName() const
Returns name of class to which the object belongs.
Definition: TObject.cxx:128
@ Generation
Definition: RooGlobalFunc.h:67