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