Utilities/src/NumericalIntegration.cc

0001 #include "PhysicsTools/Utilities/interface/NumericalIntegration.h"
0002 #include "FWCore/Utilities/interface/EDMException.h"
0003 #include <cmath>
0004
0005 funct::GaussLegendreIntegrator::GaussLegendreIntegrator(unsigned int samples, double epsilon) : samples_(samples) {
0006   if (samples <= 0)
0007     throw edm::Exception(edm::errors::Configuration) << "gauss_legendre_integral: number of samples must be positive\n";
0008   if (epsilon <= 0)
0009     throw edm::Exception(edm::errors::Configuration)
0010         << "gauss_legendre_integral: numerical precision must be positive\n";
0011
0012   x.resize(samples);
0013   w.resize(samples);
0014   const unsigned int m = (samples + 1) / 2;
0015
0016   double z, zSqr, pp, p1, p2, p3;
0017
0018   for (unsigned int i = 0; i < m; ++i) {
0019     z = std::cos(3.14159265358979323846 * (i + 0.75) / (samples + 0.5));
0020     zSqr = z * z;
0021     do {
0022       p1 = 1.0;
0023       p2 = 0.0;
0024       for (unsigned int j = 0; j < samples; ++j) {
0025         p3 = p2;
0026         p2 = p1;
0027         p1 = ((2.0 * j + 1.0) * z * p2 - j * p3) / (j + 1.0);
0028       }
0029       pp = samples * (z * p1 - p2) / (zSqr - 1.0);
0030       z -= p1 / pp;
0031     } while (std::fabs(p1 / pp) > epsilon);
0032
0033     x[i] = -z;
0034     x[samples - i - 1] = z;
0035     w[i] = 2.0 / ((1.0 - zSqr) * pp * pp);
0036     w[samples - i - 1] = w[i];
0037   }
0038 }
0039
0040 const double funct::GaussIntegrator::x[12] = {0.96028985649753623,
0041                                               0.79666647741362674,
0042                                               0.52553240991632899,
0043                                               0.18343464249564980,
0044                                               0.98940093499164993,
0045                                               0.94457502307323258,
0046                                               0.86563120238783174,
0047                                               0.75540440835500303,
0048                                               0.61787624440264375,
0049                                               0.45801677765722739,
0050                                               0.28160355077925891,
0051                                               0.09501250983763744};
0052
0053 const double funct::GaussIntegrator::w[12] = {0.10122853629037626,
0054                                               0.22238103445337447,
0055                                               0.31370664587788729,
0056                                               0.36268378337836198,
0057                                               0.02715245941175409,
0058                                               0.06225352393864789,
0059                                               0.09515851168249278,
0060                                               0.12462897125553387,
0061                                               0.14959598881657673,
0062                                               0.16915651939500254,
0063                                               0.18260341504492359,
0064                                               0.18945061045506850};
0065
0066 const double funct::GaussIntegrator::kHF = 0.5;
0067 const double funct::GaussIntegrator::kCST = 5. / 1000;