Geant4  10.02.p02
RandPoisson.cc
Go to the documentation of this file.
1 // $Id:$
2 // -*- C++ -*-
3 //
4 // -----------------------------------------------------------------------
5 // HEP Random
6 // --- RandPoisson ---
7 // class implementation file
8 // -----------------------------------------------------------------------
9 // This file is part of Geant4 (simulation toolkit for HEP).
10 
11 // =======================================================================
12 // Gabriele Cosmo - Created: 5th September 1995
13 // - Added not static Shoot() method: 17th May 1996
14 // - Algorithm now operates on doubles: 31st Oct 1996
15 // - Added methods to shoot arrays: 28th July 1997
16 // - Added check in case xm=-1: 4th February 1998
17 // J.Marraffino - Added default mean as attribute and
18 // operator() with mean: 16th Feb 1998
19 // Gabriele Cosmo - Relocated static data from HepRandom: 5th Jan 1999
20 // M Fischler - put and get to/from streams 12/15/04
21 // M Fischler - put/get to/from streams uses pairs of ulongs when
22 // + storing doubles avoid problems with precision
23 // 4/14/05
24 // Mark Fischler - Repaired BUG - when mean > 2 billion, was returning
25 // mean instead of the proper value. 01/13/06
26 // =======================================================================
27 
28 #include "CLHEP/Random/RandPoisson.h"
29 #include "CLHEP/Units/PhysicalConstants.h"
30 #include "CLHEP/Random/DoubConv.h"
31 #include <cmath> // for std::floor()
32 
33 namespace CLHEP {
34 
35 std::string RandPoisson::name() const {return "RandPoisson";}
36 HepRandomEngine & RandPoisson::engine() {return *localEngine;}
37 
38 // Initialisation of static data
39 CLHEP_THREAD_LOCAL double RandPoisson::status_st[3] = {0., 0., 0.};
40 CLHEP_THREAD_LOCAL double RandPoisson::oldm_st = -1.0;
41 const double RandPoisson::meanMax_st = 2.0E9;
42 
43 RandPoisson::~RandPoisson() {
44 }
45 
46 double RandPoisson::operator()() {
47  return double(fire( defaultMean ));
48 }
49 
50 double RandPoisson::operator()( double mean ) {
51  return double(fire( mean ));
52 }
53 
54 double gammln(double xx) {
55 
56 // Returns the value ln(Gamma(xx) for xx > 0. Full accuracy is obtained for
57 // xx > 1. For 0 < xx < 1. the reflection formula (6.1.4) can be used first.
58 // (Adapted from Numerical Recipes in C)
59 
60  static const double cof[6] = {76.18009172947146,-86.50532032941677,
61  24.01409824083091, -1.231739572450155,
62  0.1208650973866179e-2, -0.5395239384953e-5};
63  int j;
64  double x = xx - 1.0;
65  double tmp = x + 5.5;
66  tmp -= (x + 0.5) * std::log(tmp);
67  double ser = 1.000000000190015;
68 
69  for ( j = 0; j <= 5; j++ ) {
70  x += 1.0;
71  ser += cof[j]/x;
72  }
73  return -tmp + std::log(2.5066282746310005*ser);
74 }
75 
76 static
77 double normal (HepRandomEngine* eptr) // mf 1/13/06
78 {
79  double r;
80  double v1,v2,fac;
81  do {
82  v1 = 2.0 * eptr->flat() - 1.0;
83  v2 = 2.0 * eptr->flat() - 1.0;
84  r = v1*v1 + v2*v2;
85  } while ( r > 1.0 );
86 
87  fac = std::sqrt(-2.0*std::log(r)/r);
88  return v2*fac;
89 }
90 
91 long RandPoisson::shoot(double xm) {
92 
93 // Returns as a floating-point number an integer value that is a random
94 // deviation drawn from a Poisson distribution of mean xm, using flat()
95 // as a source of uniform random numbers.
96 // (Adapted from Numerical Recipes in C)
97 
98  double em, t, y;
99  double sq, alxm, g1;
100  double om = getOldMean();
101  HepRandomEngine* anEngine = HepRandom::getTheEngine();
102 
103  double* status = getPStatus();
104  sq = status[0];
105  alxm = status[1];
106  g1 = status[2];
107 
108  if( xm == -1 ) return 0;
109  if( xm < 12.0 ) {
110  if( xm != om ) {
111  setOldMean(xm);
112  g1 = std::exp(-xm);
113  }
114  em = -1;
115  t = 1.0;
116  do {
117  em += 1.0;
118  t *= anEngine->flat();
119  } while( t > g1 );
120  }
121  else if ( xm < getMaxMean() ) {
122  if ( xm != om ) {
123  setOldMean(xm);
124  sq = std::sqrt(2.0*xm);
125  alxm = std::log(xm);
126  g1 = xm*alxm - gammln(xm + 1.0);
127  }
128  do {
129  do {
130  y = std::tan(CLHEP::pi*anEngine->flat());
131  em = sq*y + xm;
132  } while( em < 0.0 );
133  em = std::floor(em);
134  t = 0.9*(1.0 + y*y)* std::exp(em*alxm - gammln(em + 1.0) - g1);
135  } while( anEngine->flat() > t );
136  }
137  else {
138  em = xm + std::sqrt(xm) * normal (anEngine); // mf 1/13/06
139  if ( static_cast<long>(em) < 0 )
140  em = static_cast<long>(xm) >= 0 ? xm : getMaxMean();
141  }
142  setPStatus(sq,alxm,g1);
143  return long(em);
144 }
145 
146 void RandPoisson::shootArray(const int size, long* vect, double m1)
147 {
148  for( long* v = vect; v != vect + size; ++v )
149  *v = shoot(m1);
150 }
151 
152 long RandPoisson::shoot(HepRandomEngine* anEngine, double xm) {
153 
154 // Returns as a floating-point number an integer value that is a random
155 // deviation drawn from a Poisson distribution of mean xm, using flat()
156 // of a given Random Engine as a source of uniform random numbers.
157 // (Adapted from Numerical Recipes in C)
158 
159  double em, t, y;
160  double sq, alxm, g1;
161  double om = getOldMean();
162 
163  double* status = getPStatus();
164  sq = status[0];
165  alxm = status[1];
166  g1 = status[2];
167 
168  if( xm == -1 ) return 0;
169  if( xm < 12.0 ) {
170  if( xm != om ) {
171  setOldMean(xm);
172  g1 = std::exp(-xm);
173  }
174  em = -1;
175  t = 1.0;
176  do {
177  em += 1.0;
178  t *= anEngine->flat();
179  } while( t > g1 );
180  }
181  else if ( xm < getMaxMean() ) {
182  if ( xm != om ) {
183  setOldMean(xm);
184  sq = std::sqrt(2.0*xm);
185  alxm = std::log(xm);
186  g1 = xm*alxm - gammln(xm + 1.0);
187  }
188  do {
189  do {
190  y = std::tan(CLHEP::pi*anEngine->flat());
191  em = sq*y + xm;
192  } while( em < 0.0 );
193  em = std::floor(em);
194  t = 0.9*(1.0 + y*y)* std::exp(em*alxm - gammln(em + 1.0) - g1);
195  } while( anEngine->flat() > t );
196  }
197  else {
198  em = xm + std::sqrt(xm) * normal (anEngine); // mf 1/13/06
199  if ( static_cast<long>(em) < 0 )
200  em = static_cast<long>(xm) >= 0 ? xm : getMaxMean();
201  }
202  setPStatus(sq,alxm,g1);
203  return long(em);
204 }
205 
206 void RandPoisson::shootArray(HepRandomEngine* anEngine, const int size,
207  long* vect, double m1)
208 {
209  for( long* v = vect; v != vect + size; ++v )
210  *v = shoot(anEngine,m1);
211 }
212 
213 long RandPoisson::fire() {
214  return long(fire( defaultMean ));
215 }
216 
217 long RandPoisson::fire(double xm) {
218 
219 // Returns as a floating-point number an integer value that is a random
220 // deviation drawn from a Poisson distribution of mean xm, using flat()
221 // as a source of uniform random numbers.
222 // (Adapted from Numerical Recipes in C)
223 
224  double em, t, y;
225  double sq, alxm, g1;
226 
227  sq = status[0];
228  alxm = status[1];
229  g1 = status[2];
230 
231  if( xm == -1 ) return 0;
232  if( xm < 12.0 ) {
233  if( xm != oldm ) {
234  oldm = xm;
235  g1 = std::exp(-xm);
236  }
237  em = -1;
238  t = 1.0;
239  do {
240  em += 1.0;
241  t *= localEngine->flat();
242  } while( t > g1 );
243  }
244  else if ( xm < meanMax ) {
245  if ( xm != oldm ) {
246  oldm = xm;
247  sq = std::sqrt(2.0*xm);
248  alxm = std::log(xm);
249  g1 = xm*alxm - gammln(xm + 1.0);
250  }
251  do {
252  do {
253  y = std::tan(CLHEP::pi*localEngine->flat());
254  em = sq*y + xm;
255  } while( em < 0.0 );
256  em = std::floor(em);
257  t = 0.9*(1.0 + y*y)* std::exp(em*alxm - gammln(em + 1.0) - g1);
258  } while( localEngine->flat() > t );
259  }
260  else {
261  em = xm + std::sqrt(xm) * normal (localEngine.get()); // mf 1/13/06
262  if ( static_cast<long>(em) < 0 )
263  em = static_cast<long>(xm) >= 0 ? xm : getMaxMean();
264  }
265  status[0] = sq; status[1] = alxm; status[2] = g1;
266  return long(em);
267 }
268 
269 void RandPoisson::fireArray(const int size, long* vect )
270 {
271  for( long* v = vect; v != vect + size; ++v )
272  *v = fire( defaultMean );
273 }
274 
275 void RandPoisson::fireArray(const int size, long* vect, double m1)
276 {
277  for( long* v = vect; v != vect + size; ++v )
278  *v = fire( m1 );
279 }
280 
281 std::ostream & RandPoisson::put ( std::ostream & os ) const {
282  int pr=os.precision(20);
283  std::vector<unsigned long> t(2);
284  os << " " << name() << "\n";
285  os << "Uvec" << "\n";
286  t = DoubConv::dto2longs(meanMax);
287  os << meanMax << " " << t[0] << " " << t[1] << "\n";
288  t = DoubConv::dto2longs(defaultMean);
289  os << defaultMean << " " << t[0] << " " << t[1] << "\n";
290  t = DoubConv::dto2longs(status[0]);
291  os << status[0] << " " << t[0] << " " << t[1] << "\n";
292  t = DoubConv::dto2longs(status[1]);
293  os << status[1] << " " << t[0] << " " << t[1] << "\n";
294  t = DoubConv::dto2longs(status[2]);
295  os << status[2] << " " << t[0] << " " << t[1] << "\n";
296  t = DoubConv::dto2longs(oldm);
297  os << oldm << " " << t[0] << " " << t[1] << "\n";
298  os.precision(pr);
299  return os;
300 }
301 
302 std::istream & RandPoisson::get ( std::istream & is ) {
303  std::string inName;
304  is >> inName;
305  if (inName != name()) {
306  is.clear(std::ios::badbit | is.rdstate());
307  std::cerr << "Mismatch when expecting to read state of a "
308  << name() << " distribution\n"
309  << "Name found was " << inName
310  << "\nistream is left in the badbit state\n";
311  return is;
312  }
313  if (possibleKeywordInput(is, "Uvec", meanMax)) {
314  std::vector<unsigned long> t(2);
315  is >> meanMax >> t[0] >> t[1]; meanMax = DoubConv::longs2double(t);
316  is >> defaultMean >> t[0] >> t[1]; defaultMean = DoubConv::longs2double(t);
317  is >> status[0] >> t[0] >> t[1]; status[0] = DoubConv::longs2double(t);
318  is >> status[1] >> t[0] >> t[1]; status[1] = DoubConv::longs2double(t);
319  is >> status[2] >> t[0] >> t[1]; status[2] = DoubConv::longs2double(t);
320  is >> oldm >> t[0] >> t[1]; oldm = DoubConv::longs2double(t);
321  return is;
322  }
323  // is >> meanMax encompassed by possibleKeywordInput
324  is >> defaultMean >> status[0] >> status[1] >> status[2];
325  return is;
326 }
327 
328 } // namespace CLHEP
329 
ThreeVector shoot(const G4int Ap, const G4int Af)
double gammln(double xx)
Definition: RandPoisson.cc:54
G4String name
Definition: TRTMaterials.hh:40
static int engine(pchar, pchar, double &, pchar &, const dic_type &)
Definition: Evaluator.cc:358
static double normal(HepRandomEngine *eptr)
Definition: RandPoisson.cc:77
static const double pi
Definition: G4SIunits.hh:74
const G4double x[NPOINTSGL]
static const G4double fac