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
34namespace CLHEP {
35
36std::string RandPoisson::name() const {return "RandPoisson";}
37HepRandomEngine & RandPoisson::engine() {return *localEngine;}
38
39// Initialisation of static data
40double RandPoisson::status_st[3] = {0., 0., 0.};
41double RandPoisson::oldm_st = -1.0;
42const double RandPoisson::meanMax_st = 2.0E9;
43
45}
46
48 return double(fire( defaultMean ));
49}
50
51double RandPoisson::operator()( double mean ) {
52 return double(fire( mean ));
53}
54
55double 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
77static
78double 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
92long 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
147void 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
153long 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
207void 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
218long 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
270void RandPoisson::fireArray(const int size, long* vect )
271{
272 for( long* v = vect; v != vect + size; ++v )
273 *v = fire( defaultMean );
274}
275
276void 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
282std::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
311std::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
static double longs2double(const std::vector< unsigned long > &v)
Definition: DoubConv.cc:106
static std::vector< unsigned long > dto2longs(double d)
Definition: DoubConv.cc:90
virtual double flat()=0
static HepRandomEngine * getTheEngine()
Definition: Random.cc:166
std::string name() const
Definition: RandPoisson.cc:36
std::ostream & put(std::ostream &os) const
Definition: RandPoisson.cc:282
virtual ~RandPoisson()
Definition: RandPoisson.cc:44
static long shoot(double m=1.0)
Definition: RandPoisson.cc:92
static void setOldMean(double val)
static void setPStatus(double sq, double alxm, double g1)
static void shootArray(const int size, long *vect, double m=1.0)
Definition: RandPoisson.cc:147
void fireArray(const int size, long *vect)
Definition: RandPoisson.cc:270
std::istream & get(std::istream &is)
Definition: RandPoisson.cc:311
HepRandomEngine & engine()
Definition: RandPoisson.cc:37
#define double(obj)
Definition: excDblThrow.cc:32
bool possibleKeywordInput(IS &is, const std::string &key, T &t)
double gammln(double xx)
Definition: RandPoisson.cc:55