| Geant4 Cross Reference |
1 // -*- C++ -*- 1
2 //
3 // -------------------------------------------
4 // HEP Random
5 // --- RandPoissonQ --
6 // class implementation f
7 // -------------------------------------------
8
9 // ===========================================
10 // M. Fischler - Implemented new, much fast
11 // applicable for mu < 100
12 // - Implemented "quick()" methods, shich
13 // new methods for mu < 100 and are a s
14 // approximation for large mu.
15 // M. Fischler - Removed mean=100 from the
16 // uses a value just off the end of the
17 // M Fischler - put and get to/from stream
18 // M Fischler - Utilize RandGaussQ rather
19 // intended by the inclusion of RandGau
20 // introduces a subtle trap in that the
21 // can never be properly captured witho
22 // state of RandGauss! RandGaussQ is,
23 // stateless except for the engine used
24 // M Fischler - Modified use of wrong engine
25 // is called. This flaw was preventing
26 // saving and restoring in the instance
27 // M Fischler - fireArray using defaultMea
28 // M Fischler - put/get to/from streams
29 // + storing doubles avoid problems with
30 // 4/14/05
31 // M Fischler - Modified use of shoot (mean)
32 // shoot(getLocalEngine(), mean) when f
33 // This flaw was causing bad "cross-tal
34 // in CMS, where one used its own engin
35 // used the static generator. 10/18/07
36 //
37 // ===========================================
38
39 #include "CLHEP/Random/RandPoissonQ.h"
40 #include "CLHEP/Random/RandGaussQ.h"
41 #include "CLHEP/Random/DoubConv.h"
42 #include "CLHEP/Random/Stat.h"
43 #include "CLHEP/Utility/thread_local.h"
44 #include <cmath> // for std::pow()
45 #include <iostream>
46 #include <string>
47 #include <vector>
48
49 namespace CLHEP {
50
51 std::string RandPoissonQ::name() const {return
52 HepRandomEngine & RandPoissonQ::engine() {retu
53
54 // Initialization of static data: Note that t
55 // so that saveEngineStatus properly saves all
56
57 // The following MUST MATCH the correspondin
58 // poissonTables.cc) when poissonTables.cdat
59
60 const double RandPoissonQ::FIRST_MU = 10;// lo
61 const double RandPoissonQ::LAST_MU = 95;// hi
62 const double RandPoissonQ::S = 5; // Sp
63 const int RandPoissonQ::BELOW = 30; // St
64 const int RandPoissonQ::ENTRIES = 51; // Nu
65
66 const double RandPoissonQ::MAXIMUM_POISSON_DEV
67 // Careful -- this is NOT the maximum number
68 // a long. It actually should be some large
69 // that.
70
71 // Here comes the big (9K bytes) table, kept
72 // ENTRIES * (FIRST_MU - LAST_MU + 1)/S doub
73
74 static const double poissonTables [ 51 * ( (95
75 #include "CLHEP/Random/poissonTables.cdat"
76 };
77
78 //
79 // Constructors and destructors:
80 //
81
82 RandPoissonQ::~RandPoissonQ() {
83 }
84
85 void RandPoissonQ::setupForDefaultMu() {
86
87 // The following are useful for quick approx
88
89 double sig2 = defaultMean * (.9998654 - .083
90 sigma = std::sqrt(sig2);
91 // sigma for the Guassian which approximates
92 // sqrt (defaultMean).
93 //
94 // The multiplier corrects for fact that dis
95 // [gaussian+.5] increases the second moment
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 wh
104 // The coeffeicients are chosen to match the
105 // true Poisson distribution.
106 //
107 // Actually, if the correction for discretiz
108 // a2 could be taken one order higher by add
109 // the discretization correction is not perf
110 // on the order to 1/mu**2, so adding a thir
111
112 } // setupForDefaultMu()
113
114
115 //
116 // fire, quick, operator(), and shoot methods:
117 //
118
119 long RandPoissonQ::shoot(double xm) {
120 return shoot(getTheEngine(), xm);
121 }
122
123 double RandPoissonQ::operator()() {
124 return (double) fire();
125 }
126
127 double RandPoissonQ::operator()( double mean )
128 return (double) fire(mean);
129 }
130
131 long RandPoissonQ::fire(double mean) {
132 return shoot(getLocalEngine(), mean);
133 }
134
135 long RandPoissonQ::fire() {
136 if ( defaultMean < LAST_MU + S ) {
137 return poissonDeviateSmall ( getLocalEngin
138 } else {
139 return poissonDeviateQuick ( getLocalEngin
140 }
141 } // fire()
142
143 long RandPoissonQ::shoot(HepRandomEngine* anEn
144
145 // The following variables, static to this m
146 // last time a large mean was supplied; they
147 // if consecutive calls use the same mean.
148
149 static CLHEP_THREAD_LOCAL double lastLargeMe
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, mea
158 } else {
159 if ( mean != lastLargeMean ) {
160 // Compute the coefficients defining the
161 // Gaussian to a Poisson for this mean.
162 double sig2 = mean * (.9998654 - .08346/
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*si
167 lastA0 = mean + .5 - sig2 * lastA2;
168 }
169 return poissonDeviateQuick ( anEngine, las
170 }
171
172 } // shoot (anEngine, mean)
173
174 void RandPoissonQ::shootArray(const int size,
175 for( long* v = vect; v != vect + size; ++v )
176 *v = shoot(m);
177 // Note: We could test for m > 100, and i
178 // and sigma and call the appropriate for
179 // But since those are cached anyway, not
180 }
181
182 void RandPoissonQ::shootArray(HepRandomEngine*
183 long* vect, doubl
184 for( long* v = vect; v != vect + size; ++v )
185 *v = shoot(anEngine,m1);
186 }
187
188 void RandPoissonQ::fireArray(const int size, l
189 for( long* v = vect; v != vect + size; ++v )
190 *v = fire( m );
191 }
192
193 void RandPoissonQ::fireArray(const int size, l
194 for( long* v = vect; v != vect + size; ++v )
195 *v = fire( defaultMean );
196 }
197
198
199 // Quick Poisson deviate algorithm used by qui
200
201 long RandPoissonQ::poissonDeviateQuick ( HepRa
202
203 // Compute the coefficients defining the qua
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 dis
209 // [gaussian+.5] increases the second moment
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
218 // The coeffeicients are chosen to match the
219 // true Poisson distribution.
220
221 return poissonDeviateQuick ( e, sa0, sa1, sa
222 }
223
224
225 long RandPoissonQ::poissonDeviateQuick ( HepRa
226 double A0, double A1, double A2, double si
227 //
228 // Quick Poisson deviate algorithm used by qui
229 //
230 // The principle: For very large mu, a poisso
231 // by a gaussian: return the integer part of m
232 // normal. However, this yelds a miserable ap
233 // "large" as 100. The primary problem is tha
234 // supposed to have a skew of 1/mu**2, and the
235 // leads to errors of order as big as 1/mu**2.
236 //
237 // We substitute for the gaussian a quadratic
238 // The expression looks very nearly like mu +
239 // The small positive quadratic term causes th
240 // a positive skew; the -1/6 constant term is
241 // in the mean. By adjusting these two and th
242 // first three moments to high accuracy in 1/m
243 //
244 // The sigma used is not precisely sqrt(mu) si
245 // has a second moment which is slightly large
246 // To compensate, sig is multiplied by a facto
247
248 // double g = RandGauss::shootQuick( e );
249 double g = RandGaussQ::shoot( e ); // Unit
250 g *= sig;
251 double p = A2*g*g + A1*g + A0;
252 if ( p < 0 ) return 0; // Shouldn't ever po
253 // mean should not be less than 100, b
254 // we check due to paranoia.
255 if ( p > MAXIMUM_POISSON_DEVIATE ) p = MAXIM
256
257 return long(p);
258
259 } // poissonDeviateQuick ()
260
261
262
263 long RandPoissonQ::poissonDeviateSmall (HepRan
264 long N1;
265 long N2;
266 // The following are for later use to form a
267 double rRange; // This will hold the in
268 // computed N1 and cdf for N1+1.
269 double rRemainder = 0; // This will hold the
270
271 // Coming in, mean should not be more than L
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
280 }
281
282 // >>> 1 <<<
283 // Generate the first random, which we alwa
284
285 double r = e->flat();
286
287 // >>> 2 <<<
288 // For small mean, below the start of the t
289 // do the series for cdf directly.
290
291 // In this case, since we know the series wi
292 // almost alwaye use a precomputed 1/N array
293
294 static const double oneOverN[50] =
295 { 0, 1., 1/2., 1/3., 1/4., 1/5.,
296 1/10., 1/11., 1/12., 1/13., 1/14., 1/15.,
297 1/20., 1/21., 1/22., 1/23., 1/24., 1/25.,
298 1/30., 1/31., 1/32., 1/33., 1/34., 1/35.,
299 1/40., 1/41., 1/42., 1/43., 1/44., 1/45.,
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
313 // before the 30-th term in the series,
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
327 // off the end of the table of 1/N. Sin
328 // ignore the table and do the identical
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
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
345 // which is no greater than our actual mean
346
347 int rowNumber = int((mean - FIRST_MU)/S);
348 const double * cdfs = &poissonTables [rowNum
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
358 // generate the deviate directly from the s
359
360 if ( r < cdfs[0] ) {
361
362 // In this case, we are tempted to use the
363 // generate a second deviate to account fo
364 // That would be an error, generating a di
365 // at Nmin + (mean-mu)/2 to be detectable
366
367 // Since this case is very rare (never mor
368 // and can happen where N will be large (u
369 // we use explicit division so as to avoid
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
383 }
384 N1 = N;
385 // std::cout << r << " " << N << " ";
386 // DBG_small = true;
387 rRange = 0; // In this case there is alw
388
389 } // end of small-r case
390
391
392 // >>> 5 <<<
393 // Assuming r lies within the scope of the
394 // largest entry not greater than r. N1 is
395 // index a.
396
397 else if ( r < cdfs[ENTRIES-1] ) { // r i
398
399 //
400 // **** This is the normal code path ****
401 //
402
403 int a = 0; // Highest val
404 // is known NOT to be greater than r.
405 int b = ENTRIES - 1; // Lowest value of i
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 t
430 // from that cdf, and use the series to com
431 // exceeded.
432
433 else { // if ( r >= cdfs[ENTRIES-1] ) {
434
435 // Here, division must be done explicitly,
436 // roundoff preventing termination.
437
438 //
439 //+++ cdfs[ENTRIES-1] is exp(-mu) sum (mu**m
440 //+++ (where Nmax = mu - BELOW + ENTRIES - 1
441 //+++ cdfs[ENTRIES-1]-cdfs[ENTRIES-2] is exp
442 //+++ If the sum up to k-1 <= r < sum up to
443 //+++ Consider k = Nmax in the above stateme
444 //+++ If cdfs[ENTRIES-2] <= r < cdfs[ENTRIES
445 //+++ But here r >= cdfs[ENTRIES-1] so N >=
446 //
447
448 // Erroneous:
449 //+++ cdfs[ENTRIES-1] is exp(-mu) sum (mu**m
450 //+++ cdfs[ENTRIES-1]-cdfs[ENTRIES-2] is exp
451 //+++ If a sum up to k-1 <= r < sum up to k,
452 //+++ So if cdfs[ENTRIES-1] were > r, N woul
453 //+++ But here r >= cdfs[ENTRIES-1] so N >=
454 //
455
456 // std::cout << "r = " << r << " mu = " <
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 s
469 // terminate using that value of N sin
470 // would never reach r.
471 }
472 N1 = N;
473 rRange = 0; // We can't validly omit the
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
486
487 } // end of large-r case
488
489
490
491 // >>> 7 <<<
492 // Form a second random, s, based on the po
493 // of this table entry to the next entry.
494
495 // However, if this range is very small, the
496 // randomness. In that situation, we genera
497
498 double s;
499
500 static const double MINRANGE = .01; // Sac
501 // randomness when using r to
502 // a second random s. Leads to up t
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 neede
512 }
513
514 // >>> 8 <<<
515 // Use the direct summation method to form
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 devia
543
544 // if (DBG_small) {
545 // std::cout << N2 << " " << N1+N2 <<
546 // DBG_small = false;
547 // }
548
549 return N1 + N2;
550
551 } // poissonDeviate()
552
553 std::ostream & RandPoissonQ::put ( std::ostrea
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] << "\
560 t = DoubConv::dto2longs(a1);
561 os << a1 << " " << t[0] << " " << t[1] << "\
562 t = DoubConv::dto2longs(a2);
563 os << a2 << " " << t[0] << " " << t[1] << "\
564 t = DoubConv::dto2longs(sigma);
565 os << sigma << " " << t[0] << " " << t[1] <<
566 RandPoisson::put(os);
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";
574 RandPoisson::put(os);
575 os.precision(pr);
576 return os;
577 #endif
578 }
579
580 std::istream & RandPoissonQ::get ( std::istrea
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 r
586 << name() << " distribution\n"
587 << "Name found was " << inName
588 << "\nistream is left in the badbit st
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::l
594 is >> a1 >> t[0] >> t[1]; a1 = DoubConv::l
595 is >> a2 >> t[0] >> t[1]; a2 = DoubConv::l
596 is >> sigma >> t[0] >> t[1]; sigma = DoubC
597 RandPoisson::get(is);
598 return is;
599 }
600 // is >> a0 encompassed by possibleKeywordIn
601 is >> a1 >> a2 >> sigma;
602 RandPoisson::get(is);
603 return is;
604 }
605
606 } // namespace CLHEP
607
608