CLHEP VERSION Reference Documentation
   
CLHEP Home Page     CLHEP Documentation     CLHEP Bug Reports

RandGamma.cc
Go to the documentation of this file.
1 // $Id: RandGamma.cc,v 1.6 2010/06/16 17:24:53 garren Exp $
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/defs.h"
19 #include "CLHEP/Random/RandGamma.h"
20 #include "CLHEP/Random/DoubConv.hh"
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 
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 ) {
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 
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 double aa = -1.0, aaa = -1.0, b, c, d, e, r, s, si, ss, q0,
81  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 #ifdef REMOVED
235  int pr=os.precision(20);
236  os << " " << name() << "\n";
237  os << defaultK << " " << defaultLambda << "\n";
238  os.precision(pr);
239  return os;
240 #endif
241 }
242 
243 std::istream & RandGamma::get ( std::istream & is ) {
244  std::string inName;
245  is >> inName;
246  if (inName != name()) {
247  is.clear(std::ios::badbit | is.rdstate());
248  std::cerr << "Mismatch when expecting to read state of a "
249  << name() << " distribution\n"
250  << "Name found was " << inName
251  << "\nistream is left in the badbit state\n";
252  return is;
253  }
254  if (possibleKeywordInput(is, "Uvec", defaultK)) {
255  std::vector<unsigned long> t(2);
256  is >> defaultK >> t[0] >> t[1]; defaultK = DoubConv::longs2double(t);
257  is >> defaultLambda>>t[0]>>t[1]; defaultLambda = DoubConv::longs2double(t);
258  return is;
259  }
260  // is >> defaultK encompassed by possibleKeywordInput
261  is >> defaultLambda;
262  return is;
263 }
264 
265 } // namespace CLHEP
266 
a
@ a
Definition: testCategories.cc:125
CLHEP::HepRandomEngine
Definition: Matrix/CLHEP/Random/RandomEngine.h:55
b
@ b
Definition: testCategories.cc:125
CLHEP::RandGamma::fire
double fire()
CLHEP::RandGamma::put
std::ostream & put(std::ostream &os) const
Definition: RandGamma.cc:223
is
HepRotation and so forth isNear() norm2() rectify() static Rotation row1 row4(To avoid bloat in the code pulled in for programs which don 't use all these features, we split the implementation .cc files. Only isNear() goes into the original Rotation.cc) --------------------------------------- HepAxisAngle and HepEulerAngles classes --------------------------------------- These classes are very useful and simple structures for holding the result of a nice intuituve decomposition of a rotation there is no longer much content in the distinct ZOOM PhysicsVectors library The only content left in the library is the object files representing the various Exception objects When we build the CLHEP classes for the ZOOM we will set up so as to use ZOOM SpaceVector is(but we can disable namespace usage and most of our users do so at this point). What I do is leave Hep3Vector in the global namespace
CLHEP::RandGamma::shootArray
static void shootArray(const int size, double *vect, double k=1.0, double lambda=1.0)
Definition: RandGamma.cc:45
CLHEP::RandGamma::~RandGamma
virtual ~RandGamma()
Definition: RandGamma.cc:28
CLHEP::RandGamma::name
std::string name() const
Definition: RandGamma.cc:25
CLHEP::DoubConv::longs2double
static double longs2double(const std::vector< unsigned long > &v)
Definition: DoubConv.cc:106
size
user code seldom needs to call this function directly ZMerrno whether or not they are still recorded ZMerrno size() Return the(integer) number of ZMthrow 'n exceptions currently recorded. 5) ZMerrno.clear() Set an internal counter to zero. This counter is available(see next function) to user code to track ZMthrow 'n exceptions that have occurred during any arbitrary time interval. 6) ZMerrno.countSinceCleared() Return the(integer) number of ZMthrow 'n exceptions that have been recorded via ZMerrno.write()
CLHEP::RandGamma::engine
HepRandomEngine & engine()
Definition: RandGamma.cc:26
CLHEP
Definition: ClhepVersion.h:13
v
they are gone ZOOM Features Discontinued The following features of the ZOOM package were felt to be extreme overkill These have been after checking that no existing user code was utilizing as in SpaceVector v
Definition: keyMergeIssues.doc:324
CLHEP::RandGamma::fireArray
void fireArray(const int size, double *vect)
Definition: RandGamma.cc:60
CLHEP::possibleKeywordInput
bool possibleKeywordInput(IS &is, const std::string &key, T &t)
Definition: Matrix/CLHEP/Random/RandomEngine.h:168
CLHEP::HepRandomEngine::flat
virtual double flat()=0
s
Methods applicble to containers of as in std::list< LorentzVector > s
Definition: keyMergeIssues.doc:328
CLHEP::DoubConv::dto2longs
static std::vector< unsigned long > dto2longs(double d)
Definition: DoubConv.cc:90
CLHEP::RandGamma::get
std::istream & get(std::istream &is)
Definition: RandGamma.cc:243
CLHEP::HepRandom::getTheEngine
static HepRandomEngine * getTheEngine()
Definition: Random.cc:166
x
any side effects of that construction would occur twice The semantics of throw x
Definition: whyZMthrowRethrows.txt:37
k
long k
Definition: JamesRandomSeeding.txt:29
CLHEP::RandGamma::shoot
static double shoot()