Erfc.cxx
Go to the documentation of this file.
1 
13 #ifdef _MSC_VER
14 #include "msdevstudio/MSconfig.h"
15 #endif
16 
17 #include "Erfc.h"
18 #include "FunctionHelper.h"
19 
20 #include <cmath>
21 
22 #include <cassert>
23 
24 #ifdef ITERATOR_MEMBER_DEFECT
25 using namespace std;
26 #else
27 using std::exp;
28 using std::vector;
29 #endif
30 
31 namespace hippodraw {
32 
33 Erfc::Erfc ( )
34 {
35  initialize ();
36 }
37 
38 Erfc::Erfc(double m, double s, double n)
39 {
40  initialize();
41 
42  m_parms[MEAN] = m;
43  m_parms[SIGMA] = s;
44  m_parms[NORM] = n;
45 }
46 
47 void Erfc::initialize ()
48 {
49  m_name = "Erfc";
50 
51  m_parm_names.push_back ( "Mean" );
52  m_parm_names.push_back ( "Sigma" );
53  m_parm_names.push_back ( "Normalization" );
54 
55  resize ();
56 }
57 
59 {
60  return new Erfc ( *this );
61 }
62 
65 double Erfc::operator () (double value) const
66 {
67  double result = 1; // The return value
68 
69  double x = calcRed(value); //reduced variable
70 
71  result = calcErfc(x); //this method is a copy of ROOT method in TMath.cxx
72 
73  result *= m_parms[NORM];
74 
75  return result;
76 }
77 
78 
79 /* virtual */
80 void Erfc::initialParameters ( const FunctionHelper * helper )
81 {
82  m_parms[MEAN] = helper->meanCoord();
83  m_parms[SIGMA] = 1;
84  m_parms[NORM] = helper->maxValue () / 2;
85 }
86 
87 
88 double Erfc::derivByParm ( int i, double value ) const
89 {
90  double result = 0;
91 
92  switch ( i ) {
93  case MEAN :
94  result = -1.0/m_parms[SIGMA] * m_parms[NORM] * derivByRed ( calcRed(value) );
95  return result;
96  break;
97 
98  case SIGMA :
99  result = - calcRed(value)/m_parms[SIGMA] * m_parms[NORM] * derivByRed ( calcRed(value) );
100  return result;
101  break;
102 
103  case NORM :
104  result = operator () ( value ) / m_parms[NORM];
105  return result;
106  break;
107 
108  default :
109  assert ( false );
110  break;
111  }
112 
113  return 0.0;
114 }
115 
116 
117 //Derivative wrt reduced is gaussian!
118 double Erfc::derivByRed( const double x ) const
119 {
120  //I don't know if this is platform independant: std::M_2_SQRTPI from math.h
121  static const double m_2_sqrtpi = 1.12837916709551257390; // 2/sqrt(pi)
122 
123  double result = - m_2_sqrtpi; // minus sign comes from the fact that x is lower
124  // bound of Erfc definition
125  result *= exp(-x*x);
126  return result;
127 }
128 
129 
130 //this is the copy from ROOT TMath.cxx
131 double Erfc::calcErfc(double x) const
132 {
133  // Compute the complementary error function erfc(x).
134  // Erfc(x) = (2/sqrt(pi)) Integral(exp(-t^2))dt between x and infinity
135  //
136  //--- Nve 14-nov-1998 UU-SAP Utrecht
137 
138  // The parameters of the Chebyshev fit
139  const double
140  a1 = -1.26551223, a2 = 1.00002368,
141  a3 = 0.37409196, a4 = 0.09678418,
142  a5 = -0.18628806, a6 = 0.27886807,
143  a7 = -1.13520398, a8 = 1.48851587,
144  a9 = -0.82215223, a10 = 0.17087277;
145 
146  double v = 1; // The return value
147 
148  double z = x>=0 ? x : -x; //absolute value
149 
150  if (z <= 0) return v; // erfc(0)=1
151 
152  double t = 1/(1+0.5*z);
153 
154  v = t*exp((-z*z) +a1+t*(a2+t*(a3+t*(a4+t*(a5+t*(a6+t*(a7+t*(a8+t*(a9+t*a10)))))))));
155 
156  if (x < 0) v = 2-v; // erfc(-x)=2-erfc(x)
157 
158  return v;
159 }
160 
161 } // namespace hippodraw

Generated for HippoDraw Class Library by doxygen