| Geant4 Cross Reference |
1 // -*- C++ -*- 1
2 //
3 // -------------------------------------------
4 // HEP Random
5 // --- flatToGaussian
6 // class implementation f
7 // -------------------------------------------
8
9 // Contains the methods that depend on the 30K
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"
21 #include "CLHEP/Units/PhysicalConstants.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)
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
40 // info is establised a priori and not at each
41
42 // The main data is of course the gaussTables
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+Tab
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 + Table1si
73 #define Table3offset (2*(Table0size + Table1si
74 #define Table4offset (2*(Table0size + Table1si
75
76 static const int Toffsets[5] = { Table0offset,
77 Table1offset,
78 Table2offset,
79 Table3offset,
80 Table4offset };
81
82 // Here comes the big (30K bytes) table, kep
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 n
91 // sigmas. For r > 0 we will multiply
92 // sign = -1 to return a positive numb
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 entrie
108 // of the way in the relevant bin dx and the
109
110 // Optimize for the case of table 4 by testi
111 // By removing that case from the for loop b
112 // several table lookups, but also several i
113 // now become (compile-time) constants.
114 //
115 // Past the case of table 4, we need not be
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 tp
123 // table is applicable:
124
125 int index;
126
127 if ( r >= Table4step ) {
128
129 index = int((Table4size<<1) * r); // 1 to
130 if (index <= 0) index = 1; // in case
131 if (index >= Table4size) index = Table4siz
132 dx = (Table4size<<1) * r - index; // f
133 h = Table4step;
134 #ifdef Trace2
135 std::cout << "index = " << index << " dx = " <
136 #endif
137 index = (index<<1) + (Table4offset-2);
138 // at r = table4step+eps, index refers to th
139 // and at r = .5 - eps, index refers to the
140 // remember, the table has two numbers per a
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 app
150 return (sign * transformSmall (r));
151
152 } else {
153
154 for ( int tableN = 3; tableN >= 0; tableN-
155 if ( r < Tsteps[tableN] ) {continue;}
156 #ifdef Trace2
157 std::cout << "Using table " << tableN << "\n";
158 #endif
159 double step = Tsteps[tableN];
160 index = int(r/step); // 1 to TableN
161 // The following two tests should prob
162 // roundoff might cause index to be ou
163 // In such a case, the interpolation s
164 // need to take care that tptr points
165 if (index == 0) index = 1;
166 if (index >= Tsizes[tableN]) index = Tsi
167 dx = r/step - index; // fraction
168 h = step;
169 #ifdef Trace2
170 std::cout << "index = " << index << " dx = " <
171 #endif
172 index = (index<<1) + Toffsets[tableN] -
173 tptr = &gaussTables [index];
174 break;
175 } // end of the for loop which finds tptr,
176
177 // The code can only get to here by a brea
178 // It not get to here without going throug
179 // because before the for loop we test for
180
181 } // End of the big if block.
182
183 // At the end of this if block, we have tptr
184 // than the smallest step, we have already r
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 <<
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 <<
206 #endif
207
208 double answer = f0 * y0 + f1 * y1 + g0 * d0
209
210 #ifdef Trace1
211 std::cout << "variate is: " << sign*answer <<
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
223 // ------------ * (1 - ---- + ---- - -
224 // v*sqrt(2*pi) v**2 v**4 v
225
226 // The value of r (=errInt(-v)) supplied is
227 // which is such that v < -7.25. Since the
228 // to an absolute error of 1E-16 (double pre
229 // which on the high side could be of the fo
230 // v to more than 3-4 digits of accuracy is
231 // smoothness with the table generator (whic
232 // also use terms up to 1*3*5* ... *13/v**14
233 // solution at the level of 1.0e-7.
234
235 // This routine is called less than one time
236 // speed is of no concern. As a matter of t
237 // iterations in case they would be infinite
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*v
246 s1 += 11*9*7*5*3 * vn2*vn2*vn2*
247 s1 += -9*7*5*3 * vn2*vn2*vn2*
248 s1 += 7*5*3 * vn2*vn2*vn2*
249 s1 += -5*3 * vn2*vn2*vn2;
250 s1 += 3 * vn2*vn2 -
251 v = std::sqrt ( 2.0 * std::log ( s1 / (r*g
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)
263
264 return std::sqrt(0.5) * flatToGaussian(.5*(t
265
266 }
267
268 double HepStat::erf (double x) {
269
270 // Math:
271 //
272 // For any given x we can "quickly" find t0 =
273 //
274 // Then we can find x1 = inverseErf (t0) = inv
275 //
276 // Expanding in the small epsion,
277 //
278 // x1 = x + epsilon * [deriv(inverseErf(u),u)
279 //
280 // so epsilon is (x1-x) / [deriv(inverseErf(u)
281 //
282 // But the derivative of an inverse function i
283 // function, so
284 // epsilon = (x1-x) * [deriv(erf(v),v) at v =
285 //
286 // And the definition of the erf integral make
287 // Ultimately,
288 //
289 // erf(x) = erfQ(x) - (inverseErf(erfQ(x))-x)
290 //
291
292 double t0 = erfQ(x);
293 double deriv = std::exp(-x*x) * (2.0 / std::
294
295 return t0 - (inverseErf (t0) - x) * deriv;
296
297 }
298
299
300 } // namespace CLHEP
301
302