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.

References erfQ(), inverseErf(), and test::x.

Referenced by G4QMDMeanField::Cal2BodyQuantities(), and G4FPYSamplingOps::ShiftParameters().

                             {

// 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) * exp(-x**2) * 2/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;

}

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

static

Definition at line 24 of file erfQ.cc.

References z.

Referenced by erf().

                              {
//
// 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;

}

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

static

Definition at line 89 of file flatToGaussian.cc.

References int(), G4INCL::Math::sign(), Table4offset, Table4size, Table4step, and CLHEP::transformSmall().

Referenced by inverseErf().

                                        {

  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:

  register 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

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

static

Definition at line 19 of file gammln.cc.

References test::x.

                                {

// 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 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);
}

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

static

Definition at line 261 of file flatToGaussian.cc.

References flatToGaussian().

Referenced by erf().

                                    {

  // This uses erf(x) = 2*gaussCDF(sqrt(2)*x) - 1

  return std::sqrt(0.5) * flatToGaussian(.5*(t+1));

}

---

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

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