Utilities/interface/NumericalIntegration.h

0001 #ifndef PhysicsTools_Utilities_NumericalIntegration_h
0002 #define PhysicsTools_Utilities_NumericalIntegration_h
0003 /*
0004  * Numerical integration utilities
0005  *
0006  * \author: Luca Lista
0007  * Gauss Legendre and Gauss algorithms based on ROOT implementation
0008  *
0009  */
0010 #include "Math/AllIntegrationTypes.h"
0011 #include "Math/Functor.h"
0012 #include "Math/Integrator.h"
0013 #include <cmath>
0014 #include <memory>
0015 #include <vector>
0016
0017 namespace funct {
0018
0019   template <typename F>
0020   double trapezoid_integral(const F& f, double min, double max, unsigned int samples) {
0021     const double l = max - min, delta = l / samples;
0022     double sum = 0;
0023     for (unsigned int i = 0; i < samples; ++i) {
0024       sum += f(min + (i + 0.5) * delta);
0025     }
0026     return sum * delta;
0027   }
0028
0029   class TrapezoidIntegrator {
0030   public:
0031     TrapezoidIntegrator() : samples_(0) {}
0032     explicit TrapezoidIntegrator(unsigned int samples) : samples_(samples) {}
0033     template <typename F>
0034     double operator()(const F& f, double min, double max) const {
0035       return trapezoid_integral(f, min, max, samples_);
0036     }
0037
0038   private:
0039     unsigned int samples_;
0040   };
0041
0042   class GaussLegendreIntegrator {
0043   public:
0044     GaussLegendreIntegrator() : samples_(0) {}
0045     GaussLegendreIntegrator(unsigned int samples, double epsilon);
0046     template <typename F>
0047     double operator()(const F& f, double min, double max) const {
0048       a0 = 0.5 * (max + min);
0049       b0 = 0.5 * (max - min);
0050       result = 0.0;
0051       for (i = 0; i < samples_; ++i) {
0052         result += w[i] * f(a0 + b0 * x[i]);
0053       }
0054
0055       return result * b0;
0056     }
0057
0058   private:
0059     unsigned int samples_;
0060     std::vector<double> x, w;
0061     mutable double a0, b0, result;
0062     mutable unsigned int i;
0063   };
0064
0065   class GaussIntegrator {
0066   public:
0067     GaussIntegrator() {}
0068     GaussIntegrator(double epsilon) : epsilon_(epsilon) {}
0069     template <typename F>
0070     double operator()(const F& f, double a, double b) const {
0071       h = 0;
0072       if (b == a)
0073         return h;
0074       aconst = kCST / std::abs(b - a);
0075       bb = a;
0076     CASE1:
0077       aa = bb;
0078       bb = b;
0079     CASE2:
0080       c1 = kHF * (bb + aa);
0081       c2 = kHF * (bb - aa);
0082       s8 = 0;
0083       for (i = 0; i < 4; ++i) {
0084         u = c2 * x[i];
0085         xx = c1 + u;
0086         f1 = f(xx);
0087         xx = c1 - u;
0088         f2 = f(xx);
0089         s8 += w[i] * (f1 + f2);
0090       }
0091       s16 = 0;
0092       for (i = 4; i < 12; ++i) {
0093         u = c2 * x[i];
0094         xx = c1 + u;
0095         f1 = f(xx);
0096         xx = c1 - u;
0097         f2 = f(xx);
0098         s16 += w[i] * (f1 + f2);
0099       }
0100       s16 = c2 * s16;
0101       if (std::abs(s16 - c2 * s8) <= epsilon_ * (1. + std::abs(s16))) {
0102         h += s16;
0103         if (bb != b)
0104           goto CASE1;
0105       } else {
0106         bb = c1;
0107         if (1. + aconst * std::abs(c2) != 1)
0108           goto CASE2;
0109         h = s8;
0110       }
0111
0112       error_ = std::abs(s16 - c2 * s8);
0113       return h;
0114     }
0115     double error() const { return error_; }
0116
0117   private:
0118     mutable double error_;
0119     double epsilon_;
0120     static const double x[12], w[12];
0121     static const double kHF;
0122     static const double kCST;
0123     mutable double h, aconst, bb, aa, c1, c2, u, s8, s16, f1, f2, xx;
0124     mutable unsigned int i;
0125   };
0126
0127   struct RootIntegrator {
0128     RootIntegrator(ROOT::Math::IntegrationOneDim::Type type = ROOT::Math::IntegrationOneDim::kADAPTIVE,
0129                    double absTol = 1e-9,
0130                    double relTol = 1e-6,
0131                    unsigned int size = 1000,
0132                    unsigned int rule = 3)
0133         : type_(type),
0134           absTol_(absTol),
0135           relTol_(relTol),
0136           size_(size),
0137           rule_(rule),
0138           integrator_(new ROOT::Math::Integrator(type, absTol, relTol, size, rule)) {}
0139     RootIntegrator(const RootIntegrator& o) {
0140       type_ = o.type_;
0141       absTol_ = o.absTol_;
0142       relTol_ = o.relTol_;
0143       size_ = o.size_;
0144       rule_ = o.rule_;
0145       integrator_ = std::make_unique<ROOT::Math::Integrator>(type_, absTol_, relTol_, size_, rule_);
0146     }
0147     RootIntegrator& operator=(const RootIntegrator& o) {
0148       type_ = o.type_;
0149       absTol_ = o.absTol_;
0150       relTol_ = o.relTol_;
0151       size_ = o.size_;
0152       rule_ = o.rule_;
0153       integrator_ = std::make_unique<ROOT::Math::Integrator>(type_, absTol_, relTol_, size_, rule_);
0154       return *this;
0155     }
0156     template <typename F>
0157     double operator()(const F& f, double a, double b) const {
0158       ROOT::Math::Functor1D wrapper(&f, &F::operator());
0159       integrator_->SetFunction(wrapper);
0160       return integrator_->Integral(a, b);
0161     }
0162
0163   private:
0164     ROOT::Math::IntegrationOneDim::Type type_;
0165     double absTol_, relTol_;
0166     unsigned int size_, rule_;
0167     mutable std::unique_ptr<ROOT::Math::Integrator> integrator_;
0168   };
0169 }  // namespace funct
0170
0171 #endif