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 
22 Class RooAcceptReject is a generic toy monte carlo generator implement
23 the accept/reject sampling technique on any positively valued function.
24 The RooAcceptReject generator is used by the various generator context
25 classes to take care of generation of observables for which p.d.fs
26 do 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 
52 using 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 
80 RooAcceptReject::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 
178 const 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) {
191  addEventToCache();
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--) {
228  addEventToCache();
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) {
251  addEventToCache() ;
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?
275  if(r*_maxFuncVal > _funcValPtr->getVal()) {
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
300  _nextCatVar->Reset();
301  RooCategory *cat = 0;
302  while((cat= (RooCategory*)_nextCatVar->Next())) cat->randomize();
303 
304  // randomize each real argument
305  _nextRealVar->Reset();
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) {
337  addEventToCache();
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 
RooAbsCategoryLValue::randomize
virtual void randomize(const char *rangeName=0)
Randomize current value.
Definition: RooAbsCategoryLValue.cxx:150
RooRealVar::setVal
virtual void setVal(Double_t value)
Set value of variable to 'value'.
Definition: RooRealVar.cxx:226
RooAbsReal.h
RooAcceptReject::generateEvent
const RooArgSet * generateEvent(UInt_t remaining, Double_t &resampleRatio)
Return a pointer to a generated event.
Definition: RooAcceptReject.cxx:178
RooMsgService.h
RooFit.h
RooArgSet::getRealValue
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:474
ClassImp
#define ClassImp(name)
Definition: Rtypes.h:364
RooAbsNumGenerator::_verbose
Bool_t _verbose
Definition: RooAbsNumGenerator.h:77
r
ROOT::R::TRInterface & r
Definition: Object.C:4
RooAbsData::fill
virtual void fill()
Definition: RooAbsData.cxx:300
Long64_t
long long Long64_t
Definition: RtypesCore.h:73
coutE
#define coutE(a)
Definition: RooMsgService.h:33
coutW
#define coutW(a)
Definition: RooMsgService.h:32
RooAcceptReject::_nextCatVar
TIterator * _nextCatVar
Definition: RooAcceptReject.h:63
RooAcceptReject::registerSampler
static void registerSampler(RooNumGenFactory &fact)
Register RooIntegrator1D, is parameters and capabilities with RooNumIntFactory.
Definition: RooAcceptReject.cxx:61
RooAbsReal::getVal
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:91
RooNumGenFactory.h
RooAbsNumGenerator::_cache
RooDataSet * _cache
Definition: RooAbsNumGenerator.h:80
RooAcceptReject::_eventsUsed
UInt_t _eventsUsed
Definition: RooAcceptReject.h:62
TNamed::fName
TString fName
Definition: TNamed.h:38
coutI
#define coutI(a)
Definition: RooMsgService.h:30
TClass.h
RooAbsReal
Definition: RooAbsReal.h:61
RooAbsNumGenerator::_catVars
RooArgSet _catVars
Definition: RooAbsNumGenerator.h:76
RooNumGenFactory::storeProtoSampler
Bool_t storeProtoSampler(RooAbsNumGenerator *proto, const RooArgSet &defConfig)
Method accepting registration of a prototype numeric integrator along with a RooArgSet of its default...
Definition: RooNumGenFactory.cxx:120
RooDataSet.h
TString.h
RooNumGenFactory
Definition: RooNumGenFactory.h:30
TFoam.h
RooAbsCategory::numTypes
Int_t numTypes(const char *=0) const
Return number of types defined (in range named rangeName if rangeName!=0)
Definition: RooAbsCategory.h:133
bool
TIterator
Definition: TIterator.h:30
RooDataSet::get
virtual const RooArgSet * get(Int_t index) const override
Return RooArgSet with coordinates of event 'index'.
Definition: RooDataSet.cxx:1038
RooAcceptReject::_totalEvents
UInt_t _totalEvents
Definition: RooAcceptReject.h:61
RooAbsNumGenerator::_funcClone
RooAbsReal * _funcClone
Definition: RooAbsNumGenerator.h:74
RooAbsCategory
Definition: RooAbsCategory.h:38
RooAcceptReject::_maxFuncVal
Double_t _maxFuncVal
Definition: RooAcceptReject.h:58
RooAbsNumGenerator::_funcValPtr
RooRealVar * _funcValPtr
Definition: RooAbsNumGenerator.h:78
RooAcceptReject::getFuncMax
Double_t getFuncMax()
Definition: RooAcceptReject.cxx:329
ccoutI
#define ccoutI(a)
Definition: RooMsgService.h:38
RooPrintable::kName
@ kName
Definition: RooPrintable.h:33
RooAcceptReject.h
TIterator.h
RooAcceptReject::~RooAcceptReject
virtual ~RooAcceptReject()
Destructor.
Definition: RooAcceptReject.cxx:165
RooFit::Generation
@ Generation
Definition: RooGlobalFunc.h:67
RooRandom.h
RooAbsCollection::createIterator
TIterator * createIterator(Bool_t dir=kIterForward) const
TIterator-style iteration over contained elements.
Definition: RooAbsCollection.h:118
RooAbsNumGenerator
Definition: RooAbsNumGenerator.h:30
RooAcceptReject
Definition: RooAcceptReject.h:29
RooAbsData::numEntries
virtual Int_t numEntries() const
Definition: RooAbsData.cxx:307
RooCategory.h
RooRealVar.h
RooNumGenConfig::getConfigSection
const RooArgSet & getConfigSection(const char *name) const
Retrieve configuration information specific to integrator with given name.
Definition: RooNumGenConfig.cxx:295
TIterator::Next
virtual TObject * Next()=0
unsigned int
RooAcceptReject::_funcSum
Double_t _funcSum
Definition: RooAcceptReject.h:58
TIterator::Reset
virtual void Reset()=0
proto
const char * proto
Definition: civetweb.c:16604
Double_t
double Double_t
Definition: RtypesCore.h:59
RooNumGenConfig
Definition: RooNumGenConfig.h:25
RooAbsNumGenerator::_realVars
RooArgSet _realVars
Definition: RooAbsNumGenerator.h:76
RooCategory
Definition: RooCategory.h:27
RooAbsNumGenerator::_funcMaxVal
const RooAbsReal * _funcMaxVal
Definition: RooAbsNumGenerator.h:75
RooAcceptReject::_minTrialsArray
UInt_t _minTrialsArray[4]
Definition: RooAcceptReject.h:66
RooNumGenConfig.h
RooAbsRealLValue::randomize
virtual void randomize(const char *rangeName=0)
Set a new value sampled from a uniform distribution over the fit range.
Definition: RooAbsRealLValue.cxx:413
RooPrintable::printStream
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,...
Definition: RooPrintable.cxx:75
genreflex::verbose
bool verbose
Definition: rootcling_impl.cxx:133
RooAcceptReject::_minTrials
UInt_t _minTrials
Definition: RooAcceptReject.h:60
RooPrintable::kSingleLine
@ kSingleLine
Definition: RooPrintable.h:34
RooErrorHandler.h
RooRealBinding.h
RooAcceptReject::RooAcceptReject
RooAcceptReject()
Definition: RooAcceptReject.h:31
RooAcceptReject::_realSampleDim
UInt_t _realSampleDim
Definition: RooAcceptReject.h:59
RooRealVar
Definition: RooRealVar.h:35
Riostream.h
RooRandom::uniform
static Double_t uniform(TRandom *generator=randomGenerator())
Return a number uniformly distributed from (0,1)
Definition: RooRandom.cxx:83
TObject::ClassName
virtual const char * ClassName() const
Returns name of class to which the object belongs.
Definition: TObject.cxx:130
RooAcceptReject::addEventToCache
void addEventToCache()
Add a trial event to our cache and update our estimates of the function maximum value and integral.
Definition: RooAcceptReject.cxx:297
RooAcceptReject::_catSampleMult
UInt_t _catSampleMult
Definition: RooAcceptReject.h:59
RooAbsCollection::getSize
Int_t getSize() const
Definition: RooAbsCollection.h:171
RooAcceptReject::nextAcceptedEvent
const RooArgSet * nextAcceptedEvent()
Scan through events in the cache which have not been used yet, looking for the first accepted one whi...
Definition: RooAcceptReject.cxx:268
RooArgSet
Definition: RooArgSet.h:28
cxcoutD
#define cxcoutD(a)
Definition: RooMsgService.h:81
int
RooAbsData::reset
virtual void reset()
Definition: RooAbsData.cxx:314
RooAcceptReject::_nextRealVar
TIterator * _nextRealVar
Definition: RooAcceptReject.h:64