Integral.h

Go to the documentation of this file.

#ifndef PhysicsTools_Utilities_Integral_h
#define PhysicsTools_Utilities_Integral_h
#include "PhysicsTools/Utilities/interface/Primitive.h"
#include "PhysicsTools/Utilities/interface/NumericalIntegration.h"
 
namespace funct {
 
  struct no_var;
 
  template <typename F, typename X = no_var>
  struct IntegralStruct {
    IntegralStruct(const F& f) : p(primitive<X>(f)) {}
    double operator()(double min, double max) const {
      X::set(min);
      double pMin = p();
      X::set(max);
      double pMax = p();
      return pMax - pMin;
    }
 
  private:
    typename Primitive<F, X>::type p;
  };
 
  template <typename F>
  struct IntegralStruct<F> {
    IntegralStruct(const F& f) : p(primitive(f)) {}
    double operator()(double min, double max) const { return p(max) - p(min); }
 
  private:
    typename Primitive<F>::type p;
  };
 
  template <typename Integrator, typename F, typename X = no_var>
  struct NumericalIntegral {
    NumericalIntegral(const F& f, const Integrator& integrator) : f_(f), integrator_(integrator) {}
    inline double operator()(double min, double max) const { return integrator_(f_, min, max); }
 
  private:
    struct function {
      function(const F& f) : f_(f) {}
      inline double operator()(double x) const {
        X::set(x);
        return f_();
      }
 
    private:
      F f_;
    };
    function f_;
    Integrator integrator_;
  };
 
  template <typename Integrator, typename F>
  struct NumericalIntegral<Integrator, F, no_var> {
    NumericalIntegral(const F& f, const Integrator& integrator) : f_(f), integrator_(integrator) {}
    double operator()(double min, double max) const { return integrator_(f_, min, max); }
 
  private:
    F f_;
    Integrator integrator_;
  };
 
  template <typename F, typename X = no_var>
  struct Integral {
    typedef IntegralStruct<F, X> type;
  };
 
  template <typename X, typename F>
  typename Integral<F, X>::type integral(const F& f) {
    return typename Integral<F, X>::type(f);
  }
 
  template <typename X, typename F, typename Integrator>
  typename Integral<F, X>::type integral(const F& f, const Integrator& integrator) {
    return typename Integral<F, X>::type(f, integrator);
  }
 
  template <typename F, typename Integrator>
  typename Integral<F>::type integral_f(const F& f, const Integrator& integrator) {
    return typename Integral<F>::type(f, integrator);
  }
 
  template <typename X, typename F>
  double integral(const F& f, double min, double max) {
    return integral<X>(f)(min, max);
  }
 
  template <typename X, typename F, typename Integrator>
  double integral(const F& f, double min, double max, const Integrator& integrator) {
    return integral<X>(f, integrator)(min, max);
  }
 
  template <typename F>
  typename Integral<F>::type integral_f(const F& f) {
    return typename Integral<F>::type(f);
  }
 
  template <typename F>
  double integral_f(const F& f, double min, double max) {
    return integral_f(f)(min, max);
  }
 
  template <typename F, typename Integrator>
  double integral_f(const F& f, double min, double max, const Integrator& integrator) {
    return integral_f(f, integrator)(min, max);
  }
 
  template <typename F, typename MIN, typename MAX, typename Integrator = no_var, typename X = no_var>
  struct DefIntegral {
    DefIntegral(const F& f, const MIN& min, const MAX& max, const Integrator& integrator)
        : f_(f), min_(min), max_(max), integrator_(integrator) {}
    double operator()() const { return integral<X>(f_, min_(), max_(), integrator_); }
 
  private:
    F f_;
    MIN min_;
    MAX max_;
    Integrator integrator_;
  };
 
  template <typename F, typename MIN, typename MAX, typename X>
  struct DefIntegral<F, MIN, MAX, no_var, X> {
    DefIntegral(const F& f, const MIN& min, const MAX& max) : f_(f), min_(min), max_(max) {}
    double operator()(double x) const { return integral<X>(f_, min_(x), max_(x)); }
 
  private:
    F f_;
    MIN min_;
    MAX max_;
  };
 
  template <typename F, typename MIN, typename MAX, typename Integrator>
  struct DefIntegral<F, MIN, MAX, Integrator, no_var> {
    DefIntegral(const F& f, const MIN& min, const MAX& max, const Integrator& integrator)
        : f_(f), min_(min), max_(max), integrator_(integrator) {}
    double operator()(double x) const { return integral_f(f_, min_(x), max_(x), integrator_); }
 
  private:
    F f_;
    MIN min_;
    MAX max_;
    Integrator integrator_;
  };
 
  template <typename F, typename MIN, typename MAX>
  struct DefIntegral<F, MIN, MAX, no_var, no_var> {
    DefIntegral(const F& f, const MIN& min, const MAX& max) : f_(f), min_(min), max_(max) {}
    double operator()(double x) const { return integral_f(f_, min_(x), max_(x)); }
 
  private:
    F f_;
    MIN min_;
    MAX max_;
  };
 
}  // namespace funct
 
#define NUMERICAL_INTEGRAL(X, F, INTEGRATOR)            \
  namespace funct {                                     \
    template <typename X>                               \
    struct Integral<F, X> {                             \
      typedef NumericalIntegral<INTEGRATOR, F, X> type; \
    };                                                  \
  }                                                     \
  struct __useless_ignoreme
 
#define NUMERICAL_FUNCT_INTEGRAL(F, INTEGRATOR)      \
  namespace funct {                                  \
    template <>                                      \
    struct Integral<F, no_var> {                     \
      typedef NumericalIntegral<INTEGRATOR, F> type; \
    };                                               \
  }                                                  \
  struct __useless_ignoreme
 
#endif