Convolution.h

Go to the documentation of this file.

#ifndef PhysicsTools_Utilities_Convolution_h
#define PhysicsTools_Utilities_Convolution_h
#include "FWCore/Utilities/interface/EDMException.h"
#include "PhysicsTools/Utilities/interface/NumericalIntegration.h"

namespace funct {
  template <typename A, typename B, typename Integrator>
  class ConvolutionStruct {
  public:
    // min and max are defined in the domain of b
    ConvolutionStruct(const A& a, const B& b, double min, double max, const Integrator& integrator)
        : f_(a, b), min_(min), max_(max), integrator_(integrator) {
      if (max < min)
        throw edm::Exception(edm::errors::Configuration) << "Convolution: min must be smaller than max\n";
    }
    double operator()(double x) const {
      f_.setX(x);
      return integrator_(f_, x - max_, x - min_);
    }

  private:
    struct function {
      function(const A& a, const B& b) : _1(a), _2(b) {}
      void setX(double x) const { x_ = x; }
      double operator()(double y) const { return _1(y) * _2(x_ - y); }

    private:
      A _1;
      B _2;
      mutable double x_;
    };
    function f_;
    double min_, max_, delta_;
    Integrator integrator_;
  };

  template <typename A, typename B, typename Integrator>
  struct Convolution {
    typedef ConvolutionStruct<A, B, Integrator> type;
    static type compose(const A& a, const B& b, double min, double max, const Integrator& i) {
      return type(a, b, min, max, i);
    }
  };

  template <typename A, typename B, typename Integrator>
  inline typename funct::Convolution<A, B, Integrator>::type conv(
      const A& a, const B& b, double min, double max, const Integrator& i) {
    return funct::Convolution<A, B, Integrator>::compose(a, b, min, max, i);
  }

}  // namespace funct

#endif