Erfc.cxx

Go to the documentation of this file.
00001 
00013 #ifdef _MSC_VER
00014 #include "msdevstudio/MSconfig.h"
00015 #endif
00016 
00017 #include "Erfc.h"
00018 #include "FunctionHelper.h"
00019 
00020 #include <cmath>
00021 
00022 #include <cassert>
00023 
00024 #ifdef ITERATOR_MEMBER_DEFECT
00025 using namespace std;
00026 #else
00027 using std::exp;
00028 using std::vector;
00029 #endif
00030 
00031 namespace hippodraw {
00032 
00033 Erfc::Erfc ( )
00034 {
00035   initialize ();
00036 }
00037 
00038 Erfc::Erfc(double m, double s, double n)
00039 {
00040   initialize();
00041   
00042   m_parms[MEAN]  = m;
00043   m_parms[SIGMA] = s;
00044   m_parms[NORM]  = n;  
00045 }
00046 
00047 void Erfc::initialize ()
00048 {
00049   m_name = "Erfc";
00050 
00051   m_parm_names.push_back ( "Mean" );
00052   m_parm_names.push_back ( "Sigma" );
00053   m_parm_names.push_back ( "Normalization" );
00054 
00055   resize ();
00056 }
00057 
00058 FunctionBase * Erfc::clone () const
00059 {
00060   return new Erfc ( *this );
00061 }
00062 
00065 double Erfc::operator () (double value) const
00066 {
00067   double result = 1; // The return value
00068 
00069   double x = calcRed(value); //reduced variable
00070 
00071   result = calcErfc(x);  //this method is a copy of ROOT method in TMath.cxx
00072   
00073   result *= m_parms[NORM];
00074 
00075   return result;
00076 }
00077 
00078 
00079 /* virtual */
00080 void  Erfc::initialParameters ( const FunctionHelper * helper )
00081 {
00082    m_parms[MEAN] = helper->meanCoord();
00083    m_parms[SIGMA] = 1;
00084    m_parms[NORM] = helper->maxValue () / 2;
00085 }
00086 
00087 
00088 double Erfc::derivByParm ( int i, double value ) const
00089 {
00090    double result = 0;
00091 
00092   switch ( i ) {
00093   case MEAN :
00094     result = -1.0/m_parms[SIGMA] * m_parms[NORM] * derivByRed ( calcRed(value) );
00095     return result; 
00096     break;
00097 
00098   case SIGMA :
00099     result = - calcRed(value)/m_parms[SIGMA] * m_parms[NORM] * derivByRed ( calcRed(value) );
00100     return result; 
00101     break;
00102 
00103   case NORM :
00104    result = operator () ( value ) / m_parms[NORM];
00105    return result; 
00106    break;
00107 
00108   default :
00109     assert ( false );
00110     break;
00111   }
00112 
00113   return 0.0; 
00114 }
00115 
00116 
00117 //Derivative wrt reduced is gaussian!
00118 double Erfc::derivByRed( const double x ) const
00119 { 
00120   //I don't know if this is platform independant: std::M_2_SQRTPI from math.h
00121   static const double m_2_sqrtpi = 1.12837916709551257390; // 2/sqrt(pi)
00122 
00123   double result = - m_2_sqrtpi; // minus sign comes from the fact that x is lower
00124                                 // bound of Erfc definition
00125   result *= exp(-x*x);
00126   return result;
00127 }
00128 
00129 
00130 //this is the copy from ROOT TMath.cxx
00131 double Erfc::calcErfc(double x) const
00132 {
00133   // Compute the complementary error function erfc(x).
00134   // Erfc(x) = (2/sqrt(pi)) Integral(exp(-t^2))dt between x and infinity
00135   //
00136   //--- Nve 14-nov-1998 UU-SAP Utrecht
00137   
00138   // The parameters of the Chebyshev fit
00139   const double 
00140     a1 = -1.26551223,   a2 = 1.00002368,
00141     a3 =  0.37409196,   a4 = 0.09678418,
00142     a5 = -0.18628806,   a6 = 0.27886807,
00143     a7 = -1.13520398,   a8 = 1.48851587,
00144     a9 = -0.82215223,  a10 = 0.17087277;
00145 
00146    double v = 1; // The return value
00147 
00148    double z = x>=0 ? x : -x; //absolute value
00149 
00150    if (z <= 0) return v; // erfc(0)=1
00151 
00152    double t = 1/(1+0.5*z);
00153 
00154    v = t*exp((-z*z) +a1+t*(a2+t*(a3+t*(a4+t*(a5+t*(a6+t*(a7+t*(a8+t*(a9+t*a10)))))))));
00155 
00156    if (x < 0) v = 2-v; // erfc(-x)=2-erfc(x)
00157 
00158    return v;
00159 }
00160 
00161 } // namespace hippodraw

Generated for HippoDraw Class Library by doxygen