NumericalIntegration.h

Go to the documentation of this file.

#ifndef PhysicsTools_Utilities_NumericalIntegration_h
#define PhysicsTools_Utilities_NumericalIntegration_h
/* 
 * Numerical integration utilities
 *
 * \author: Luca Lista
 * Gauss Legendre and Gauss algorithms based on ROOT implementation
 *
 */
#include "Math/Functor.h"
#include "Math/Integrator.h"
#include "Math/AllIntegrationTypes.h"
#include <vector>
#include <cmath>
#include <memory>

namespace funct {

  template<typename F>
  double trapezoid_integral(const F& f, double min, double max, unsigned int samples) { 
      const double l = max - min, delta = l / samples;
      double sum = 0;
      for(unsigned int i = 0; i < samples; ++i) {
    sum += f(min + (i + 0.5) * delta);
      }
      return sum * delta;
    }
  
  class TrapezoidIntegrator {
  public:
    TrapezoidIntegrator() : samples_(0) { }
    explicit TrapezoidIntegrator(unsigned int samples) : samples_(samples) { }
    template<typename F>
    double operator()(const F& f, double min, double max) const { 
      return trapezoid_integral(f, min, max, samples_);
    }
  private:
    unsigned int samples_;
  };

  class GaussLegendreIntegrator {
  public:
    GaussLegendreIntegrator() : samples_(0) { }
    GaussLegendreIntegrator(unsigned int samples, double epsilon);
    template<typename F>
    double operator()(const F& f, double min, double max) const {
      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;      
    }
  private:
    unsigned int samples_;
    std::vector<double> x, w;
    mutable double a0, b0, result;
    mutable unsigned int i;
  };

  class GaussIntegrator {
  public:
    GaussIntegrator() { }
    GaussIntegrator(double epsilon) : epsilon_(epsilon) { }
    template<typename F>
    double operator()(const F& f, double a, double b) const {
      h = 0;
      if (b == a) return h;
      aconst = kCST/std::abs(b - a);
      bb = a;
    CASE1:
      aa = bb;
      bb = b;
    CASE2:
      c1 = kHF*(bb + aa);
      c2 = kHF*(bb - aa);
      s8 = 0;
      for(i = 0; i < 4; ++i) {
    u = c2*x[i];
    xx = c1 + u;
    f1 = f(xx);
    xx = c1-u;
    f2 = f(xx);
    s8 += w[i]*(f1 + f2);
      }
      s16 = 0;
      for(i = 4; i < 12; ++i) {
    u = c2*x[i];
    xx = c1+u;
    f1 = f(xx);
    xx = c1-u;
    f2 = f(xx);
    s16 += w[i]*(f1 + f2);
      }
      s16 = c2*s16;
      if (std::abs(s16 - c2*s8) <= epsilon_*(1. + std::abs(s16))) {
    h += s16;
    if(bb != b) goto CASE1;
      } else {
    bb = c1;
    if(1. + aconst*std::abs(c2) != 1) goto CASE2;
    h = s8;
      }
      
      error_ = std::abs(s16 - c2*s8);
      return h;   
    }    
    double error() const { return error_; }
  private:
    mutable double error_;
    double epsilon_;
    static const double x[12], w[12];
    static const double kHF;
    static const double kCST;
    mutable double h, aconst, bb, aa, c1, c2, u, s8, s16, f1, f2, xx;
    mutable unsigned int i;
  };

  struct RootIntegrator {
    RootIntegrator(ROOT::Math::IntegrationOneDim::Type type = ROOT::Math::IntegrationOneDim::kADAPTIVE, 
           double absTol = 1e-9, double relTol = 1e-6, unsigned int size = 1000, unsigned int rule = 3) :
      type_(type), absTol_(absTol), relTol_(relTol), size_(size), rule_(rule),
      integrator_(new ROOT::Math::Integrator(type, absTol, relTol, size, rule)) { }
    RootIntegrator(const RootIntegrator & o) { 
      type_ = o.type_; 
      absTol_ = o.absTol_; relTol_ = o.relTol_;
      size_ = o.size_; rule_ = o.rule_;
      integrator_.reset(new ROOT::Math::Integrator(type_, absTol_, relTol_, size_, rule_));
    }
    RootIntegrator & operator=(const RootIntegrator & o) {
      type_ = o.type_; 
      absTol_ = o.absTol_; relTol_ = o.relTol_;
      size_ = o.size_; rule_ = o.rule_;
      integrator_.reset(new ROOT::Math::Integrator(type_, absTol_, relTol_, size_, rule_));
      return * this;
    }
    template<typename F>
    double operator()(const F& f, double a, double b) const {
      ROOT::Math::Functor1D wrapper(&f, &F::operator());
      integrator_->SetFunction(wrapper);
      return integrator_->Integral(a, b);
    }
  private:
    ROOT::Math::IntegrationOneDim::Type type_;
    double absTol_, relTol_;
    unsigned int size_, rule_;
    mutable std::unique_ptr<ROOT::Math::Integrator> integrator_;
  };    
}

#endif