CLHEP VERSION 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 AbsFunction(right), 00025 _peak(right._peak), 00026 _width(right._width) 00027 { 00028 } 00029 00030 double Landau::operator() (double x) const { 00031 double s = _width.getValue(); 00032 double x0 = _peak.getValue(); 00033 double xs = x0 + 0.222782*s; 00034 return _denlan((x-xs)/s)/s; 00035 } 00036 00037 Parameter & Landau::peak() { 00038 return _peak; 00039 } 00040 00041 Parameter & Landau::width() { 00042 return _width; 00043 } 00044 00045 const Parameter & Landau::peak() const { 00046 return _peak; 00047 } 00048 00049 const Parameter & Landau::width() const { 00050 return _width; 00051 } 00052 00053 double Landau::_denlan(double x) const { 00054 /* Initialized data */ 00055 00056 static float p1[5] = { (float).4259894875,(float)-.124976255,(float) 00057 .039842437,(float)-.006298287635,(float).001511162253 }; 00058 static float q5[5] = { (float)1.,(float)156.9424537,(float)3745.310488,( 00059 float)9834.698876,(float)66924.28357 }; 00060 static float p6[5] = { (float)1.000827619,(float)664.9143136,(float) 00061 62972.92665,(float)475554.6998,(float)-5743609.109 }; 00062 static float q6[5] = { (float)1.,(float)651.4101098,(float)56974.73333,( 00063 float)165917.4725,(float)-2815759.939 }; 00064 static float a1[3] = { (float).04166666667,(float)-.01996527778,(float) 00065 .02709538966 }; 00066 static float a2[2] = { (float)-1.84556867,(float)-4.284640743 }; 00067 static float q1[5] = { (float)1.,(float)-.3388260629,(float).09594393323,( 00068 float)-.01608042283,(float).003778942063 }; 00069 static float p2[5] = { (float).1788541609,(float).1173957403,(float) 00070 .01488850518,(float)-.001394989411,(float)1.283617211e-4 }; 00071 static float q2[5] = { (float)1.,(float).7428795082,(float).3153932961,( 00072 float).06694219548,(float).008790609714 }; 00073 static float p3[5] = { (float).1788544503,(float).09359161662,(float) 00074 .006325387654,(float)6.611667319e-5,(float)-2.031049101e-6 }; 00075 static float q3[5] = { (float)1.,(float).6097809921,(float).2560616665,( 00076 float).04746722384,(float).006957301675 }; 00077 static float p4[5] = { (float).9874054407,(float)118.6723273,(float) 00078 849.279436,(float)-743.7792444,(float)427.0262186 }; 00079 static float q4[5] = { (float)1.,(float)106.8615961,(float)337.6496214,( 00080 float)2016.712389,(float)1597.063511 }; 00081 static float p5[5] = { (float)1.003675074,(float)167.5702434,(float) 00082 4789.711289,(float)21217.86767,(float)-22324.9491 }; 00083 00084 /* System generated locals */ 00085 float ret_val, r__1; 00086 00087 /* Local variables */ 00088 static float u, v; 00089 v = x; 00090 if (v < (float)-5.5) { 00091 u = exp(v + (float)1.); 00092 ret_val = exp(-1 / u) / sqrt(u) * (float).3989422803 * ((a1[0] + (a1[ 00093 1] + a1[2] * u) * u) * u + 1); 00094 } else if (v < (float)-1.) { 00095 u = exp(-v - 1); 00096 ret_val = exp(-u) * sqrt(u) * (p1[0] + (p1[1] + (p1[2] + (p1[3] + p1[ 00097 4] * v) * v) * v) * v) / (q1[0] + (q1[1] + (q1[2] + (q1[3] + 00098 q1[4] * v) * v) * v) * v); 00099 } else if (v < (float)1.) { 00100 ret_val = (p2[0] + (p2[1] + (p2[2] + (p2[3] + p2[4] * v) * v) * v) * 00101 v) / (q2[0] + (q2[1] + (q2[2] + (q2[3] + q2[4] * v) * v) * v) 00102 * v); 00103 } else if (v < (float)5.) { 00104 ret_val = (p3[0] + (p3[1] + (p3[2] + (p3[3] + p3[4] * v) * v) * v) * 00105 v) / (q3[0] + (q3[1] + (q3[2] + (q3[3] + q3[4] * v) * v) * v) 00106 * v); 00107 } else if (v < (float)12.) { 00108 u = 1 / v; 00109 /* Computing 2nd power */ 00110 r__1 = u; 00111 ret_val = r__1 * r__1 * (p4[0] + (p4[1] + (p4[2] + (p4[3] + p4[4] * u) 00112 * u) * u) * u) / (q4[0] + (q4[1] + (q4[2] + (q4[3] + q4[4] * 00113 u) * u) * u) * u); 00114 } else if (v < (float)50.) { 00115 u = 1 / v; 00116 /* Computing 2nd power */ 00117 r__1 = u; 00118 ret_val = r__1 * r__1 * (p5[0] + (p5[1] + (p5[2] + (p5[3] + p5[4] * u) 00119 * u) * u) * u) / (q5[0] + (q5[1] + (q5[2] + (q5[3] + q5[4] * 00120 u) * u) * u) * u); 00121 } else if (v < (float)300.) { 00122 u = 1 / v; 00123 /* Computing 2nd power */ 00124 r__1 = u; 00125 ret_val = r__1 * r__1 * (p6[0] + (p6[1] + (p6[2] + (p6[3] + p6[4] * u) 00126 * u) * u) * u) / (q6[0] + (q6[1] + (q6[2] + (q6[3] + q6[4] * 00127 u) * u) * u) * u); 00128 } else { 00129 u = 1 / (v - v * log(v) / (v + 1)); 00130 /* Computing 2nd power */ 00131 r__1 = u; 00132 ret_val = r__1 * r__1 * ((a2[0] + a2[1] * u) * u + 1); 00133 } 00134 return ret_val; 00135 00136 } 00137 00138 } // namespace Genfun