Geant4  10.01.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 "CLHEP/Utility/thread_local.h"
21 #include <cmath> // for std::log()
22 
23 namespace CLHEP {
24 
25 std::string RandGamma::name() const {return "RandGamma";}
26 HepRandomEngine & RandGamma::engine() {return *localEngine;}
27 
28 RandGamma::~RandGamma() {
29 }
30 
31 double RandGamma::shoot( HepRandomEngine *anEngine, double k,
32  double lambda ) {
33  return genGamma( anEngine, k, lambda );
34 }
35 
36 double RandGamma::shoot( double k, double lambda ) {
37  HepRandomEngine *anEngine = HepRandom::getTheEngine();
38  return genGamma( anEngine, k, lambda );
39 }
40 
41 double RandGamma::fire( double k, double lambda ) {
42  return genGamma( localEngine.get(), k, lambda );
43 }
44 
45 void RandGamma::shootArray( const int size, double* vect,
46  double k, double lambda )
47 {
48  for( double* v = vect; v != vect + size; ++v )
49  *v = shoot(k,lambda);
50 }
51 
52 void RandGamma::shootArray( HepRandomEngine* anEngine,
53  const int size, double* vect,
54  double k, double lambda )
55 {
56  for( double* v = vect; v != vect + size; ++v )
57  *v = shoot(anEngine,k,lambda);
58 }
59 
60 void RandGamma::fireArray( const int size, double* vect)
61 {
62  for( double* v = vect; v != vect + size; ++v )
63  *v = fire(defaultK,defaultLambda);
64 }
65 
66 void RandGamma::fireArray( const int size, double* vect,
67  double k, double lambda )
68 {
69  for( double* v = vect; v != vect + size; ++v )
70  *v = fire(k,lambda);
71 }
72 
73 double RandGamma::genGamma( HepRandomEngine *anEngine,
74  double a, double lambda ) {
75 /*************************************************************************
76  * Gamma Distribution - Rejection algorithm gs combined with *
77  * Acceptance complement method gd *
78  *************************************************************************/
79 
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,
89  e7 = 0.000247453;
90 
91 double gds,p,q,t,sign_u,u,v,w,x;
92 double v1,v2,v12;
93 
94 // Check for invalid input values
95 
96  if( a <= 0.0 ) return (-1.0);
97  if( lambda <= 0.0 ) return (-1.0);
98 
99  if (a < 1.0)
100  { // CASE A: Acceptance rejection algorithm gs
101  b = 1.0 + 0.36788794412 * a; // Step 1
102  for(;;)
103  {
104  p = b * anEngine->flat();
105  if (p <= 1.0)
106  { // Step 2. Case gds <= 1
107  gds = std::exp(std::log(p) / a);
108  if (std::log(anEngine->flat()) <= -gds) return(gds/lambda);
109  }
110  else
111  { // Step 3. Case gds > 1
112  gds = - std::log ((b - p) / a);
113  if (std::log(anEngine->flat()) <= ((a - 1.0) * std::log(gds))) return(gds/lambda);
114  }
115  }
116  }
117  else
118  { // CASE B: Acceptance complement algorithm gd
119  if (a != aa)
120  { // Step 1. Preparations
121  aa = a;
122  ss = a - 0.5;
123  s = std::sqrt(ss);
124  d = 5.656854249 - 12.0 * s;
125  }
126  // Step 2. Normal deviate
127  do {
128  v1 = 2.0 * anEngine->flat() - 1.0;
129  v2 = 2.0 * anEngine->flat() - 1.0;
130  v12 = v1*v1 + v2*v2;
131  } while ( v12 > 1.0 );
132  t = v1*std::sqrt(-2.0*std::log(v12)/v12);
133  x = s + 0.5 * t;
134  gds = x * x;
135  if (t >= 0.0) return(gds/lambda); // Immediate acceptance
136 
137  u = anEngine->flat(); // Step 3. Uniform random number
138  if (d * u <= t * t * t) return(gds/lambda); // Squeeze acceptance
139 
140  if (a != aaa)
141  { // Step 4. Set-up for hat case
142  aaa = a;
143  r = 1.0 / a;
144  q0 = ((((((((q9 * r + q8) * r + q7) * r + q6) * r + q5) * r + q4) *
145  r + q3) * r + q2) * r + q1) * r;
146  if (a > 3.686)
147  {
148  if (a > 13.022)
149  {
150  b = 1.77;
151  si = 0.75;
152  c = 0.1515 / s;
153  }
154  else
155  {
156  b = 1.654 + 0.0076 * ss;
157  si = 1.68 / s + 0.275;
158  c = 0.062 / s + 0.024;
159  }
160  }
161  else
162  {
163  b = 0.463 + s - 0.178 * ss;
164  si = 1.235;
165  c = 0.195 / s - 0.079 + 0.016 * s;
166  }
167  }
168  if (x > 0.0) // Step 5. Calculation of q
169  {
170  v = t / (s + s); // Step 6.
171  if (std::fabs(v) > 0.25)
172  {
173  q = q0 - s * t + 0.25 * t * t + (ss + ss) * std::log(1.0 + v);
174  }
175  else
176  {
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;
179  } // Step 7. Quotient acceptance
180  if (std::log(1.0 - u) <= q) return(gds/lambda);
181  }
182 
183  for(;;)
184  { // Step 8. Double exponential deviate t
185  do
186  {
187  e = -std::log(anEngine->flat());
188  u = anEngine->flat();
189  u = u + u - 1.0;
190  sign_u = (u > 0)? 1.0 : -1.0;
191  t = b + (e * si) * sign_u;
192  }
193  while (t <= -0.71874483771719); // Step 9. Rejection of t
194  v = t / (s + s); // Step 10. New q(t)
195  if (std::fabs(v) > 0.25)
196  {
197  q = q0 - s * t + 0.25 * t * t + (ss + ss) * std::log(1.0 + v);
198  }
199  else
200  {
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;
203  }
204  if (q <= 0.0) continue; // Step 11.
205  if (q > 0.5)
206  {
207  w = std::exp(q) - 1.0;
208  }
209  else
210  {
211  w = ((((((e7 * q + e6) * q + e5) * q + e4) * q + e3) * q + e2) *
212  q + e1) * q;
213  } // Step 12. Hat acceptance
214  if ( c * u * sign_u <= w * std::exp(e - 0.5 * t * t))
215  {
216  x = s + 0.5 * t;
217  return(x*x/lambda);
218  }
219  }
220  }
221 }
222 
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";
232  os.precision(pr);
233  return os;
234 }
235 
236 std::istream & RandGamma::get ( std::istream & is ) {
237  std::string inName;
238  is >> inName;
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";
245  return is;
246  }
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);
251  return is;
252  }
253  // is >> defaultK encompassed by possibleKeywordInput
254  is >> defaultLambda;
255  return is;
256 }
257 
258 } // namespace CLHEP
259 
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