Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/externals/clhep/src/RandGamma.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

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