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