NumericalIntegration.cc

Go to the documentation of this file.

#include "PhysicsTools/Utilities/interface/NumericalIntegration.h"
#include "FWCore/Utilities/interface/EDMException.h"
#include <cmath>

funct::GaussLegendreIntegrator::GaussLegendreIntegrator(unsigned int samples, double epsilon) : 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];
  }
}

const double funct::GaussIntegrator::x[12] = {0.96028985649753623,
                                              0.79666647741362674,
                                              0.52553240991632899,
                                              0.18343464249564980,
                                              0.98940093499164993,
                                              0.94457502307323258,
                                              0.86563120238783174,
                                              0.75540440835500303,
                                              0.61787624440264375,
                                              0.45801677765722739,
                                              0.28160355077925891,
                                              0.09501250983763744};

const double funct::GaussIntegrator::w[12] = {0.10122853629037626,
                                              0.22238103445337447,
                                              0.31370664587788729,
                                              0.36268378337836198,
                                              0.02715245941175409,
                                              0.06225352393864789,
                                              0.09515851168249278,
                                              0.12462897125553387,
                                              0.14959598881657673,
                                              0.16915651939500254,
                                              0.18260341504492359,
                                              0.18945061045506850};

const double funct::GaussIntegrator::kHF = 0.5;
const double funct::GaussIntegrator::kCST = 5. / 1000;