Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/externals/clhep/src/RandPoissonQ.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 
 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 "RandPoissonQ";}
 52 HepRandomEngine & RandPoissonQ::engine() {return RandPoisson::engine();}
 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 
 60 const double RandPoissonQ::FIRST_MU = 10;// lowest mu value in table
 61 const double RandPoissonQ::LAST_MU =  95;// highest mu value
 62 const double RandPoissonQ::S = 5;        // Spacing between mu values
 63 const int RandPoissonQ::BELOW = 30;      // Starting point for N is at mu - BELOW
 64 const int RandPoissonQ::ENTRIES = 51;    // Number of entries in each mu row
 65 
 66 const 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 
 74 static const double poissonTables [ 51 * ( (95-10)/5 + 1 ) ] = {
 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 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 
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 ( getLocalEngine(), defaultMean );
138   } else {
139     return poissonDeviateQuick ( getLocalEngine(), a0, a1, a2, sigma );
140   }
141 } // fire()
142 
143 long 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 
174 void 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 
182 void 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 
188 void RandPoissonQ::fireArray(const int size, long* vect, double m) {
189   for( long* v = vect; v != vect + size; ++v )
190     *v = fire( m );
191 }
192 
193 void 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 
201 long 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 
225 long 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.
255   if ( p > MAXIMUM_POISSON_DEVIATE ) p = MAXIMUM_POISSON_DEVIATE;
256 
257   return long(p);
258 
259 } // poissonDeviateQuick ()
260 
261 
262 
263 long 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 
553 std::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";
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::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); 
597     RandPoisson::get(is);
598     return is;
599   }
600   // is >> a0 encompassed by possibleKeywordInput
601   is >> a1 >> a2 >> sigma;
602   RandPoisson::get(is);
603   return is;
604 }
605 
606 }  // namespace CLHEP
607 
608