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