Geant4  10.00.p03
RandGamma.cc
Go to the documentation of this file.
1 // $Id:$
2 // -*- C++ -*-
3 //
4 // -----------------------------------------------------------------------
5 // HEP Random
6 // --- RandGamma ---
7 // class implementation file
8 // -----------------------------------------------------------------------
9 
10 // =======================================================================
11 // John Marraffino - Created: 12th May 1998
12 // M Fischler - put and get to/from streams 12/13/04
13 // M Fischler - put/get to/from streams uses pairs of ulongs when
14 // + storing doubles avoid problems with precision
15 // 4/14/05
16 // =======================================================================
17 
18 #include "CLHEP/Random/RandGamma.h"
19 #include "CLHEP/Random/DoubConv.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 double aa = -1.0, aaa = -1.0, b, c, d, e, r, s, si, ss, q0,
80  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 
ThreeVector shoot(const G4int Ap, const G4int Af)
G4String name
Definition: TRTMaterials.hh:40
static const G4double a1
static const G4double e2
static int engine(pchar, pchar, double &, pchar &, const dic_type &)
Definition: Evaluator.cc:358
static const G4double a4
G4double a
Definition: TRTMaterials.hh:39
static const G4double e4
static const double s
Definition: G4SIunits.hh:150
static const G4double a3
static const G4double e1
static const G4double a5
static const G4double e3
static const G4double a2