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 53 of file Stat.h.

Member Function Documentation

erf()

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

static

Definition at line 275 of file flatToGaussian.cc.

                             {
 
// Math:
//
// For any given x we can "quickly" find t0 = erfQ (x) = erf(x) + epsilon.
//
// Then we can find x1 = inverseErf (t0) = inverseErf (erf(x)+epsilon)
//
// Expanding in the small epsion, 
// 
//  x1 = x + epsilon * [deriv(inverseErf(u),u) at u = t0] + O(epsilon**2)
//
// so epsilon is (x1-x) / [deriv(inverseErf(u),u) at u = t0] + O(epsilon**2)
//
// But the derivative of an inverse function is one over the derivative of the
// function, so 
// epsilon  = (x1-x) * [deriv(erf(v),v) at v = x] + O(epsilon**2)
//
// And the definition of the erf integral makes that derivative trivial.
// Ultimately,
//
// erf(x) = erfQ(x) - (inverseErf(erfQ(x))-x) * std::exp(-x**2) * 2/std::sqrt(CLHEP::pi)
//
 
  double t0 = erfQ(x);
  double deriv = std::exp(-x*x) * (2.0 / std::sqrt(CLHEP::pi));
 
  return t0 - (inverseErf (t0) - x) * deriv;
 
}

erfQ()

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

static

Definition at line 25 of file erfQ.cc.

                              {
//
// erfQ is accurate to 7 places.
// See Numerical Recipes P 221
//
 
  double t, z, erfc; 
 
  z = std::abs(x);
  t = 1.0/(1.0+.5*z);
 
  erfc= t*std::exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+t*(0.09678418+
    t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+
    t*(-0.82215223+t*0.17087277 ))) ))) )));
 
  // (The derivation of this formula should be obvious.)
 
  if ( x < 0 ) erfc = 2.0 - erfc;
 
  return 1 - erfc;
 
}

Referenced by erf().

flatToGaussian()

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

static

Definition at line 92 of file flatToGaussian.cc.

                                        {
 
  double sign = +1.0;   // We always compute a negative number of 
                // sigmas.  For r > 0 we will multiply by
                // sign = -1 to return a positive number.
#ifdef Trace1
std::cout << "r = " << r << "\n";
#endif
 
  if ( r > .5 ) {
    r = 1-r;
    sign = -1.0;
#ifdef Trace1
std::cout << "r = " << r << "sign negative \n";
#endif
  } else if ( r == .5 ) {
    return 0.0;
  }  
 
  // Find a pointer to the proper table entries, along with the fraction 
  // of the way in the relevant bin dx and the bin size h.
  
  // Optimize for the case of table 4 by testing for that first.  
  // By removing that case from the for loop below, we save not only
  // several table lookups, but also several index calculations that
  // now become (compile-time) constants.
  //
  // Past the case of table 4, we need not be as concerned about speed since
  // this will happen only .1% of the time.
 
  const double* tptr = 0;
  double  dx = 0;
  double  h = 0;
 
  // The following big if block will locate tptr, h, and dx from whichever
  // table is applicable:
 
  int index;
 
  if ( r >= Table4step ) {
 
    index = int((Table4size<<1) * r);   // 1 to Table4size-1 
    if (index <= 0) index = 1;          // in case of rounding problem
    if (index >= Table4size) index = Table4size-1;
    dx = (Table4size<<1) * r - index;       // fraction of way to next bin
    h = Table4step;
#ifdef Trace2 
std::cout << "index = " << index << " dx = " << dx << " h = " << h << "\n";
#endif
    index = (index<<1) + (Table4offset-2);  
    // at r = table4step+eps, index refers to the start of table 4 
    // and at r = .5 - eps, index refers to the next-to-last entry:
    // remember, the table has two numbers per actual entry.
#ifdef Trace2 
std::cout << "offset index = " << index << "\n";
#endif
 
    tptr = &gaussTables [index];
    
  } else if (r < Tsteps[0])  {
 
    // If r is so small none of the tables apply, use the asymptotic formula.
    return (sign * transformSmall (r));
 
  } else {
    
    for ( int tableN = 3; tableN >= 0; tableN-- ) {
      if ( r < Tsteps[tableN] ) {continue;}     // can't happen when tableN==0
#ifdef Trace2 
std::cout << "Using table " << tableN << "\n";
#endif
      double step = Tsteps[tableN];
      index = int(r/step);          // 1 to TableNsize-1 
        // The following two tests should probably never be true, but
        // roundoff might cause index to be outside its proper range.
        // In such a case, the interpolation still makes sense, but we
        // need to take care that tptr points to valid data in the right table.
      if (index == 0) index = 1;            
      if (index >= Tsizes[tableN]) index = Tsizes[tableN] - 1;
      dx =  r/step - index;             // fraction of way to next bin
      h  =  step;
#ifdef Trace2 
std::cout << "index = " << index << " dx = " << dx << " h = " << h << "\n";
#endif
      index = (index<<1) + Toffsets[tableN] - 2;
      tptr = &gaussTables [index];
      break;
    } // end of the for loop which finds tptr, dx and h in one of the tables
 
    // The code can only get to here by a break statement, having set dx etc.
    // It not get to here without going through one of the breaks,
    // because before the for loop we test for the case of r < Tsteps[0].
 
  } // End of the big if block.
 
  // At the end of this if block, we have tptr, dx and h -- and if r is less 
  // than the smallest step, we have already returned the proper answer.  
 
  double  y0 = *tptr++;
  double  d0 = *tptr++;
  double  y1 = *tptr++;
  double  d1 = *tptr;
 
#ifdef Trace3
std::cout << "y0: " << y0 << " y1: " << y1 << " d0: " << d0 << " d1: " << d1 << "\n";
#endif
 
  double  x2 = dx * dx;
  double  oneMinusX = 1 - dx;
  double  oneMinusX2 = oneMinusX * oneMinusX;
 
  double  f0 = (2. * dx + 1.) * oneMinusX2;
  double  f1 = (3. - 2. * dx) * x2;
  double  g0 =  h * dx * oneMinusX2;
  double  g1 =  - h * oneMinusX * x2;
 
#ifdef Trace3
std::cout << "f0: " << f0 << " f1: " << f1 << " g0: " << g0 << " g1: " << g1 << "\n";
#endif
 
  double answer = f0 * y0 + f1 * y1 + g0 * d0 + g1 * d1;
 
#ifdef Trace1
std::cout << "variate is: " << sign*answer << "\n";
#endif
 
  return sign * answer;
 
} // flatToGaussian

Referenced by inverseErf(), and CLHEP::RandGaussT::operator()().

gammln()

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

static

Definition at line 20 of file gammln.cc.

                                {
 
// Returns the value ln(Gamma(xx) for xx > 0.  Full accuracy is obtained for
// xx > 1. For 0 < xx < 1. the reflection formula (6.1.4) can be used first.
// (Adapted from Numerical Recipes in C.  Relative to that routine, this 
// subtracts one from x at the very start, and in exchange does not have to 
// divide ser by x at the end.  The results are formally equal, and practically
// indistinguishable.)
 
  static const double cof[6] = {76.18009172947146,-86.50532032941677,
                             24.01409824083091, -1.231739572450155,
                             0.1208650973866179e-2, -0.5395239384953e-5};
  int j;
  double x = xx - 1.0;
  double tmp = x + 5.5;
  tmp -= (x + 0.5) * std::log(tmp);
  double ser = 1.000000000190015;
 
  for ( j = 0; j <= 5; j++ ) {
    x += 1.0;
    ser += cof[j]/x;
  }
  return -tmp + std::log(2.5066282746310005*ser);
}

inverseErf()

double CLHEP::HepStat::inverseErf ( double  t )

static

Definition at line 266 of file flatToGaussian.cc.

                                    {
 
  // This uses erf(x) = 2*gaussCDF(std::sqrt(2)*x) - 1
 
  return std::sqrt(0.5) * flatToGaussian(.5*(t+1));
 
}

Referenced by erf().

---

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

• Stat.h
• erfQ.cc
• flatToGaussian.cc
• gammln.cc