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