Geant4 Cross Reference |
1 // -*- C++ -*- 2 // 3 // ----------------------------------------------------------------------- 4 // HEP Random 5 // --- RandGaussQ --- 6 // class implementation file 7 // ----------------------------------------------------------------------- 8 9 // ======================================================================= 10 // M Fischler - Created 24 Jan 2000 11 // M Fischler - put and get to/from streams 12/13/04 12 // ======================================================================= 13 14 #include "CLHEP/Random/RandGaussQ.h" 15 #include "CLHEP/Units/PhysicalConstants.h" 16 #include <iostream> 17 #include <cmath> // for std::log() 18 19 namespace CLHEP { 20 21 std::string RandGaussQ::name() const {return "RandGaussQ";} 22 HepRandomEngine & RandGaussQ::engine() {return RandGauss::engine();} 23 24 RandGaussQ::~RandGaussQ() { 25 } 26 27 double RandGaussQ::operator()() { 28 return transformQuick(localEngine->flat()) * defaultStdDev + defaultMean; 29 } 30 31 double RandGaussQ::operator()( double mean, double stdDev ) { 32 return transformQuick(localEngine->flat()) * stdDev + mean; 33 } 34 35 void RandGaussQ::shootArray( const int size, double* vect, 36 double mean, double stdDev ) 37 { 38 for( double* v = vect; v != vect + size; ++v ) 39 *v = shoot(mean,stdDev); 40 } 41 42 void RandGaussQ::shootArray( HepRandomEngine* anEngine, 43 const int size, double* vect, 44 double mean, double stdDev ) 45 { 46 for( double* v = vect; v != vect + size; ++v ) 47 *v = shoot(anEngine,mean,stdDev); 48 } 49 50 void RandGaussQ::fireArray( const int size, double* vect) 51 { 52 for( double* v = vect; v != vect + size; ++v ) 53 *v = fire( defaultMean, defaultStdDev ); 54 } 55 56 void RandGaussQ::fireArray( const int size, double* vect, 57 double mean, double stdDev ) 58 { 59 for( double* v = vect; v != vect + size; ++v ) 60 *v = fire( mean, stdDev ); 61 } 62 63 // 64 // Table of errInts, for use with transform(r) and quickTransform(r) 65 // 66 67 // Since all these are this is static to this compilation unit only, the 68 // info is establised a priori and not at each invocation. 69 70 // The main data is of course the gaussQTables table; the rest is all 71 // bookkeeping to know what the tables mean. 72 73 #define Table0size 250 74 #define Table1size 1000 75 #define TableSize (Table0size+Table1size) 76 77 #define Table0step (2.0E-6) 78 #define Table1step (5.0E-4) 79 80 #define Table0scale (1.0/Table1step) 81 82 #define Table0offset 0 83 #define Table1offset (Table0size) 84 85 // Here comes the big (5K bytes) table, kept in a file --- 86 87 static const float gaussTables [TableSize] = { 88 #include "CLHEP/Random/gaussQTables.cdat" 89 }; 90 91 double RandGaussQ::transformQuick (double r) { 92 double sign = +1.0; // We always compute a negative number of 93 // sigmas. For r > 0 we will multiply by 94 // sign = -1 to return a positive number. 95 if ( r > .5 ) { 96 r = 1-r; 97 sign = -1.0; 98 } 99 100 int index; 101 double dx; 102 103 if ( r >= Table1step ) { 104 index = int((Table1size<<1) * r); // 1 to Table1size 105 if (index == Table1size) return 0.0; 106 dx = (Table1size<<1) * r - index; // fraction of way to next bin 107 index += Table1offset-1; 108 } else if ( r > Table0step ) { 109 double rr = r * Table0scale; 110 index = int(Table0size * rr); // 1 to Table0size 111 dx = Table0size * rr - index; // fraction of way to next bin 112 index += Table0offset-1; 113 } else { // r <= Table0step - not in tables 114 return sign*transformSmall(r); 115 } 116 117 double y0 = gaussTables [index++]; 118 double y1 = gaussTables [index]; 119 120 return (float) (sign * ( y1 * dx + y0 * (1.0-dx) )); 121 122 } // transformQuick() 123 124 double RandGaussQ::transformSmall (double r) { 125 126 // Solve for -v in the asymtotic formula 127 // 128 // errInt (-v) = exp(-v*v/2) 1 1*3 1*3*5 129 // ------------ * (1 - ---- + ---- - ----- + ... ) 130 // v*sqrt(2*pi) v**2 v**4 v**6 131 132 // The value of r (=errInt(-v)) supplied is going to less than 2.0E-13, 133 // which is such that v < -7.25. Since the value of r is meaningful only 134 // to an absolute error of 1E-16 (double precision accuracy for a number 135 // which on the high side could be of the form 1-epsilon), computing 136 // v to more than 3-4 digits of accuracy is suspect; however, to ensure 137 // smoothness with the table generator (which uses quite a few terms) we 138 // also use terms up to 1*3*5* ... *13/v**14, and insist on accuracy of 139 // solution at the level of 1.0e-7. 140 141 // This routine is called less than one time in a million firings, so 142 // speed is of no concern. As a matter of technique, we terminate the 143 // iterations in case they would be infinite, but this should not happen. 144 145 double eps = 1.0e-7; 146 double guess = 7.5; 147 double v; 148 149 for ( int i = 1; i < 50; i++ ) { 150 double vn2 = 1.0/(guess*guess); 151 double s1 = -13*11*9*7*5*3 * vn2*vn2*vn2*vn2*vn2*vn2*vn2; 152 s1 += 11*9*7*5*3 * vn2*vn2*vn2*vn2*vn2*vn2; 153 s1 += -9*7*5*3 * vn2*vn2*vn2*vn2*vn2; 154 s1 += 7*5*3 * vn2*vn2*vn2*vn2; 155 s1 += -5*3 * vn2*vn2*vn2; 156 s1 += 3 * vn2*vn2 - vn2 + 1.0; 157 v = std::sqrt ( 2.0 * std::log ( s1 / (r*guess*std::sqrt(CLHEP::twopi)) ) ); 158 if ( std::fabs(v-guess) < eps ) break; 159 guess = v; 160 } 161 return -v; 162 163 } // transformSmall() 164 165 std::ostream & RandGaussQ::put ( std::ostream & os ) const { 166 long pr=os.precision(20); 167 os << " " << name() << "\n"; 168 RandGauss::put(os); 169 os.precision(pr); 170 return os; 171 } 172 173 std::istream & RandGaussQ::get ( std::istream & is ) { 174 std::string inName; 175 is >> inName; 176 if (inName != name()) { 177 is.clear(std::ios::badbit | is.rdstate()); 178 std::cerr << "Mismatch when expecting to read state of a " 179 << name() << " distribution\n" 180 << "Name found was " << inName 181 << "\nistream is left in the badbit state\n"; 182 return is; 183 } 184 RandGauss::get(is); 185 return is; 186 } 187 188 } // namespace CLHEP 189