00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #include "CLHEP/Random/Stat.h"
00022 #include "CLHEP/Units/PhysicalConstants.h"
00023 #include <iostream>
00024 #include <cmath>
00025
00026 namespace CLHEP {
00027
00028 double transformSmall (double r);
00029
00030
00031
00032
00033
00034 #ifdef Traces
00035 #define Trace1
00036 #define Trace2
00037 #define Trace3
00038 #endif
00039
00040
00041
00042
00043
00044
00045
00046 #define Table0size 200
00047 #define Table1size 250
00048 #define Table2size 200
00049 #define Table3size 250
00050 #define Table4size 1000
00051 #define TableSize (Table0size+Table1size+Table2size+Table3size+Table4size)
00052
00053 static const int Tsizes[5] = { Table0size,
00054 Table1size,
00055 Table2size,
00056 Table3size,
00057 Table4size };
00058
00059 #define Table0step (2.0E-13)
00060 #define Table1step (4.0E-11)
00061 #define Table2step (1.0E-8)
00062 #define Table3step (2.0E-6)
00063 #define Table4step (5.0E-4)
00064
00065 static const double Tsteps[5] = { Table0step,
00066 Table1step,
00067 Table2step,
00068 Table3step,
00069 Table4step };
00070
00071 #define Table0offset 0
00072 #define Table1offset (2*(Table0size))
00073 #define Table2offset (2*(Table0size + Table1size))
00074 #define Table3offset (2*(Table0size + Table1size + Table2size))
00075 #define Table4offset (2*(Table0size + Table1size + Table2size + Table3size))
00076
00077 static const int Toffsets[5] = { Table0offset,
00078 Table1offset,
00079 Table2offset,
00080 Table3offset,
00081 Table4offset };
00082
00083
00084
00085 static const double gaussTables [2*TableSize] = {
00086 #include "CLHEP/Random/gaussTables.cdat"
00087 };
00088
00089 double HepStat::flatToGaussian (double r) {
00090
00091 double sign = +1.0;
00092
00093
00094 #ifdef Trace1
00095 std::cout << "r = " << r << "\n";
00096 #endif
00097
00098 if ( r > .5 ) {
00099 r = 1-r;
00100 sign = -1.0;
00101 #ifdef Trace1
00102 std::cout << "r = " << r << "sign negative \n";
00103 #endif
00104 } else if ( r == .5 ) {
00105 return 0.0;
00106 }
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119 const double* tptr = 0;
00120 double dx = 0;
00121 double h = 0;
00122
00123
00124
00125
00126 register int index;
00127
00128 if ( r >= Table4step ) {
00129
00130 index = int((Table4size<<1) * r);
00131 if (index <= 0) index = 1;
00132 if (index >= Table4size) index = Table4size-1;
00133 dx = (Table4size<<1) * r - index;
00134 h = Table4step;
00135 #ifdef Trace2
00136 std::cout << "index = " << index << " dx = " << dx << " h = " << h << "\n";
00137 #endif
00138 index = (index<<1) + (Table4offset-2);
00139
00140
00141
00142 #ifdef Trace2
00143 std::cout << "offset index = " << index << "\n";
00144 #endif
00145
00146 tptr = &gaussTables [index];
00147
00148 } else if (r < Tsteps[0]) {
00149
00150
00151 return (sign * transformSmall (r));
00152
00153 } else {
00154
00155 for ( int tableN = 3; tableN >= 0; tableN-- ) {
00156 if ( r < Tsteps[tableN] ) {continue;}
00157 #ifdef Trace2
00158 std::cout << "Using table " << tableN << "\n";
00159 #endif
00160 double step = Tsteps[tableN];
00161 index = int(r/step);
00162
00163
00164
00165
00166 if (index == 0) index = 1;
00167 if (index >= Tsizes[tableN]) index = Tsizes[tableN] - 1;
00168 dx = r/step - index;
00169 h = step;
00170 #ifdef Trace2
00171 std::cout << "index = " << index << " dx = " << dx << " h = " << h << "\n";
00172 #endif
00173 index = (index<<1) + Toffsets[tableN] - 2;
00174 tptr = &gaussTables [index];
00175 break;
00176 }
00177
00178
00179
00180
00181
00182 }
00183
00184
00185
00186
00187 double y0 = *tptr++;
00188 double d0 = *tptr++;
00189 double y1 = *tptr++;
00190 double d1 = *tptr;
00191
00192 #ifdef Trace3
00193 std::cout << "y0: " << y0 << " y1: " << y1 << " d0: " << d0 << " d1: " << d1 << "\n";
00194 #endif
00195
00196 double x2 = dx * dx;
00197 double oneMinusX = 1 - dx;
00198 double oneMinusX2 = oneMinusX * oneMinusX;
00199
00200 double f0 = (2. * dx + 1.) * oneMinusX2;
00201 double f1 = (3. - 2. * dx) * x2;
00202 double g0 = h * dx * oneMinusX2;
00203 double g1 = - h * oneMinusX * x2;
00204
00205 #ifdef Trace3
00206 std::cout << "f0: " << f0 << " f1: " << f1 << " g0: " << g0 << " g1: " << g1 << "\n";
00207 #endif
00208
00209 double answer = f0 * y0 + f1 * y1 + g0 * d0 + g1 * d1;
00210
00211 #ifdef Trace1
00212 std::cout << "variate is: " << sign*answer << "\n";
00213 #endif
00214
00215 return sign * answer;
00216
00217 }
00218
00219 double transformSmall (double r) {
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240 double eps = 1.0e-7;
00241 double guess = 7.5;
00242 double v;
00243
00244 for ( int i = 1; i < 50; i++ ) {
00245 double vn2 = 1.0/(guess*guess);
00246 double s1 = -13*11*9*7*5*3 * vn2*vn2*vn2*vn2*vn2*vn2*vn2;
00247 s1 += 11*9*7*5*3 * vn2*vn2*vn2*vn2*vn2*vn2;
00248 s1 += -9*7*5*3 * vn2*vn2*vn2*vn2*vn2;
00249 s1 += 7*5*3 * vn2*vn2*vn2*vn2;
00250 s1 += -5*3 * vn2*vn2*vn2;
00251 s1 += 3 * vn2*vn2 - vn2 + 1.0;
00252 v = std::sqrt ( 2.0 * std::log ( s1 / (r*guess*std::sqrt(CLHEP::twopi)) ) );
00253 if ( std::abs(v-guess) < eps ) break;
00254 guess = v;
00255 }
00256
00257 return -v;
00258
00259 }
00260
00261 double HepStat::inverseErf (double t) {
00262
00263
00264
00265 return std::sqrt(0.5) * flatToGaussian(.5*(t+1));
00266
00267 }
00268
00269 double HepStat::erf (double x) {
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293 double t0 = erfQ(x);
00294 double deriv = std::exp(-x*x) * (2.0 / std::sqrt(CLHEP::pi));
00295
00296 return t0 - (inverseErf (t0) - x) * deriv;
00297
00298 }
00299
00300
00301 }
00302