Geant4 Cross Reference |
1 // -*- C++ -*- 1 // -*- C++ -*- 2 // 2 // 3 // ------------------------------------------- 3 // ----------------------------------------------------------------------- 4 // HEP Random 4 // HEP Random 5 // --- RandGamma --- 5 // --- RandGamma --- 6 // class implementation f 6 // class implementation file 7 // ------------------------------------------- 7 // ----------------------------------------------------------------------- 8 8 9 // =========================================== 9 // ======================================================================= 10 // John Marraffino - Created: 12th May 1998 10 // John Marraffino - Created: 12th May 1998 11 // M Fischler - put and get to/from strea 11 // M Fischler - put and get to/from streams 12/13/04 12 // M Fischler - put/get to/from streams 12 // M Fischler - put/get to/from streams uses pairs of ulongs when 13 // + storing doubles avoid problems with 13 // + storing doubles avoid problems with precision 14 // 4/14/05 14 // 4/14/05 15 // =========================================== 15 // ======================================================================= 16 16 17 #include "CLHEP/Random/RandGamma.h" 17 #include "CLHEP/Random/RandGamma.h" 18 #include "CLHEP/Random/DoubConv.h" 18 #include "CLHEP/Random/DoubConv.h" 19 #include "CLHEP/Utility/thread_local.h" 19 #include "CLHEP/Utility/thread_local.h" 20 #include <cmath> // for std::log() 20 #include <cmath> // for std::log() 21 #include <iostream> 21 #include <iostream> 22 #include <string> 22 #include <string> 23 #include <vector> 23 #include <vector> 24 24 25 namespace CLHEP { 25 namespace CLHEP { 26 26 27 std::string RandGamma::name() const {return "R 27 std::string RandGamma::name() const {return "RandGamma";} 28 HepRandomEngine & RandGamma::engine() {return 28 HepRandomEngine & RandGamma::engine() {return *localEngine;} 29 29 30 RandGamma::~RandGamma() { 30 RandGamma::~RandGamma() { 31 } 31 } 32 32 33 double RandGamma::shoot( HepRandomEngine *anEn 33 double RandGamma::shoot( HepRandomEngine *anEngine, double k, 34 34 double lambda ) { 35 return genGamma( anEngine, k, lambda ); 35 return genGamma( anEngine, k, lambda ); 36 } 36 } 37 37 38 double RandGamma::shoot( double k, double lamb 38 double RandGamma::shoot( double k, double lambda ) { 39 HepRandomEngine *anEngine = HepRandom::getTh 39 HepRandomEngine *anEngine = HepRandom::getTheEngine(); 40 return genGamma( anEngine, k, lambda ); 40 return genGamma( anEngine, k, lambda ); 41 } 41 } 42 42 43 double RandGamma::fire( double k, double lambd 43 double RandGamma::fire( double k, double lambda ) { 44 return genGamma( localEngine.get(), k, lambd 44 return genGamma( localEngine.get(), k, lambda ); 45 } 45 } 46 46 47 void RandGamma::shootArray( const int size, do 47 void RandGamma::shootArray( const int size, double* vect, 48 double k, double l 48 double k, double lambda ) 49 { 49 { 50 for( double* v = vect; v != vect + size; ++v 50 for( double* v = vect; v != vect + size; ++v ) 51 *v = shoot(k,lambda); 51 *v = shoot(k,lambda); 52 } 52 } 53 53 54 void RandGamma::shootArray( HepRandomEngine* a 54 void RandGamma::shootArray( HepRandomEngine* anEngine, 55 const int size, do 55 const int size, double* vect, 56 double k, double l 56 double k, double lambda ) 57 { 57 { 58 for( double* v = vect; v != vect + size; ++v 58 for( double* v = vect; v != vect + size; ++v ) 59 *v = shoot(anEngine,k,lambda); 59 *v = shoot(anEngine,k,lambda); 60 } 60 } 61 61 62 void RandGamma::fireArray( const int size, dou 62 void RandGamma::fireArray( const int size, double* vect) 63 { 63 { 64 for( double* v = vect; v != vect + size; ++v 64 for( double* v = vect; v != vect + size; ++v ) 65 *v = fire(defaultK,defaultLambda); 65 *v = fire(defaultK,defaultLambda); 66 } 66 } 67 67 68 void RandGamma::fireArray( const int size, dou 68 void RandGamma::fireArray( const int size, double* vect, 69 double k, double la 69 double k, double lambda ) 70 { 70 { 71 for( double* v = vect; v != vect + size; ++v 71 for( double* v = vect; v != vect + size; ++v ) 72 *v = fire(k,lambda); 72 *v = fire(k,lambda); 73 } 73 } 74 74 75 double RandGamma::genGamma( HepRandomEngine *a 75 double RandGamma::genGamma( HepRandomEngine *anEngine, 76 double a, doubl 76 double a, double lambda ) { 77 /********************************************* 77 /************************************************************************* 78 * Gamma Distribution - Rejection algo 78 * Gamma Distribution - Rejection algorithm gs combined with * 79 * Acceptance com 79 * Acceptance complement method gd * 80 ********************************************* 80 *************************************************************************/ 81 81 82 double aa = -1.0, aaa = -1.0, b{0.}, c{0.}, << 82 static CLHEP_THREAD_LOCAL double aa = -1.0, aaa = -1.0, b, c, d, e, r, s, si, ss, q0; 83 constexpr double q1 = 0.0416666664, q2 = 0. << 83 static const double q1 = 0.0416666664, q2 = 0.0208333723, q3 = 0.0079849875, 84 q4 = 0.0015746717, q5 = -0.0003349403, 84 q4 = 0.0015746717, q5 = -0.0003349403, q6 = 0.0003340332, 85 q7 = 0.0006053049, q8 = -0.0004701849, 85 q7 = 0.0006053049, q8 = -0.0004701849, q9 = 0.0001710320, 86 a1 = 0.333333333, a2 = -0.249999949, 86 a1 = 0.333333333, a2 = -0.249999949, a3 = 0.199999867, 87 a4 =-0.166677482, a5 = 0.142873973, 87 a4 =-0.166677482, a5 = 0.142873973, a6 =-0.124385581, 88 a7 = 0.110368310, a8 = -0.112750886, 88 a7 = 0.110368310, a8 = -0.112750886, a9 = 0.104089866, 89 e1 = 1.000000000, e2 = 0.499999994, 89 e1 = 1.000000000, e2 = 0.499999994, e3 = 0.166666848, 90 e4 = 0.041664508, e5 = 0.008345522, 90 e4 = 0.041664508, e5 = 0.008345522, e6 = 0.001353826, 91 e7 = 0.000247453; 91 e7 = 0.000247453; 92 92 93 double gds{0.},p{0.},q{0.},t{0.},sign_u{0.},u << 93 double gds,p,q,t,sign_u,u,v,w,x; 94 double v1{0.},v2{0.},v12{0.}; << 94 double v1,v2,v12; 95 95 96 // Check for invalid input values 96 // Check for invalid input values 97 97 98 if( a <= 0.0 ) return (-1.0); 98 if( a <= 0.0 ) return (-1.0); 99 if( lambda <= 0.0 ) return (-1.0); 99 if( lambda <= 0.0 ) return (-1.0); 100 100 101 if (a < 1.0) 101 if (a < 1.0) 102 { // CASE A: Acceptance rejection 102 { // CASE A: Acceptance rejection algorithm gs 103 b = 1.0 + 0.36788794412 * a; // Step 103 b = 1.0 + 0.36788794412 * a; // Step 1 104 for(;;) 104 for(;;) 105 { 105 { 106 p = b * anEngine->flat(); 106 p = b * anEngine->flat(); 107 if (p <= 1.0) 107 if (p <= 1.0) 108 { // Step 2. Ca 108 { // Step 2. Case gds <= 1 109 gds = std::exp(std::log(p) / a); 109 gds = std::exp(std::log(p) / a); 110 if (std::log(anEngine->flat()) <= -gds) r 110 if (std::log(anEngine->flat()) <= -gds) return(gds/lambda); 111 } 111 } 112 else 112 else 113 { // Step 3. Ca 113 { // Step 3. Case gds > 1 114 gds = - std::log ((b - p) / a); 114 gds = - std::log ((b - p) / a); 115 if (std::log(anEngine->flat()) <= ((a - 1 115 if (std::log(anEngine->flat()) <= ((a - 1.0) * std::log(gds))) return(gds/lambda); 116 } 116 } 117 } 117 } 118 } 118 } 119 else 119 else 120 { // CASE B: Acceptance complement 120 { // CASE B: Acceptance complement algorithm gd 121 if (a != aa) 121 if (a != aa) 122 { // Step 122 { // Step 1. Preparations 123 aa = a; 123 aa = a; 124 ss = a - 0.5; 124 ss = a - 0.5; 125 s = std::sqrt(ss); 125 s = std::sqrt(ss); 126 d = 5.656854249 - 12.0 * s; 126 d = 5.656854249 - 12.0 * s; 127 } 127 } 128 128 // Step 2. Normal deviate 129 do { 129 do { 130 v1 = 2.0 * anEngine->flat() - 1.0; 130 v1 = 2.0 * anEngine->flat() - 1.0; 131 v2 = 2.0 * anEngine->flat() - 1.0; 131 v2 = 2.0 * anEngine->flat() - 1.0; 132 v12 = v1*v1 + v2*v2; 132 v12 = v1*v1 + v2*v2; 133 } while ( v12 > 1.0 ); 133 } while ( v12 > 1.0 ); 134 t = v1*std::sqrt(-2.0*std::log(v12)/v12); 134 t = v1*std::sqrt(-2.0*std::log(v12)/v12); 135 x = s + 0.5 * t; 135 x = s + 0.5 * t; 136 gds = x * x; 136 gds = x * x; 137 if (t >= 0.0) return(gds/lambda); 137 if (t >= 0.0) return(gds/lambda); // Immediate acceptance 138 138 139 u = anEngine->flat(); // Step 3 139 u = anEngine->flat(); // Step 3. Uniform random number 140 if (d * u <= t * t * t) return(gds/lambda) 140 if (d * u <= t * t * t) return(gds/lambda); // Squeeze acceptance 141 141 142 if (a != aaa) 142 if (a != aaa) 143 { // Step 143 { // Step 4. Set-up for hat case 144 aaa = a; 144 aaa = a; 145 r = 1.0 / a; 145 r = 1.0 / a; 146 q0 = ((((((((q9 * r + q8) * r + q7) * r + q6 146 q0 = ((((((((q9 * r + q8) * r + q7) * r + q6) * r + q5) * r + q4) * 147 r + q3) * r + q2) * r + q1) * r; 147 r + q3) * r + q2) * r + q1) * r; 148 if (a > 3.686) 148 if (a > 3.686) 149 { 149 { 150 if (a > 13.022) 150 if (a > 13.022) 151 { 151 { 152 b = 1.77; 152 b = 1.77; 153 si = 0.75; 153 si = 0.75; 154 c = 0.1515 / s; 154 c = 0.1515 / s; 155 } 155 } 156 else 156 else 157 { 157 { 158 b = 1.654 + 0.0076 * ss; 158 b = 1.654 + 0.0076 * ss; 159 si = 1.68 / s + 0.275; 159 si = 1.68 / s + 0.275; 160 c = 0.062 / s + 0.024; 160 c = 0.062 / s + 0.024; 161 } 161 } 162 } 162 } 163 else 163 else 164 { 164 { 165 b = 0.463 + s - 0.178 * ss; 165 b = 0.463 + s - 0.178 * ss; 166 si = 1.235; 166 si = 1.235; 167 c = 0.195 / s - 0.079 + 0.016 * s; 167 c = 0.195 / s - 0.079 + 0.016 * s; 168 } 168 } 169 } 169 } 170 if (x > 0.0) // Step 170 if (x > 0.0) // Step 5. Calculation of q 171 { 171 { 172 v = t / (s + s); // Step 6. 172 v = t / (s + s); // Step 6. 173 if (std::fabs(v) > 0.25) 173 if (std::fabs(v) > 0.25) 174 { 174 { 175 q = q0 - s * t + 0.25 * t * t + (ss + ss 175 q = q0 - s * t + 0.25 * t * t + (ss + ss) * std::log(1.0 + v); 176 } 176 } 177 else 177 else 178 { 178 { 179 q = q0 + 0.5 * t * t * ((((((((a9 * v + 179 q = q0 + 0.5 * t * t * ((((((((a9 * v + a8) * v + a7) * v + a6) * 180 v + a5) * v + a4) * v + a3) * v + a2 180 v + a5) * v + a4) * v + a3) * v + a2) * v + a1) * v; 181 } // Step 7. Quotient acce 181 } // Step 7. Quotient acceptance 182 if (std::log(1.0 - u) <= q) return(gds/lambd 182 if (std::log(1.0 - u) <= q) return(gds/lambda); 183 } 183 } 184 184 185 for(;;) 185 for(;;) 186 { // Step 8. Double 186 { // Step 8. Double exponential deviate t 187 do 187 do 188 { 188 { 189 e = -std::log(anEngine->flat()); 189 e = -std::log(anEngine->flat()); 190 u = anEngine->flat(); 190 u = anEngine->flat(); 191 u = u + u - 1.0; 191 u = u + u - 1.0; 192 sign_u = (u > 0)? 1.0 : -1.0; 192 sign_u = (u > 0)? 1.0 : -1.0; 193 t = b + (e * si) * sign_u; 193 t = b + (e * si) * sign_u; 194 } 194 } 195 while (t <= -0.71874483771719); // Step 9. 195 while (t <= -0.71874483771719); // Step 9. Rejection of t 196 v = t / (s + s); // Step 10 196 v = t / (s + s); // Step 10. New q(t) 197 if (std::fabs(v) > 0.25) 197 if (std::fabs(v) > 0.25) 198 { 198 { 199 q = q0 - s * t + 0.25 * t * t + (ss + ss 199 q = q0 - s * t + 0.25 * t * t + (ss + ss) * std::log(1.0 + v); 200 } 200 } 201 else 201 else 202 { 202 { 203 q = q0 + 0.5 * t * t * ((((((((a9 * v + 203 q = q0 + 0.5 * t * t * ((((((((a9 * v + a8) * v + a7) * v + a6) * 204 v + a5) * v + a4) * v + a3) * v + a2 204 v + a5) * v + a4) * v + a3) * v + a2) * v + a1) * v; 205 } 205 } 206 if (q <= 0.0) continue; // Step 11 206 if (q <= 0.0) continue; // Step 11. 207 if (q > 0.5) 207 if (q > 0.5) 208 { 208 { 209 w = std::exp(q) - 1.0; 209 w = std::exp(q) - 1.0; 210 } 210 } 211 else 211 else 212 { 212 { 213 w = ((((((e7 * q + e6) * q + e5) * q + e 213 w = ((((((e7 * q + e6) * q + e5) * q + e4) * q + e3) * q + e2) * 214 q + e1) * q; 214 q + e1) * q; 215 } // Step 12. Hat acce 215 } // Step 12. Hat acceptance 216 if ( c * u * sign_u <= w * std::exp(e - 0.5 216 if ( c * u * sign_u <= w * std::exp(e - 0.5 * t * t)) 217 { 217 { 218 x = s + 0.5 * t; 218 x = s + 0.5 * t; 219 return(x*x/lambda); 219 return(x*x/lambda); 220 } 220 } 221 } 221 } 222 } 222 } 223 } 223 } 224 224 225 std::ostream & RandGamma::put ( std::ostream & 225 std::ostream & RandGamma::put ( std::ostream & os ) const { 226 long pr=os.precision(20); << 226 int pr=os.precision(20); 227 std::vector<unsigned long> t(2); 227 std::vector<unsigned long> t(2); 228 os << " " << name() << "\n"; 228 os << " " << name() << "\n"; 229 os << "Uvec" << "\n"; 229 os << "Uvec" << "\n"; 230 t = DoubConv::dto2longs(defaultK); 230 t = DoubConv::dto2longs(defaultK); 231 os << defaultK << " " << t[0] << " " << t[1] 231 os << defaultK << " " << t[0] << " " << t[1] << "\n"; 232 t = DoubConv::dto2longs(defaultLambda); 232 t = DoubConv::dto2longs(defaultLambda); 233 os << defaultLambda << " " << t[0] << " " << 233 os << defaultLambda << " " << t[0] << " " << t[1] << "\n"; 234 os.precision(pr); 234 os.precision(pr); 235 return os; 235 return os; 236 } 236 } 237 237 238 std::istream & RandGamma::get ( std::istream & 238 std::istream & RandGamma::get ( std::istream & is ) { 239 std::string inName; 239 std::string inName; 240 is >> inName; 240 is >> inName; 241 if (inName != name()) { 241 if (inName != name()) { 242 is.clear(std::ios::badbit | is.rdstate()); 242 is.clear(std::ios::badbit | is.rdstate()); 243 std::cerr << "Mismatch when expecting to r 243 std::cerr << "Mismatch when expecting to read state of a " 244 << name() << " distribution\n" 244 << name() << " distribution\n" 245 << "Name found was " << inName 245 << "Name found was " << inName 246 << "\nistream is left in the badbit st 246 << "\nistream is left in the badbit state\n"; 247 return is; 247 return is; 248 } 248 } 249 if (possibleKeywordInput(is, "Uvec", default 249 if (possibleKeywordInput(is, "Uvec", defaultK)) { 250 std::vector<unsigned long> t(2); 250 std::vector<unsigned long> t(2); 251 is >> defaultK >> t[0] >> t[1]; defaultK = 251 is >> defaultK >> t[0] >> t[1]; defaultK = DoubConv::longs2double(t); 252 is >> defaultLambda>>t[0]>>t[1]; defaultLa 252 is >> defaultLambda>>t[0]>>t[1]; defaultLambda = DoubConv::longs2double(t); 253 return is; 253 return is; 254 } 254 } 255 // is >> defaultK encompassed by possibleKey 255 // is >> defaultK encompassed by possibleKeywordInput 256 is >> defaultLambda; 256 is >> defaultLambda; 257 return is; 257 return is; 258 } 258 } 259 259 260 } // namespace CLHEP 260 } // namespace CLHEP 261 261 262 262