ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
flatToGaussian.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file flatToGaussian.cc
1 // -*- C++ -*-
2 //
3 // -----------------------------------------------------------------------
4 // HEP Random
5 // --- flatToGaussian ---
6 // class implementation file
7 // -----------------------------------------------------------------------
8 
9 // Contains the methods that depend on the 30K-footpring gaussTables.cdat.
10 //
11 // flatToGaussian (double x)
12 // inverseErf (double x)
13 // erf (double x)
14 
15 // =======================================================================
16 // M Fischler - Created 1/25/00.
17 //
18 // =======================================================================
19 
20 #include "CLHEP/Random/Stat.h"
22 #include <iostream>
23 #include <cmath>
24 
25 namespace CLHEP {
26 
27 double transformSmall (double r);
28 
29 //
30 // Table of errInts, for use with transform(r) and quickTransform(r)
31 //
32 
33 #ifdef Traces
34 #define Trace1
35 #define Trace2
36 #define Trace3
37 #endif
38 
39 // Since all these are this is static to this compilation unit only, the
40 // info is establised a priori and not at each invocation.
41 
42 // The main data is of course the gaussTables table; the rest is all
43 // bookkeeping to know what the tables mean.
44 
45 #define Table0size 200
46 #define Table1size 250
47 #define Table2size 200
48 #define Table3size 250
49 #define Table4size 1000
50 #define TableSize (Table0size+Table1size+Table2size+Table3size+Table4size)
51 
52 static const int Tsizes[5] = { Table0size,
53  Table1size,
54  Table2size,
55  Table3size,
56  Table4size };
57 
58 #define Table0step (2.0E-13)
59 #define Table1step (4.0E-11)
60 #define Table2step (1.0E-8)
61 #define Table3step (2.0E-6)
62 #define Table4step (5.0E-4)
63 
64 static const double Tsteps[5] = { Table0step,
65  Table1step,
66  Table2step,
67  Table3step,
68  Table4step };
69 
70 #define Table0offset 0
71 #define Table1offset (2*(Table0size))
72 #define Table2offset (2*(Table0size + Table1size))
73 #define Table3offset (2*(Table0size + Table1size + Table2size))
74 #define Table4offset (2*(Table0size + Table1size + Table2size + Table3size))
75 
76 static const int Toffsets[5] = { Table0offset,
80  Table4offset };
81 
82  // Here comes the big (30K bytes) table, kept in a file ---
83 
84 static const double gaussTables [2*TableSize] = {
85 #include "CLHEP/Random/gaussTables.cdat"
86 };
87 
88 double HepStat::flatToGaussian (double r) {
89 
90  double sign = +1.0; // We always compute a negative number of
91  // sigmas. For r > 0 we will multiply by
92  // sign = -1 to return a positive number.
93 #ifdef Trace1
94 std::cout << "r = " << r << "\n";
95 #endif
96 
97  if ( r > .5 ) {
98  r = 1-r;
99  sign = -1.0;
100 #ifdef Trace1
101 std::cout << "r = " << r << "sign negative \n";
102 #endif
103  } else if ( r == .5 ) {
104  return 0.0;
105  }
106 
107  // Find a pointer to the proper table entries, along with the fraction
108  // of the way in the relevant bin dx and the bin size h.
109 
110  // Optimize for the case of table 4 by testing for that first.
111  // By removing that case from the for loop below, we save not only
112  // several table lookups, but also several index calculations that
113  // now become (compile-time) constants.
114  //
115  // Past the case of table 4, we need not be as concerned about speed since
116  // this will happen only .1% of the time.
117 
118  const double* tptr = 0;
119  double dx = 0;
120  double h = 0;
121 
122  // The following big if block will locate tptr, h, and dx from whichever
123  // table is applicable:
124 
125  int index;
126 
127  if ( r >= Table4step ) {
128 
129  index = int((Table4size<<1) * r); // 1 to Table4size-1
130  if (index <= 0) index = 1; // in case of rounding problem
131  if (index >= Table4size) index = Table4size-1;
132  dx = (Table4size<<1) * r - index; // fraction of way to next bin
133  h = Table4step;
134 #ifdef Trace2
135 std::cout << "index = " << index << " dx = " << dx << " h = " << h << "\n";
136 #endif
137  index = (index<<1) + (Table4offset-2);
138  // at r = table4step+eps, index refers to the start of table 4
139  // and at r = .5 - eps, index refers to the next-to-last entry:
140  // remember, the table has two numbers per actual entry.
141 #ifdef Trace2
142 std::cout << "offset index = " << index << "\n";
143 #endif
144 
145  tptr = &gaussTables [index];
146 
147  } else if (r < Tsteps[0]) {
148 
149  // If r is so small none of the tables apply, use the asymptotic formula.
150  return (sign * transformSmall (r));
151 
152  } else {
153 
154  for ( int tableN = 3; tableN >= 0; tableN-- ) {
155  if ( r < Tsteps[tableN] ) {continue;} // can't happen when tableN==0
156 #ifdef Trace2
157 std::cout << "Using table " << tableN << "\n";
158 #endif
159  double step = Tsteps[tableN];
160  index = int(r/step); // 1 to TableNsize-1
161  // The following two tests should probably never be true, but
162  // roundoff might cause index to be outside its proper range.
163  // In such a case, the interpolation still makes sense, but we
164  // need to take care that tptr points to valid data in the right table.
165  if (index == 0) index = 1;
166  if (index >= Tsizes[tableN]) index = Tsizes[tableN] - 1;
167  dx = r/step - index; // fraction of way to next bin
168  h = step;
169 #ifdef Trace2
170 std::cout << "index = " << index << " dx = " << dx << " h = " << h << "\n";
171 #endif
172  index = (index<<1) + Toffsets[tableN] - 2;
173  tptr = &gaussTables [index];
174  break;
175  } // end of the for loop which finds tptr, dx and h in one of the tables
176 
177  // The code can only get to here by a break statement, having set dx etc.
178  // It not get to here without going through one of the breaks,
179  // because before the for loop we test for the case of r < Tsteps[0].
180 
181  } // End of the big if block.
182 
183  // At the end of this if block, we have tptr, dx and h -- and if r is less
184  // than the smallest step, we have already returned the proper answer.
185 
186  double y0 = *tptr++;
187  double d0 = *tptr++;
188  double y1 = *tptr++;
189  double d1 = *tptr;
190 
191 #ifdef Trace3
192 std::cout << "y0: " << y0 << " y1: " << y1 << " d0: " << d0 << " d1: " << d1 << "\n";
193 #endif
194 
195  double x2 = dx * dx;
196  double oneMinusX = 1 - dx;
197  double oneMinusX2 = oneMinusX * oneMinusX;
198 
199  double f0 = (2. * dx + 1.) * oneMinusX2;
200  double f1 = (3. - 2. * dx) * x2;
201  double g0 = h * dx * oneMinusX2;
202  double g1 = - h * oneMinusX * x2;
203 
204 #ifdef Trace3
205 std::cout << "f0: " << f0 << " f1: " << f1 << " g0: " << g0 << " g1: " << g1 << "\n";
206 #endif
207 
208  double answer = f0 * y0 + f1 * y1 + g0 * d0 + g1 * d1;
209 
210 #ifdef Trace1
211 std::cout << "variate is: " << sign*answer << "\n";
212 #endif
213 
214  return sign * answer;
215 
216 } // flatToGaussian
217 
218 double transformSmall (double r) {
219 
220  // Solve for -v in the asymtotic formula
221  //
222  // errInt (-v) = exp(-v*v/2) 1 1*3 1*3*5
223  // ------------ * (1 - ---- + ---- - ----- + ... )
224  // v*sqrt(2*pi) v**2 v**4 v**6
225 
226  // The value of r (=errInt(-v)) supplied is going to less than 2.0E-13,
227  // which is such that v < -7.25. Since the value of r is meaningful only
228  // to an absolute error of 1E-16 (double precision accuracy for a number
229  // which on the high side could be of the form 1-epsilon), computing
230  // v to more than 3-4 digits of accuracy is suspect; however, to ensure
231  // smoothness with the table generator (which uses quite a few terms) we
232  // also use terms up to 1*3*5* ... *13/v**14, and insist on accuracy of
233  // solution at the level of 1.0e-7.
234 
235  // This routine is called less than one time in a trillion firings, so
236  // speed is of no concern. As a matter of technique, we terminate the
237  // iterations in case they would be infinite, but this should not happen.
238 
239  double eps = 1.0e-7;
240  double guess = 7.5;
241  double v;
242 
243  for ( int i = 1; i < 50; i++ ) {
244  double vn2 = 1.0/(guess*guess);
245  double s1 = -13*11*9*7*5*3 * vn2*vn2*vn2*vn2*vn2*vn2*vn2;
246  s1 += 11*9*7*5*3 * vn2*vn2*vn2*vn2*vn2*vn2;
247  s1 += -9*7*5*3 * vn2*vn2*vn2*vn2*vn2;
248  s1 += 7*5*3 * vn2*vn2*vn2*vn2;
249  s1 += -5*3 * vn2*vn2*vn2;
250  s1 += 3 * vn2*vn2 - vn2 + 1.0;
251  v = std::sqrt ( 2.0 * std::log ( s1 / (r*guess*std::sqrt(CLHEP::twopi)) ) );
252  if ( std::abs(v-guess) < eps ) break;
253  guess = v;
254  }
255 
256  return -v;
257 
258 } // transformSmall()
259 
260 double HepStat::inverseErf (double t) {
261 
262  // This uses erf(x) = 2*gaussCDF(sqrt(2)*x) - 1
263 
264  return std::sqrt(0.5) * flatToGaussian(.5*(t+1));
265 
266 }
267 
268 double HepStat::erf (double x) {
269 
270 // Math:
271 //
272 // For any given x we can "quickly" find t0 = erfQ (x) = erf(x) + epsilon.
273 //
274 // Then we can find x1 = inverseErf (t0) = inverseErf (erf(x)+epsilon)
275 //
276 // Expanding in the small epsion,
277 //
278 // x1 = x + epsilon * [deriv(inverseErf(u),u) at u = t0] + O(epsilon**2)
279 //
280 // so epsilon is (x1-x) / [deriv(inverseErf(u),u) at u = t0] + O(epsilon**2)
281 //
282 // But the derivative of an inverse function is one over the derivative of the
283 // function, so
284 // epsilon = (x1-x) * [deriv(erf(v),v) at v = x] + O(epsilon**2)
285 //
286 // And the definition of the erf integral makes that derivative trivial.
287 // Ultimately,
288 //
289 // erf(x) = erfQ(x) - (inverseErf(erfQ(x))-x) * exp(-x**2) * 2/sqrt(CLHEP::pi)
290 //
291 
292  double t0 = erfQ(x);
293  double deriv = std::exp(-x*x) * (2.0 / std::sqrt(CLHEP::pi));
294 
295  return t0 - (inverseErf (t0) - x) * deriv;
296 
297 }
298 
299 
300 } // namespace CLHEP
301