Geant4  10.01.p02
G4MTRandGamma.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // $Id:$
28 //
29 #if __clang__
30  #if ((defined(G4MULTITHREADED) && !defined(G4USE_STD11)) || \
31  !__has_feature(cxx_thread_local))
32  #define CLANG_NOSTDTLS
33  #endif
34 #endif
35 
36 #if (defined(G4MULTITHREADED) && \
37  (!defined(G4USE_STD11) || (defined(CLANG_NOSTDTLS))))
38 
39 #include <cmath> // for log()
40 
41 #include "G4MTRandGamma.hh"
42 
44  if ( deleteEngine ) delete localEngine;
45 }
46 
47 G4double G4MTRandGamma::shoot( CLHEP::HepRandomEngine *anEngine, G4double k,
48  G4double lambda ) {
49  return genGamma( anEngine, k, lambda );
50 }
51 
53  CLHEP::HepRandomEngine *anEngine = G4MTHepRandom::getTheEngine();
54  return genGamma( anEngine, k, lambda );
55 }
56 
58  return genGamma( localEngine, k, lambda );
59 }
60 
61 void G4MTRandGamma::shootArray( const G4int size, G4double* vect,
63 {
64  G4int i;
65 
66  for (i=0; i<size; ++i)
67  vect[i] = shoot(k,lambda);
68 }
69 
70 void G4MTRandGamma::shootArray( CLHEP::HepRandomEngine* anEngine,
71  const G4int size, G4double* vect,
72  G4double k, G4double lambda )
73 {
74  G4int i;
75 
76  for (i=0; i<size; ++i)
77  vect[i] = shoot(anEngine,k,lambda);
78 }
79 
80 void G4MTRandGamma::fireArray( const G4int size, G4double* vect)
81 {
82  G4int i;
83 
84  for (i=0; i<size; ++i)
85  vect[i] = fire(defaultK,defaultLambda);
86 }
87 
88 void G4MTRandGamma::fireArray( const G4int size, G4double* vect,
89  G4double k, G4double lambda )
90 {
91  G4int i;
92 
93  for (i=0; i<size; ++i)
94  vect[i] = fire(k,lambda);
95 }
96 
97 G4double G4MTRandGamma::genGamma( CLHEP::HepRandomEngine *anEngine,
98  G4double a, G4double lambda ) {
99 /*************************************************************************
100  * Gamma Distribution - Rejection algorithm gs combined with *
101  * Acceptance complement method gd *
102  *************************************************************************/
103 
104 static G4ThreadLocal G4double aa = -1.0, aaa = -1.0, b, c, d, e, r, s, si, ss, q0,
105  q1 = 0.0416666664, q2 = 0.0208333723, q3 = 0.0079849875,
106  q4 = 0.0015746717, q5 = -0.0003349403, q6 = 0.0003340332,
107  q7 = 0.0006053049, q8 = -0.0004701849, q9 = 0.0001710320,
108  a1 = 0.333333333, a2 = -0.249999949, a3 = 0.199999867,
109  a4 =-0.166677482, a5 = 0.142873973, a6 =-0.124385581,
110  a7 = 0.110368310, a8 = -0.112750886, a9 = 0.104089866,
111  e1 = 1.000000000, e2 = 0.499999994, e3 = 0.166666848,
112  e4 = 0.041664508, e5 = 0.008345522, e6 = 0.001353826,
113  e7 = 0.000247453;
114 
115 G4double gds,p,q,t,sign_u,u,v,w,x;
116 G4double v1,v2,v12;
117 
118 // Check for invalid input values
119 
120  if( a <= 0.0 ) return (-1.0);
121  if( lambda <= 0.0 ) return (-1.0);
122 
123  if (a < 1.0)
124  { // CASE A: Acceptance rejection algorithm gs
125  b = 1.0 + 0.36788794412 * a; // Step 1
126  for(;;)
127  {
128  p = b * anEngine->flat();
129  if (p <= 1.0)
130  { // Step 2. Case gds <= 1
131  gds = std::exp(std::log(p) / a);
132  if (std::log(anEngine->flat()) <= -gds) return(gds/lambda);
133  }
134  else
135  { // Step 3. Case gds > 1
136  gds = - std::log ((b - p) / a);
137  if (std::log(anEngine->flat()) <= ((a - 1.0) * std::log(gds))) return(gds/lambda);
138  }
139  }
140  }
141  else
142  { // CASE B: Acceptance complement algorithm gd
143  if (a != aa)
144  { // Step 1. Preparations
145  aa = a;
146  ss = a - 0.5;
147  s = std::sqrt(ss);
148  d = 5.656854249 - 12.0 * s;
149  }
150  // Step 2. Normal deviate
151  do {
152  v1 = 2.0 * anEngine->flat() - 1.0;
153  v2 = 2.0 * anEngine->flat() - 1.0;
154  v12 = v1*v1 + v2*v2;
155  } while ( v12 > 1.0 );
156  t = v1*std::sqrt(-2.0*std::log(v12)/v12);
157  x = s + 0.5 * t;
158  gds = x * x;
159  if (t >= 0.0) return(gds/lambda); // Immediate acceptance
160 
161  u = anEngine->flat(); // Step 3. Uniform random number
162  if (d * u <= t * t * t) return(gds/lambda); // Squeeze acceptance
163 
164  if (a != aaa)
165  { // Step 4. Set-up for hat case
166  aaa = a;
167  r = 1.0 / a;
168  q0 = ((((((((q9 * r + q8) * r + q7) * r + q6) * r + q5) * r + q4) *
169  r + q3) * r + q2) * r + q1) * r;
170  if (a > 3.686)
171  {
172  if (a > 13.022)
173  {
174  b = 1.77;
175  si = 0.75;
176  c = 0.1515 / s;
177  }
178  else
179  {
180  b = 1.654 + 0.0076 * ss;
181  si = 1.68 / s + 0.275;
182  c = 0.062 / s + 0.024;
183  }
184  }
185  else
186  {
187  b = 0.463 + s - 0.178 * ss;
188  si = 1.235;
189  c = 0.195 / s - 0.079 + 0.016 * s;
190  }
191  }
192  if (x > 0.0) // Step 5. Calculation of q
193  {
194  v = t / (s + s); // Step 6.
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  } // Step 7. Quotient acceptance
204  if (std::log(1.0 - u) <= q) return(gds/lambda);
205  }
206 
207  for(;;)
208  { // Step 8. Double exponential deviate t
209  do
210  {
211  e = -std::log(anEngine->flat());
212  u = anEngine->flat();
213  u = u + u - 1.0;
214  sign_u = (u > 0)? 1.0 : -1.0;
215  t = b + (e * si) * sign_u;
216  }
217  while (t <= -0.71874483771719); // Step 9. Rejection of t
218  v = t / (s + s); // Step 10. New q(t)
219  if (std::fabs(v) > 0.25)
220  {
221  q = q0 - s * t + 0.25 * t * t + (ss + ss) * std::log(1.0 + v);
222  }
223  else
224  {
225  q = q0 + 0.5 * t * t * ((((((((a9 * v + a8) * v + a7) * v + a6) *
226  v + a5) * v + a4) * v + a3) * v + a2) * v + a1) * v;
227  }
228  if (q <= 0.0) continue; // Step 11.
229  if (q > 0.5)
230  {
231  w = std::exp(q) - 1.0;
232  }
233  else
234  {
235  w = ((((((e7 * q + e6) * q + e5) * q + e4) * q + e3) * q + e2) *
236  q + e1) * q;
237  } // Step 12. Hat acceptance
238  if ( c * u * sign_u <= w * std::exp(e - 0.5 * t * t))
239  {
240  x = s + 0.5 * t;
241  return(x*x/lambda);
242  }
243  }
244  }
245 }
246 
247 #endif
static CLHEP::HepRandomEngine * getTheEngine()
void fireArray(const G4int size, G4double *vect)
static const G4double a1
G4double fire()
static G4double shoot()
G4double defaultK
static const G4double e2
static const G4double a4
G4double a
Definition: TRTMaterials.hh:39
#define G4ThreadLocal
Definition: tls.hh:89
static const G4double e4
int G4int
Definition: G4Types.hh:78
CLHEP::HepRandomEngine * localEngine
static const double s
Definition: G4SIunits.hh:150
static const G4double a3
static const G4double e1
G4double defaultLambda
static void shootArray(const G4int size, G4double *vect, G4double k=1.0, G4double lambda=1.0)
virtual ~G4MTRandGamma()
static const G4double a5
double G4double
Definition: G4Types.hh:76
static const G4double e3
static const G4double a2
static G4double genGamma(CLHEP::HepRandomEngine *anEngine, G4double k, G4double lambda)