CLHEP 2.0.4.7 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() 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 _p0(right._p0), 00033 _p1(right._p1), 00034 _p2(right._p2), 00035 _p3(right._p3), 00036 _p4(right._p4), 00037 _p5(right._p5) 00038 { 00039 } 00040 00041 double PtRelFcn::operator() (double x) const { 00042 00043 double p0 = _p0.getValue(); 00044 double p1 = _p1.getValue(); 00045 double p2 = _p2.getValue(); 00046 double p3 = _p3.getValue(); 00047 double p4 = _p4.getValue(); 00048 double p5 = _p5.getValue(); 00049 00050 //assert ((p0>=0.0) && (p0<=1.0)); 00051 if (p0<0.0) p0=FLT_MIN; 00052 if (p0>1.0) p0=1.0-FLT_MIN; 00053 00054 if (x<=0.0) return 1.0E-10; 00055 00056 double n = (1+p1)/p3; 00057 double a = (1/p3)*pow(p2,-n); 00058 00059 double norm = 1.0/(a*exp(_logGamma(n))); 00060 static const double s2 = sqrt(2.0); 00061 double retVal= 00062 norm*p0*pow(x,p1)*exp(-p2*pow(x,p3)) + 00063 (2.0/(1+_erf(p5/p4/s2))*(1.0-p0)/(sqrt(2*M_PI)*p4))*exp(-(x-p5)*(x-p5)/(2.0*p4*p4)); 00064 00065 //if (!finite(retVal)) return 1.0E-10; 00066 00067 return std::max(retVal,1.0E-10); 00068 } 00069 00070 Parameter & PtRelFcn::P0() { 00071 return _p0; 00072 } 00073 00074 const Parameter & PtRelFcn::P0() const { 00075 return _p0; 00076 } 00077 00078 Parameter & PtRelFcn::P1() { 00079 return _p1; 00080 } 00081 00082 const Parameter & PtRelFcn::P1() const { 00083 return _p1; 00084 } 00085 00086 Parameter & PtRelFcn::P2() { 00087 return _p2; 00088 } 00089 00090 const Parameter & PtRelFcn::P2() const { 00091 return _p2; 00092 } 00093 00094 Parameter & PtRelFcn::P3() { 00095 return _p3; 00096 } 00097 00098 const Parameter & PtRelFcn::P3() const { 00099 return _p3; 00100 } 00101 00102 Parameter & PtRelFcn::P4() { 00103 return _p4; 00104 } 00105 00106 const Parameter & PtRelFcn::P4() const { 00107 return _p4; 00108 } 00109 00110 Parameter & PtRelFcn::P5() { 00111 return _p5; 00112 } 00113 00114 const Parameter & PtRelFcn::P5() const { 00115 return _p5; 00116 } 00117 00118 00119 00120 00121 00122 } // namespace Genfun