33 return genGamma( anEngine, k, lambda );
38 return genGamma( anEngine, k, lambda );
48 for(
double* v = vect; v != vect + size; ++v )
53 const int size,
double* vect,
56 for(
double* v = vect; v != vect + size; ++v )
57 *v =
shoot(anEngine,k,lambda);
62 for(
double* v = vect; v != vect + size; ++v )
69 for(
double* v = vect; v != vect + size; ++v )
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))
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";
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";
248 std::vector<unsigned long> t(2);
bool possibleKeywordInput(IS &is, const std::string &key, T &t)
std::vector< ExP01TrackerHit * > a
static HepRandomEngine * getTheEngine()
std::istream & get(std::istream &is)
HepRandomEngine & engine()
static double longs2double(const std::vector< unsigned long > &v)
void fireArray(const int size, double *vect)
std::ostream & put(std::ostream &os) const
static double genGamma(HepRandomEngine *anEngine, double k, double lambda)
static constexpr double s
static std::vector< unsigned long > dto2longs(double d)
#define CLHEP_THREAD_LOCAL
std::shared_ptr< HepRandomEngine > localEngine
static void shootArray(const int size, double *vect, double k=1.0, double lambda=1.0)