00001 #ifndef PhysicsTools_Utilities_Integral_h 00002 #define PhysicsTools_Utilities_Integral_h 00003 #include "PhysicsTools/Utilities/interface/Primitive.h" 00004 #include "PhysicsTools/Utilities/interface/NumericalIntegration.h" 00005 00006 namespace funct { 00007 00008 struct no_var; 00009 00010 template<typename F, typename X = no_var> 00011 struct IntegralStruct { 00012 IntegralStruct(const F& f) : p(primitive<X>(f)) { } 00013 double operator()(double min, double max) const { 00014 X::set(min); double pMin = p(); 00015 X::set(max); double pMax = p(); 00016 return pMax - pMin; 00017 } 00018 private: 00019 typename Primitive<F, X>::type p; 00020 }; 00021 00022 template<typename F> 00023 struct IntegralStruct<F> { 00024 IntegralStruct(const F& f) : p(primitive(f)) { } 00025 double operator()(double min, double max) const { 00026 return p(max) - p(min); 00027 } 00028 private: 00029 typename Primitive<F>::type p; 00030 }; 00031 00032 template<typename Integrator, typename F, typename X = no_var> 00033 struct NumericalIntegral { 00034 NumericalIntegral(const F& f, const Integrator & integrator) : 00035 f_(f), integrator_(integrator) { } 00036 inline double operator()(double min, double max) const { 00037 return integrator_(f_, min, max); 00038 } 00039 private: 00040 struct function { 00041 function(const F & f) : f_(f) { } 00042 inline double operator()(double x) const { 00043 X::set(x); return f_(); 00044 } 00045 private: 00046 F f_; 00047 }; 00048 function f_; 00049 Integrator integrator_; 00050 }; 00051 00052 template<typename Integrator, typename F> 00053 struct NumericalIntegral<Integrator, F, no_var> { 00054 NumericalIntegral(const F& f, const Integrator & integrator) : 00055 f_(f), integrator_(integrator) { } 00056 double operator()(double min, double max) const { 00057 return integrator_(f_, min, max); 00058 } 00059 private: 00060 F f_; 00061 Integrator integrator_; 00062 }; 00063 00064 template<typename F, typename X = no_var> struct Integral { 00065 typedef IntegralStruct<F, X> type; 00066 }; 00067 00068 template<typename X, typename F> 00069 typename Integral<F, X>::type integral(const F& f) { 00070 return typename Integral<F, X>::type(f); 00071 } 00072 00073 template<typename X, typename F, typename Integrator> 00074 typename Integral<F, X>::type integral(const F& f, const Integrator & integrator) { 00075 return typename Integral<F, X>::type(f, integrator); 00076 } 00077 00078 template<typename F, typename Integrator> 00079 typename Integral<F>::type integral_f(const F& f, const Integrator & integrator) { 00080 return typename Integral<F>::type(f, integrator); 00081 } 00082 00083 template<typename X, typename F> 00084 double integral(const F& f, double min, double max) { 00085 return integral<X>(f)(min, max); 00086 } 00087 00088 template<typename X, typename F, typename Integrator> 00089 double integral(const F& f, double min, double max, const Integrator & integrator) { 00090 return integral<X>(f, integrator)(min, max); 00091 } 00092 00093 template<typename F> 00094 typename Integral<F>::type integral_f(const F& f) { 00095 return typename Integral<F>::type(f); 00096 } 00097 00098 template<typename F> 00099 double integral_f(const F& f, double min, double max) { 00100 return integral_f(f)(min, max); 00101 } 00102 00103 template<typename F, typename Integrator> 00104 double integral_f(const F& f, double min, double max, const Integrator & integrator) { 00105 return integral_f(f, integrator)(min, max); 00106 } 00107 00108 template<typename F, typename MIN, typename MAX, typename Integrator = no_var, typename X = no_var> 00109 struct DefIntegral { 00110 DefIntegral(const F & f, const MIN & min, const MAX & max, const Integrator & integrator) : 00111 f_(f), min_(min), max_(max), integrator_(integrator) { } 00112 double operator()() const { 00113 return integral<X>(f_, min_(), max_(), integrator_); 00114 } 00115 private: 00116 F f_; 00117 MIN min_; 00118 MAX max_; 00119 Integrator integrator_; 00120 }; 00121 00122 template<typename F, typename MIN, typename MAX, typename X> 00123 struct DefIntegral<F, MIN, MAX, no_var, X> { 00124 DefIntegral(const F & f, const MIN & min, const MAX & max) : 00125 f_(f), min_(min), max_(max) { } 00126 double operator()(double x) const { 00127 return integral<X>(f_, min_(x), max_(x)); 00128 } 00129 private: 00130 F f_; 00131 MIN min_; 00132 MAX max_; 00133 }; 00134 00135 template<typename F, typename MIN, typename MAX, typename Integrator> 00136 struct DefIntegral<F, MIN, MAX, Integrator, no_var> { 00137 DefIntegral(const F & f, const MIN & min, const MAX & max, const Integrator & integrator) : 00138 f_(f), min_(min), max_(max), integrator_(integrator) { } 00139 double operator()(double x) const { 00140 return integral_f(f_, min_(x), max_(x), integrator_); 00141 } 00142 private: 00143 F f_; 00144 MIN min_; 00145 MAX max_; 00146 Integrator integrator_; 00147 }; 00148 00149 template<typename F, typename MIN, typename MAX> 00150 struct DefIntegral<F, MIN, MAX, no_var, no_var> { 00151 DefIntegral(const F & f, const MIN & min, const MAX & max) : f_(f), min_(min), max_(max) { } 00152 double operator()(double x) const { 00153 return integral_f(f_, min_(x), max_(x)); 00154 } 00155 private: 00156 F f_; 00157 MIN min_; 00158 MAX max_; 00159 }; 00160 00161 } 00162 00163 #define NUMERICAL_INTEGRAL(X, F, INTEGRATOR) \ 00164 namespace funct { \ 00165 template<typename X> struct Integral<F, X> { \ 00166 typedef NumericalIntegral<INTEGRATOR, F, X> type; \ 00167 }; \ 00168 } \ 00169 struct __useless_ignoreme 00170 00171 #define NUMERICAL_FUNCT_INTEGRAL(F, INTEGRATOR) \ 00172 namespace funct { \ 00173 template<> struct Integral<F, no_var> { \ 00174 typedef NumericalIntegral<INTEGRATOR, F> type; \ 00175 }; \ 00176 } \ 00177 struct __useless_ignoreme 00178 00179 #endif