Logo ROOT   6.14/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 
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 "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 
53 using 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 
81 RooAcceptReject::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 
179 const 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) {
192  addEventToCache();
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--) {
229  addEventToCache();
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) {
252  addEventToCache() ;
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?
276  if(r*_maxFuncVal > _funcValPtr->getVal()) {
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
301  _nextCatVar->Reset();
302  RooCategory *cat = 0;
303  while((cat= (RooCategory*)_nextCatVar->Next())) cat->randomize();
304 
305  // randomize each real argument
306  _nextRealVar->Reset();
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) {
338  addEventToCache();
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 
TIterator * createIterator(Bool_t dir=kIterForward) const
virtual ~RooAcceptReject()
Destructor.
#define coutE(a)
Definition: RooMsgService.h:34
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, which is interpreted as an OR of &#39;enum ContentsOptions&#39; values and in the style given by &#39;enum StyleOption&#39;.
const RooArgSet * nextAcceptedEvent()
Scan through events in the cache which have not been used yet, looking for the first accepted one whi...
long long Long64_t
Definition: RtypesCore.h:69
virtual void Reset()=0
virtual void randomize(const char *rangeName=0)
Set a new value sampled from a uniform distribution over the fit range.
#define coutI(a)
Definition: RooMsgService.h:31
Double_t getVal(const RooArgSet *set=0) const
Definition: RooAbsReal.h:64
#define cxcoutD(a)
Definition: RooMsgService.h:79
int Int_t
Definition: RtypesCore.h:41
const RooAbsReal * _funcMaxVal
bool Bool_t
Definition: RtypesCore.h:59
STL namespace.
#define coutW(a)
Definition: RooMsgService.h:33
void addEventToCache()
Add a trial event to our cache and update our estimates of the function maximum value and integral...
#define ccoutI(a)
Definition: RooMsgService.h:38
Iterator abstract base class.
Definition: TIterator.h:30
virtual void randomize(const char *rangeName=0)
Randomize current value.
virtual const char * ClassName() const
Returns name of class to which the object belongs.
Definition: TObject.cxx:128
Double_t _maxFuncVal
virtual void reset()
Definition: RooAbsData.cxx:292
Int_t numTypes(const char *=0) const
RooRealVar represents a fundamental (non-derived) real valued object.
Definition: RooRealVar.h:36
virtual void setVal(Double_t value)
Set value of variable to &#39;value&#39;.
Definition: RooRealVar.cxx:204
Int_t getSize() const
const RooArgSet * generateEvent(UInt_t remaining, Double_t &resampleRatio)
Return a pointer to a generated event.
const RooArgSet & getConfigSection(const char *name) const
Retrieve configuration information specific to integrator with given name.
ROOT::R::TRInterface & r
Definition: Object.C:4
Class RooAcceptReject is a generic toy monte carlo generator implement the accept/reject sampling tec...
RooRealVar * _funcValPtr
unsigned int UInt_t
Definition: RtypesCore.h:42
Class RooAbsNumGenerator is the abstract base class for MC event generator implementations like RooAc...
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 fill()
Definition: RooAbsData.cxx:278
TString fName
Definition: TNamed.h:32
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 &#39;index&#39;.
UInt_t _minTrialsArray[4]
#define ClassImp(name)
Definition: Rtypes.h:359
static Double_t uniform(TRandom *generator=randomGenerator())
Return a number uniformly distributed from (0,1)
Definition: RooRandom.cxx:84
double Double_t
Definition: RtypesCore.h:55
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:53
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:526
virtual TObject * Next()=0
RooAbsCategory is the common abstract base class for objects that represent a discrete value with a f...
const char * proto
Definition: civetweb.c:15049
TIterator * _nextRealVar
TIterator * _nextCatVar
RooNumGenConfig holds the configuration parameters of the various numeric integrators used by RooReal...
static void registerSampler(RooNumGenFactory &fact)
Register RooIntegrator1D, is parameters and capabilities with RooNumIntFactory.
RooNumGenFactory is a factory to instantiate numeric integrators from a given function binding and a ...
virtual Int_t numEntries() const
Definition: RooAbsData.cxx:285