CLHEP VERSION Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // -*- C++ -*- 00002 // $Id: 00003 #include "CLHEP/GenericFunctions/defs.h" 00004 #include "CLHEP/GenericFunctions/PtRelFcn.hh" 00005 #include "CLHEP/GenericFunctions/Variable.hh" 00006 #include <assert.h> 00007 #include <cmath> // for pow() and exp() and isfinite() 00008 #include <float.h> 00009 00010 #if (defined __STRICT_ANSI__) || (defined _WIN32) 00011 #ifndef M_PI 00012 #define M_PI 3.14159265358979323846 00013 #endif // M_PI 00014 #endif // __STRICT_ANSI__ 00015 00016 namespace Genfun { 00017 FUNCTION_OBJECT_IMP(PtRelFcn) 00018 00019 PtRelFcn::PtRelFcn(): 00020 _p0("P0", 0, 0, 1), 00021 _p1("P1", 0, 0, 2), 00022 _p2("P2", 1, 0, 10), 00023 _p3("P3", 0, 0, 10), 00024 _p4("P4", 1.0, 0.1, 5.0), 00025 _p5("P5", 0.0, 0, 50) 00026 {} 00027 00028 PtRelFcn::~PtRelFcn() { 00029 } 00030 00031 PtRelFcn::PtRelFcn(const PtRelFcn & right): 00032 AbsFunction(right), 00033 _p0(right._p0), 00034 _p1(right._p1), 00035 _p2(right._p2), 00036 _p3(right._p3), 00037 _p4(right._p4), 00038 _p5(right._p5) 00039 { 00040 } 00041 00042 double PtRelFcn::operator() (double x) const { 00043 00044 double p0 = _p0.getValue(); 00045 double p1 = _p1.getValue(); 00046 double p2 = _p2.getValue(); 00047 double p3 = _p3.getValue(); 00048 double p4 = _p4.getValue(); 00049 double p5 = _p5.getValue(); 00050 00051 //assert ((p0>=0.0) && (p0<=1.0)); 00052 if (p0<0.0) p0=FLT_MIN; 00053 if (p0>1.0) p0=1.0-FLT_MIN; 00054 00055 if (x<=0.0) return 1.0E-10; 00056 00057 double n = (1+p1)/p3; 00058 double a = (1/p3)*std::pow(p2,-n); 00059 00060 double norm = 1.0/(a*exp(_logGamma(n))); 00061 static const double s2 = sqrt(2.0); 00062 double retVal= 00063 norm*p0*std::pow(x,p1)*exp(-p2*std::pow(x,p3)) + 00064 (2.0/(1+_erf(p5/p4/s2))*(1.0-p0)/(sqrt(2*M_PI)*p4))*exp(-(x-p5)*(x-p5)/(2.0*p4*p4)); 00065 00066 //if (!std::isfinite(retVal)) return 1.0E-10; 00067 00068 return std::max(retVal,1.0E-10); 00069 } 00070 00071 Parameter & PtRelFcn::P0() { 00072 return _p0; 00073 } 00074 00075 const Parameter & PtRelFcn::P0() const { 00076 return _p0; 00077 } 00078 00079 Parameter & PtRelFcn::P1() { 00080 return _p1; 00081 } 00082 00083 const Parameter & PtRelFcn::P1() const { 00084 return _p1; 00085 } 00086 00087 Parameter & PtRelFcn::P2() { 00088 return _p2; 00089 } 00090 00091 const Parameter & PtRelFcn::P2() const { 00092 return _p2; 00093 } 00094 00095 Parameter & PtRelFcn::P3() { 00096 return _p3; 00097 } 00098 00099 const Parameter & PtRelFcn::P3() const { 00100 return _p3; 00101 } 00102 00103 Parameter & PtRelFcn::P4() { 00104 return _p4; 00105 } 00106 00107 const Parameter & PtRelFcn::P4() const { 00108 return _p4; 00109 } 00110 00111 Parameter & PtRelFcn::P5() { 00112 return _p5; 00113 } 00114 00115 const Parameter & PtRelFcn::P5() const { 00116 return _p5; 00117 } 00118 00119 00120 00121 00122 00123 } // namespace Genfun