72 #include "CLHEP/Random/poissonTables.cdat"
87 sigma = std::sqrt(sig2);
97 a1 = std::sqrt (1-2*
a2*
a2*sig2);
121 return (
double)
fire();
125 return (
double)
fire(mean);
156 if ( mean != lastLargeMean ) {
159 double sig2 = mean * (.9998654 - .08346/mean);
160 lastSigma = std::sqrt(sig2);
162 lastA2 = t*(1./6.) + t*t*(1./324.);
163 lastA1 = std::sqrt (1-2*lastA2*lastA2*sig2);
164 lastA0 = mean + .5 - sig2 * lastA2;
172 for(
long*
v = vect;
v != vect + size; ++
v )
180 for(
long*
v = vect;
v != vect + size; ++
v )
185 for(
long*
v = vect;
v != vect + size; ++
v )
197 double sig2 = mu * (.9998654 - .08346/mu);
198 double sig = std::sqrt(sig2);
204 double sa2 = t*(1./6.) + t*t*(1./324.);
205 double sa1 = std::sqrt (1-2*sa2*sa2*sig2);
206 double sa0 = mu + .5 - sig2 * sa2;
217 double A0,
double A1,
double A2,
double sig) {
242 double p = A2*g*g + A1*g + A0;
243 if ( p < 0 )
return 0;
260 double rRemainder = 0;
276 double r = e->
flat();
285 static const double oneOverN[50] =
286 { 0, 1., 1/2., 1/3., 1/4., 1/5., 1/6., 1/7., 1/8., 1/9.,
287 1/10., 1/11., 1/12., 1/13., 1/14., 1/15., 1/16., 1/17., 1/18., 1/19.,
288 1/20., 1/21., 1/22., 1/23., 1/24., 1/25., 1/26., 1/27., 1/28., 1/29.,
289 1/30., 1/31., 1/32., 1/33., 1/34., 1/35., 1/36., 1/37., 1/38., 1/39.,
290 1/40., 1/41., 1/42., 1/43., 1/44., 1/45., 1/46., 1/47., 1/48., 1/49. };
296 double term = std::exp(-mean);
299 if ( r < (1 - 1.0
E-9) ) {
305 const double* oneOverNptr = oneOverN;
309 term *= ( mean * (*oneOverNptr) );
324 term *= ( mean /
N );
327 if (cdf == cdf0)
break;
341 double deltaMu = mean - mu;
343 if (Nmin < 1) Nmin = 1;
344 int Nmax = Nmin + (
ENTRIES - 1);
364 double term = std::exp(-mu);
373 if (cdf == cdf0)
break;
388 else if ( r < cdfs[
ENTRIES-1] ) {
399 while (b != (a+1) ) {
409 rRange = cdfs[a+1] - cdfs[
a];
410 rRemainder = r - cdfs[
a];
450 double term = cdf - cdfs[
ENTRIES-2];
459 if (cdf == cdf0)
break;
491 static const double MINRANGE = .01;
496 if ( rRange > MINRANGE ) {
500 s = rRemainder / rRange;
510 double term = std::exp(-deltaMu);
513 if ( s < (1 - 1.0
E-10) ) {
517 const double* oneOverNptr = oneOverN;
521 term *= ( deltaMu * (*oneOverNptr) );
527 term *= ( deltaMu / N2 );
545 int pr=os.precision(20);
546 std::vector<unsigned long>
t(2);
547 os <<
" " <<
name() <<
"\n";
548 os <<
"Uvec" <<
"\n";
550 os <<
a0 <<
" " << t[0] <<
" " << t[1] <<
"\n";
552 os <<
a1 <<
" " << t[0] <<
" " << t[1] <<
"\n";
554 os <<
a2 <<
" " << t[0] <<
" " << t[1] <<
"\n";
556 os <<
sigma <<
" " << t[0] <<
" " << t[1] <<
"\n";
561 int pr=os.precision(20);
562 os <<
" " <<
name() <<
"\n";
563 os <<
a0 <<
" " <<
a1 <<
" " <<
a2 <<
"\n";
574 if (inName !=
name()) {
575 is.clear(std::ios::badbit | is.rdstate());
576 std::cerr <<
"Mismatch when expecting to read state of a "
577 <<
name() <<
" distribution\n"
578 <<
"Name found was " << inName
579 <<
"\nistream is left in the badbit state\n";
583 std::vector<unsigned long>
t(2);