#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 41 of file NumericalIntegration.h.

Constructor & Destructor Documentation

funct::GaussLegendreIntegrator::GaussLegendreIntegrator ( ) [inline]

Definition at line 43 of file NumericalIntegration.h.

: samples_(0) { }

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

Definition at line 5 of file NumericalIntegration.cc.

References edm::errors::Configuration, funct::cos(), Exception, i, j, funct::m, p1, p2, p3, createTree::pp, w, x, and z.

                                                                                          :
  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];
  }
}

Member Function Documentation

template<typename F >

double funct::GaussLegendreIntegrator::operator()	(	const F &	f,
		double	min,
		double	max
	)		const `[inline]`

Definition at line 46 of file NumericalIntegration.h.

References a0, b0, f, i, min, result, samples_, w, and x.

                                                                {
      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;      
    }