00001 #ifndef PhysicsTools_Utilities_Convolution_h
00002 #define PhysicsTools_Utilities_Convolution_h
00003 #include "FWCore/Utilities/interface/EDMException.h"
00004 #include "PhysicsTools/Utilities/interface/NumericalIntegration.h"
00005
00006 namespace funct {
00007 template<typename A, typename B, typename Integrator>
00008 class ConvolutionStruct {
00009 public:
00010
00011 ConvolutionStruct(const A& a, const B& b,
00012 double min, double max, const Integrator & integrator) :
00013 f_(a, b), min_(min), max_(max), integrator_(integrator) {
00014 if(max < min)
00015 throw edm::Exception(edm::errors::Configuration)
00016 << "Convolution: min must be smaller than max\n";
00017 }
00018 double operator()(double x) const {
00019 f_.setX(x);
00020 return integrator_(f_, x - max_, x - min_);
00021 }
00022 private:
00023 struct function {
00024 function(const A& a, const B& b) : _1(a), _2(b) { }
00025 void setX(double x) const { x_ = x; }
00026 double operator()(double y) const {
00027 return _1(y) * _2(x_ - y);
00028 }
00029 private:
00030 A _1;
00031 B _2;
00032 mutable double x_;
00033 };
00034 function f_;
00035 double min_, max_, delta_;
00036 Integrator integrator_;
00037 };
00038
00039 template<typename A, typename B, typename Integrator>
00040 struct Convolution {
00041 typedef ConvolutionStruct<A, B, Integrator> type;
00042 static type compose(const A& a, const B& b, double min, double max, const Integrator& i) {
00043 return type(a, b, min, max, i);
00044 }
00045 };
00046
00047 template<typename A, typename B, typename Integrator>
00048 inline typename funct::Convolution<A, B, Integrator>::type conv(const A& a, const B& b, double min, double max, const Integrator& i) {
00049 return funct::Convolution<A, B, Integrator>::compose(a, b, min, max, i);
00050 }
00051
00052 }
00053
00054
00055
00056 #endif