Geant4  10.02.p01
flatToGaussian.cc
Go to the documentation of this file.
1 // $Id:$
2 // -*- C++ -*-
3 //
4 // -----------------------------------------------------------------------
5 // HEP Random
6 // --- flatToGaussian ---
7 // class implementation file
8 // -----------------------------------------------------------------------
9 
10 // Contains the methods that depend on the 30K-footpring gaussTables.cdat.
11 //
12 // flatToGaussian (double x)
13 // inverseErf (double x)
14 // erf (double x)
15 
16 // =======================================================================
17 // M Fischler - Created 1/25/00.
18 //
19 // =======================================================================
20 
21 #include "CLHEP/Random/Stat.h"
22 #include "CLHEP/Units/PhysicalConstants.h"
23 #include <iostream>
24 #include <cmath>
25 
26 namespace CLHEP {
27 
28 double transformSmall (double r);
29 
30 //
31 // Table of errInts, for use with transform(r) and quickTransform(r)
32 //
33 
34 #ifdef Traces
35 #define Trace1
36 #define Trace2
37 #define Trace3
38 #endif
39 
40 // Since all these are this is static to this compilation unit only, the
41 // info is establised a priori and not at each invocation.
42 
43 // The main data is of course the gaussTables table; the rest is all
44 // bookkeeping to know what the tables mean.
45 
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)
52 
53 static const int Tsizes[5] = { Table0size,
54  Table1size,
55  Table2size,
56  Table3size,
57  Table4size };
58 
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)
64 
65 static const double Tsteps[5] = { Table0step,
66  Table1step,
67  Table2step,
68  Table3step,
69  Table4step };
70 
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))
76 
77 static const int Toffsets[5] = { Table0offset,
81  Table4offset };
82 
83  // Here comes the big (30K bytes) table, kept in a file ---
84 
85 static const double gaussTables [2*TableSize] = {
86 #include "CLHEP/Random/gaussTables.cdat"
87 };
88 
89 double HepStat::flatToGaussian (double r) {
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
218 
219 double transformSmall (double r) {
220 
221  // Solve for -v in the asymtotic formula
222  //
223  // errInt (-v) = exp(-v*v/2) 1 1*3 1*3*5
224  // ------------ * (1 - ---- + ---- - ----- + ... )
225  // v*sqrt(2*pi) v**2 v**4 v**6
226 
227  // The value of r (=errInt(-v)) supplied is going to less than 2.0E-13,
228  // which is such that v < -7.25. Since the value of r is meaningful only
229  // to an absolute error of 1E-16 (double precision accuracy for a number
230  // which on the high side could be of the form 1-epsilon), computing
231  // v to more than 3-4 digits of accuracy is suspect; however, to ensure
232  // smoothness with the table generator (which uses quite a few terms) we
233  // also use terms up to 1*3*5* ... *13/v**14, and insist on accuracy of
234  // solution at the level of 1.0e-7.
235 
236  // This routine is called less than one time in a trillion firings, so
237  // speed is of no concern. As a matter of technique, we terminate the
238  // iterations in case they would be infinite, but this should not happen.
239 
240  double eps = 1.0e-7;
241  double guess = 7.5;
242  double v;
243 
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;
254  guess = v;
255  }
256 
257  return -v;
258 
259 } // transformSmall()
260 
261 double HepStat::inverseErf (double t) {
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 }
268 
269 double HepStat::erf (double x) {
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 }
299 
300 
301 } // namespace CLHEP
302 
#define Table2step
static const G4double d1
#define Table2offset
#define Table0step
#define Table1step
static const double Tsteps[5]
#define Table2size
static const G4double eps
#define Table1size
#define TableSize
#define Table3step
#define Table3size
static const double twopi
Definition: G4SIunits.hh:75
#define Table4step
#define Table0offset
static const int Toffsets[5]
static const G4double f1
static const double gaussTables[2 *TableSize]
#define Table0size
static const int Tsizes[5]
static const double pi
Definition: G4SIunits.hh:74
#define Table4size
#define Table1offset
const G4double x[NPOINTSGL]
double transformSmall(double r)
#define Table3offset
#define Table4offset
G4int sign(const T t)
A simple sign function that allows us to port fortran code to c++ more easily.