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

RandPoisson.cc
Go to the documentation of this file.
1 // $Id: RandPoisson.cc,v 1.7 2010/06/16 17:24:53 garren Exp $
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/defs.h"
29 #include "CLHEP/Random/RandPoisson.h"
30 #include "CLHEP/Units/PhysicalConstants.h"
31 #include "CLHEP/Random/DoubConv.hh"
32 #include <cmath> // for std::floor()
33 
34 namespace CLHEP {
35 
36 std::string RandPoisson::name() const {return "RandPoisson";}
37 HepRandomEngine & RandPoisson::engine() {return *localEngine;}
38 
39 // Initialisation of static data
40 double RandPoisson::status_st[3] = {0., 0., 0.};
41 double RandPoisson::oldm_st = -1.0;
42 const double RandPoisson::meanMax_st = 2.0E9;
43 
45 }
46 
48  return double(fire( defaultMean ));
49 }
50 
51 double RandPoisson::operator()( double mean ) {
52  return double(fire( mean ));
53 }
54 
55 double gammln(double xx) {
56 
57 // Returns the value ln(Gamma(xx) for xx > 0. Full accuracy is obtained for
58 // xx > 1. For 0 < xx < 1. the reflection formula (6.1.4) can be used first.
59 // (Adapted from Numerical Recipes in C)
60 
61  static double cof[6] = {76.18009172947146,-86.50532032941677,
62  24.01409824083091, -1.231739572450155,
63  0.1208650973866179e-2, -0.5395239384953e-5};
64  int j;
65  double x = xx - 1.0;
66  double tmp = x + 5.5;
67  tmp -= (x + 0.5) * std::log(tmp);
68  double ser = 1.000000000190015;
69 
70  for ( j = 0; j <= 5; j++ ) {
71  x += 1.0;
72  ser += cof[j]/x;
73  }
74  return -tmp + std::log(2.5066282746310005*ser);
75 }
76 
77 static
78 double normal (HepRandomEngine* eptr) // mf 1/13/06
79 {
80  double r;
81  double v1,v2,fac;
82  do {
83  v1 = 2.0 * eptr->flat() - 1.0;
84  v2 = 2.0 * eptr->flat() - 1.0;
85  r = v1*v1 + v2*v2;
86  } while ( r > 1.0 );
87 
88  fac = std::sqrt(-2.0*std::log(r)/r);
89  return v2*fac;
90 }
91 
92 long RandPoisson::shoot(double xm) {
93 
94 // Returns as a floating-point number an integer value that is a random
95 // deviation drawn from a Poisson distribution of mean xm, using flat()
96 // as a source of uniform random numbers.
97 // (Adapted from Numerical Recipes in C)
98 
99  double em, t, y;
100  double sq, alxm, g1;
101  double om = getOldMean();
103 
104  double* status = getPStatus();
105  sq = status[0];
106  alxm = status[1];
107  g1 = status[2];
108 
109  if( xm == -1 ) return 0;
110  if( xm < 12.0 ) {
111  if( xm != om ) {
112  setOldMean(xm);
113  g1 = std::exp(-xm);
114  }
115  em = -1;
116  t = 1.0;
117  do {
118  em += 1.0;
119  t *= anEngine->flat();
120  } while( t > g1 );
121  }
122  else if ( xm < getMaxMean() ) {
123  if ( xm != om ) {
124  setOldMean(xm);
125  sq = std::sqrt(2.0*xm);
126  alxm = std::log(xm);
127  g1 = xm*alxm - gammln(xm + 1.0);
128  }
129  do {
130  do {
131  y = std::tan(CLHEP::pi*anEngine->flat());
132  em = sq*y + xm;
133  } while( em < 0.0 );
134  em = std::floor(em);
135  t = 0.9*(1.0 + y*y)* std::exp(em*alxm - gammln(em + 1.0) - g1);
136  } while( anEngine->flat() > t );
137  }
138  else {
139  em = xm + std::sqrt(xm) * normal (anEngine); // mf 1/13/06
140  if ( static_cast<long>(em) < 0 )
141  em = static_cast<long>(xm) >= 0 ? xm : getMaxMean();
142  }
143  setPStatus(sq,alxm,g1);
144  return long(em);
145 }
146 
147 void RandPoisson::shootArray(const int size, long* vect, double m1)
148 {
149  for( long* v = vect; v != vect + size; ++v )
150  *v = shoot(m1);
151 }
152 
153 long RandPoisson::shoot(HepRandomEngine* anEngine, double xm) {
154 
155 // Returns as a floating-point number an integer value that is a random
156 // deviation drawn from a Poisson distribution of mean xm, using flat()
157 // of a given Random Engine as a source of uniform random numbers.
158 // (Adapted from Numerical Recipes in C)
159 
160  double em, t, y;
161  double sq, alxm, g1;
162  double om = getOldMean();
163 
164  double* status = getPStatus();
165  sq = status[0];
166  alxm = status[1];
167  g1 = status[2];
168 
169  if( xm == -1 ) return 0;
170  if( xm < 12.0 ) {
171  if( xm != om ) {
172  setOldMean(xm);
173  g1 = std::exp(-xm);
174  }
175  em = -1;
176  t = 1.0;
177  do {
178  em += 1.0;
179  t *= anEngine->flat();
180  } while( t > g1 );
181  }
182  else if ( xm < getMaxMean() ) {
183  if ( xm != om ) {
184  setOldMean(xm);
185  sq = std::sqrt(2.0*xm);
186  alxm = std::log(xm);
187  g1 = xm*alxm - gammln(xm + 1.0);
188  }
189  do {
190  do {
191  y = std::tan(CLHEP::pi*anEngine->flat());
192  em = sq*y + xm;
193  } while( em < 0.0 );
194  em = std::floor(em);
195  t = 0.9*(1.0 + y*y)* std::exp(em*alxm - gammln(em + 1.0) - g1);
196  } while( anEngine->flat() > t );
197  }
198  else {
199  em = xm + std::sqrt(xm) * normal (anEngine); // mf 1/13/06
200  if ( static_cast<long>(em) < 0 )
201  em = static_cast<long>(xm) >= 0 ? xm : getMaxMean();
202  }
203  setPStatus(sq,alxm,g1);
204  return long(em);
205 }
206 
207 void RandPoisson::shootArray(HepRandomEngine* anEngine, const int size,
208  long* vect, double m1)
209 {
210  for( long* v = vect; v != vect + size; ++v )
211  *v = shoot(anEngine,m1);
212 }
213 
215  return long(fire( defaultMean ));
216 }
217 
218 long RandPoisson::fire(double xm) {
219 
220 // Returns as a floating-point number an integer value that is a random
221 // deviation drawn from a Poisson distribution of mean xm, using flat()
222 // as a source of uniform random numbers.
223 // (Adapted from Numerical Recipes in C)
224 
225  double em, t, y;
226  double sq, alxm, g1;
227 
228  sq = status[0];
229  alxm = status[1];
230  g1 = status[2];
231 
232  if( xm == -1 ) return 0;
233  if( xm < 12.0 ) {
234  if( xm != oldm ) {
235  oldm = xm;
236  g1 = std::exp(-xm);
237  }
238  em = -1;
239  t = 1.0;
240  do {
241  em += 1.0;
242  t *= localEngine->flat();
243  } while( t > g1 );
244  }
245  else if ( xm < meanMax ) {
246  if ( xm != oldm ) {
247  oldm = xm;
248  sq = std::sqrt(2.0*xm);
249  alxm = std::log(xm);
250  g1 = xm*alxm - gammln(xm + 1.0);
251  }
252  do {
253  do {
254  y = std::tan(CLHEP::pi*localEngine->flat());
255  em = sq*y + xm;
256  } while( em < 0.0 );
257  em = std::floor(em);
258  t = 0.9*(1.0 + y*y)* std::exp(em*alxm - gammln(em + 1.0) - g1);
259  } while( localEngine->flat() > t );
260  }
261  else {
262  em = xm + std::sqrt(xm) * normal (localEngine.get()); // mf 1/13/06
263  if ( static_cast<long>(em) < 0 )
264  em = static_cast<long>(xm) >= 0 ? xm : getMaxMean();
265  }
266  status[0] = sq; status[1] = alxm; status[2] = g1;
267  return long(em);
268 }
269 
270 void RandPoisson::fireArray(const int size, long* vect )
271 {
272  for( long* v = vect; v != vect + size; ++v )
273  *v = fire( defaultMean );
274 }
275 
276 void RandPoisson::fireArray(const int size, long* vect, double m1)
277 {
278  for( long* v = vect; v != vect + size; ++v )
279  *v = fire( m1 );
280 }
281 
282 std::ostream & RandPoisson::put ( std::ostream & os ) const {
283  int pr=os.precision(20);
284  std::vector<unsigned long> t(2);
285  os << " " << name() << "\n";
286  os << "Uvec" << "\n";
288  os << meanMax << " " << t[0] << " " << t[1] << "\n";
290  os << defaultMean << " " << t[0] << " " << t[1] << "\n";
291  t = DoubConv::dto2longs(status[0]);
292  os << status[0] << " " << t[0] << " " << t[1] << "\n";
293  t = DoubConv::dto2longs(status[1]);
294  os << status[1] << " " << t[0] << " " << t[1] << "\n";
295  t = DoubConv::dto2longs(status[2]);
296  os << status[2] << " " << t[0] << " " << t[1] << "\n";
297  t = DoubConv::dto2longs(oldm);
298  os << oldm << " " << t[0] << " " << t[1] << "\n";
299  os.precision(pr);
300  return os;
301 #ifdef REMOVED
302  int pr=os.precision(20);
303  os << " " << name() << "\n";
304  os << meanMax << " " << defaultMean << "\n";
305  os << status[0] << " " << status[1] << " " << status[2] << "\n";
306  os.precision(pr);
307  return os;
308 #endif
309 }
310 
311 std::istream & RandPoisson::get ( std::istream & is ) {
312  std::string inName;
313  is >> inName;
314  if (inName != name()) {
315  is.clear(std::ios::badbit | is.rdstate());
316  std::cerr << "Mismatch when expecting to read state of a "
317  << name() << " distribution\n"
318  << "Name found was " << inName
319  << "\nistream is left in the badbit state\n";
320  return is;
321  }
322  if (possibleKeywordInput(is, "Uvec", meanMax)) {
323  std::vector<unsigned long> t(2);
324  is >> meanMax >> t[0] >> t[1]; meanMax = DoubConv::longs2double(t);
325  is >> defaultMean >> t[0] >> t[1]; defaultMean = DoubConv::longs2double(t);
326  is >> status[0] >> t[0] >> t[1]; status[0] = DoubConv::longs2double(t);
327  is >> status[1] >> t[0] >> t[1]; status[1] = DoubConv::longs2double(t);
328  is >> status[2] >> t[0] >> t[1]; status[2] = DoubConv::longs2double(t);
329  is >> oldm >> t[0] >> t[1]; oldm = DoubConv::longs2double(t);
330  return is;
331  }
332  // is >> meanMax encompassed by possibleKeywordInput
333  is >> defaultMean >> status[0] >> status[1] >> status[2];
334  return is;
335 }
336 
337 } // namespace CLHEP
338 
bool possibleKeywordInput(IS &is, const std::string &key, T &t)
#define double(obj)
Definition: excDblThrow.cc:32
static HepRandomEngine * getTheEngine()
Definition: Random.cc:166
static void shootArray(const int size, long *vect, double m=1.0)
Definition: RandPoisson.cc:147
virtual ~RandPoisson()
Definition: RandPoisson.cc:44
void fireArray(const int size, long *vect)
Definition: RandPoisson.cc:270
static long shoot(double m=1.0)
Definition: RandPoisson.cc:92
static double longs2double(const std::vector< unsigned long > &v)
Definition: DoubConv.cc:106
double gammln(double xx)
Definition: RandPoisson.cc:55
static void setPStatus(double sq, double alxm, double g1)
virtual double flat()=0
static void setOldMean(double val)
HepRandomEngine & engine()
Definition: RandPoisson.cc:37
std::istream & get(std::istream &is)
Definition: RandPoisson.cc:311
static std::vector< unsigned long > dto2longs(double d)
Definition: DoubConv.cc:90
std::ostream & put(std::ostream &os) const
Definition: RandPoisson.cc:282
std::string name() const
Definition: RandPoisson.cc:36