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
name() const36 std::string RandPoisson::name() const {return "RandPoisson";}
engine()37 HepRandomEngine & RandPoisson::engine() {return *localEngine;}
38
39 // Initialisation of static data
40 CLHEP_THREAD_LOCAL double RandPoisson::status_st[3] = {0., 0., 0.};
41 CLHEP_THREAD_LOCAL double RandPoisson::oldm_st = -1.0;
42 const double RandPoisson::meanMax_st = 2.0E9;
43
~RandPoisson()44 RandPoisson::~RandPoisson() {
45 }
46
operator ()()47 double RandPoisson::operator()() {
48 return double(fire( defaultMean ));
49 }
50
operator ()(double mean)51 double RandPoisson::operator()( double mean ) {
52 return double(fire( mean ));
53 }
54
gammln(double xx)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 const 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
normal(HepRandomEngine * eptr)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
shoot(double xm)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();
102 HepRandomEngine* anEngine = HepRandom::getTheEngine();
103
104 double* pstatus = getPStatus();
105 sq = pstatus[0];
106 alxm = pstatus[1];
107 g1 = pstatus[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
shootArray(const int size,long * vect,double m1)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
shoot(HepRandomEngine * anEngine,double xm)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* pstatus = getPStatus();
165 sq = pstatus[0];
166 alxm = pstatus[1];
167 g1 = pstatus[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
shootArray(HepRandomEngine * anEngine,const int size,long * vect,double m1)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
fire()214 long RandPoisson::fire() {
215 return long(fire( defaultMean ));
216 }
217
fire(double xm)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
fireArray(const int size,long * vect)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
fireArray(const int size,long * vect,double m1)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
put(std::ostream & os) const282 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";
287 t = DoubConv::dto2longs(meanMax);
288 os << meanMax << " " << t[0] << " " << t[1] << "\n";
289 t = DoubConv::dto2longs(defaultMean);
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
get(std::istream & is)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
339