Go to the documentation of this file.00001 #include "PhysicsTools/Utilities/interface/NumericalIntegration.h"
00002 #include "FWCore/Utilities/interface/EDMException.h"
00003 #include <cmath>
00004
00005 funct::GaussLegendreIntegrator::GaussLegendreIntegrator(unsigned int samples, double epsilon) :
00006 samples_(samples) {
00007 if (samples <= 0)
00008 throw edm::Exception(edm::errors::Configuration)
00009 << "gauss_legendre_integral: number of samples must be positive\n";
00010 if(epsilon <= 0)
00011 throw edm::Exception(edm::errors::Configuration)
00012 << "gauss_legendre_integral: numerical precision must be positive\n";
00013
00014 x.resize(samples);
00015 w.resize(samples);
00016 const unsigned int m = (samples + 1)/2;
00017
00018 double z, zSqr, pp, p1, p2, p3;
00019
00020 for (unsigned int i = 0; i < m; ++i) {
00021 z = std::cos(3.14159265358979323846 * (i + 0.75)/(samples + 0.5));
00022 zSqr = z*z;
00023 do {
00024 p1 = 1.0;
00025 p2 = 0.0;
00026 for (unsigned int j = 0; j < samples; ++j) {
00027 p3 = p2;
00028 p2 = p1;
00029 p1 = ((2.0*j + 1.0)*z*p2 - j*p3)/(j + 1.0);
00030 }
00031 pp = samples*(z*p1 - p2)/(zSqr - 1.0);
00032 z -= p1/pp;
00033 } while (std::fabs(p1/pp) > epsilon);
00034
00035 x[i] = -z;
00036 x[samples - i - 1] = z;
00037 w[i] = 2.0/((1.0 - zSqr)*pp*pp);
00038 w[samples - i -1] = w[i];
00039 }
00040 }
00041
00042
00043 const double funct::GaussIntegrator::x[12] = {
00044 0.96028985649753623, 0.79666647741362674,
00045 0.52553240991632899, 0.18343464249564980,
00046 0.98940093499164993, 0.94457502307323258,
00047 0.86563120238783174, 0.75540440835500303,
00048 0.61787624440264375, 0.45801677765722739,
00049 0.28160355077925891, 0.09501250983763744 };
00050
00051 const double funct::GaussIntegrator::w[12] = {
00052 0.10122853629037626, 0.22238103445337447,
00053 0.31370664587788729, 0.36268378337836198,
00054 0.02715245941175409, 0.06225352393864789,
00055 0.09515851168249278, 0.12462897125553387,
00056 0.14959598881657673, 0.16915651939500254,
00057 0.18260341504492359, 0.18945061045506850 };
00058
00059 const double funct::GaussIntegrator::kHF = 0.5;
00060 const double funct::GaussIntegrator::kCST = 5./1000;
1.7.3