CLHEP 2.0.4.7 Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // -*- C++ -*- 00002 // $Id: Landau.cc,v 1.4 2003/10/10 17:40:39 garren Exp $ 00003 // --------------------------------------------------------------------------- 00004 00005 #include "CLHEP/GenericFunctions/Landau.hh" 00006 #include "CLHEP/GenericFunctions/Variable.hh" 00007 #include <cmath> 00008 #include <assert.h> 00009 00010 using namespace std; 00011 00012 namespace Genfun { 00013 FUNCTION_OBJECT_IMP(Landau) 00014 00015 Landau::Landau(): 00016 _peak("Peak", 5.0 ,0, 10), 00017 _width("Width",1.0,0, 10) 00018 {} 00019 00020 Landau::~Landau() { 00021 } 00022 00023 Landau::Landau(const Landau & right): 00024 _peak(right._peak), 00025 _width(right._width) 00026 { 00027 } 00028 00029 double Landau::operator() (double x) const { 00030 double s = _width.getValue(); 00031 double x0 = _peak.getValue(); 00032 double xs = x0 + 0.222782*s; 00033 return _denlan((x-xs)/s)/s; 00034 } 00035 00036 Parameter & Landau::peak() { 00037 return _peak; 00038 } 00039 00040 Parameter & Landau::width() { 00041 return _width; 00042 } 00043 00044 const Parameter & Landau::peak() const { 00045 return _peak; 00046 } 00047 00048 const Parameter & Landau::width() const { 00049 return _width; 00050 } 00051 00052 double Landau::_denlan(double x) const { 00053 /* Initialized data */ 00054 00055 static float p1[5] = { (float).4259894875,(float)-.124976255,(float) 00056 .039842437,(float)-.006298287635,(float).001511162253 }; 00057 static float q5[5] = { (float)1.,(float)156.9424537,(float)3745.310488,( 00058 float)9834.698876,(float)66924.28357 }; 00059 static float p6[5] = { (float)1.000827619,(float)664.9143136,(float) 00060 62972.92665,(float)475554.6998,(float)-5743609.109 }; 00061 static float q6[5] = { (float)1.,(float)651.4101098,(float)56974.73333,( 00062 float)165917.4725,(float)-2815759.939 }; 00063 static float a1[3] = { (float).04166666667,(float)-.01996527778,(float) 00064 .02709538966 }; 00065 static float a2[2] = { (float)-1.84556867,(float)-4.284640743 }; 00066 static float q1[5] = { (float)1.,(float)-.3388260629,(float).09594393323,( 00067 float)-.01608042283,(float).003778942063 }; 00068 static float p2[5] = { (float).1788541609,(float).1173957403,(float) 00069 .01488850518,(float)-.001394989411,(float)1.283617211e-4 }; 00070 static float q2[5] = { (float)1.,(float).7428795082,(float).3153932961,( 00071 float).06694219548,(float).008790609714 }; 00072 static float p3[5] = { (float).1788544503,(float).09359161662,(float) 00073 .006325387654,(float)6.611667319e-5,(float)-2.031049101e-6 }; 00074 static float q3[5] = { (float)1.,(float).6097809921,(float).2560616665,( 00075 float).04746722384,(float).006957301675 }; 00076 static float p4[5] = { (float).9874054407,(float)118.6723273,(float) 00077 849.279436,(float)-743.7792444,(float)427.0262186 }; 00078 static float q4[5] = { (float)1.,(float)106.8615961,(float)337.6496214,( 00079 float)2016.712389,(float)1597.063511 }; 00080 static float p5[5] = { (float)1.003675074,(float)167.5702434,(float) 00081 4789.711289,(float)21217.86767,(float)-22324.9491 }; 00082 00083 /* System generated locals */ 00084 float ret_val, r__1; 00085 00086 /* Local variables */ 00087 static float u, v; 00088 v = x; 00089 if (v < (float)-5.5) { 00090 u = exp(v + (float)1.); 00091 ret_val = exp(-1 / u) / sqrt(u) * (float).3989422803 * ((a1[0] + (a1[ 00092 1] + a1[2] * u) * u) * u + 1); 00093 } else if (v < (float)-1.) { 00094 u = exp(-v - 1); 00095 ret_val = exp(-u) * sqrt(u) * (p1[0] + (p1[1] + (p1[2] + (p1[3] + p1[ 00096 4] * v) * v) * v) * v) / (q1[0] + (q1[1] + (q1[2] + (q1[3] + 00097 q1[4] * v) * v) * v) * v); 00098 } else if (v < (float)1.) { 00099 ret_val = (p2[0] + (p2[1] + (p2[2] + (p2[3] + p2[4] * v) * v) * v) * 00100 v) / (q2[0] + (q2[1] + (q2[2] + (q2[3] + q2[4] * v) * v) * v) 00101 * v); 00102 } else if (v < (float)5.) { 00103 ret_val = (p3[0] + (p3[1] + (p3[2] + (p3[3] + p3[4] * v) * v) * v) * 00104 v) / (q3[0] + (q3[1] + (q3[2] + (q3[3] + q3[4] * v) * v) * v) 00105 * v); 00106 } else if (v < (float)12.) { 00107 u = 1 / v; 00108 /* Computing 2nd power */ 00109 r__1 = u; 00110 ret_val = r__1 * r__1 * (p4[0] + (p4[1] + (p4[2] + (p4[3] + p4[4] * u) 00111 * u) * u) * u) / (q4[0] + (q4[1] + (q4[2] + (q4[3] + q4[4] * 00112 u) * u) * u) * u); 00113 } else if (v < (float)50.) { 00114 u = 1 / v; 00115 /* Computing 2nd power */ 00116 r__1 = u; 00117 ret_val = r__1 * r__1 * (p5[0] + (p5[1] + (p5[2] + (p5[3] + p5[4] * u) 00118 * u) * u) * u) / (q5[0] + (q5[1] + (q5[2] + (q5[3] + q5[4] * 00119 u) * u) * u) * u); 00120 } else if (v < (float)300.) { 00121 u = 1 / v; 00122 /* Computing 2nd power */ 00123 r__1 = u; 00124 ret_val = r__1 * r__1 * (p6[0] + (p6[1] + (p6[2] + (p6[3] + p6[4] * u) 00125 * u) * u) * u) / (q6[0] + (q6[1] + (q6[2] + (q6[3] + q6[4] * 00126 u) * u) * u) * u); 00127 } else { 00128 u = 1 / (v - v * log(v) / (v + 1)); 00129 /* Computing 2nd power */ 00130 r__1 = u; 00131 ret_val = r__1 * r__1 * ((a2[0] + a2[1] * u) * u + 1); 00132 } 00133 return ret_val; 00134 00135 } 00136 00137 } // namespace Genfun