CLHEP VERSION Reference Documentation
   
CLHEP Home Page     CLHEP Documentation     CLHEP Bug Reports

PtRelFcn.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 // $Id:
6 #include <assert.h>
7 #include <cmath> // for pow() and exp() and isfinite()
8 #include <float.h>
9 
10 #if (defined __STRICT_ANSI__) || (defined _WIN32)
11 #ifndef M_PI
12 #define M_PI 3.14159265358979323846
13 #endif // M_PI
14 #endif // __STRICT_ANSI__
15 
16 namespace Genfun {
17 FUNCTION_OBJECT_IMP(PtRelFcn)
18 
20  _p0("P0", 0, 0, 1),
21  _p1("P1", 0, 0, 2),
22  _p2("P2", 1, 0, 10),
23  _p3("P3", 0, 0, 10),
24  _p4("P4", 1.0, 0.1, 5.0),
25  _p5("P5", 0.0, 0, 50)
26 {}
27 
29 }
30 
32 _p0(right._p0),
33 _p1(right._p1),
34 _p2(right._p2),
35 _p3(right._p3),
36 _p4(right._p4),
37 _p5(right._p5)
38 {
39 }
40 
41 double PtRelFcn::operator() (double x) const {
42 
43  double p0 = _p0.getValue();
44  double p1 = _p1.getValue();
45  double p2 = _p2.getValue();
46  double p3 = _p3.getValue();
47  double p4 = _p4.getValue();
48  double p5 = _p5.getValue();
49 
50  //assert ((p0>=0.0) && (p0<=1.0));
51  if (p0<0.0) p0=FLT_MIN;
52  if (p0>1.0) p0=1.0-FLT_MIN;
53 
54  if (x<=0.0) return 1.0E-10;
55 
56  double n = (1+p1)/p3;
57  double a = (1/p3)*pow(p2,-n);
58 
59  double norm = 1.0/(a*exp(_logGamma(n)));
60  static const double s2 = sqrt(2.0);
61  double retVal=
62  norm*p0*pow(x,p1)*exp(-p2*pow(x,p3)) +
63  (2.0/(1+_erf(p5/p4/s2))*(1.0-p0)/(sqrt(2*M_PI)*p4))*exp(-(x-p5)*(x-p5)/(2.0*p4*p4));
64 
65  //if (!std::isfinite(retVal)) return 1.0E-10;
66 
67  return std::max(retVal,1.0E-10);
68 }
69 
71  return _p0;
72 }
73 
74 const Parameter & PtRelFcn::P0() const {
75  return _p0;
76 }
77 
79  return _p1;
80 }
81 
82 const Parameter & PtRelFcn::P1() const {
83  return _p1;
84 }
85 
87  return _p2;
88 }
89 
90 const Parameter & PtRelFcn::P2() const {
91  return _p2;
92 }
93 
95  return _p3;
96 }
97 
98 const Parameter & PtRelFcn::P3() const {
99  return _p3;
100 }
101 
103  return _p4;
104 }
105 
106 const Parameter & PtRelFcn::P4() const {
107  return _p4;
108 }
109 
111  return _p5;
112 }
113 
114 const Parameter & PtRelFcn::P5() const {
115  return _p5;
116 }
117 
118 
119 
120 
121 
122 } // namespace Genfun