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