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