Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
RooGaussKronrodIntegrator1D.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 RooGaussKronrodIntegrator1D.cxx
19\class RooGaussKronrodIntegrator1D
20\ingroup Roofitcore
21
22Implements the Gauss-Kronrod integration algorithm.
23
24An Gaussian quadrature method for numerical integration in which
25error is estimation based on evaluation at special points known as
26"Kronrod points." By suitably picking these points, abscissas from
27previous iterations can be reused as part of the new set of points,
28whereas usual Gaussian quadrature would require recomputation of
29all abscissas at each iteration.
30
31This class automatically handles (-inf,+inf) integrals by dividing
32the integration in three regions (-inf,-1), (-1,1), (1,inf) and
33calculating the 1st and 3rd term using a x -> 1/x coordinate
34transformation
35
36This class embeds the Gauss-Kronrod integrator from the GNU
37Scientific Library version 1.5 and applies the 10-, 21-, 43- and
3887-point rule in succession until the required target precision is
39reached
40**/
41
43
44#include <RooArgSet.h>
45#include <RooMsgService.h>
46#include <RooNumIntFactory.h>
47#include <RooNumber.h>
48#include <RooRealVar.h>
49
50#include <Riostream.h>
51
52#include <TMath.h>
53
54#include <gsl/gsl_integration.h>
55
56#include <cassert>
57#include <cfloat>
58#include <cmath>
59
60using std::endl;
61
62/// \cond ROOFIT_INTERNAL
63
64// register integrator class
65// create a derived class in order to call the protected method of the
66// RoodaptiveGaussKronrodIntegrator1D
67namespace RooFit_internal {
68struct Roo_internal_GKInteg1D : public RooGaussKronrodIntegrator1D {
69
70 static void registerIntegrator()
71 {
72 auto &intFactory = RooNumIntFactory::instance();
74 }
75};
76// class used to register integrator at loafing time
77struct Roo_reg_GKInteg1D {
78 Roo_reg_GKInteg1D() { Roo_internal_GKInteg1D::registerIntegrator(); }
79};
80
81static Roo_reg_GKInteg1D instance;
82} // namespace RooFit_internal
83
84/// \endcond
85
86
87////////////////////////////////////////////////////////////////////////////////
88/// Register RooGaussKronrodIntegrator1D, its parameters and capabilities with RooNumIntConfig
89
91{
92 auto creator = [](const RooAbsFunc &function, const RooNumIntConfig &config) {
93 return std::make_unique<RooGaussKronrodIntegrator1D>(function, config);
94 };
95
96 fact.registerPlugin("RooGaussKronrodIntegrator1D", creator, {},
97 /*canIntegrate1D=*/true,
98 /*canIntegrate2D=*/false,
99 /*canIntegrateND=*/false,
100 /*canIntegrateOpenEnded=*/true);
101
102 oocoutI(nullptr, Integration) << "RooGaussKronrodIntegrator1D has been registered" << std::endl;
103}
104
105
106
107////////////////////////////////////////////////////////////////////////////////
108/// Construct integral on 'function' using given configuration object. The integration
109/// range is taken from the definition in the function binding
110
112 : RooAbsIntegrator(function), _useIntegrandLimits(true), _epsAbs(config.epsRel()), _epsRel(config.epsAbs())
113{
114
116}
117
118
119
120////////////////////////////////////////////////////////////////////////////////
121/// Construct integral on 'function' using given configuration object in the given range
122
124 const RooNumIntConfig &config)
125 : RooAbsIntegrator(function),
126 _useIntegrandLimits(false),
127 _epsAbs(config.epsRel()),
128 _epsRel(config.epsAbs()),
129 _xmin(xmin),
130 _xmax(xmax)
131{
133}
134
135
136////////////////////////////////////////////////////////////////////////////////
137/// Perform one-time initialization of integrator
138
140{
141 // Allocate coordinate buffer size after number of function dimensions
142 _x.resize(_function->getDimension());
143
144 return checkLimits();
145}
146
147
148
149////////////////////////////////////////////////////////////////////////////////
150/// Change our integration limits. Return true if the new limits are
151/// ok, or otherwise false. Always returns false and does nothing
152/// if this object was constructed to always use our integrand's limits.
153
155{
157 oocoutE(nullptr,Eval) << "RooGaussKronrodIntegrator1D::setLimits: cannot override integrand's limits" << endl;
158 return false;
159 }
160 _xmin= *xmin;
161 _xmax= *xmax;
162 return checkLimits();
163}
164
165
166
167////////////////////////////////////////////////////////////////////////////////
168/// Check that our integration range is finite and otherwise return false.
169/// Update the limits from the integrand if requested.
170
172{
174 assert(nullptr != integrand() && integrand()->isValid());
177 }
178 return true ;
179}
180
181
182
184{
185 auto instance = reinterpret_cast<RooGaussKronrodIntegrator1D*>(data);
186 return instance->integrand(instance->xvec(x)) ;
187}
188
189
190
191////////////////////////////////////////////////////////////////////////////////
192/// Calculate and return integral
193
195{
196 assert(isValid());
197
198 // Copy yvec to xvec if provided
199 if (yvec) {
200 UInt_t i ; for (i=0 ; i<_function->getDimension()-1 ; i++) {
201 _x[i+1] = yvec[i] ;
202 }
203 }
204
205 // Setup glue function
208 F.params = this ;
209
210 // Return values
211 double result;
212 double error;
213 size_t neval = 0 ;
214
215 // Call GSL implementation of integeator
216 gsl_integration_qng (&F, _xmin, _xmax, _epsAbs, _epsRel, &result, &error, &neval);
217
218 return result;
219}
static Roo_reg_AGKInteg1D instance
double RooGaussKronrodIntegrator1D_GSL_GlueFunction(double x, void *data)
#define oocoutE(o, a)
#define oocoutI(o, a)
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t result
float xmin
float xmax
Abstract interface for evaluating a real-valued function of one real variable and performing numerica...
Definition RooAbsFunc.h:27
virtual double getMaxLimit(UInt_t dimension) const =0
virtual double getMinLimit(UInt_t dimension) const =0
UInt_t getDimension() const
Definition RooAbsFunc.h:33
Abstract interface for integrators of real-valued functions that implement the RooAbsFunc interface.
bool isValid() const
Is integrator in valid state.
double integrand(const double x[]) const
Return value of integrand at given observable values.
const RooAbsFunc * _function
Pointer to function binding of integrand.
const RooAbsFunc * integrand() const
Return integrand function binding.
bool _valid
Is integrator in valid state?
Implements the Gauss-Kronrod integration algorithm.
double integral(const double *yvec=nullptr) override
Calculate and return integral.
double _xmax
Lower integration bound.
RooGaussKronrodIntegrator1D(const RooAbsFunc &function, const RooNumIntConfig &config)
Construct integral on 'function' using given configuration object.
friend double RooGaussKronrodIntegrator1D_GSL_GlueFunction(double x, void *data)
bool initialize()
Perform one-time initialization of integrator.
bool checkLimits() const override
Check that our integration range is finite and otherwise return false.
bool setLimits(double *xmin, double *xmax) override
Change our integration limits.
static void registerIntegrator(RooNumIntFactory &fact)
Register RooGaussKronrodIntegrator1D, its parameters and capabilities with RooNumIntConfig.
Holds the configuration parameters of the various numeric integrators used by RooRealIntegral.
Factory to instantiate numeric integrators from a given function binding and a given configuration.
static RooNumIntFactory & instance()
Static method returning reference to singleton instance of factory.
bool registerPlugin(std::string const &name, Creator const &creator, const RooArgSet &defConfig, bool canIntegrate1D, bool canIntegrate2D, bool canIntegrateND, bool canIntegrateOpenEnded, const char *depName="")
Method accepting registration of a prototype numeric integrator along with a RooArgSet of its default...
Double_t x[n]
Definition legend1.C:17
#define F(x, y, z)
double(* function)(double x, void *params)