40#include "CLHEP/Random/defs.h"
41#include "CLHEP/Random/RandPoissonQ.h"
42#include "CLHEP/Random/RandGaussQ.h"
43#include "CLHEP/Random/DoubConv.hh"
44#include "CLHEP/Random/Stat.h"
58const double RandPoissonQ::FIRST_MU = 10;
59const double RandPoissonQ::LAST_MU = 95;
60const double RandPoissonQ::S = 5;
61const int RandPoissonQ::BELOW = 30;
62const int RandPoissonQ::ENTRIES = 51;
72static const double poissonTables [ 51 * ( (95-10)/5 + 1 ) ] = {
73#include "poissonTables.cdat"
84void RandPoissonQ::setupForDefaultMu() {
89 sigma = std::sqrt(sig2);
99 a1 = std::sqrt (1-2*a2*a2*sig2);
123 return (
double)
fire();
127 return (
double)
fire(mean);
138 return poissonDeviateQuick (
getLocalEngine(), a0, a1, a2, sigma );
148 static double lastLargeMean = -1.;
150 static double lastA0;
151 static double lastA1;
152 static double lastA2;
153 static double lastSigma;
155 if ( mean < LAST_MU + S ) {
156 return poissonDeviateSmall ( anEngine, mean );
158 if ( mean != lastLargeMean ) {
161 double sig2 = mean * (.9998654 - .08346/mean);
162 lastSigma = std::sqrt(sig2);
164 lastA2 = t*(1./6.) + t*t*(1./324.);
165 lastA1 = std::sqrt (1-2*lastA2*lastA2*sig2);
166 lastA0 = mean + .5 - sig2 * lastA2;
168 return poissonDeviateQuick ( anEngine, lastA0, lastA1, lastA2, lastSigma );
174 for(
long* v = vect; v != vect + size; ++v )
182 for(
long* v = vect; v != vect + size; ++v )
187 for(
long* v = vect; v != vect + size; ++v )
194long RandPoissonQ::poissonDeviateQuick (
HepRandomEngine *e,
double mu ) {
199 double sig2 = mu * (.9998654 - .08346/mu);
200 double sig = std::sqrt(sig2);
206 double sa2 = t*(1./6.) + t*t*(1./324.);
207 double sa1 = std::sqrt (1-2*sa2*sa2*sig2);
208 double sa0 = mu + .5 - sig2 * sa2;
214 return poissonDeviateQuick ( e, sa0, sa1, sa2, sig );
218long RandPoissonQ::poissonDeviateQuick ( HepRandomEngine *e,
219 double A0,
double A1,
double A2,
double sig) {
244 double p = A2*
g*
g + A1*
g + A0;
245 if ( p < 0 )
return 0;
256long RandPoissonQ::poissonDeviateSmall (HepRandomEngine * e,
double mean) {
262 double rRemainder = 0;
267 if ( mean > LAST_MU + S ) {
278 double r = e->flat();
287 static const double oneOverN[50] =
288 { 0, 1., 1/2., 1/3., 1/4., 1/5., 1/6., 1/7., 1/8., 1/9.,
289 1/10., 1/11., 1/12., 1/13., 1/14., 1/15., 1/16., 1/17., 1/18., 1/19.,
290 1/20., 1/21., 1/22., 1/23., 1/24., 1/25., 1/26., 1/27., 1/28., 1/29.,
291 1/30., 1/31., 1/32., 1/33., 1/34., 1/35., 1/36., 1/37., 1/38., 1/39.,
292 1/40., 1/41., 1/42., 1/43., 1/44., 1/45., 1/46., 1/47., 1/48., 1/49. };
295 if ( mean < FIRST_MU ) {
298 double term = std::exp(-mean);
301 if ( r < (1 - 1.0E-9) ) {
307 const double* oneOverNptr = oneOverN;
311 term *= ( mean * (*oneOverNptr) );
326 term *= ( mean / N );
329 if (cdf == cdf0)
break;
340 int rowNumber = int((mean - FIRST_MU)/S);
341 const double * cdfs = &poissonTables [rowNumber*ENTRIES];
342 double mu = FIRST_MU + rowNumber*S;
343 double deltaMu = mean - mu;
344 int Nmin = int(mu - BELOW);
345 if (Nmin < 1) Nmin = 1;
346 int Nmax = Nmin + (ENTRIES - 1);
366 double term = std::exp(-mu);
375 if (cdf == cdf0)
break;
390 else if ( r < cdfs[ENTRIES-1] ) {
401 while (
b != (
a+1) ) {
411 rRange = cdfs[
a+1] - cdfs[
a];
412 rRemainder = r - cdfs[
a];
451 double cdf = cdfs[ENTRIES-1];
452 double term = cdf - cdfs[ENTRIES-2];
461 if (cdf == cdf0)
break;
493 static const double MINRANGE = .01;
498 if ( rRange > MINRANGE ) {
502 s = rRemainder / rRange;
512 double term = std::exp(-deltaMu);
515 if ( s < (1 - 1.0E-10) ) {
519 const double* oneOverNptr = oneOverN;
523 term *= ( deltaMu * (*oneOverNptr) );
529 term *= ( deltaMu / N2 );
547 int pr=os.precision(20);
548 std::vector<unsigned long> t(2);
549 os <<
" " <<
name() <<
"\n";
550 os <<
"Uvec" <<
"\n";
552 os << a0 <<
" " << t[0] <<
" " << t[1] <<
"\n";
554 os << a1 <<
" " << t[0] <<
" " << t[1] <<
"\n";
556 os << a2 <<
" " << t[0] <<
" " << t[1] <<
"\n";
558 os << sigma <<
" " << t[0] <<
" " << t[1] <<
"\n";
563 int pr=os.precision(20);
564 os <<
" " <<
name() <<
"\n";
565 os << a0 <<
" " << a1 <<
" " << a2 <<
"\n";
576 if (inName !=
name()) {
577 is.clear(std::ios::badbit | is.rdstate());
578 std::cerr <<
"Mismatch when expecting to read state of a "
579 <<
name() <<
" distribution\n"
580 <<
"Name found was " << inName
581 <<
"\nistream is left in the badbit state\n";
585 std::vector<unsigned long> t(2);
594 is >> a1 >> a2 >> sigma;
static double longs2double(const std::vector< unsigned long > &v)
static std::vector< unsigned long > dto2longs(double d)
static HepRandomEngine * getTheEngine()
static long shoot(double m=1.0)
std::ostream & put(std::ostream &os) const
static const double MAXIMUM_POISSON_DEVIATE
void fireArray(const int size, long *vect)
static void shootArray(const int size, long *vect, double m=1.0)
std::istream & get(std::istream &is)
HepRandomEngine & engine()
std::ostream & put(std::ostream &os) const
HepRandomEngine * getLocalEngine()
static long shoot(double m=1.0)
std::istream & get(std::istream &is)
HepRandomEngine & engine()
bool possibleKeywordInput(IS &is, const std::string &key, T &t)