Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
CLHEP::HepStat Class Reference

#include <Stat.h>

Static Public Member Functions

static double flatToGaussian (double r)
 
static double inverseErf (double t)
 
static double erf (double x)
 
static double erfQ (double x)
 
static double gammln (double x)
 

Detailed Description

Author

Definition at line 51 of file Stat.h.

Member Function Documentation

double CLHEP::HepStat::erf ( double  x)
static

Definition at line 269 of file flatToGaussian.cc.

269  {
270 
271 // Math:
272 //
273 // For any given x we can "quickly" find t0 = erfQ (x) = erf(x) + epsilon.
274 //
275 // Then we can find x1 = inverseErf (t0) = inverseErf (erf(x)+epsilon)
276 //
277 // Expanding in the small epsion,
278 //
279 // x1 = x + epsilon * [deriv(inverseErf(u),u) at u = t0] + O(epsilon**2)
280 //
281 // so epsilon is (x1-x) / [deriv(inverseErf(u),u) at u = t0] + O(epsilon**2)
282 //
283 // But the derivative of an inverse function is one over the derivative of the
284 // function, so
285 // epsilon = (x1-x) * [deriv(erf(v),v) at v = x] + O(epsilon**2)
286 //
287 // And the definition of the erf integral makes that derivative trivial.
288 // Ultimately,
289 //
290 // erf(x) = erfQ(x) - (inverseErf(erfQ(x))-x) * exp(-x**2) * 2/sqrt(CLHEP::pi)
291 //
292 
293  double t0 = erfQ(x);
294  double deriv = std::exp(-x*x) * (2.0 / std::sqrt(CLHEP::pi));
295 
296  return t0 - (inverseErf (t0) - x) * deriv;
297 
298 }
static double inverseErf(double t)
static double erfQ(double x)
Definition: erfQ.cc:24
static constexpr double pi
Definition: SystemOfUnits.h:54

Here is the call graph for this function:

Here is the caller graph for this function:

double CLHEP::HepStat::erfQ ( double  x)
static

Definition at line 24 of file erfQ.cc.

24  {
25 //
26 // erfQ is accurate to 7 places.
27 // See Numerical Recipes P 221
28 //
29 
30  double t, z, erfc;
31 
32  z = std::abs(x);
33  t = 1.0/(1.0+.5*z);
34 
35  erfc= t*std::exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+t*(0.09678418+
36  t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+
37  t*(-0.82215223+t*0.17087277 ))) ))) )));
38 
39  // (The derivation of this formula should be obvious.)
40 
41  if ( x < 0 ) erfc = 2.0 - erfc;
42 
43  return 1 - erfc;
44 
45 }

Here is the caller graph for this function:

double CLHEP::HepStat::flatToGaussian ( double  r)
static

Definition at line 89 of file flatToGaussian.cc.

89  {
90 
91  double sign = +1.0; // We always compute a negative number of
92  // sigmas. For r > 0 we will multiply by
93  // sign = -1 to return a positive number.
94 #ifdef Trace1
95 std::cout << "r = " << r << "\n";
96 #endif
97 
98  if ( r > .5 ) {
99  r = 1-r;
100  sign = -1.0;
101 #ifdef Trace1
102 std::cout << "r = " << r << "sign negative \n";
103 #endif
104  } else if ( r == .5 ) {
105  return 0.0;
106  }
107 
108  // Find a pointer to the proper table entries, along with the fraction
109  // of the way in the relevant bin dx and the bin size h.
110 
111  // Optimize for the case of table 4 by testing for that first.
112  // By removing that case from the for loop below, we save not only
113  // several table lookups, but also several index calculations that
114  // now become (compile-time) constants.
115  //
116  // Past the case of table 4, we need not be as concerned about speed since
117  // this will happen only .1% of the time.
118 
119  const double* tptr = 0;
120  double dx = 0;
121  double h = 0;
122 
123  // The following big if block will locate tptr, h, and dx from whichever
124  // table is applicable:
125 
126  int index;
127 
128  if ( r >= Table4step ) {
129 
130  index = int((Table4size<<1) * r); // 1 to Table4size-1
131  if (index <= 0) index = 1; // in case of rounding problem
132  if (index >= Table4size) index = Table4size-1;
133  dx = (Table4size<<1) * r - index; // fraction of way to next bin
134  h = Table4step;
135 #ifdef Trace2
136 std::cout << "index = " << index << " dx = " << dx << " h = " << h << "\n";
137 #endif
138  index = (index<<1) + (Table4offset-2);
139  // at r = table4step+eps, index refers to the start of table 4
140  // and at r = .5 - eps, index refers to the next-to-last entry:
141  // remember, the table has two numbers per actual entry.
142 #ifdef Trace2
143 std::cout << "offset index = " << index << "\n";
144 #endif
145 
146  tptr = &gaussTables [index];
147 
148  } else if (r < Tsteps[0]) {
149 
150  // If r is so small none of the tables apply, use the asymptotic formula.
151  return (sign * transformSmall (r));
152 
153  } else {
154 
155  for ( int tableN = 3; tableN >= 0; tableN-- ) {
156  if ( r < Tsteps[tableN] ) {continue;} // can't happen when tableN==0
157 #ifdef Trace2
158 std::cout << "Using table " << tableN << "\n";
159 #endif
160  double step = Tsteps[tableN];
161  index = int(r/step); // 1 to TableNsize-1
162  // The following two tests should probably never be true, but
163  // roundoff might cause index to be outside its proper range.
164  // In such a case, the interpolation still makes sense, but we
165  // need to take care that tptr points to valid data in the right table.
166  if (index == 0) index = 1;
167  if (index >= Tsizes[tableN]) index = Tsizes[tableN] - 1;
168  dx = r/step - index; // fraction of way to next bin
169  h = step;
170 #ifdef Trace2
171 std::cout << "index = " << index << " dx = " << dx << " h = " << h << "\n";
172 #endif
173  index = (index<<1) + Toffsets[tableN] - 2;
174  tptr = &gaussTables [index];
175  break;
176  } // end of the for loop which finds tptr, dx and h in one of the tables
177 
178  // The code can only get to here by a break statement, having set dx etc.
179  // It not get to here without going through one of the breaks,
180  // because before the for loop we test for the case of r < Tsteps[0].
181 
182  } // End of the big if block.
183 
184  // At the end of this if block, we have tptr, dx and h -- and if r is less
185  // than the smallest step, we have already returned the proper answer.
186 
187  double y0 = *tptr++;
188  double d0 = *tptr++;
189  double y1 = *tptr++;
190  double d1 = *tptr;
191 
192 #ifdef Trace3
193 std::cout << "y0: " << y0 << " y1: " << y1 << " d0: " << d0 << " d1: " << d1 << "\n";
194 #endif
195 
196  double x2 = dx * dx;
197  double oneMinusX = 1 - dx;
198  double oneMinusX2 = oneMinusX * oneMinusX;
199 
200  double f0 = (2. * dx + 1.) * oneMinusX2;
201  double f1 = (3. - 2. * dx) * x2;
202  double g0 = h * dx * oneMinusX2;
203  double g1 = - h * oneMinusX * x2;
204 
205 #ifdef Trace3
206 std::cout << "f0: " << f0 << " f1: " << f1 << " g0: " << g0 << " g1: " << g1 << "\n";
207 #endif
208 
209  double answer = f0 * y0 + f1 * y1 + g0 * d0 + g1 * d1;
210 
211 #ifdef Trace1
212 std::cout << "variate is: " << sign*answer << "\n";
213 #endif
214 
215  return sign * answer;
216 
217 } // flatToGaussian
static const double Tsteps[5]
#define Table4step
typedef int(XMLCALL *XML_NotStandaloneHandler)(void *userData)
static const int Toffsets[5]
static const G4double d1
static const double gaussTables[2 *TableSize]
static const int Tsizes[5]
#define Table4size
double transformSmall(double r)
#define Table4offset
G4int sign(const T t)

Here is the call graph for this function:

Here is the caller graph for this function:

double CLHEP::HepStat::gammln ( double  x)
static

Definition at line 19 of file gammln.cc.

19  {
20 
21 // Returns the value ln(Gamma(xx) for xx > 0. Full accuracy is obtained for
22 // xx > 1. For 0 < xx < 1. the reflection formula (6.1.4) can be used first.
23 // (Adapted from Numerical Recipes in C. Relative to that routine, this
24 // subtracts one from x at the very start, and in exchange does not have to
25 // divide ser by x at the end. The results are formally equal, and practically
26 // indistinguishable.)
27 
28  static const double cof[6] = {76.18009172947146,-86.50532032941677,
29  24.01409824083091, -1.231739572450155,
30  0.1208650973866179e-2, -0.5395239384953e-5};
31  int j;
32  double x = xx - 1.0;
33  double tmp = x + 5.5;
34  tmp -= (x + 0.5) * std::log(tmp);
35  double ser = 1.000000000190015;
36 
37  for ( j = 0; j <= 5; j++ ) {
38  x += 1.0;
39  ser += cof[j]/x;
40  }
41  return -tmp + std::log(2.5066282746310005*ser);
42 }
double CLHEP::HepStat::inverseErf ( double  t)
static

Definition at line 261 of file flatToGaussian.cc.

261  {
262 
263  // This uses erf(x) = 2*gaussCDF(sqrt(2)*x) - 1
264 
265  return std::sqrt(0.5) * flatToGaussian(.5*(t+1));
266 
267 }
static double flatToGaussian(double r)

Here is the call graph for this function:

Here is the caller graph for this function:


The documentation for this class was generated from the following files: