00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018 #include "CLHEP/Random/RandGamma.h"
00019 #include "CLHEP/Random/DoubConv.h"
00020 #include <cmath>
00021
00022 namespace CLHEP {
00023
00024 std::string RandGamma::name() const {return "RandGamma";}
00025 HepRandomEngine & RandGamma::engine() {return *localEngine;}
00026
00027 RandGamma::~RandGamma() {
00028 }
00029
00030 double RandGamma::shoot( HepRandomEngine *anEngine, double k,
00031 double lambda ) {
00032 return genGamma( anEngine, k, lambda );
00033 }
00034
00035 double RandGamma::shoot( double k, double lambda ) {
00036 HepRandomEngine *anEngine = HepRandom::getTheEngine();
00037 return genGamma( anEngine, k, lambda );
00038 }
00039
00040 double RandGamma::fire( double k, double lambda ) {
00041 return genGamma( localEngine.get(), k, lambda );
00042 }
00043
00044 void RandGamma::shootArray( const int size, double* vect,
00045 double k, double lambda )
00046 {
00047 for( double* v = vect; v != vect + size; ++v )
00048 *v = shoot(k,lambda);
00049 }
00050
00051 void RandGamma::shootArray( HepRandomEngine* anEngine,
00052 const int size, double* vect,
00053 double k, double lambda )
00054 {
00055 for( double* v = vect; v != vect + size; ++v )
00056 *v = shoot(anEngine,k,lambda);
00057 }
00058
00059 void RandGamma::fireArray( const int size, double* vect)
00060 {
00061 for( double* v = vect; v != vect + size; ++v )
00062 *v = fire(defaultK,defaultLambda);
00063 }
00064
00065 void RandGamma::fireArray( const int size, double* vect,
00066 double k, double lambda )
00067 {
00068 for( double* v = vect; v != vect + size; ++v )
00069 *v = fire(k,lambda);
00070 }
00071
00072 double RandGamma::genGamma( HepRandomEngine *anEngine,
00073 double a, double lambda ) {
00074
00075
00076
00077
00078
00079 static double aa = -1.0, aaa = -1.0, b, c, d, e, r, s, si, ss, q0,
00080 q1 = 0.0416666664, q2 = 0.0208333723, q3 = 0.0079849875,
00081 q4 = 0.0015746717, q5 = -0.0003349403, q6 = 0.0003340332,
00082 q7 = 0.0006053049, q8 = -0.0004701849, q9 = 0.0001710320,
00083 a1 = 0.333333333, a2 = -0.249999949, a3 = 0.199999867,
00084 a4 =-0.166677482, a5 = 0.142873973, a6 =-0.124385581,
00085 a7 = 0.110368310, a8 = -0.112750886, a9 = 0.104089866,
00086 e1 = 1.000000000, e2 = 0.499999994, e3 = 0.166666848,
00087 e4 = 0.041664508, e5 = 0.008345522, e6 = 0.001353826,
00088 e7 = 0.000247453;
00089
00090 double gds,p,q,t,sign_u,u,v,w,x;
00091 double v1,v2,v12;
00092
00093
00094
00095 if( a <= 0.0 ) return (-1.0);
00096 if( lambda <= 0.0 ) return (-1.0);
00097
00098 if (a < 1.0)
00099 {
00100 b = 1.0 + 0.36788794412 * a;
00101 for(;;)
00102 {
00103 p = b * anEngine->flat();
00104 if (p <= 1.0)
00105 {
00106 gds = std::exp(std::log(p) / a);
00107 if (std::log(anEngine->flat()) <= -gds) return(gds/lambda);
00108 }
00109 else
00110 {
00111 gds = - std::log ((b - p) / a);
00112 if (std::log(anEngine->flat()) <= ((a - 1.0) * std::log(gds))) return(gds/lambda);
00113 }
00114 }
00115 }
00116 else
00117 {
00118 if (a != aa)
00119 {
00120 aa = a;
00121 ss = a - 0.5;
00122 s = std::sqrt(ss);
00123 d = 5.656854249 - 12.0 * s;
00124 }
00125
00126 do {
00127 v1 = 2.0 * anEngine->flat() - 1.0;
00128 v2 = 2.0 * anEngine->flat() - 1.0;
00129 v12 = v1*v1 + v2*v2;
00130 } while ( v12 > 1.0 );
00131 t = v1*std::sqrt(-2.0*std::log(v12)/v12);
00132 x = s + 0.5 * t;
00133 gds = x * x;
00134 if (t >= 0.0) return(gds/lambda);
00135
00136 u = anEngine->flat();
00137 if (d * u <= t * t * t) return(gds/lambda);
00138
00139 if (a != aaa)
00140 {
00141 aaa = a;
00142 r = 1.0 / a;
00143 q0 = ((((((((q9 * r + q8) * r + q7) * r + q6) * r + q5) * r + q4) *
00144 r + q3) * r + q2) * r + q1) * r;
00145 if (a > 3.686)
00146 {
00147 if (a > 13.022)
00148 {
00149 b = 1.77;
00150 si = 0.75;
00151 c = 0.1515 / s;
00152 }
00153 else
00154 {
00155 b = 1.654 + 0.0076 * ss;
00156 si = 1.68 / s + 0.275;
00157 c = 0.062 / s + 0.024;
00158 }
00159 }
00160 else
00161 {
00162 b = 0.463 + s - 0.178 * ss;
00163 si = 1.235;
00164 c = 0.195 / s - 0.079 + 0.016 * s;
00165 }
00166 }
00167 if (x > 0.0)
00168 {
00169 v = t / (s + s);
00170 if (std::fabs(v) > 0.25)
00171 {
00172 q = q0 - s * t + 0.25 * t * t + (ss + ss) * std::log(1.0 + v);
00173 }
00174 else
00175 {
00176 q = q0 + 0.5 * t * t * ((((((((a9 * v + a8) * v + a7) * v + a6) *
00177 v + a5) * v + a4) * v + a3) * v + a2) * v + a1) * v;
00178 }
00179 if (std::log(1.0 - u) <= q) return(gds/lambda);
00180 }
00181
00182 for(;;)
00183 {
00184 do
00185 {
00186 e = -std::log(anEngine->flat());
00187 u = anEngine->flat();
00188 u = u + u - 1.0;
00189 sign_u = (u > 0)? 1.0 : -1.0;
00190 t = b + (e * si) * sign_u;
00191 }
00192 while (t <= -0.71874483771719);
00193 v = t / (s + s);
00194 if (std::fabs(v) > 0.25)
00195 {
00196 q = q0 - s * t + 0.25 * t * t + (ss + ss) * std::log(1.0 + v);
00197 }
00198 else
00199 {
00200 q = q0 + 0.5 * t * t * ((((((((a9 * v + a8) * v + a7) * v + a6) *
00201 v + a5) * v + a4) * v + a3) * v + a2) * v + a1) * v;
00202 }
00203 if (q <= 0.0) continue;
00204 if (q > 0.5)
00205 {
00206 w = std::exp(q) - 1.0;
00207 }
00208 else
00209 {
00210 w = ((((((e7 * q + e6) * q + e5) * q + e4) * q + e3) * q + e2) *
00211 q + e1) * q;
00212 }
00213 if ( c * u * sign_u <= w * std::exp(e - 0.5 * t * t))
00214 {
00215 x = s + 0.5 * t;
00216 return(x*x/lambda);
00217 }
00218 }
00219 }
00220 }
00221
00222 std::ostream & RandGamma::put ( std::ostream & os ) const {
00223 int pr=os.precision(20);
00224 std::vector<unsigned long> t(2);
00225 os << " " << name() << "\n";
00226 os << "Uvec" << "\n";
00227 t = DoubConv::dto2longs(defaultK);
00228 os << defaultK << " " << t[0] << " " << t[1] << "\n";
00229 t = DoubConv::dto2longs(defaultLambda);
00230 os << defaultLambda << " " << t[0] << " " << t[1] << "\n";
00231 os.precision(pr);
00232 return os;
00233 }
00234
00235 std::istream & RandGamma::get ( std::istream & is ) {
00236 std::string inName;
00237 is >> inName;
00238 if (inName != name()) {
00239 is.clear(std::ios::badbit | is.rdstate());
00240 std::cerr << "Mismatch when expecting to read state of a "
00241 << name() << " distribution\n"
00242 << "Name found was " << inName
00243 << "\nistream is left in the badbit state\n";
00244 return is;
00245 }
00246 if (possibleKeywordInput(is, "Uvec", defaultK)) {
00247 std::vector<unsigned long> t(2);
00248 is >> defaultK >> t[0] >> t[1]; defaultK = DoubConv::longs2double(t);
00249 is >> defaultLambda>>t[0]>>t[1]; defaultLambda = DoubConv::longs2double(t);
00250 return is;
00251 }
00252
00253 is >> defaultLambda;
00254 return is;
00255 }
00256
00257 }
00258