Logo ROOT   6.10/09
Reference Guide
RooIntegrator1D.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 RooIntegrator1D.cxx
19 \class RooIntegrator1D
20 \ingroup Roofitcore
21 
22 RooIntegrator1D implements an adaptive one-dimensional
23 numerical integration algorithm.
24 **/
25 
26 
27 #include "RooFit.h"
28 #include "Riostream.h"
29 
30 #include "TClass.h"
31 #include "RooIntegrator1D.h"
32 #include "RooArgSet.h"
33 #include "RooRealVar.h"
34 #include "RooNumber.h"
35 #include "RooIntegratorBinding.h"
36 #include "RooNumIntConfig.h"
37 #include "RooNumIntFactory.h"
38 #include "RooMsgService.h"
39 
40 #include <assert.h>
41 
42 
43 
44 using namespace std;
45 
47 ;
48 
49 // Register this class with RooNumIntConfig
50 
51 ////////////////////////////////////////////////////////////////////////////////
52 /// Register RooIntegrator1D, is parameters and capabilities with RooNumIntFactory
53 
55 {
56  RooCategory sumRule("sumRule","Summation Rule") ;
57  sumRule.defineType("Trapezoid",RooIntegrator1D::Trapezoid) ;
58  sumRule.defineType("Midpoint",RooIntegrator1D::Midpoint) ;
59  sumRule.setLabel("Trapezoid") ;
60  RooCategory extrap("extrapolation","Extrapolation procedure") ;
61  extrap.defineType("None",0) ;
62  extrap.defineType("Wynn-Epsilon",1) ;
63  extrap.setLabel("Wynn-Epsilon") ;
64  RooRealVar maxSteps("maxSteps","Maximum number of steps",20) ;
65  RooRealVar minSteps("minSteps","Minimum number of steps",999) ;
66  RooRealVar fixSteps("fixSteps","Fixed number of steps",0) ;
67 
69  fact.storeProtoIntegrator(proto,RooArgSet(sumRule,extrap,maxSteps,minSteps,fixSteps)) ;
71 }
72 
73 
74 
75 ////////////////////////////////////////////////////////////////////////////////
76 /// coverity[UNINIT_CTOR]
77 /// Default constructor
78 
80  _h(0), _s(0), _c(0), _d(0), _x(0)
81 {
82 }
83 
84 
85 ////////////////////////////////////////////////////////////////////////////////
86 /// Construct integrator on given function binding, using specified summation
87 /// rule, maximum number of steps and conversion tolerance. The integration
88 /// limits are taken from the function binding
89 
91  Int_t maxSteps, Double_t eps) :
93 {
95  _valid= initialize();
96 }
97 
98 
99 ////////////////////////////////////////////////////////////////////////////////
100 /// Construct integrator on given function binding for given range,
101 /// using specified summation rule, maximum number of steps and
102 /// conversion tolerance. The integration limits are taken from the
103 /// function binding
104 
106  SummationRule rule, Int_t maxSteps, Double_t eps) :
108  _rule(rule),
109  _maxSteps(maxSteps),
110  _minStepsZero(999),
111  _fixSteps(0),
112  _epsAbs(eps),
113  _epsRel(eps),
115 {
117  _xmin= xmin;
118  _xmax= xmax;
119  _valid= initialize();
120 }
121 
122 
123 ////////////////////////////////////////////////////////////////////////////////
124 /// Construct integrator on given function binding, using specified
125 /// configuration object. The integration limits are taken from the
126 /// function binding
127 
130  _epsAbs(config.epsAbs()),
131  _epsRel(config.epsRel())
132 {
133  // Extract parameters from config object
134  const RooArgSet& configSet = config.getConfigSection(IsA()->GetName()) ;
135  _rule = (SummationRule) configSet.getCatIndex("sumRule",Trapezoid) ;
136  _maxSteps = (Int_t) configSet.getRealValue("maxSteps",20) ;
137  _minStepsZero = (Int_t) configSet.getRealValue("minSteps",999) ;
138  _fixSteps = (Int_t) configSet.getRealValue("fixSteps",0) ;
139  _doExtrap = (Bool_t) configSet.getCatIndex("extrapolation",1) ;
140 
141  if (_fixSteps>_maxSteps) {
142  oocoutE((TObject*)0,Integration) << "RooIntegrator1D::ctor() ERROR: fixSteps>maxSteps, fixSteps set to maxSteps" << endl ;
143  _fixSteps = _maxSteps ;
144  }
145 
147  _valid= initialize();
148 }
149 
150 
151 
152 ////////////////////////////////////////////////////////////////////////////////
153 /// Construct integrator on given function binding, using specified
154 /// configuration object and integration range
155 
157  const RooNumIntConfig& config) :
159  _epsAbs(config.epsAbs()),
160  _epsRel(config.epsRel())
161 {
162  // Extract parameters from config object
163  const RooArgSet& configSet = config.getConfigSection(IsA()->GetName()) ;
164  _rule = (SummationRule) configSet.getCatIndex("sumRule",Trapezoid) ;
165  _maxSteps = (Int_t) configSet.getRealValue("maxSteps",20) ;
166  _minStepsZero = (Int_t) configSet.getRealValue("minSteps",999) ;
167  _fixSteps = (Int_t) configSet.getRealValue("fixSteps",0) ;
168  _doExtrap = (Bool_t) configSet.getCatIndex("extrapolation",1) ;
169 
171  _xmin= xmin;
172  _xmax= xmax;
173  _valid= initialize();
174 }
175 
176 
177 
178 ////////////////////////////////////////////////////////////////////////////////
179 /// Clone integrator with new function binding and configuration. Needed by RooNumIntFactory
180 
182 {
183  return new RooIntegrator1D(function,config) ;
184 }
185 
186 
187 
188 ////////////////////////////////////////////////////////////////////////////////
189 /// Initialize the integrator
190 
192 {
193  // apply defaults if necessary
194  if(_maxSteps <= 0) {
195  _maxSteps= (_rule == Trapezoid) ? 20 : 14;
196  }
197 
198  if(_epsRel <= 0) _epsRel= 1e-6;
199  if(_epsAbs <= 0) _epsAbs= 1e-6;
200 
201  // check that the integrand is a valid function
202  if(!isValid()) {
203  oocoutE((TObject*)0,Integration) << "RooIntegrator1D::initialize: cannot integrate invalid function" << endl;
204  return kFALSE;
205  }
206 
207  // Allocate coordinate buffer size after number of function dimensions
208  _x = new Double_t[_function->getDimension()] ;
209 
210 
211  // Allocate workspace for numerical integration engine
212  _h= new Double_t[_maxSteps + 2];
213  _s= new Double_t[_maxSteps + 2];
214  _c= new Double_t[_nPoints + 1];
215  _d= new Double_t[_nPoints + 1];
216 
217  return checkLimits();
218 }
219 
220 
221 ////////////////////////////////////////////////////////////////////////////////
222 /// Destructor
223 
225 {
226  // Release integrator workspace
227  if(_h) delete[] _h;
228  if(_s) delete[] _s;
229  if(_c) delete[] _c;
230  if(_d) delete[] _d;
231  if(_x) delete[] _x;
232 }
233 
234 
235 ////////////////////////////////////////////////////////////////////////////////
236 /// Change our integration limits. Return kTRUE if the new limits are
237 /// ok, or otherwise kFALSE. Always returns kFALSE and does nothing
238 /// if this object was constructed to always use our integrand's limits.
239 
241 {
242  if(_useIntegrandLimits) {
243  oocoutE((TObject*)0,Integration) << "RooIntegrator1D::setLimits: cannot override integrand's limits" << endl;
244  return kFALSE;
245  }
246  _xmin= *xmin;
247  _xmax= *xmax;
248  return checkLimits();
249 }
250 
251 
252 ////////////////////////////////////////////////////////////////////////////////
253 /// Check that our integration range is finite and otherwise return kFALSE.
254 /// Update the limits from the integrand if requested.
255 
257 {
258  if(_useIntegrandLimits) {
259  assert(0 != integrand() && integrand()->isValid());
260  _xmin= integrand()->getMinLimit(0);
261  _xmax= integrand()->getMaxLimit(0);
262  }
263  _range= _xmax - _xmin;
264  if(_range < 0) {
265  oocoutE((TObject*)0,Integration) << "RooIntegrator1D::checkLimits: bad range with min >= max (_xmin = " << _xmin << " _xmax = " << _xmax << ")" << endl;
266  return kFALSE;
267  }
269 }
270 
271 
272 ////////////////////////////////////////////////////////////////////////////////
273 /// Calculate numeric integral at given set of function binding parameters
274 
276 {
277  assert(isValid());
278 
279  // Copy yvec to xvec if provided
280  if (yvec) {
281  UInt_t i ; for (i=0 ; i<_function->getDimension()-1 ; i++) {
282  _x[i+1] = yvec[i] ;
283  }
284  }
285 
286  Int_t j;
287  _h[1]=1.0;
288  Double_t zeroThresh = _epsAbs/_range ;
289  for(j= 1; j<=_maxSteps; j++) {
290  // refine our estimate using the appropriate summation rule
291  _s[j]= (_rule == Trapezoid) ? addTrapezoids(j) : addMidpoints(j);
292 
293  if (j >= _minStepsZero) {
294  Bool_t allZero(kTRUE) ;
295  Int_t jj ; for (jj=0 ; jj<=j ; jj++) {
296  if (_s[j]>=zeroThresh) {
297  allZero=kFALSE ;
298  }
299  }
300  if (allZero) {
301  //cout << "Roo1DIntegrator(" << this << "): zero convergence at step " << j << ", value = " << 0 << endl ;
302  return 0;
303  }
304  }
305 
306  if (_fixSteps>0) {
307 
308  // Fixed step mode, return result after fixed number of steps
309  if (j==_fixSteps) {
310  //cout << "returning result at fixed step " << j << endl ;
311  return _s[j];
312  }
313 
314  } else if(j >= _nPoints) {
315 
316  // extrapolate the results of recent refinements and check for a stable result
317  if (_doExtrap) {
318  extrapolate(j);
319  } else {
320  _extrapValue = _s[j] ;
321  _extrapError = _s[j]-_s[j-1] ;
322  }
323 
325  return _extrapValue;
326  }
327  if(fabs(_extrapError) <= _epsAbs) {
328  return _extrapValue ;
329  }
330 
331  }
332  // update the step size for the next refinement of the summation
333  _h[j+1]= (_rule == Trapezoid) ? _h[j]/4. : _h[j]/9.;
334  }
335 
336  oocoutW((TObject*)0,Integration) << "RooIntegrator1D::integral: integral of " << _function->getName() << " over range (" << _xmin << "," << _xmax << ") did not converge after "
337  << _maxSteps << " steps" << endl;
338  for(j= 1; j <= _maxSteps; j++) {
339  ooccoutW((TObject*)0,Integration) << " [" << j << "] h = " << _h[j] << " , s = " << _s[j] << endl;
340  }
341 
342  return _s[_maxSteps] ;
343 
344 }
345 
346 
347 ////////////////////////////////////////////////////////////////////////////////
348 /// Calculate the n-th stage of refinement of the Second Euler-Maclaurin
349 /// summation rule which has the useful property of not evaluating the
350 /// integrand at either of its endpoints but requires more function
351 /// evaluations than the trapezoidal rule. This rule can be used with
352 /// a suitable change of variables to estimate improper integrals.
353 
355 {
356  Double_t x,tnm,sum,del,ddel;
357  Int_t it,j;
358 
359  if(n == 1) {
360  Double_t xmid= 0.5*(_xmin + _xmax);
361  return (_savedResult= _range*integrand(xvec(xmid)));
362  }
363  else {
364  for(it=1, j=1; j < n-1; j++) it*= 3;
365  tnm= it;
366  del= _range/(3.*tnm);
367  ddel= del+del;
368  x= _xmin + 0.5*del;
369  for(sum= 0, j= 1; j <= it; j++) {
370  sum+= integrand(xvec(x));
371  x+= ddel;
372  sum+= integrand(xvec(x));
373  x+= del;
374  }
375  return (_savedResult= (_savedResult + _range*sum/tnm)/3.);
376  }
377 }
378 
379 
380 ////////////////////////////////////////////////////////////////////////////////
381 /// Calculate the n-th stage of refinement of the extended trapezoidal
382 /// summation rule. This is the most efficient rule for a well behaved
383 /// integrand that can be evaluated over its entire range, including the
384 /// endpoints.
385 
387 {
388  Double_t x,tnm,sum,del;
389  Int_t it,j;
390 
391  if (n == 1) {
392  // use a single trapezoid to cover the full range
393  return (_savedResult= 0.5*_range*(integrand(xvec(_xmin)) + integrand(xvec(_xmax))));
394  }
395  else {
396  // break the range down into several trapezoids using 2**(n-2)
397  // equally-spaced interior points
398  for(it=1, j=1; j < n-1; j++) it <<= 1;
399  tnm= it;
400  del= _range/tnm;
401  x= _xmin + 0.5*del;
402  for(sum=0.0, j=1; j<=it; j++, x+=del) sum += integrand(xvec(x));
403  return (_savedResult= 0.5*(_savedResult + _range*sum/tnm));
404  }
405 }
406 
407 
408 
409 ////////////////////////////////////////////////////////////////////////////////
410 /// Extrapolate result to final value
411 
413 {
414  Double_t *xa= &_h[n-_nPoints];
415  Double_t *ya= &_s[n-_nPoints];
416  Int_t i,m,ns=1;
417  Double_t den,dif,dift,ho,hp,w;
418 
419  dif=fabs(xa[1]);
420  for (i= 1; i <= _nPoints; i++) {
421  if ((dift=fabs(xa[i])) < dif) {
422  ns=i;
423  dif=dift;
424  }
425  _c[i]=ya[i];
426  _d[i]=ya[i];
427  }
428  _extrapValue= ya[ns--];
429  for(m= 1; m < _nPoints; m++) {
430  for(i= 1; i <= _nPoints-m; i++) {
431  ho=xa[i];
432  hp=xa[i+m];
433  w=_c[i+1]-_d[i];
434  if((den=ho-hp) == 0.0) {
435  oocoutE((TObject*)0,Integration) << "RooIntegrator1D::extrapolate: internal error" << endl;
436  }
437  den=w/den;
438  _d[i]=hp*den;
439  _c[i]=ho*den;
440  }
441  _extrapValue += (_extrapError=(2*ns < (_nPoints-m) ? _c[ns+1] : _d[ns--]));
442  }
443 }
static RooNumIntConfig & defaultConfig()
Return reference to instance of default numeric integrator configuration object.
RooAbsIntegrator is the abstract interface for integrators of real-valued functions that implement th...
Int_t getCatIndex(const char *name, Int_t defVal=0, Bool_t verbose=kFALSE) const
Get index value of a RooAbsCategory stored in set with given name.
Definition: RooArgSet.cxx:613
static long int sum(long int i)
Definition: Factory.cxx:2162
Double_t * _s
Integrator workspace.
float xmin
Definition: THbookFile.cxx:93
RooNumIntConfig holds the configuration parameters of the various numeric integrators used by RooReal...
Double_t _extrapError
Extrapolated value.
RooIntegrator1D()
coverity[UNINIT_CTOR] Default constructor
Double_t * _c
Integrator workspace.
Double_t * xvec(Double_t &xx)
Integrator workspace.
RooNumIntFactory is a factory to instantiate numeric integrators from a given function binding and a ...
Bool_t initialize()
Initialize the integrator.
static void registerIntegrator(RooNumIntFactory &fact)
Register RooIntegrator1D, is parameters and capabilities with RooNumIntFactory.
Double_t _extrapValue
Size of integration range.
Bool_t printEvalCounter() const
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
Double_t addMidpoints(Int_t n)
Calculate the n-th stage of refinement of the Second Euler-Maclaurin summation rule which has the use...
const RooArgSet & getConfigSection(const char *name) const
Retrieve configuration information specific to integrator with given name.
static Int_t isInfinite(Double_t x)
Return true if x is infinite by RooNumBer internal specification.
Definition: RooNumber.cxx:58
STL namespace.
virtual ~RooIntegrator1D()
Destructor.
Bool_t isValid() const
Double_t addTrapezoids(Int_t n)
Calculate the n-th stage of refinement of the extended trapezoidal summation rule.
#define ooccoutW(o, a)
Definition: RooMsgService.h:53
virtual Bool_t setLabel(const char *label, Bool_t printError=kTRUE)
Set value by specifying the name of the desired state If printError is set, a message will be printed...
Double_t _savedResult
Integrator workspace.
RooCategory & method1D()
Double_t x[n]
Definition: legend1.C:17
#define oocoutE(o, a)
Definition: RooMsgService.h:47
void function(const Char_t *name_, T fun, const Char_t *docstring=0)
Definition: RExports.h:146
Double_t _range
Upper integration bound.
virtual const char * getName() const
Definition: RooAbsFunc.h:59
Double_t * _d
Integrator workspace.
RooRealVar represents a fundamental (non-derived) real valued object.
Definition: RooRealVar.h:36
virtual Bool_t checkLimits() const
Check that our integration range is finite and otherwise return kFALSE.
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
virtual RooAbsIntegrator * clone(const RooAbsFunc &function, const RooNumIntConfig &config) const
Clone integrator with new function binding and configuration. Needed by RooNumIntFactory.
virtual Double_t integral(const Double_t *yvec=0)
Calculate numeric integral at given set of function binding parameters.
UInt_t getDimension() const
Definition: RooAbsFunc.h:29
unsigned int UInt_t
Definition: RtypesCore.h:42
TMarker * m
Definition: textangle.C:8
const RooAbsFunc * _function
virtual Double_t getMinLimit(UInt_t dimension) const =0
float xmax
Definition: THbookFile.cxx:93
Bool_t setLimits(Double_t *xmin, Double_t *xmax)
Change our integration limits.
RooCategory represents a fundamental (non-derived) discrete value object.
Definition: RooCategory.h:24
const Bool_t kFALSE
Definition: RtypesCore.h:92
const RooAbsFunc * integrand() const
RooIntegrator1D implements an adaptive one-dimensional numerical integration algorithm.
#define ClassImp(name)
Definition: Rtypes.h:336
virtual Double_t getMaxLimit(UInt_t dimension) const =0
double Double_t
Definition: RtypesCore.h:55
you should not use this method at all Int_t Int_t Double_t Double_t Double_t e
Definition: TRolke.cxx:630
Double_t * _h
Error on extrapolated value.
#define oocoutW(o, a)
Definition: RooMsgService.h:46
Mother of all ROOT objects.
Definition: TObject.h:37
Bool_t _useIntegrandLimits
Double_t _xmax
Lower integration bound.
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:527
const char * proto
Definition: civetweb.c:11652
Bool_t defineType(const char *label)
Define a state with given name, the lowest available positive integer is assigned as index...
void extrapolate(Int_t n)
Extrapolate result to final value.
virtual const char * GetName() const
Returns name of object.
Definition: TObject.cxx:364
SummationRule _rule
const Bool_t kTRUE
Definition: RtypesCore.h:91
Abstract interface for evaluating a real-valued function of one real variable and performing numerica...
Definition: RooAbsFunc.h:23
const Int_t n
Definition: legend1.C:16
Bool_t storeProtoIntegrator(RooAbsIntegrator *proto, const RooArgSet &defConfig, const char *depName="")
Method accepting registration of a prototype numeric integrator along with a RooArgSet of its default...