Geant4 11.2.2
Toolkit for the simulation of the passage of particles through matter
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
RandPoissonQ.cc
Go to the documentation of this file.
1// -*- C++ -*-
2//
3// -----------------------------------------------------------------------
4// HEP Random
5// --- RandPoissonQ ---
6// class implementation file
7// -----------------------------------------------------------------------
8
9// =======================================================================
10// M. Fischler - Implemented new, much faster table-driven algorithm
11// applicable for mu < 100
12// - Implemented "quick()" methods, shich are the same as the
13// new methods for mu < 100 and are a skew-corrected gaussian
14// approximation for large mu.
15// M. Fischler - Removed mean=100 from the table-driven set, since it
16// uses a value just off the end of the table. (April 2004)
17// M Fischler - put and get to/from streams 12/15/04
18// M Fischler - Utilize RandGaussQ rather than RandGauss, as clearly
19// intended by the inclusion of RandGaussQ.h. Using RandGauss
20// introduces a subtle trap in that the state of RandPoissonQ
21// can never be properly captured without also saveing the
22// state of RandGauss! RandGaussQ is, on the other hand,
23// stateless except for the engine used.
24// M Fischler - Modified use of wrong engine when shoot (anEngine, mean)
25// is called. This flaw was preventing any hope of proper
26// saving and restoring in the instance cases.
27// M Fischler - fireArray using defaultMean 2/10/05
28// M Fischler - put/get to/from streams uses pairs of ulongs when
29// + storing doubles avoid problems with precision
30// 4/14/05
31// M Fischler - Modified use of shoot (mean) instead of
32// shoot(getLocalEngine(), mean) when fire(mean) is called.
33// This flaw was causing bad "cross-talk" between modules
34// in CMS, where one used its own engine, and the other
35// used the static generator. 10/18/07
36//
37// =======================================================================
38
42#include "CLHEP/Random/Stat.h"
44#include <cmath> // for std::pow()
45#include <iostream>
46#include <string>
47#include <vector>
48
49namespace CLHEP {
50
51std::string RandPoissonQ::name() const {return "RandPoissonQ";}
53
54// Initialization of static data: Note that this is all const static data,
55// so that saveEngineStatus properly saves all needed information.
56
57 // The following MUST MATCH the corresponding values used (in
58 // poissonTables.cc) when poissonTables.cdat was created.
59
60const double RandPoissonQ::FIRST_MU = 10;// lowest mu value in table
61const double RandPoissonQ::LAST_MU = 95;// highest mu value
62const double RandPoissonQ::S = 5; // Spacing between mu values
63const int RandPoissonQ::BELOW = 30; // Starting point for N is at mu - BELOW
64const int RandPoissonQ::ENTRIES = 51; // Number of entries in each mu row
65
66const double RandPoissonQ::MAXIMUM_POISSON_DEVIATE = 2.0E9;
67 // Careful -- this is NOT the maximum number that can be held in
68 // a long. It actually should be some large number of sigma below
69 // that.
70
71 // Here comes the big (9K bytes) table, kept in a file of
72 // ENTRIES * (FIRST_MU - LAST_MU + 1)/S doubles
73
74static const double poissonTables [ 51 * ( (95-10)/5 + 1 ) ] = {
75#include "CLHEP/Random/poissonTables.cdat"
76};
77
78//
79// Constructors and destructors:
80//
81
84
85void RandPoissonQ::setupForDefaultMu() {
86
87 // The following are useful for quick approximation, for large mu
88
89 double sig2 = defaultMean * (.9998654 - .08346/defaultMean);
90 sigma = std::sqrt(sig2);
91 // sigma for the Guassian which approximates the Poisson -- naively
92 // sqrt (defaultMean).
93 //
94 // The multiplier corrects for fact that discretization of the form
95 // [gaussian+.5] increases the second moment by a small amount.
96
97 double t = 1./(sig2);
98
99 a2 = t/6 + t*t/324;
100 a1 = std::sqrt (1-2*a2*a2*sig2);
101 a0 = defaultMean + .5 - sig2 * a2;
102
103 // The formula will be a0 + a1*x + a2*x*x where x has 2nd moment of sigma.
104 // The coeffeicients are chosen to match the first THREE moments of the
105 // true Poisson distribution.
106 //
107 // Actually, if the correction for discretization were not needed, then
108 // a2 could be taken one order higher by adding t*t*t/5832. However,
109 // the discretization correction is not perfect, leading to inaccuracy
110 // on the order to 1/mu**2, so adding a third term is overkill.
111
112} // setupForDefaultMu()
113
114
115//
116// fire, quick, operator(), and shoot methods:
117//
118
119long RandPoissonQ::shoot(double xm) {
120 return shoot(getTheEngine(), xm);
121}
122
124 return (double) fire();
125}
126
127double RandPoissonQ::operator()( double mean ) {
128 return (double) fire(mean);
129}
130
131long RandPoissonQ::fire(double mean) {
132 return shoot(getLocalEngine(), mean);
133}
134
136 if ( defaultMean < LAST_MU + S ) {
137 return poissonDeviateSmall ( getLocalEngine(), defaultMean );
138 } else {
139 return poissonDeviateQuick ( getLocalEngine(), a0, a1, a2, sigma );
140 }
141} // fire()
142
143long RandPoissonQ::shoot(HepRandomEngine* anEngine, double mean) {
144
145 // The following variables, static to this method, apply to the
146 // last time a large mean was supplied; they obviate certain calculations
147 // if consecutive calls use the same mean.
148
149 static CLHEP_THREAD_LOCAL double lastLargeMean = -1.; // Mean from previous shoot
150 // requiring poissonDeviateQuick()
151 static CLHEP_THREAD_LOCAL double lastA0;
152 static CLHEP_THREAD_LOCAL double lastA1;
153 static CLHEP_THREAD_LOCAL double lastA2;
154 static CLHEP_THREAD_LOCAL double lastSigma;
155
156 if ( mean < LAST_MU + S ) {
157 return poissonDeviateSmall ( anEngine, mean );
158 } else {
159 if ( mean != lastLargeMean ) {
160 // Compute the coefficients defining the quadratic transformation from a
161 // Gaussian to a Poisson for this mean. Also save these for next time.
162 double sig2 = mean * (.9998654 - .08346/mean);
163 lastSigma = std::sqrt(sig2);
164 double t = 1./sig2;
165 lastA2 = t*(1./6.) + t*t*(1./324.);
166 lastA1 = std::sqrt (1-2*lastA2*lastA2*sig2);
167 lastA0 = mean + .5 - sig2 * lastA2;
168 }
169 return poissonDeviateQuick ( anEngine, lastA0, lastA1, lastA2, lastSigma );
170 }
171
172} // shoot (anEngine, mean)
173
174void RandPoissonQ::shootArray(const int size, long* vect, double m) {
175 for( long* v = vect; v != vect + size; ++v )
176 *v = shoot(m);
177 // Note: We could test for m > 100, and if it is, precompute a0, a1, a2,
178 // and sigma and call the appropriate form of poissonDeviateQuick.
179 // But since those are cached anyway, not much time would be saved.
180}
181
182void RandPoissonQ::shootArray(HepRandomEngine* anEngine, const int size,
183 long* vect, double m1) {
184 for( long* v = vect; v != vect + size; ++v )
185 *v = shoot(anEngine,m1);
186}
187
188void RandPoissonQ::fireArray(const int size, long* vect, double m) {
189 for( long* v = vect; v != vect + size; ++v )
190 *v = fire( m );
191}
192
193void RandPoissonQ::fireArray(const int size, long* vect) {
194 for( long* v = vect; v != vect + size; ++v )
195 *v = fire( defaultMean );
196}
197
198
199// Quick Poisson deviate algorithm used by quick for large mu:
200
201long RandPoissonQ::poissonDeviateQuick ( HepRandomEngine *e, double mu ) {
202
203 // Compute the coefficients defining the quadratic transformation from a
204 // Gaussian to a Poisson:
205
206 double sig2 = mu * (.9998654 - .08346/mu);
207 double sig = std::sqrt(sig2);
208 // The multiplier corrects for fact that discretization of the form
209 // [gaussian+.5] increases the second moment by a small amount.
210
211 double t = 1./sig2;
212
213 double sa2 = t*(1./6.) + t*t*(1./324.);
214 double sa1 = std::sqrt (1-2*sa2*sa2*sig2);
215 double sa0 = mu + .5 - sig2 * sa2;
216
217 // The formula will be sa0 + sa1*x + sa2*x*x where x has sigma of sq.
218 // The coeffeicients are chosen to match the first THREE moments of the
219 // true Poisson distribution.
220
221 return poissonDeviateQuick ( e, sa0, sa1, sa2, sig );
222}
223
224
225long RandPoissonQ::poissonDeviateQuick ( HepRandomEngine *e,
226 double A0, double A1, double A2, double sig) {
227//
228// Quick Poisson deviate algorithm used by quick for large mu:
229//
230// The principle: For very large mu, a poisson distribution can be approximated
231// by a gaussian: return the integer part of mu + .5 + g where g is a unit
232// normal. However, this yelds a miserable approximation at values as
233// "large" as 100. The primary problem is that the poisson distribution is
234// supposed to have a skew of 1/mu**2, and the zero skew of the Guassian
235// leads to errors of order as big as 1/mu**2.
236//
237// We substitute for the gaussian a quadratic function of that gaussian random.
238// The expression looks very nearly like mu + .5 - 1/6 + g + g**2/(6*mu).
239// The small positive quadratic term causes the resulting variate to have
240// a positive skew; the -1/6 constant term is there to correct for this bias
241// in the mean. By adjusting these two and the linear term, we can match the
242// first three moments to high accuracy in 1/mu.
243//
244// The sigma used is not precisely sqrt(mu) since a rounded-off Gaussian
245// has a second moment which is slightly larger than that of the Gaussian.
246// To compensate, sig is multiplied by a factor which is slightly less than 1.
247
248 // double g = RandGauss::shootQuick( e ); // TEMPORARY MOD:
249 double g = RandGaussQ::shoot( e ); // Unit normal
250 g *= sig;
251 double p = A2*g*g + A1*g + A0;
252 if ( p < 0 ) return 0; // Shouldn't ever possibly happen since
253 // mean should not be less than 100, but
254 // we check due to paranoia.
256
257 return long(p);
258
259} // poissonDeviateQuick ()
260
261
262
263long RandPoissonQ::poissonDeviateSmall (HepRandomEngine * e, double mean) {
264 long N1;
265 long N2;
266 // The following are for later use to form a secondary random s:
267 double rRange; // This will hold the interval between cdf for the
268 // computed N1 and cdf for N1+1.
269 double rRemainder = 0; // This will hold the length into that interval.
270
271 // Coming in, mean should not be more than LAST_MU + S. However, we will
272 // be paranoid and test for this:
273
274 if ( mean > LAST_MU + S ) {
275 return RandPoisson::shoot(e, mean);
276 }
277
278 if (mean <= 0) {
279 return 0; // Perhaps we ought to balk harder here!
280 }
281
282 // >>> 1 <<<
283 // Generate the first random, which we always will need.
284
285 double r = e->flat();
286
287 // >>> 2 <<<
288 // For small mean, below the start of the tables,
289 // do the series for cdf directly.
290
291 // In this case, since we know the series will terminate relatively quickly,
292 // almost alwaye use a precomputed 1/N array without fear of overrunning it.
293
294 static const double oneOverN[50] =
295 { 0, 1., 1/2., 1/3., 1/4., 1/5., 1/6., 1/7., 1/8., 1/9.,
296 1/10., 1/11., 1/12., 1/13., 1/14., 1/15., 1/16., 1/17., 1/18., 1/19.,
297 1/20., 1/21., 1/22., 1/23., 1/24., 1/25., 1/26., 1/27., 1/28., 1/29.,
298 1/30., 1/31., 1/32., 1/33., 1/34., 1/35., 1/36., 1/37., 1/38., 1/39.,
299 1/40., 1/41., 1/42., 1/43., 1/44., 1/45., 1/46., 1/47., 1/48., 1/49. };
300
301
302 if ( mean < FIRST_MU ) {
303
304 long N = 0;
305 double term = std::exp(-mean);
306 double cdf = term;
307
308 if ( r < (1 - 1.0E-9) ) {
309 //
310 // **** This is a normal path: ****
311 //
312 // Except when r is very close to 1, it is certain that we will exceed r
313 // before the 30-th term in the series, so a simple while loop is OK.
314 const double* oneOverNptr = oneOverN;
315 while( cdf <= r ) {
316 N++ ;
317 oneOverNptr++;
318 term *= ( mean * (*oneOverNptr) );
319 cdf += term;
320 }
321 return N;
322 //
323 // **** ****
324 //
325 } else { // r is almost 1...
326 // For r very near to 1 we would have to check that we don't fall
327 // off the end of the table of 1/N. Since this is very rare, we just
328 // ignore the table and do the identical while loop, using explicit
329 // division.
330 double cdf0;
331 while ( cdf <= r ) {
332 N++ ;
333 term *= ( mean / N );
334 cdf0 = cdf;
335 cdf += term;
336 if (cdf == cdf0) break; // Can't happen, but just in case...
337 }
338 return N;
339 } // end of if ( r compared to (1 - 1.0E-9) )
340
341 } // End of the code for mean < FIRST_MU
342
343 // >>> 3 <<<
344 // Find the row of the tables corresponding to the highest tabulated mu
345 // which is no greater than our actual mean.
346
347 int rowNumber = int((mean - FIRST_MU)/S);
348 const double * cdfs = &poissonTables [rowNumber*ENTRIES];
349 double mu = FIRST_MU + rowNumber*S;
350 double deltaMu = mean - mu;
351 int Nmin = int(mu - BELOW);
352 if (Nmin < 1) Nmin = 1;
353 int Nmax = Nmin + (ENTRIES - 1);
354
355
356 // >>> 4 <<<
357 // If r is less that the smallest entry in the row, then
358 // generate the deviate directly from the series.
359
360 if ( r < cdfs[0] ) {
361
362 // In this case, we are tempted to use the actual mean, and not
363 // generate a second deviate to account for the leftover part mean - mu.
364 // That would be an error, generating a distribution with enough excess
365 // at Nmin + (mean-mu)/2 to be detectable in 4,000,000 trials.
366
367 // Since this case is very rare (never more than .2% of the r values)
368 // and can happen where N will be large (up to 65 for the mu=95 row)
369 // we use explicit division so as to avoid having to worry about running
370 // out of oneOverN table.
371
372 long N = 0;
373 double term = std::exp(-mu);
374 double cdf = term;
375 double cdf0;
376
377 while(cdf <= r) {
378 N++ ;
379 term *= ( mu / N );
380 cdf0 = cdf;
381 cdf += term;
382 if (cdf == cdf0) break; // Can't happen, but just in case...
383 }
384 N1 = N;
385 // std::cout << r << " " << N << " ";
386 // DBG_small = true;
387 rRange = 0; // In this case there is always a second r needed
388
389 } // end of small-r case
390
391
392 // >>> 5 <<<
393 // Assuming r lies within the scope of the row for this mu, find the
394 // largest entry not greater than r. N1 is the N corresponding to the
395 // index a.
396
397 else if ( r < cdfs[ENTRIES-1] ) { // r is also >= cdfs[0]
398
399 //
400 // **** This is the normal code path ****
401 //
402
403 int a = 0; // Highest value of index such that cdfs[a]
404 // is known NOT to be greater than r.
405 int b = ENTRIES - 1; // Lowest value of index such that cdfs[b] is
406 // known to exeed r.
407
408 while (b != (a+1) ) {
409 int c = (a+b+1)>>1;
410 if (r > cdfs[c]) {
411 a = c;
412 } else {
413 b = c;
414 }
415 }
416
417 N1 = Nmin + a;
418 rRange = cdfs[a+1] - cdfs[a];
419 rRemainder = r - cdfs[a];
420
421 //
422 // **** ****
423 //
424
425 } // end of medium-r (normal) case
426
427
428 // >>> 6 <<<
429 // If r exceeds the greatest entry in the table for this mu, then start
430 // from that cdf, and use the series to compute from there until r is
431 // exceeded.
432
433 else { // if ( r >= cdfs[ENTRIES-1] ) {
434
435 // Here, division must be done explicitly, and we must also protect against
436 // roundoff preventing termination.
437
438 //
439 //+++ cdfs[ENTRIES-1] is exp(-mu) sum (mu**m/m! , m=0 to Nmax)
440 //+++ (where Nmax = mu - BELOW + ENTRIES - 1)
441 //+++ cdfs[ENTRIES-1]-cdfs[ENTRIES-2] is exp(-mu) mu**(Nmax)/(Nmax)!
442 //+++ If the sum up to k-1 <= r < sum up to k, then N = k-1
443 //+++ Consider k = Nmax in the above statement:
444 //+++ If cdfs[ENTRIES-2] <= r < cdfs[ENTRIES-1], N would be Nmax-1
445 //+++ But here r >= cdfs[ENTRIES-1] so N >= Nmax
446 //
447
448 // Erroneous:
449 //+++ cdfs[ENTRIES-1] is exp(-mu) sum (mu**m/m! , m=0 to Nmax-1)
450 //+++ cdfs[ENTRIES-1]-cdfs[ENTRIES-2] is exp(-mu) mu**(Nmax-1)/(Nmax-1)!
451 //+++ If a sum up to k-1 <= r < sum up to k, then N = k-1
452 //+++ So if cdfs[ENTRIES-1] were > r, N would be Nmax-1 (or less)
453 //+++ But here r >= cdfs[ENTRIES-1] so N >= Nmax
454 //
455
456 // std::cout << "r = " << r << " mu = " << mu << "\n";
457 long N = Nmax -1;
458 double cdf = cdfs[ENTRIES-1];
459 double term = cdf - cdfs[ENTRIES-2];
460 double cdf0;
461 while(cdf <= r) {
462 N++ ;
463 // std::cout << " N " << N << " term " <<
464 // term << " cdf " << cdf << "\n";
465 term *= ( mu / N );
466 cdf0 = cdf;
467 cdf += term;
468 if (cdf == cdf0) break; // If term gets so small cdf stops increasing,
469 // terminate using that value of N since we
470 // would never reach r.
471 }
472 N1 = N;
473 rRange = 0; // We can't validly omit the second true random
474
475 // N = Nmax -1;
476 // cdf = cdfs[ENTRIES-1];
477 // term = cdf - cdfs[ENTRIES-2];
478 // for (int isxz=0; isxz < 100; isxz++) {
479 // N++ ;
480 // term *= ( mu / N );
481 // cdf0 = cdf;
482 // cdf += term;
483 // }
484 // std::cout.precision(20);
485 // std::cout << "Final sum is " << cdf << "\n";
486
487 } // end of large-r case
488
489
490
491 // >>> 7 <<<
492 // Form a second random, s, based on the position of r within the range
493 // of this table entry to the next entry.
494
495 // However, if this range is very small, then we lose too many bits of
496 // randomness. In that situation, we generate a second random for s.
497
498 double s;
499
500 static const double MINRANGE = .01; // Sacrifice up to two digits of
501 // randomness when using r to produce
502 // a second random s. Leads to up to
503 // .09 extra randoms each time.
504
505 if ( rRange > MINRANGE ) {
506 //
507 // **** This path taken 90% of the time ****
508 //
509 s = rRemainder / rRange;
510 } else {
511 s = e->flat(); // extra true random needed about one time in 10.
512 }
513
514 // >>> 8 <<<
515 // Use the direct summation method to form a second poisson deviate N2
516 // from deltaMu and s.
517
518 N2 = 0;
519 double term = std::exp(-deltaMu);
520 double cdf = term;
521
522 if ( s < (1 - 1.0E-10) ) {
523 //
524 // This is the normal path:
525 //
526 const double* oneOverNptr = oneOverN;
527 while( cdf <= s ) {
528 N2++ ;
529 oneOverNptr++;
530 term *= ( deltaMu * (*oneOverNptr) );
531 cdf += term;
532 }
533 } else { // s is almost 1...
534 while( cdf <= s ) {
535 N2++ ;
536 term *= ( deltaMu / N2 );
537 cdf += term;
538 }
539 } // end of if ( s compared to (1 - 1.0E-10) )
540
541 // >>> 9 <<<
542 // The result is the sum of those two deviates
543
544 // if (DBG_small) {
545 // std::cout << N2 << " " << N1+N2 << "\n";
546 // DBG_small = false;
547 // }
548
549 return N1 + N2;
550
551} // poissonDeviate()
552
553std::ostream & RandPoissonQ::put ( std::ostream & os ) const {
554 long pr=os.precision(20);
555 std::vector<unsigned long> t(2);
556 os << " " << name() << "\n";
557 os << "Uvec" << "\n";
558 t = DoubConv::dto2longs(a0);
559 os << a0 << " " << t[0] << " " << t[1] << "\n";
560 t = DoubConv::dto2longs(a1);
561 os << a1 << " " << t[0] << " " << t[1] << "\n";
562 t = DoubConv::dto2longs(a2);
563 os << a2 << " " << t[0] << " " << t[1] << "\n";
564 t = DoubConv::dto2longs(sigma);
565 os << sigma << " " << t[0] << " " << t[1] << "\n";
567 os.precision(pr);
568 return os;
569#ifdef REMOVED
570 int pr=os.precision(20);
571 os << " " << name() << "\n";
572 os << a0 << " " << a1 << " " << a2 << "\n";
573 os << sigma << "\n";
575 os.precision(pr);
576 return os;
577#endif
578}
579
580std::istream & RandPoissonQ::get ( std::istream & is ) {
581 std::string inName;
582 is >> inName;
583 if (inName != name()) {
584 is.clear(std::ios::badbit | is.rdstate());
585 std::cerr << "Mismatch when expecting to read state of a "
586 << name() << " distribution\n"
587 << "Name found was " << inName
588 << "\nistream is left in the badbit state\n";
589 return is;
590 }
591 if (possibleKeywordInput(is, "Uvec", a0)) {
592 std::vector<unsigned long> t(2);
593 is >> a0 >> t[0] >> t[1]; a0 = DoubConv::longs2double(t);
594 is >> a1 >> t[0] >> t[1]; a1 = DoubConv::longs2double(t);
595 is >> a2 >> t[0] >> t[1]; a2 = DoubConv::longs2double(t);
596 is >> sigma >> t[0] >> t[1]; sigma = DoubConv::longs2double(t);
598 return is;
599 }
600 // is >> a0 encompassed by possibleKeywordInput
601 is >> a1 >> a2 >> sigma;
603 return is;
604}
605
606} // namespace CLHEP
607
static double longs2double(const std::vector< unsigned long > &v)
Definition DoubConv.cc:110
static std::vector< unsigned long > dto2longs(double d)
Definition DoubConv.cc:94
static HepRandomEngine * getTheEngine()
Definition Random.cc:268
static double shoot()
static void shootArray(const int size, long *vect, double mean=1.0)
std::ostream & put(std::ostream &os) const
void fireArray(const int size, long *vect)
std::istream & get(std::istream &is)
HepRandomEngine & engine()
std::string name() const
static const double MAXIMUM_POISSON_DEVIATE
static long shoot(double mean=1.0)
std::ostream & put(std::ostream &os) const
HepRandomEngine * getLocalEngine()
static long shoot(double mean=1.0)
std::istream & get(std::istream &is)
HepRandomEngine & engine()
#define N
Definition crc32.c:57
bool possibleKeywordInput(IS &is, const std::string &key, T &t)
#define CLHEP_THREAD_LOCAL