18 #include "CLHEP/Random/RandGamma.h"
19 #include "CLHEP/Random/DoubConv.h"
20 #include "CLHEP/Utility/thread_local.h"
28 RandGamma::~RandGamma() {
33 return genGamma( anEngine, k, lambda );
37 HepRandomEngine *anEngine = HepRandom::getTheEngine();
38 return genGamma( anEngine, k, lambda );
41 double RandGamma::fire(
double k,
double lambda ) {
42 return genGamma( localEngine.get(), k,
lambda );
45 void RandGamma::shootArray(
const int size,
double* vect,
46 double k,
double lambda )
48 for(
double* v = vect; v != vect + size; ++v )
52 void RandGamma::shootArray( HepRandomEngine* anEngine,
53 const int size,
double* vect,
54 double k,
double lambda )
56 for(
double* v = vect; v != vect + size; ++v )
57 *v =
shoot(anEngine,k,lambda);
60 void RandGamma::fireArray(
const int size,
double* vect)
62 for(
double* v = vect; v != vect + size; ++v )
63 *v = fire(defaultK,defaultLambda);
66 void RandGamma::fireArray(
const int size,
double* vect,
67 double k,
double lambda )
69 for(
double* v = vect; v != vect + size; ++v )
73 double RandGamma::genGamma( HepRandomEngine *anEngine,
74 double a,
double lambda ) {
80 static CLHEP_THREAD_LOCAL
double aa = -1.0, aaa = -1.0, b, c, d, e, r,
s, si, ss, q0;
81 static const double 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))
223 std::ostream & RandGamma::put ( std::ostream & os )
const {
224 int pr=os.precision(20);
225 std::vector<unsigned long> t(2);
226 os <<
" " <<
name() <<
"\n";
227 os <<
"Uvec" <<
"\n";
228 t = DoubConv::dto2longs(defaultK);
229 os << defaultK <<
" " << t[0] <<
" " << t[1] <<
"\n";
230 t = DoubConv::dto2longs(defaultLambda);
231 os << defaultLambda <<
" " << t[0] <<
" " << t[1] <<
"\n";
236 std::istream & RandGamma::get ( std::istream & is ) {
239 if (inName !=
name()) {
240 is.clear(std::ios::badbit | is.rdstate());
241 std::cerr <<
"Mismatch when expecting to read state of a "
242 <<
name() <<
" distribution\n"
243 <<
"Name found was " << inName
244 <<
"\nistream is left in the badbit state\n";
247 if (possibleKeywordInput(is,
"Uvec", defaultK)) {
248 std::vector<unsigned long> t(2);
249 is >> defaultK >> t[0] >> t[1]; defaultK = DoubConv::longs2double(t);
250 is >> defaultLambda>>t[0]>>t[1]; defaultLambda = DoubConv::longs2double(t);
ThreeVector shoot(const G4int Ap, const G4int Af)
const G4double w[NPOINTSGL]
static int engine(pchar, pchar, double &, pchar &, const dic_type &)
const G4double x[NPOINTSGL]