73 #include "CLHEP/Random/poissonTables.cdat"
88 sigma = std::sqrt(sig2);
98 a1 = std::sqrt (1-2*
a2*
a2*sig2);
122 return (
double)
fire();
126 return (
double)
fire(mean);
157 if ( mean != lastLargeMean ) {
160 double sig2 = mean * (.9998654 - .08346/mean);
161 lastSigma = std::sqrt(sig2);
163 lastA2 = t*(1./6.) + t*t*(1./324.);
164 lastA1 = std::sqrt (1-2*lastA2*lastA2*sig2);
165 lastA0 = mean + .5 - sig2 * lastA2;
173 for(
long* v = vect; v != vect + size; ++v )
181 for(
long* v = vect; v != vect + size; ++v )
186 for(
long* v = vect; v != vect + size; ++v )
198 double sig2 = mu * (.9998654 - .08346/mu);
199 double sig = std::sqrt(sig2);
205 double sa2 = t*(1./6.) + t*t*(1./324.);
206 double sa1 = std::sqrt (1-2*sa2*sa2*sig2);
207 double sa0 = mu + .5 - sig2 * sa2;
218 double A0,
double A1,
double A2,
double sig) {
243 double p = A2*g*g + A1*g + A0;
244 if ( p < 0 )
return 0;
261 double rRemainder = 0;
277 double r = e->
flat();
286 static const double oneOverN[50] =
287 { 0, 1., 1/2., 1/3., 1/4., 1/5., 1/6., 1/7., 1/8., 1/9.,
288 1/10., 1/11., 1/12., 1/13., 1/14., 1/15., 1/16., 1/17., 1/18., 1/19.,
289 1/20., 1/21., 1/22., 1/23., 1/24., 1/25., 1/26., 1/27., 1/28., 1/29.,
290 1/30., 1/31., 1/32., 1/33., 1/34., 1/35., 1/36., 1/37., 1/38., 1/39.,
291 1/40., 1/41., 1/42., 1/43., 1/44., 1/45., 1/46., 1/47., 1/48., 1/49. };
297 double term = std::exp(-mean);
300 if ( r < (1 - 1.0
E-9) ) {
306 const double* oneOverNptr = oneOverN;
310 term *= ( mean * (*oneOverNptr) );
325 term *= ( mean /
N );
328 if (cdf == cdf0)
break;
342 double deltaMu = mean - mu;
344 if (Nmin < 1) Nmin = 1;
345 int Nmax = Nmin + (
ENTRIES - 1);
365 double term = std::exp(-mu);
374 if (cdf == cdf0)
break;
389 else if ( r < cdfs[
ENTRIES-1] ) {
400 while (b != (a+1) ) {
410 rRange = cdfs[a+1] - cdfs[
a];
411 rRemainder = r - cdfs[
a];
451 double term = cdf - cdfs[
ENTRIES-2];
460 if (cdf == cdf0)
break;
492 static const double MINRANGE = .01;
497 if ( rRange > MINRANGE ) {
501 s = rRemainder / rRange;
511 double term = std::exp(-deltaMu);
514 if ( s < (1 - 1.0
E-10) ) {
518 const double* oneOverNptr = oneOverN;
522 term *= ( deltaMu * (*oneOverNptr) );
528 term *= ( deltaMu / N2 );
546 int pr=os.precision(20);
547 std::vector<unsigned long> t(2);
548 os <<
" " <<
name() <<
"\n";
549 os <<
"Uvec" <<
"\n";
551 os <<
a0 <<
" " << t[0] <<
" " << t[1] <<
"\n";
553 os <<
a1 <<
" " << t[0] <<
" " << t[1] <<
"\n";
555 os <<
a2 <<
" " << t[0] <<
" " << t[1] <<
"\n";
557 os <<
sigma <<
" " << t[0] <<
" " << t[1] <<
"\n";
562 int pr=os.precision(20);
563 os <<
" " <<
name() <<
"\n";
564 os <<
a0 <<
" " <<
a1 <<
" " <<
a2 <<
"\n";
575 if (inName !=
name()) {
576 is.clear(std::ios::badbit | is.rdstate());
577 std::cerr <<
"Mismatch when expecting to read state of a "
578 <<
name() <<
" distribution\n"
579 <<
"Name found was " << inName
580 <<
"\nistream is left in the badbit state\n";
584 std::vector<unsigned long> t(2);
bool possibleKeywordInput(IS &is, const std::string &key, T &t)
std::vector< ExP01TrackerHit * > a
static HepRandomEngine * getTheEngine()
static const double LAST_MU
std::istream & get(std::istream &is)
static constexpr double g
static long poissonDeviateSmall(HepRandomEngine *e, double mean)
HepRandomEngine * getLocalEngine()
static const double FIRST_MU
std::ostream & put(std::ostream &os) const
static long shoot(double m=1.0)
static double longs2double(const std::vector< unsigned long > &v)
HepRandomEngine & engine()
typedef int(XMLCALL *XML_NotStandaloneHandler)(void *userData)
static constexpr double s
std::istream & get(std::istream &is)
static const double poissonTables[51 *((95-10)/5+1)]
static std::vector< unsigned long > dto2longs(double d)
#define CLHEP_THREAD_LOCAL
static constexpr double m
std::ostream & put(std::ostream &os) const
static long poissonDeviateQuick(HepRandomEngine *e, double mean)
static const double MAXIMUM_POISSON_DEVIATE
static void shootArray(const int size, long *vect, double m=1.0)
void fireArray(const int size, long *vect)
static long shoot(double m=1.0)
HepRandomEngine & engine()