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