Landau.cxx

Go to the documentation of this file.
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 //                            REAL FUNCTION FITLAND(X)
00089 
00090 //                             DOUBLE PRECISION FITEMP
00091 
00092 //                             COMMON/PAWPAR/ PAR(3)
00093 
00094 
00095 //                             PI=3.1415926
00096 
00097 //                            Y=(X-PAR(1))/PAR(3)
00098 
00099 //            FITEMP=DBLE(PAR(2)*EXP(-0.5*(Y+EXP(-1.*Y)))/SQRT(2.*PI))
00100 
00101 //                               FITLAND=REAL(FITEMP)
00102 
00103 //                            END
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    //Moyal old formula
00145   //  double t = exp ( -0.5 * ( y + exp ( -1.0 * y ) ) ) 
00146   //    / sqrt ( 2.0 * M_PI );
00147    //  return m_parms[norm] * t;
00148 }
00149 
00150 /* virtual */
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 // double Landau::derivByParm ( int i, double x ) const
00161 // {
00162 //   switch ( i ) {
00163 //   case peak :
00164 //     return derivByPeak ( x );
00165 //     break;
00166 
00167 //   case norm :
00168 //     return derivByNorm ( x );
00169 //     break;
00170 
00171 //   case sigma :
00172 //     return derivBySigma ( x );
00173 //     break;
00174 
00175 //   default :
00176 //     assert ( false );
00177 //     break;
00178 //   }
00179 //   return 0.0;
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 } // namespace hippodraw
00209 

Generated for HippoDraw Class Library by doxygen