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

Landau.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 // $Id: Landau.cc,v 1.4 2003/10/10 17:40:39 garren Exp $
3 // ---------------------------------------------------------------------------
4 
7 #include <cmath>
8 #include <assert.h>
9 
10 using namespace std;
11 
12 namespace Genfun {
13 FUNCTION_OBJECT_IMP(Landau)
14 
16  _peak("Peak", 5.0 ,0, 10),
17  _width("Width",1.0,0, 10)
18 {}
19 
20 Landau::~Landau() {
21 }
22 
23 Landau::Landau(const Landau & right):
24 AbsFunction(right),
25 _peak(right._peak),
26 _width(right._width)
27 {
28 }
29 
30 double Landau::operator() (double x) const {
31  double s = _width.getValue();
32  double x0 = _peak.getValue();
33  double xs = x0 + 0.222782*s;
34  return _denlan((x-xs)/s)/s;
35 }
36 
38  return _peak;
39 }
40 
42  return _width;
43 }
44 
45 const Parameter & Landau::peak() const {
46  return _peak;
47 }
48 
49 const Parameter & Landau::width() const {
50  return _width;
51 }
52 
53 double Landau::_denlan(double x) const {
54  /* Initialized data */
55 
56  static float p1[5] = { (float).4259894875,(float)-.124976255,(float)
57  .039842437,(float)-.006298287635,(float).001511162253 };
58  static float q5[5] = { (float)1.,(float)156.9424537,(float)3745.310488,(
59  float)9834.698876,(float)66924.28357 };
60  static float p6[5] = { (float)1.000827619,(float)664.9143136,(float)
61  62972.92665,(float)475554.6998,(float)-5743609.109 };
62  static float q6[5] = { (float)1.,(float)651.4101098,(float)56974.73333,(
63  float)165917.4725,(float)-2815759.939 };
64  static float a1[3] = { (float).04166666667,(float)-.01996527778,(float)
65  .02709538966 };
66  static float a2[2] = { (float)-1.84556867,(float)-4.284640743 };
67  static float q1[5] = { (float)1.,(float)-.3388260629,(float).09594393323,(
68  float)-.01608042283,(float).003778942063 };
69  static float p2[5] = { (float).1788541609,(float).1173957403,(float)
70  .01488850518,(float)-.001394989411,(float)1.283617211e-4 };
71  static float q2[5] = { (float)1.,(float).7428795082,(float).3153932961,(
72  float).06694219548,(float).008790609714 };
73  static float p3[5] = { (float).1788544503,(float).09359161662,(float)
74  .006325387654,(float)6.611667319e-5,(float)-2.031049101e-6 };
75  static float q3[5] = { (float)1.,(float).6097809921,(float).2560616665,(
76  float).04746722384,(float).006957301675 };
77  static float p4[5] = { (float).9874054407,(float)118.6723273,(float)
78  849.279436,(float)-743.7792444,(float)427.0262186 };
79  static float q4[5] = { (float)1.,(float)106.8615961,(float)337.6496214,(
80  float)2016.712389,(float)1597.063511 };
81  static float p5[5] = { (float)1.003675074,(float)167.5702434,(float)
82  4789.711289,(float)21217.86767,(float)-22324.9491 };
83 
84  /* System generated locals */
85  float ret_val, r__1;
86 
87  /* Local variables */
88  static float u, v;
89  v = x;
90  if (v < (float)-5.5) {
91  u = exp(v + (float)1.);
92  ret_val = exp(-1 / u) / sqrt(u) * (float).3989422803 * ((a1[0] + (a1[
93  1] + a1[2] * u) * u) * u + 1);
94  } else if (v < (float)-1.) {
95  u = exp(-v - 1);
96  ret_val = exp(-u) * sqrt(u) * (p1[0] + (p1[1] + (p1[2] + (p1[3] + p1[
97  4] * v) * v) * v) * v) / (q1[0] + (q1[1] + (q1[2] + (q1[3] +
98  q1[4] * v) * v) * v) * v);
99  } else if (v < (float)1.) {
100  ret_val = (p2[0] + (p2[1] + (p2[2] + (p2[3] + p2[4] * v) * v) * v) *
101  v) / (q2[0] + (q2[1] + (q2[2] + (q2[3] + q2[4] * v) * v) * v)
102  * v);
103  } else if (v < (float)5.) {
104  ret_val = (p3[0] + (p3[1] + (p3[2] + (p3[3] + p3[4] * v) * v) * v) *
105  v) / (q3[0] + (q3[1] + (q3[2] + (q3[3] + q3[4] * v) * v) * v)
106  * v);
107  } else if (v < (float)12.) {
108  u = 1 / v;
109 /* Computing 2nd power */
110  r__1 = u;
111  ret_val = r__1 * r__1 * (p4[0] + (p4[1] + (p4[2] + (p4[3] + p4[4] * u)
112  * u) * u) * u) / (q4[0] + (q4[1] + (q4[2] + (q4[3] + q4[4] *
113  u) * u) * u) * u);
114  } else if (v < (float)50.) {
115  u = 1 / v;
116 /* Computing 2nd power */
117  r__1 = u;
118  ret_val = r__1 * r__1 * (p5[0] + (p5[1] + (p5[2] + (p5[3] + p5[4] * u)
119  * u) * u) * u) / (q5[0] + (q5[1] + (q5[2] + (q5[3] + q5[4] *
120  u) * u) * u) * u);
121  } else if (v < (float)300.) {
122  u = 1 / v;
123 /* Computing 2nd power */
124  r__1 = u;
125  ret_val = r__1 * r__1 * (p6[0] + (p6[1] + (p6[2] + (p6[3] + p6[4] * u)
126  * u) * u) * u) / (q6[0] + (q6[1] + (q6[2] + (q6[3] + q6[4] *
127  u) * u) * u) * u);
128  } else {
129  u = 1 / (v - v * log(v) / (v + 1));
130 /* Computing 2nd power */
131  r__1 = u;
132  ret_val = r__1 * r__1 * ((a2[0] + a2[1] * u) * u + 1);
133  }
134  return ret_val;
135 
136 }
137 
138 } // namespace Genfun
virtual double operator()(double argument) const
Definition: Landau.cc:30
Parameter & peak()
Definition: Landau.cc:37
#define FUNCTION_OBJECT_IMP(classname)
Parameter & width()
Definition: Landau.cc:41
virtual double getValue() const
Definition: Parameter.cc:27