21 #include "CLHEP/Random/Stat.h"
22 #include "CLHEP/Units/PhysicalConstants.h"
46 #define Table0size 200
47 #define Table1size 250
48 #define Table2size 200
49 #define Table3size 250
50 #define Table4size 1000
51 #define TableSize (Table0size+Table1size+Table2size+Table3size+Table4size)
59 #define Table0step (2.0E-13)
60 #define Table1step (4.0E-11)
61 #define Table2step (1.0E-8)
62 #define Table3step (2.0E-6)
63 #define Table4step (5.0E-4)
71 #define Table0offset 0
72 #define Table1offset (2*(Table0size))
73 #define Table2offset (2*(Table0size + Table1size))
74 #define Table3offset (2*(Table0size + Table1size + Table2size))
75 #define Table4offset (2*(Table0size + Table1size + Table2size + Table3size))
86 #include "CLHEP/Random/gaussTables.cdat"
89 double HepStat::flatToGaussian (
double r) {
95 std::cout <<
"r = " << r <<
"\n";
102 std::cout <<
"r = " << r <<
"sign negative \n";
104 }
else if ( r == .5 ) {
119 const double* tptr = 0;
131 if (index <= 0) index = 1;
136 std::cout <<
"index = " << index <<
" dx = " << dx <<
" h = " << h <<
"\n";
143 std::cout <<
"offset index = " << index <<
"\n";
146 tptr = &gaussTables [index];
148 }
else if (r < Tsteps[0]) {
155 for (
int tableN = 3; tableN >= 0; tableN-- ) {
156 if ( r < Tsteps[tableN] ) {
continue;}
158 std::cout <<
"Using table " << tableN <<
"\n";
160 double step = Tsteps[tableN];
166 if (index == 0) index = 1;
167 if (index >= Tsizes[tableN]) index = Tsizes[tableN] - 1;
171 std::cout <<
"index = " << index <<
" dx = " << dx <<
" h = " << h <<
"\n";
173 index = (index<<1) + Toffsets[tableN] - 2;
174 tptr = &gaussTables [index];
193 std::cout <<
"y0: " << y0 <<
" y1: " << y1 <<
" d0: " << d0 <<
" d1: " << d1 <<
"\n";
197 double oneMinusX = 1 - dx;
198 double oneMinusX2 = oneMinusX * oneMinusX;
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;
206 std::cout <<
"f0: " << f0 <<
" f1: " << f1 <<
" g0: " << g0 <<
" g1: " << g1 <<
"\n";
209 double answer = f0 * y0 + f1 * y1 + g0 * d0 + g1 *
d1;
212 std::cout <<
"variate is: " << sign*answer <<
"\n";
215 return sign * answer;
244 for (
int i = 1; i < 50; i++ ) {
245 double vn2 = 1.0/(guess*guess);
246 double s1 = -13*11*9*7*5*3 * vn2*vn2*vn2*vn2*vn2*vn2*vn2;
247 s1 += 11*9*7*5*3 * vn2*vn2*vn2*vn2*vn2*vn2;
248 s1 += -9*7*5*3 * vn2*vn2*vn2*vn2*vn2;
249 s1 += 7*5*3 * vn2*vn2*vn2*vn2;
250 s1 += -5*3 * vn2*vn2*vn2;
251 s1 += 3 * vn2*vn2 - vn2 + 1.0;
252 v = std::sqrt ( 2.0 * std::log ( s1 / (r*guess*std::sqrt(
CLHEP::twopi)) ) );
253 if ( std::abs(v-guess) < eps )
break;
261 double HepStat::inverseErf (
double t) {
265 return std::sqrt(0.5) * flatToGaussian(.5*(t+1));
269 double HepStat::erf (
double x) {
294 double deriv = std::exp(-x*x) * (2.0 / std::sqrt(
CLHEP::pi));
296 return t0 - (inverseErf (t0) -
x) * deriv;
static const double Tsteps[5]
static const G4double eps
static const double twopi
static const int Toffsets[5]
static const double gaussTables[2 *TableSize]
static const int Tsizes[5]
const G4double x[NPOINTSGL]
double transformSmall(double r)
G4int sign(const T t)
A simple sign function that allows us to port fortran code to c++ more easily.