CLHEP 2.0.4.7 Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // $Id: RandGaussQ.cc,v 1.4.4.1 2005/03/18 22:26:48 garren Exp $ 00002 // -*- C++ -*- 00003 // 00004 // ----------------------------------------------------------------------- 00005 // HEP Random 00006 // --- RandGaussQ --- 00007 // class implementation file 00008 // ----------------------------------------------------------------------- 00009 00010 // ======================================================================= 00011 // M Fischler - Created 24 Jan 2000 00012 // M Fischler - put and get to/from streams 12/13/04 00013 // ======================================================================= 00014 00015 #include "CLHEP/Random/defs.h" 00016 #include "CLHEP/Random/RandGaussQ.h" 00017 #include "CLHEP/Units/PhysicalConstants.h" 00018 #include <iostream> 00019 #include <cmath> // for log() 00020 00021 namespace CLHEP { 00022 00023 std::string RandGaussQ::name() const {return "RandGaussQ";} 00024 HepRandomEngine & RandGaussQ::engine() {return RandGauss::engine();} 00025 00026 RandGaussQ::~RandGaussQ() { 00027 } 00028 00029 RandGaussQ::RandGaussQ(const RandGaussQ& right) : RandGauss(right) { 00030 } 00031 00032 double RandGaussQ::operator()() { 00033 return transformQuick(localEngine->flat()) * defaultStdDev + defaultMean; 00034 } 00035 00036 double RandGaussQ::operator()( double mean, double stdDev ) { 00037 return transformQuick(localEngine->flat()) * stdDev + mean; 00038 } 00039 00040 void RandGaussQ::shootArray( const int size, double* vect, 00041 double mean, double stdDev ) 00042 { 00043 int i; 00044 00045 for (i=0; i<size; ++i) { 00046 vect[i] = shoot(mean,stdDev); 00047 } 00048 } 00049 00050 void RandGaussQ::shootArray( HepRandomEngine* anEngine, 00051 const int size, double* vect, 00052 double mean, double stdDev ) 00053 { 00054 int i; 00055 00056 for (i=0; i<size; ++i) { 00057 vect[i] = shoot(anEngine,mean,stdDev); 00058 } 00059 } 00060 00061 void RandGaussQ::fireArray( const int size, double* vect) 00062 { 00063 int i; 00064 00065 for (i=0; i<size; ++i) { 00066 vect[i] = fire( defaultMean, defaultStdDev ); 00067 } 00068 } 00069 00070 void RandGaussQ::fireArray( const int size, double* vect, 00071 double mean, double stdDev ) 00072 { 00073 int i; 00074 00075 for (i=0; i<size; ++i) { 00076 vect[i] = fire( mean, stdDev ); 00077 } 00078 } 00079 00080 00081 // 00082 // Table of errInts, for use with transform(r) and quickTransform(r) 00083 // 00084 00085 // Since all these are this is static to this compilation unit only, the 00086 // info is establised a priori and not at each invocation. 00087 00088 // The main data is of course the gaussQTables table; the rest is all 00089 // bookkeeping to know what the tables mean. 00090 00091 #define Table0size 250 00092 #define Table1size 1000 00093 #define TableSize (Table0size+Table1size) 00094 00095 #define Table0step (2.0E-6) 00096 #define Table1step (5.0E-4) 00097 00098 #define Table0scale (1.0/Table1step) 00099 00100 #define Table0offset 0 00101 #define Table1offset (Table0size) 00102 00103 // Here comes the big (5K bytes) table, kept in a file --- 00104 00105 static const float gaussTables [TableSize] = { 00106 #include "gaussQTables.cdat" 00107 }; 00108 00109 00110 double RandGaussQ::transformQuick (double r) { 00111 double sign = +1.0; // We always compute a negative number of 00112 // sigmas. For r > 0 we will multiply by 00113 // sign = -1 to return a positive number. 00114 if ( r > .5 ) { 00115 r = 1-r; 00116 sign = -1.0; 00117 } 00118 00119 register int index; 00120 double dx; 00121 00122 if ( r >= Table1step ) { 00123 index = int((Table1size<<1) * r); // 1 to Table1size 00124 if (index == Table1size) return 0.0; 00125 dx = (Table1size<<1) * r - index; // fraction of way to next bin 00126 index += Table1offset-1; 00127 } else if ( r > Table0step ) { 00128 double rr = r * Table0scale; 00129 index = int(Table0size * rr); // 1 to Table0size 00130 dx = Table0size * rr - index; // fraction of way to next bin 00131 index += Table0offset-1; 00132 } else { // r <= Table0step - not in tables 00133 return sign*transformSmall(r); 00134 } 00135 00136 double y0 = gaussTables [index++]; 00137 double y1 = gaussTables [index]; 00138 00139 return (float) (sign * ( y1 * dx + y0 * (1.0-dx) )); 00140 00141 } // transformQuick() 00142 00143 00144 00145 double RandGaussQ::transformSmall (double r) { 00146 00147 // Solve for -v in the asymtotic formula 00148 // 00149 // errInt (-v) = exp(-v*v/2) 1 1*3 1*3*5 00150 // ------------ * (1 - ---- + ---- - ----- + ... ) 00151 // v*sqrt(2*pi) v**2 v**4 v**6 00152 00153 // The value of r (=errInt(-v)) supplied is going to less than 2.0E-13, 00154 // which is such that v < -7.25. Since the value of r is meaningful only 00155 // to an absolute error of 1E-16 (double precision accuracy for a number 00156 // which on the high side could be of the form 1-epsilon), computing 00157 // v to more than 3-4 digits of accuracy is suspect; however, to ensure 00158 // smoothness with the table generator (which uses quite a few terms) we 00159 // also use terms up to 1*3*5* ... *13/v**14, and insist on accuracy of 00160 // solution at the level of 1.0e-7. 00161 00162 // This routine is called less than one time in a million firings, so 00163 // speed is of no concern. As a matter of technique, we terminate the 00164 // iterations in case they would be infinite, but this should not happen. 00165 00166 double eps = 1.0e-7; 00167 double guess = 7.5; 00168 double v; 00169 00170 for ( int i = 1; i < 50; i++ ) { 00171 double vn2 = 1.0/(guess*guess); 00172 double s = -13*11*9*7*5*3 * vn2*vn2*vn2*vn2*vn2*vn2*vn2; 00173 s += 11*9*7*5*3 * vn2*vn2*vn2*vn2*vn2*vn2; 00174 s += -9*7*5*3 * vn2*vn2*vn2*vn2*vn2; 00175 s += 7*5*3 * vn2*vn2*vn2*vn2; 00176 s += -5*3 * vn2*vn2*vn2; 00177 s += 3 * vn2*vn2 - vn2 + 1.0; 00178 v = sqrt ( 2.0 * log ( s / (r*guess*sqrt(CLHEP::twopi)) ) ); 00179 if ( fabs(v-guess) < eps ) break; 00180 guess = v; 00181 } 00182 return -v; 00183 00184 } // transformSmall() 00185 00186 std::ostream & RandGaussQ::put ( std::ostream & os ) const { 00187 int pr=os.precision(20); 00188 os << " " << name() << "\n"; 00189 RandGauss::put(os); 00190 os.precision(pr); 00191 return os; 00192 } 00193 00194 std::istream & RandGaussQ::get ( std::istream & is ) { 00195 std::string inName; 00196 is >> inName; 00197 if (inName != name()) { 00198 is.clear(std::ios::badbit | is.rdstate()); 00199 std::cerr << "Mismatch when expecting to read state of a " 00200 << name() << " distribution\n" 00201 << "Name found was " << inName 00202 << "\nistream is left in the badbit state\n"; 00203 return is; 00204 } 00205 RandGauss::get(is); 00206 return is; 00207 } 00208 00209 } // namespace CLHEP