Landau.cxx
Go to the documentation of this file.
1 
12 #ifdef _MSC_VER
13 #include "msdevstudio/MSconfig.h"
14 #endif
15 
16 #include "Landau.h"
17 
18 #include "FunctionHelper.h"
19 
20 #include <cmath>
21 #include <cassert>
22 
23 using std::distance;
24 
25 #ifdef ITERATOR_MEMBER_DEFECT
26 using namespace std;
27 #else
28 using std::exp;
29 using std::vector;
30 #endif
31 
32 namespace hippodraw {
33 
34  static double p1[5] = {0.4259894875,-0.1249762550, 0.03984243700, -0.006298287635, 0.001511162253};
35  static double q1[5] = {1.0 ,-0.3388260629, 0.09594393323, -0.01608042283, 0.003778942063};
36 
37  static double p2[5] = {0.1788541609, 0.1173957403, 0.01488850518, -0.001394989411, 0.0001283617211};
38  static double q2[5] = {1.0 , 0.7428795082, 0.3153932961, 0.06694219548, 0.008790609714};
39 
40  static double p3[5] = {0.1788544503, 0.09359161662,0.006325387654, 0.00006611667319,-0.000002031049101};
41  static double q3[5] = {1.0 , 0.6097809921, 0.2560616665, 0.04746722384, 0.006957301675};
42 
43  static double p4[5] = {0.9874054407, 118.6723273, 849.2794360, -743.7792444, 427.0262186};
44  static double q4[5] = {1.0 , 106.8615961, 337.6496214, 2016.712389, 1597.063511};
45 
46  static double p5[5] = {1.003675074, 167.5702434, 4789.711289, 21217.86767, -22324.94910};
47  static double q5[5] = {1.0 , 156.9424537, 3745.310488, 9834.698876, 66924.28357};
48 
49  static double p6[5] = {1.000827619, 664.9143136, 62972.92665, 475554.6998, -5743609.109};
50  static double q6[5] = {1.0 , 651.4101098, 56974.73333, 165917.4725, -2815759.939};
51 
52  static double a1[3] = {0.04166666667,-0.01996527778, 0.02709538966};
53 
54  static double a2[2] = {-1.845568670,-4.284640743};
55 
56 
57 Landau::Landau ( )
58 {
59  initialize ();
60 }
61 
62 Landau::Landau ( double p, double c, double s )
63 {
64  initialize ();
65 
66  m_parms[peak] = p;
67  m_parms[norm] = c;
68  m_parms[sigma] = s;
69 }
70 
71 void Landau::initialize ()
72 {
73  m_name = "Landau";
74 
75  m_parm_names.push_back ( "Peak" );
76  m_parm_names.push_back ( "Normalization" );
77  m_parm_names.push_back ( "Sigma" );
78 
79  resize ();
80 }
81 
83 {
84  return new Landau ( *this );
85 }
86 
88 // REAL FUNCTION FITLAND(X)
89 
90 // DOUBLE PRECISION FITEMP
91 
92 // COMMON/PAWPAR/ PAR(3)
93 
94 
95 // PI=3.1415926
96 
97 // Y=(X-PAR(1))/PAR(3)
98 
99 // FITEMP=DBLE(PAR(2)*EXP(-0.5*(Y+EXP(-1.*Y)))/SQRT(2.*PI))
100 
101 // FITLAND=REAL(FITEMP)
102 
103 // END
104 double Landau::operator () ( double x ) const
105 {
106  if(m_parms[sigma]==0) return 0;
107  double v = calcY ( x );
108  double u, ue, us, den;
109  if (v < -5.5) {
110  u = exp(v+1.0);
111  if (u < 1e-10) return 0.0;
112  ue = exp(-1/u);
113  us = sqrt(u);
114  den = 0.3989422803*(ue/us)*(1+(a1[0]+(a1[1]+a1[2]*u)*u)*u);
115  } else if(v < -1) {
116  u = exp(-v-1);
117  den = exp(-u)*sqrt(u)*
118  (p1[0]+(p1[1]+(p1[2]+(p1[3]+p1[4]*v)*v)*v)*v)/
119  (q1[0]+(q1[1]+(q1[2]+(q1[3]+q1[4]*v)*v)*v)*v);
120  } else if(v < 1) {
121  den = (p2[0]+(p2[1]+(p2[2]+(p2[3]+p2[4]*v)*v)*v)*v)/
122  (q2[0]+(q2[1]+(q2[2]+(q2[3]+q2[4]*v)*v)*v)*v);
123  } else if(v < 5) {
124  den = (p3[0]+(p3[1]+(p3[2]+(p3[3]+p3[4]*v)*v)*v)*v)/
125  (q3[0]+(q3[1]+(q3[2]+(q3[3]+q3[4]*v)*v)*v)*v);
126  } else if(v < 12) {
127  u = 1/v;
128  den = u*u*(p4[0]+(p4[1]+(p4[2]+(p4[3]+p4[4]*u)*u)*u)*u)/
129  (q4[0]+(q4[1]+(q4[2]+(q4[3]+q4[4]*u)*u)*u)*u);
130  } else if(v < 50) {
131  u = 1/v;
132  den = u*u*(p5[0]+(p5[1]+(p5[2]+(p5[3]+p5[4]*u)*u)*u)*u)/
133  (q5[0]+(q5[1]+(q5[2]+(q5[3]+q5[4]*u)*u)*u)*u);
134  } else if(v < 300) {
135  u = 1/v;
136  den = u*u*(p6[0]+(p6[1]+(p6[2]+(p6[3]+p6[4]*u)*u)*u)*u)/
137  (q6[0]+(q6[1]+(q6[2]+(q6[3]+q6[4]*u)*u)*u)*u);
138  } else {
139  u = 1/(v-v*log(v)/(v+1));
140  den = u*u*(1+(a2[0]+a2[1]*u)*u);
141  }
142  return m_parms[norm] * den/sigma;
143 
144  //Moyal old formula
145  // double t = exp ( -0.5 * ( y + exp ( -1.0 * y ) ) )
146  // / sqrt ( 2.0 * M_PI );
147  // return m_parms[norm] * t;
148 }
149 
150 /* virtual */
151 void
152 Landau::
153 initialParameters ( const FunctionHelper * helper )
154 {
155  m_parms[norm] = helper->maxValue () * sqrt ( 2.0 * M_PI * M_E );
156  m_parms[peak] = helper->meanCoord ();
157  m_parms[sigma] = helper->stdCoord ();
158 }
159 
160 // double Landau::derivByParm ( int i, double x ) const
161 // {
162 // switch ( i ) {
163 // case peak :
164 // return derivByPeak ( x );
165 // break;
166 
167 // case norm :
168 // return derivByNorm ( x );
169 // break;
170 
171 // case sigma :
172 // return derivBySigma ( x );
173 // break;
174 
175 // default :
176 // assert ( false );
177 // break;
178 // }
179 // return 0.0;
180 //}
181 
182 double Landau::derivByNorm ( double x ) const
183 {
184  double norm_aux = 0.0001;
185  if(m_parms[norm] != 0) norm_aux = m_parms[norm];
186  return operator () (x) / norm_aux;
187 }
188 
189 double Landau::derivByPeak ( double x ) const
190 {
191  return operator () ( x ) * calcZ ( x ) * ( ( -1.0 ) / m_parms[sigma] );
192 }
193 
194 double Landau::derivBySigma ( double x ) const
195 {
196  return operator () ( x ) * calcZ ( x )
197  * ( - ( x - m_parms[peak] ) / ( m_parms[sigma] * m_parms[sigma] ) );
198 }
199 
200 bool
201 Landau::
202 hasDerivatives () const
203 {
204  return false;
205 }
206 
207 
208 } // namespace hippodraw
209 

Generated for HippoDraw Class Library by doxygen