Utilities/interface/Convolution.h

0001 #ifndef PhysicsTools_Utilities_Convolution_h
0002 #define PhysicsTools_Utilities_Convolution_h
0003 #include "FWCore/Utilities/interface/EDMException.h"
0004 #include "PhysicsTools/Utilities/interface/NumericalIntegration.h"
0005
0006 namespace funct {
0007   template <typename A, typename B, typename Integrator>
0008   class ConvolutionStruct {
0009   public:
0010     // min and max are defined in the domain of b
0011     ConvolutionStruct(const A& a, const B& b, double min, double max, const Integrator& integrator)
0012         : f_(a, b), min_(min), max_(max), integrator_(integrator) {
0013       if (max < min)
0014         throw edm::Exception(edm::errors::Configuration) << "Convolution: min must be smaller than max\n";
0015     }
0016     double operator()(double x) const {
0017       f_.setX(x);
0018       return integrator_(f_, x - max_, x - min_);
0019     }
0020
0021   private:
0022     struct function {
0023       function(const A& a, const B& b) : _1(a), _2(b) {}
0024       void setX(double x) const { x_ = x; }
0025       double operator()(double y) const { return _1(y) * _2(x_ - y); }
0026
0027     private:
0028       A _1;
0029       B _2;
0030       mutable double x_;
0031     };
0032     function f_;
0033     double min_, max_, delta_;
0034     Integrator integrator_;
0035   };
0036
0037   template <typename A, typename B, typename Integrator>
0038   struct Convolution {
0039     typedef ConvolutionStruct<A, B, Integrator> type;
0040     static type compose(const A& a, const B& b, double min, double max, const Integrator& i) {
0041       return type(a, b, min, max, i);
0042     }
0043   };
0044
0045   template <typename A, typename B, typename Integrator>
0046   inline typename funct::Convolution<A, B, Integrator>::type conv(
0047       const A& a, const B& b, double min, double max, const Integrator& i) {
0048     return funct::Convolution<A, B, Integrator>::compose(a, b, min, max, i);
0049   }
0050
0051 }  // namespace funct
0052
0053 #endif