#include <NumericalIntegration.h>

Public Member Functions
	GaussLegendreIntegrator ()

	GaussLegendreIntegrator (unsigned int samples, double epsilon)

template<typename F >
double	operator() (const F &f, double min, double max) const

Private Attributes
double	a0

double	b0

unsigned int	i

double	result

unsigned int	samples_

std::vector< double >	w

std::vector< double >	x

Detailed Description

Definition at line 42 of file NumericalIntegration.h.

Constructor & Destructor Documentation

◆ GaussLegendreIntegrator() [1/2]

funct::GaussLegendreIntegrator::GaussLegendreIntegrator ( )

inline

Definition at line 44 of file NumericalIntegration.h.

44 : samples_(0) {}

◆ GaussLegendreIntegrator() [2/2]

funct::GaussLegendreIntegrator::GaussLegendreIntegrator	(	unsigned int	samples,
		double	epsilon
	)

Definition at line 5 of file NumericalIntegration.cc.

                                                                                           : samples_(samples) {
   if (samples <= 0)
     throw edm::Exception(edm::errors::Configuration) << "gauss_legendre_integral: number of samples must be positive\n";
   if (epsilon <= 0)
     throw edm::Exception(edm::errors::Configuration)
         << "gauss_legendre_integral: numerical precision must be positive\n";
  
   x.resize(samples);
   w.resize(samples);
   const unsigned int m = (samples + 1) / 2;
  
   double z, zSqr, pp, p1, p2, p3;
  
   for (unsigned int i = 0; i < m; ++i) {
     z = std::cos(3.14159265358979323846 * (i + 0.75) / (samples + 0.5));
     zSqr = z * z;
     do {
       p1 = 1.0;
       p2 = 0.0;
       for (unsigned int j = 0; j < samples; ++j) {
         p3 = p2;
         p2 = p1;
         p1 = ((2.0 * j + 1.0) * z * p2 - j * p3) / (j + 1.0);
       }
       pp = samples * (z * p1 - p2) / (zSqr - 1.0);
       z -= p1 / pp;
     } while (std::fabs(p1 / pp) > epsilon);
  
     x[i] = -z;
     x[samples - i - 1] = z;
     w[i] = 2.0 / ((1.0 - zSqr) * pp * pp);
     w[samples - i - 1] = w[i];
   }
 }

References edm::errors::Configuration, funct::cos(), geometryDiff::epsilon, Exception, i, dqmiolumiharvest::j, funct::m, p1, p2, p3, createTree::pp, EgammaValidation_cff::samples, w, and x.

Member Function Documentation

◆ operator()()

template<typename F >

double funct::GaussLegendreIntegrator::operator()	(	const F &	f,
		double	min,
		double	max
	)		const

inline

Definition at line 47 of file NumericalIntegration.h.

                                                                 {
       a0 = 0.5 * (max + min);
       b0 = 0.5 * (max - min);
       result = 0.0;
       for (i = 0; i < samples_; ++i) {
         result += w[i] * f(a0 + b0 * x[i]);
       }
  
       return result * b0;
     }