00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015 #include "CLHEP/Random/RandGaussQ.h"
00016 #include "CLHEP/Units/PhysicalConstants.h"
00017 #include <iostream>
00018 #include <cmath>
00019
00020 namespace CLHEP {
00021
00022 std::string RandGaussQ::name() const {return "RandGaussQ";}
00023 HepRandomEngine & RandGaussQ::engine() {return RandGauss::engine();}
00024
00025 RandGaussQ::~RandGaussQ() {
00026 }
00027
00028 double RandGaussQ::operator()() {
00029 return transformQuick(localEngine->flat()) * defaultStdDev + defaultMean;
00030 }
00031
00032 double RandGaussQ::operator()( double mean, double stdDev ) {
00033 return transformQuick(localEngine->flat()) * stdDev + mean;
00034 }
00035
00036 void RandGaussQ::shootArray( const int size, double* vect,
00037 double mean, double stdDev )
00038 {
00039 for( double* v = vect; v != vect + size; ++v )
00040 *v = shoot(mean,stdDev);
00041 }
00042
00043 void RandGaussQ::shootArray( HepRandomEngine* anEngine,
00044 const int size, double* vect,
00045 double mean, double stdDev )
00046 {
00047 for( double* v = vect; v != vect + size; ++v )
00048 *v = shoot(anEngine,mean,stdDev);
00049 }
00050
00051 void RandGaussQ::fireArray( const int size, double* vect)
00052 {
00053 for( double* v = vect; v != vect + size; ++v )
00054 *v = fire( defaultMean, defaultStdDev );
00055 }
00056
00057 void RandGaussQ::fireArray( const int size, double* vect,
00058 double mean, double stdDev )
00059 {
00060 for( double* v = vect; v != vect + size; ++v )
00061 *v = fire( mean, stdDev );
00062 }
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074 #define Table0size 250
00075 #define Table1size 1000
00076 #define TableSize (Table0size+Table1size)
00077
00078 #define Table0step (2.0E-6)
00079 #define Table1step (5.0E-4)
00080
00081 #define Table0scale (1.0/Table1step)
00082
00083 #define Table0offset 0
00084 #define Table1offset (Table0size)
00085
00086
00087
00088 static const float gaussTables [TableSize] = {
00089 #include "CLHEP/Random/gaussQTables.cdat"
00090 };
00091
00092 double RandGaussQ::transformQuick (double r) {
00093 double sign = +1.0;
00094
00095
00096 if ( r > .5 ) {
00097 r = 1-r;
00098 sign = -1.0;
00099 }
00100
00101 register int index;
00102 double dx;
00103
00104 if ( r >= Table1step ) {
00105 index = int((Table1size<<1) * r);
00106 if (index == Table1size) return 0.0;
00107 dx = (Table1size<<1) * r - index;
00108 index += Table1offset-1;
00109 } else if ( r > Table0step ) {
00110 double rr = r * Table0scale;
00111 index = int(Table0size * rr);
00112 dx = Table0size * rr - index;
00113 index += Table0offset-1;
00114 } else {
00115 return sign*transformSmall(r);
00116 }
00117
00118 double y0 = gaussTables [index++];
00119 double y1 = gaussTables [index];
00120
00121 return (float) (sign * ( y1 * dx + y0 * (1.0-dx) ));
00122
00123 }
00124
00125 double RandGaussQ::transformSmall (double r) {
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146 double eps = 1.0e-7;
00147 double guess = 7.5;
00148 double v;
00149
00150 for ( int i = 1; i < 50; i++ ) {
00151 double vn2 = 1.0/(guess*guess);
00152 double s1 = -13*11*9*7*5*3 * vn2*vn2*vn2*vn2*vn2*vn2*vn2;
00153 s1 += 11*9*7*5*3 * vn2*vn2*vn2*vn2*vn2*vn2;
00154 s1 += -9*7*5*3 * vn2*vn2*vn2*vn2*vn2;
00155 s1 += 7*5*3 * vn2*vn2*vn2*vn2;
00156 s1 += -5*3 * vn2*vn2*vn2;
00157 s1 += 3 * vn2*vn2 - vn2 + 1.0;
00158 v = std::sqrt ( 2.0 * std::log ( s1 / (r*guess*std::sqrt(CLHEP::twopi)) ) );
00159 if ( std::fabs(v-guess) < eps ) break;
00160 guess = v;
00161 }
00162 return -v;
00163
00164 }
00165
00166 std::ostream & RandGaussQ::put ( std::ostream & os ) const {
00167 int pr=os.precision(20);
00168 os << " " << name() << "\n";
00169 RandGauss::put(os);
00170 os.precision(pr);
00171 return os;
00172 }
00173
00174 std::istream & RandGaussQ::get ( std::istream & is ) {
00175 std::string inName;
00176 is >> inName;
00177 if (inName != name()) {
00178 is.clear(std::ios::badbit | is.rdstate());
00179 std::cerr << "Mismatch when expecting to read state of a "
00180 << name() << " distribution\n"
00181 << "Name found was " << inName
00182 << "\nistream is left in the badbit state\n";
00183 return is;
00184 }
00185 RandGauss::get(is);
00186 return is;
00187 }
00188
00189 }