CLHEP VERSION Reference Documentation
   
CLHEP Home Page     CLHEP Documentation     CLHEP Bug Reports

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