CLHEP 2.0.4.7 Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // $Id: flatToGaussian.cc,v 1.4 2003/08/13 20:00:12 garren Exp $ 00002 // -*- C++ -*- 00003 // 00004 // ----------------------------------------------------------------------- 00005 // HEP Random 00006 // --- flatToGaussian --- 00007 // class implementation file 00008 // ----------------------------------------------------------------------- 00009 00010 // Contains the methods that depend on the 30K-footpring gaussTables.cdat. 00011 // 00012 // flatToGaussian (double x) 00013 // inverseErf (double x) 00014 // erf (double x) 00015 00016 // ======================================================================= 00017 // M Fischler - Created 1/25/00. 00018 // 00019 // ======================================================================= 00020 00021 #include "CLHEP/Random/Stat.h" 00022 #include "CLHEP/Random/defs.h" 00023 #include "CLHEP/Units/PhysicalConstants.h" 00024 #include <iostream> 00025 #include <cmath> 00026 00027 using namespace std; 00028 00029 namespace CLHEP { 00030 00031 double transformSmall (double r); 00032 00033 00034 // 00035 // Table of errInts, for use with transform(r) and quickTransform(r) 00036 // 00037 00038 #ifdef Traces 00039 #define Trace1 00040 #define Trace2 00041 #define Trace3 00042 #endif 00043 00044 // Since all these are this is static to this compilation unit only, the 00045 // info is establised a priori and not at each invocation. 00046 00047 // The main data is of course the gaussTables table; the rest is all 00048 // bookkeeping to know what the tables mean. 00049 00050 #define Table0size 200 00051 #define Table1size 250 00052 #define Table2size 200 00053 #define Table3size 250 00054 #define Table4size 1000 00055 #define TableSize (Table0size+Table1size+Table2size+Table3size+Table4size) 00056 00057 static const int Tsizes[5] = { Table0size, 00058 Table1size, 00059 Table2size, 00060 Table3size, 00061 Table4size }; 00062 00063 #define Table0step (2.0E-13) 00064 #define Table1step (4.0E-11) 00065 #define Table2step (1.0E-8) 00066 #define Table3step (2.0E-6) 00067 #define Table4step (5.0E-4) 00068 00069 static const double Tsteps[5] = { Table0step, 00070 Table1step, 00071 Table2step, 00072 Table3step, 00073 Table4step }; 00074 00075 #define Table0offset 0 00076 #define Table1offset (2*(Table0size)) 00077 #define Table2offset (2*(Table0size + Table1size)) 00078 #define Table3offset (2*(Table0size + Table1size + Table2size)) 00079 #define Table4offset (2*(Table0size + Table1size + Table2size + Table3size)) 00080 00081 static const int Toffsets[5] = { Table0offset, 00082 Table1offset, 00083 Table2offset, 00084 Table3offset, 00085 Table4offset }; 00086 00087 // Here comes the big (30K bytes) table, kept in a file --- 00088 00089 static const double gaussTables [2*TableSize] = { 00090 #include "gaussTables.cdat" 00091 }; 00092 00093 00094 double HepStat::flatToGaussian (double r) { 00095 00096 double sign = +1.0; // We always compute a negative number of 00097 // sigmas. For r > 0 we will multiply by 00098 // sign = -1 to return a positive number. 00099 #ifdef Trace1 00100 std::cout << "r = " << r << "\n"; 00101 #endif 00102 00103 if ( r > .5 ) { 00104 r = 1-r; 00105 sign = -1.0; 00106 #ifdef Trace1 00107 std::cout << "r = " << r << "sign negative \n"; 00108 #endif 00109 } else if ( r == .5 ) { 00110 return 0.0; 00111 } 00112 00113 // Find a pointer to the proper table entries, along with the fraction 00114 // of the way in the relevant bin dx and the bin size h. 00115 00116 // Optimize for the case of table 4 by testing for that first. 00117 // By removing that case from the for loop below, we save not only 00118 // several table lookups, but also several index calculations that 00119 // now become (compile-time) constants. 00120 // 00121 // Past the case of table 4, we need not be as concerned about speed since 00122 // this will happen only .1% of the time. 00123 00124 const double* tptr = 0; 00125 double dx = 0; 00126 double h = 0; 00127 00128 // The following big if block will locate tptr, h, and dx from whichever 00129 // table is applicable: 00130 00131 register int index; 00132 00133 if ( r >= Table4step ) { 00134 00135 index = int((Table4size<<1) * r); // 1 to Table4size-1 00136 if (index <= 0) index = 1; // in case of rounding problem 00137 if (index >= Table4size) index = Table4size-1; 00138 dx = (Table4size<<1) * r - index; // fraction of way to next bin 00139 h = Table4step; 00140 #ifdef Trace2 00141 std::cout << "index = " << index << " dx = " << dx << " h = " << h << "\n"; 00142 #endif 00143 index = (index<<1) + (Table4offset-2); 00144 // at r = table4step+eps, index refers to the start of table 4 00145 // and at r = .5 - eps, index refers to the next-to-last entry: 00146 // remember, the table has two numbers per actual entry. 00147 #ifdef Trace2 00148 std::cout << "offset index = " << index << "\n"; 00149 #endif 00150 00151 tptr = &gaussTables [index]; 00152 00153 } else if (r < Tsteps[0]) { 00154 00155 // If r is so small none of the tables apply, use the asymptotic formula. 00156 return (sign * transformSmall (r)); 00157 00158 } else { 00159 00160 for ( int tableN = 3; tableN >= 0; tableN-- ) { 00161 if ( r < Tsteps[tableN] ) {continue;} // can't happen when tableN==0 00162 #ifdef Trace2 00163 std::cout << "Using table " << tableN << "\n"; 00164 #endif 00165 double step = Tsteps[tableN]; 00166 index = int(r/step); // 1 to TableNsize-1 00167 // The following two tests should probably never be true, but 00168 // roundoff might cause index to be outside its proper range. 00169 // In such a case, the interpolation still makes sense, but we 00170 // need to take care that tptr points to valid data in the right table. 00171 if (index == 0) index = 1; 00172 if (index >= Tsizes[tableN]) index = Tsizes[tableN] - 1; 00173 dx = r/step - index; // fraction of way to next bin 00174 h = step; 00175 #ifdef Trace2 00176 std::cout << "index = " << index << " dx = " << dx << " h = " << h << "\n"; 00177 #endif 00178 index = (index<<1) + Toffsets[tableN] - 2; 00179 tptr = &gaussTables [index]; 00180 break; 00181 } // end of the for loop which finds tptr, dx and h in one of the tables 00182 00183 // The code can only get to here by a break statement, having set dx etc. 00184 // It not get to here without going through one of the breaks, 00185 // because before the for loop we test for the case of r < Tsteps[0]. 00186 00187 } // End of the big if block. 00188 00189 // At the end of this if block, we have tptr, dx and h -- and if r is less 00190 // than the smallest step, we have already returned the proper answer. 00191 00192 double y0 = *tptr++; 00193 double d0 = *tptr++; 00194 double y1 = *tptr++; 00195 double d1 = *tptr; 00196 00197 #ifdef Trace3 00198 std::cout << "y0: " << y0 << " y1: " << y1 << " d0: " << d0 << " d1: " << d1 << "\n"; 00199 #endif 00200 00201 double x2 = dx * dx; 00202 double oneMinusX = 1 - dx; 00203 double oneMinusX2 = oneMinusX * oneMinusX; 00204 00205 double f0 = (2. * dx + 1.) * oneMinusX2; 00206 double f1 = (3. - 2. * dx) * x2; 00207 double g0 = h * dx * oneMinusX2; 00208 double g1 = - h * oneMinusX * x2; 00209 00210 #ifdef Trace3 00211 std::cout << "f0: " << f0 << " f1: " << f1 << " g0: " << g0 << " g1: " << g1 << "\n"; 00212 #endif 00213 00214 double answer = f0 * y0 + f1 * y1 + g0 * d0 + g1 * d1; 00215 00216 #ifdef Trace1 00217 std::cout << "variate is: " << sign*answer << "\n"; 00218 #endif 00219 00220 return sign * answer; 00221 00222 } // flatToGaussian 00223 00224 00225 double transformSmall (double r) { 00226 00227 // Solve for -v in the asymtotic formula 00228 // 00229 // errInt (-v) = exp(-v*v/2) 1 1*3 1*3*5 00230 // ------------ * (1 - ---- + ---- - ----- + ... ) 00231 // v*sqrt(2*pi) v**2 v**4 v**6 00232 00233 // The value of r (=errInt(-v)) supplied is going to less than 2.0E-13, 00234 // which is such that v < -7.25. Since the value of r is meaningful only 00235 // to an absolute error of 1E-16 (double precision accuracy for a number 00236 // which on the high side could be of the form 1-epsilon), computing 00237 // v to more than 3-4 digits of accuracy is suspect; however, to ensure 00238 // smoothness with the table generator (which uses quite a few terms) we 00239 // also use terms up to 1*3*5* ... *13/v**14, and insist on accuracy of 00240 // solution at the level of 1.0e-7. 00241 00242 // This routine is called less than one time in a trillion firings, so 00243 // speed is of no concern. As a matter of technique, we terminate the 00244 // iterations in case they would be infinite, but this should not happen. 00245 00246 double eps = 1.0e-7; 00247 double guess = 7.5; 00248 double v; 00249 00250 for ( int i = 1; i < 50; i++ ) { 00251 double vn2 = 1.0/(guess*guess); 00252 double s = -13*11*9*7*5*3 * vn2*vn2*vn2*vn2*vn2*vn2*vn2; 00253 s += 11*9*7*5*3 * vn2*vn2*vn2*vn2*vn2*vn2; 00254 s += -9*7*5*3 * vn2*vn2*vn2*vn2*vn2; 00255 s += 7*5*3 * vn2*vn2*vn2*vn2; 00256 s += -5*3 * vn2*vn2*vn2; 00257 s += 3 * vn2*vn2 - vn2 + 1.0; 00258 v = sqrt ( 2.0 * log ( s / (r*guess*sqrt(CLHEP::twopi)) ) ); 00259 if ( abs(v-guess) < eps ) break; 00260 guess = v; 00261 } 00262 00263 return -v; 00264 00265 } // transformSmall() 00266 00267 00268 double HepStat::inverseErf (double t) { 00269 00270 // This uses erf(x) = 2*gaussCDF(sqrt(2)*x) - 1 00271 00272 return sqrt(0.5) * flatToGaussian(.5*(t+1)); 00273 00274 } 00275 00276 00277 double HepStat::erf (double x) { 00278 00279 // Math: 00280 // 00281 // For any given x we can "quickly" find t0 = erfQ (x) = erf(x) + epsilon. 00282 // 00283 // Then we can find x1 = inverseErf (t0) = inverseErf (erf(x)+epsilon) 00284 // 00285 // Expanding in the small epsion, 00286 // 00287 // x1 = x + epsilon * [deriv(inverseErf(u),u) at u = t0] + O(epsilon**2) 00288 // 00289 // so epsilon is (x1-x) / [deriv(inverseErf(u),u) at u = t0] + O(epsilon**2) 00290 // 00291 // But the derivative of an inverse function is one over the derivative of the 00292 // function, so 00293 // epsilon = (x1-x) * [deriv(erf(v),v) at v = x] + O(epsilon**2) 00294 // 00295 // And the definition of the erf integral makes that derivative trivial. 00296 // Ultimately, 00297 // 00298 // erf(x) = erfQ(x) - (inverseErf(erfQ(x))-x) * exp(-x**2) * 2/sqrt(CLHEP::pi) 00299 // 00300 00301 double t0 = erfQ(x); 00302 double deriv = exp(-x*x) * (2.0 / sqrt(CLHEP::pi)); 00303 00304 return t0 - (inverseErf (t0) - x) * deriv; 00305 00306 } 00307 00308 00309 } // namespace CLHEP 00310