Logo ROOT   6.18/05
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
22RooIntegrator1D implements an adaptive one-dimensional
23numerical 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"
36#include "RooNumIntConfig.h"
37#include "RooNumIntFactory.h"
38#include "RooMsgService.h"
39
40#include <assert.h>
41
42
43
44using 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) :
92 RooAbsIntegrator(function), _rule(rule), _maxSteps(maxSteps), _minStepsZero(999), _fixSteps(0), _epsAbs(eps), _epsRel(eps), _doExtrap(kTRUE)
93{
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),
114 _doExtrap(kTRUE)
115{
117 _xmin= xmin;
118 _xmax= xmax;
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
129 RooAbsIntegrator(function,config.printEvalCounter()),
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 ;
144 }
145
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) :
158 RooAbsIntegrator(function,config.printEvalCounter()),
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;
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
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{
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{
259 assert(0 != integrand() && integrand()->isValid());
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
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}
#define e(i)
Definition: RSha256.hxx:103
#define oocoutW(o, a)
Definition: RooMsgService.h:46
#define oocoutE(o, a)
Definition: RooMsgService.h:47
#define ooccoutW(o, a)
Definition: RooMsgService.h:53
int Int_t
Definition: RtypesCore.h:41
unsigned int UInt_t
Definition: RtypesCore.h:42
const Bool_t kFALSE
Definition: RtypesCore.h:88
bool Bool_t
Definition: RtypesCore.h:59
double Double_t
Definition: RtypesCore.h:55
const Bool_t kTRUE
Definition: RtypesCore.h:87
#define ClassImp(name)
Definition: Rtypes.h:365
float xmin
Definition: THbookFile.cxx:93
float xmax
Definition: THbookFile.cxx:93
const char * proto
Definition: civetweb.c:16604
Abstract interface for evaluating a real-valued function of one real variable and performing numerica...
Definition: RooAbsFunc.h:23
virtual Double_t getMinLimit(UInt_t dimension) const =0
virtual Double_t getMaxLimit(UInt_t dimension) const =0
UInt_t getDimension() const
Definition: RooAbsFunc.h:29
virtual const char * getName() const
Definition: RooAbsFunc.h:59
RooAbsIntegrator is the abstract interface for integrators of real-valued functions that implement th...
const RooAbsFunc * _function
const RooAbsFunc * integrand() const
Bool_t isValid() const
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
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:558
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 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...
Bool_t defineType(const char *label)
Define a state with given name, the lowest available positive integer is assigned as index.
RooIntegrator1D implements an adaptive one-dimensional numerical integration algorithm.
Double_t addTrapezoids(Int_t n)
Calculate the n-th stage of refinement of the extended trapezoidal summation rule.
Double_t _range
Upper integration bound.
virtual Bool_t checkLimits() const
Check that our integration range is finite and otherwise return kFALSE.
virtual RooAbsIntegrator * clone(const RooAbsFunc &function, const RooNumIntConfig &config) const
Clone integrator with new function binding and configuration. Needed by RooNumIntFactory.
void extrapolate(Int_t n)
Extrapolate result to final value.
Double_t * xvec(Double_t &xx)
Integrator workspace.
Bool_t setLimits(Double_t *xmin, Double_t *xmax)
Change our integration limits.
SummationRule _rule
Double_t _extrapValue
Size of integration range.
virtual ~RooIntegrator1D()
Destructor.
Bool_t initialize()
Initialize the integrator.
virtual Double_t integral(const Double_t *yvec=0)
Calculate numeric integral at given set of function binding parameters.
Double_t addMidpoints(Int_t n)
Calculate the n-th stage of refinement of the Second Euler-Maclaurin summation rule which has the use...
Double_t * _d
Integrator workspace.
Bool_t _useIntegrandLimits
static void registerIntegrator(RooNumIntFactory &fact)
Register RooIntegrator1D, is parameters and capabilities with RooNumIntFactory.
Double_t * _h
Error on extrapolated value.
Double_t _savedResult
Integrator workspace.
Double_t _extrapError
Extrapolated value.
RooIntegrator1D()
coverity[UNINIT_CTOR] Default constructor
Double_t * _c
Integrator workspace.
Double_t * _s
Integrator workspace.
Double_t _xmax
Lower integration bound.
RooNumIntConfig 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.
RooCategory & method1D()
static RooNumIntConfig & defaultConfig()
Return reference to instance of default numeric integrator configuration object.
RooNumIntFactory is a factory to instantiate numeric integrators from a given function binding and a ...
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...
static Int_t isInfinite(Double_t x)
Return true if x is infinite by RooNumBer internal specification.
Definition: RooNumber.cxx:58
RooRealVar represents a fundamental (non-derived) real valued object.
Definition: RooRealVar.h:36
Mother of all ROOT objects.
Definition: TObject.h:37
virtual const char * GetName() const
Returns name of object.
Definition: TObject.cxx:357
Double_t x[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
void function(const Char_t *name_, T fun, const Char_t *docstring=0)
Definition: RExports.h:151
@ Integration
Definition: RooGlobalFunc.h:57
static constexpr double ns
auto * m
Definition: textangle.C:8
static long int sum(long int i)
Definition: Factory.cxx:2258