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 
double
#define double(obj)
Definition: excDblThrow.cc:32
CLHEP::RandPoisson::getPStatus
static double * getPStatus()
Definition: Matrix/CLHEP/Random/RandPoisson.h:108
CLHEP::RandPoisson::get
std::istream & get(std::istream &is)
Definition: RandPoisson.cc:311
CLHEP::HepRandomEngine
Definition: Matrix/CLHEP/Random/RandomEngine.h:55
CLHEP::RandPoisson::put
std::ostream & put(std::ostream &os) const
Definition: RandPoisson.cc:282
CLHEP::RandPoisson::engine
HepRandomEngine & engine()
Definition: RandPoisson.cc:37
CLHEP::RandPoisson::fire
long fire()
Definition: RandPoisson.cc:214
CLHEP::RandPoisson::setOldMean
static void setOldMean(double val)
Definition: Matrix/CLHEP/Random/RandPoisson.h:106
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::DoubConv::longs2double
static double longs2double(const std::vector< unsigned long > &v)
Definition: DoubConv.cc:106
CLHEP::gammln
double gammln(double xx)
Definition: RandPoisson.cc:55
CLHEP::RandPoisson::shootArray
static void shootArray(const int size, long *vect, double m=1.0)
Definition: RandPoisson.cc:147
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::RandPoisson::getMaxMean
static double getMaxMean()
Definition: Matrix/CLHEP/Random/RandPoisson.h:104
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::RandPoisson::defaultMean
double defaultMean
Definition: Matrix/CLHEP/Random/RandPoisson.h:100
CLHEP::RandPoisson::meanMax
double meanMax
Definition: Matrix/CLHEP/Random/RandPoisson.h:99
CLHEP::RandPoisson::getOldMean
static double getOldMean()
Definition: Matrix/CLHEP/Random/RandPoisson.h:102
CLHEP::possibleKeywordInput
bool possibleKeywordInput(IS &is, const std::string &key, T &t)
Definition: Matrix/CLHEP/Random/RandomEngine.h:168
CLHEP::RandPoisson::operator()
double operator()()
Definition: RandPoisson.cc:47
j
long j
Definition: JamesRandomSeeding.txt:28
CLHEP::RandPoisson::~RandPoisson
virtual ~RandPoisson()
Definition: RandPoisson.cc:44
CLHEP::HepRandomEngine::flat
virtual double flat()=0
CLHEP::DoubConv::dto2longs
static std::vector< unsigned long > dto2longs(double d)
Definition: DoubConv.cc:90
CLHEP::RandPoisson::fireArray
void fireArray(const int size, long *vect)
Definition: RandPoisson.cc:270
CLHEP::HepRandom::getTheEngine
static HepRandomEngine * getTheEngine()
Definition: Random.cc:166
CLHEP::RandPoisson::setPStatus
static void setPStatus(double sq, double alxm, double g1)
Definition: Matrix/CLHEP/Random/RandPoisson.h:110
x
any side effects of that construction would occur twice The semantics of throw x
Definition: whyZMthrowRethrows.txt:37
CLHEP::RandPoisson::shoot
static long shoot(double m=1.0)
Definition: RandPoisson.cc:92
CLHEP::RandPoisson::name
std::string name() const
Definition: RandPoisson.cc:36