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