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