00001
00012 #ifdef _MSC_VER
00013 #include "msdevstudio/MSconfig.h"
00014 #endif
00015
00016 #include "Landau.h"
00017
00018 #include "FunctionHelper.h"
00019
00020 #include <cmath>
00021 #include <cassert>
00022
00023 using std::distance;
00024
00025 #ifdef ITERATOR_MEMBER_DEFECT
00026 using namespace std;
00027 #else
00028 using std::exp;
00029 using std::vector;
00030 #endif
00031
00032 namespace hippodraw {
00033
00034 static double p1[5] = {0.4259894875,-0.1249762550, 0.03984243700, -0.006298287635, 0.001511162253};
00035 static double q1[5] = {1.0 ,-0.3388260629, 0.09594393323, -0.01608042283, 0.003778942063};
00036
00037 static double p2[5] = {0.1788541609, 0.1173957403, 0.01488850518, -0.001394989411, 0.0001283617211};
00038 static double q2[5] = {1.0 , 0.7428795082, 0.3153932961, 0.06694219548, 0.008790609714};
00039
00040 static double p3[5] = {0.1788544503, 0.09359161662,0.006325387654, 0.00006611667319,-0.000002031049101};
00041 static double q3[5] = {1.0 , 0.6097809921, 0.2560616665, 0.04746722384, 0.006957301675};
00042
00043 static double p4[5] = {0.9874054407, 118.6723273, 849.2794360, -743.7792444, 427.0262186};
00044 static double q4[5] = {1.0 , 106.8615961, 337.6496214, 2016.712389, 1597.063511};
00045
00046 static double p5[5] = {1.003675074, 167.5702434, 4789.711289, 21217.86767, -22324.94910};
00047 static double q5[5] = {1.0 , 156.9424537, 3745.310488, 9834.698876, 66924.28357};
00048
00049 static double p6[5] = {1.000827619, 664.9143136, 62972.92665, 475554.6998, -5743609.109};
00050 static double q6[5] = {1.0 , 651.4101098, 56974.73333, 165917.4725, -2815759.939};
00051
00052 static double a1[3] = {0.04166666667,-0.01996527778, 0.02709538966};
00053
00054 static double a2[2] = {-1.845568670,-4.284640743};
00055
00056
00057 Landau::Landau ( )
00058 {
00059 initialize ();
00060 }
00061
00062 Landau::Landau ( double p, double c, double s )
00063 {
00064 initialize ();
00065
00066 m_parms[peak] = p;
00067 m_parms[norm] = c;
00068 m_parms[sigma] = s;
00069 }
00070
00071 void Landau::initialize ()
00072 {
00073 m_name = "Landau";
00074
00075 m_parm_names.push_back ( "Peak" );
00076 m_parm_names.push_back ( "Normalization" );
00077 m_parm_names.push_back ( "Sigma" );
00078
00079 resize ();
00080 }
00081
00082 FunctionBase * Landau::clone () const
00083 {
00084 return new Landau ( *this );
00085 }
00086
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104 double Landau::operator () ( double x ) const
00105 {
00106 if(m_parms[sigma]==0) return 0;
00107 double v = calcY ( x );
00108 double u, ue, us, den;
00109 if (v < -5.5) {
00110 u = exp(v+1.0);
00111 if (u < 1e-10) return 0.0;
00112 ue = exp(-1/u);
00113 us = sqrt(u);
00114 den = 0.3989422803*(ue/us)*(1+(a1[0]+(a1[1]+a1[2]*u)*u)*u);
00115 } else if(v < -1) {
00116 u = exp(-v-1);
00117 den = exp(-u)*sqrt(u)*
00118 (p1[0]+(p1[1]+(p1[2]+(p1[3]+p1[4]*v)*v)*v)*v)/
00119 (q1[0]+(q1[1]+(q1[2]+(q1[3]+q1[4]*v)*v)*v)*v);
00120 } else if(v < 1) {
00121 den = (p2[0]+(p2[1]+(p2[2]+(p2[3]+p2[4]*v)*v)*v)*v)/
00122 (q2[0]+(q2[1]+(q2[2]+(q2[3]+q2[4]*v)*v)*v)*v);
00123 } else if(v < 5) {
00124 den = (p3[0]+(p3[1]+(p3[2]+(p3[3]+p3[4]*v)*v)*v)*v)/
00125 (q3[0]+(q3[1]+(q3[2]+(q3[3]+q3[4]*v)*v)*v)*v);
00126 } else if(v < 12) {
00127 u = 1/v;
00128 den = u*u*(p4[0]+(p4[1]+(p4[2]+(p4[3]+p4[4]*u)*u)*u)*u)/
00129 (q4[0]+(q4[1]+(q4[2]+(q4[3]+q4[4]*u)*u)*u)*u);
00130 } else if(v < 50) {
00131 u = 1/v;
00132 den = u*u*(p5[0]+(p5[1]+(p5[2]+(p5[3]+p5[4]*u)*u)*u)*u)/
00133 (q5[0]+(q5[1]+(q5[2]+(q5[3]+q5[4]*u)*u)*u)*u);
00134 } else if(v < 300) {
00135 u = 1/v;
00136 den = u*u*(p6[0]+(p6[1]+(p6[2]+(p6[3]+p6[4]*u)*u)*u)*u)/
00137 (q6[0]+(q6[1]+(q6[2]+(q6[3]+q6[4]*u)*u)*u)*u);
00138 } else {
00139 u = 1/(v-v*log(v)/(v+1));
00140 den = u*u*(1+(a2[0]+a2[1]*u)*u);
00141 }
00142 return m_parms[norm] * den/sigma;
00143
00144
00145
00146
00147
00148 }
00149
00150
00151 void
00152 Landau::
00153 initialParameters ( const FunctionHelper * helper )
00154 {
00155 m_parms[norm] = helper->maxValue () * sqrt ( 2.0 * M_PI * M_E );
00156 m_parms[peak] = helper->meanCoord ();
00157 m_parms[sigma] = helper->stdCoord ();
00158 }
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182 double Landau::derivByNorm ( double x ) const
00183 {
00184 double norm_aux = 0.0001;
00185 if(m_parms[norm] != 0) norm_aux = m_parms[norm];
00186 return operator () (x) / norm_aux;
00187 }
00188
00189 double Landau::derivByPeak ( double x ) const
00190 {
00191 return operator () ( x ) * calcZ ( x ) * ( ( -1.0 ) / m_parms[sigma] );
00192 }
00193
00194 double Landau::derivBySigma ( double x ) const
00195 {
00196 return operator () ( x ) * calcZ ( x )
00197 * ( - ( x - m_parms[peak] ) / ( m_parms[sigma] * m_parms[sigma] ) );
00198 }
00199
00200 bool
00201 Landau::
00202 hasDerivatives () const
00203 {
00204 return false;
00205 }
00206
00207
00208 }
00209