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