18 #include "CLHEP/Random/defs.h"
19 #include "CLHEP/Random/RandGamma.h"
20 #include "CLHEP/Random/DoubConv.hh"
33 return genGamma( anEngine,
k, lambda );
38 return genGamma( anEngine,
k, lambda );
42 return genGamma( localEngine.get(),
k, lambda );
46 double k,
double lambda )
48 for(
double*
v = vect;
v != vect +
size; ++
v )
53 const int size,
double* vect,
54 double k,
double lambda )
56 for(
double*
v = vect;
v != vect +
size; ++
v )
62 for(
double*
v = vect;
v != vect +
size; ++
v )
63 *
v =
fire(defaultK,defaultLambda);
67 double k,
double lambda )
69 for(
double*
v = vect;
v != vect +
size; ++
v )
74 double a,
double lambda ) {
80 static double aa = -1.0, aaa = -1.0,
b, c, d, e, r,
s, si, ss, q0,
81 q1 = 0.0416666664, q2 = 0.0208333723, q3 = 0.0079849875,
82 q4 = 0.0015746717, q5 = -0.0003349403, q6 = 0.0003340332,
83 q7 = 0.0006053049, q8 = -0.0004701849, q9 = 0.0001710320,
84 a1 = 0.333333333, a2 = -0.249999949, a3 = 0.199999867,
85 a4 =-0.166677482, a5 = 0.142873973, a6 =-0.124385581,
86 a7 = 0.110368310, a8 = -0.112750886, a9 = 0.104089866,
87 e1 = 1.000000000, e2 = 0.499999994, e3 = 0.166666848,
88 e4 = 0.041664508, e5 = 0.008345522, e6 = 0.001353826,
91 double gds,p,q,t,sign_u,u,
v,w,
x;
96 if(
a <= 0.0 )
return (-1.0);
97 if( lambda <= 0.0 )
return (-1.0);
101 b = 1.0 + 0.36788794412 *
a;
104 p =
b * anEngine->
flat();
107 gds = std::exp(std::log(p) /
a);
108 if (std::log(anEngine->
flat()) <= -gds)
return(gds/lambda);
112 gds = - std::log ((
b - p) /
a);
113 if (std::log(anEngine->
flat()) <= ((
a - 1.0) * std::log(gds)))
return(gds/lambda);
124 d = 5.656854249 - 12.0 *
s;
128 v1 = 2.0 * anEngine->
flat() - 1.0;
129 v2 = 2.0 * anEngine->
flat() - 1.0;
131 }
while ( v12 > 1.0 );
132 t = v1*std::sqrt(-2.0*std::log(v12)/v12);
135 if (t >= 0.0)
return(gds/lambda);
137 u = anEngine->
flat();
138 if (d * u <= t * t * t)
return(gds/lambda);
144 q0 = ((((((((q9 * r + q8) * r + q7) * r + q6) * r + q5) * r + q4) *
145 r + q3) * r + q2) * r + q1) * r;
156 b = 1.654 + 0.0076 * ss;
157 si = 1.68 /
s + 0.275;
158 c = 0.062 /
s + 0.024;
163 b = 0.463 +
s - 0.178 * ss;
165 c = 0.195 /
s - 0.079 + 0.016 *
s;
171 if (std::fabs(
v) > 0.25)
173 q = q0 -
s * t + 0.25 * t * t + (ss + ss) * std::log(1.0 +
v);
177 q = q0 + 0.5 * t * t * ((((((((a9 *
v + a8) *
v + a7) *
v + a6) *
178 v + a5) *
v + a4) *
v + a3) *
v + a2) *
v + a1) *
v;
180 if (std::log(1.0 - u) <= q)
return(gds/lambda);
187 e = -std::log(anEngine->
flat());
188 u = anEngine->
flat();
190 sign_u = (u > 0)? 1.0 : -1.0;
191 t =
b + (e * si) * sign_u;
193 while (t <= -0.71874483771719);
195 if (std::fabs(
v) > 0.25)
197 q = q0 -
s * t + 0.25 * t * t + (ss + ss) * std::log(1.0 +
v);
201 q = q0 + 0.5 * t * t * ((((((((a9 *
v + a8) *
v + a7) *
v + a6) *
202 v + a5) *
v + a4) *
v + a3) *
v + a2) *
v + a1) *
v;
204 if (q <= 0.0)
continue;
207 w = std::exp(q) - 1.0;
211 w = ((((((e7 * q + e6) * q + e5) * q + e4) * q + e3) * q + e2) *
214 if ( c * u * sign_u <= w * std::exp(e - 0.5 * t * t))
224 int pr=os.precision(20);
225 std::vector<unsigned long> t(2);
226 os <<
" " <<
name() <<
"\n";
227 os <<
"Uvec" <<
"\n";
229 os << defaultK <<
" " << t[0] <<
" " << t[1] <<
"\n";
231 os << defaultLambda <<
" " << t[0] <<
" " << t[1] <<
"\n";
235 int pr=os.precision(20);
236 os <<
" " <<
name() <<
"\n";
237 os << defaultK <<
" " << defaultLambda <<
"\n";
246 if (inName !=
name()) {
247 is.clear(std::ios::badbit |
is.rdstate());
248 std::cerr <<
"Mismatch when expecting to read state of a "
249 <<
name() <<
" distribution\n"
250 <<
"Name found was " << inName
251 <<
"\nistream is left in the badbit state\n";
255 std::vector<unsigned long> t(2);