40 #include "CLHEP/Random/RandPoissonQ.h"
41 #include "CLHEP/Random/RandGaussQ.h"
42 #include "CLHEP/Random/DoubConv.h"
43 #include "CLHEP/Random/Stat.h"
57 const double RandPoissonQ::FIRST_MU = 10;
58 const double RandPoissonQ::LAST_MU = 95;
59 const double RandPoissonQ::S = 5;
60 const int RandPoissonQ::BELOW = 30;
61 const int RandPoissonQ::ENTRIES = 51;
63 const double RandPoissonQ::MAXIMUM_POISSON_DEVIATE = 2.0E9;
72 #include "CLHEP/Random/poissonTables.cdat"
79 RandPoissonQ::~RandPoissonQ() {
82 void RandPoissonQ::setupForDefaultMu() {
86 double sig2 = defaultMean * (.9998654 - .08346/defaultMean);
87 sigma = std::sqrt(sig2);
97 a1 = std::sqrt (1-2*
a2*
a2*sig2);
98 a0 = defaultMean + .5 - sig2 *
a2;
117 return shoot(getTheEngine(), xm);
120 double RandPoissonQ::operator()() {
121 return (
double) fire();
124 double RandPoissonQ::operator()(
double mean ) {
125 return (
double) fire(mean);
128 long RandPoissonQ::fire(
double mean) {
129 return shoot(getLocalEngine(), mean);
132 long RandPoissonQ::fire() {
133 if ( defaultMean < LAST_MU + S ) {
134 return poissonDeviateSmall ( getLocalEngine(), defaultMean );
136 return poissonDeviateQuick ( getLocalEngine(),
a0,
a1,
a2, sigma );
146 static double lastLargeMean = -1.;
148 static double lastA0;
149 static double lastA1;
150 static double lastA2;
151 static double lastSigma;
153 if ( mean < LAST_MU + S ) {
154 return poissonDeviateSmall ( anEngine, 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;
166 return poissonDeviateQuick ( anEngine, lastA0, lastA1, lastA2, lastSigma );
171 void RandPoissonQ::shootArray(
const int size,
long* vect,
double m) {
172 for(
long* v = vect; v != vect + size; ++v )
179 void RandPoissonQ::fireArray(
const int size,
long* vect,
double m) {
180 for(
long* v = vect; v != vect + size; ++v )
184 void RandPoissonQ::fireArray(
const int size,
long* vect) {
185 for(
long* v = vect; v != vect + size; ++v )
186 *v = fire( defaultMean );
192 long RandPoissonQ::poissonDeviateQuick ( HepRandomEngine *e,
double mu ) {
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;
212 return poissonDeviateQuick ( e, sa0, sa1, sa2, sig );
216 long RandPoissonQ::poissonDeviateQuick ( HepRandomEngine *e,
217 double A0,
double A1,
double A2,
double sig) {
242 double p = A2*g*g + A1*g + A0;
243 if ( p < 0 )
return 0;
246 if ( p > MAXIMUM_POISSON_DEVIATE ) p = MAXIMUM_POISSON_DEVIATE;
254 long RandPoissonQ::poissonDeviateSmall (HepRandomEngine * e,
double mean) {
260 double rRemainder = 0;
265 if ( mean > LAST_MU + S ) {
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. };
293 if ( mean < FIRST_MU ) {
296 double term = std::exp(-mean);
299 if ( r < (1 - 1.0E-9) ) {
305 const double* oneOverNptr = oneOverN;
309 term *= ( mean * (*oneOverNptr) );
324 term *= ( mean / N );
327 if (cdf == cdf0)
break;
338 int rowNumber = int((mean - FIRST_MU)/S);
340 double mu = FIRST_MU + rowNumber*S;
341 double deltaMu = mean - mu;
342 int Nmin = int(mu - BELOW);
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];
449 double cdf = cdfs[ENTRIES-1];
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.0E-10) ) {
517 const double* oneOverNptr = oneOverN;
521 term *= ( deltaMu * (*oneOverNptr) );
527 term *= ( deltaMu / N2 );
544 std::ostream & RandPoissonQ::put ( std::ostream & os )
const {
545 int pr=os.precision(20);
546 std::vector<unsigned long> t(2);
547 os <<
" " <<
name() <<
"\n";
548 os <<
"Uvec" <<
"\n";
549 t = DoubConv::dto2longs(
a0);
550 os <<
a0 <<
" " << t[0] <<
" " << t[1] <<
"\n";
551 t = DoubConv::dto2longs(
a1);
552 os <<
a1 <<
" " << t[0] <<
" " << t[1] <<
"\n";
553 t = DoubConv::dto2longs(
a2);
554 os <<
a2 <<
" " << t[0] <<
" " << t[1] <<
"\n";
555 t = DoubConv::dto2longs(sigma);
556 os << sigma <<
" " << t[0] <<
" " << t[1] <<
"\n";
557 RandPoisson::put(os);
561 int pr=os.precision(20);
562 os <<
" " <<
name() <<
"\n";
563 os <<
a0 <<
" " <<
a1 <<
" " <<
a2 <<
"\n";
565 RandPoisson::put(os);
571 std::istream & RandPoissonQ::get ( std::istream & is ) {
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";
582 if (possibleKeywordInput(is,
"Uvec",
a0)) {
583 std::vector<unsigned long> t(2);
584 is >>
a0 >> t[0] >> t[1];
a0 = DoubConv::longs2double(t);
585 is >>
a1 >> t[0] >> t[1];
a1 = DoubConv::longs2double(t);
586 is >> a2 >> t[0] >> t[1]; a2 = DoubConv::longs2double(t);
587 is >> sigma >> t[0] >> t[1]; sigma = DoubConv::longs2double(t);
588 RandPoisson::get(is);
592 is >>
a1 >> a2 >> sigma;
593 RandPoisson::get(is);
ThreeVector shoot(const G4int Ap, const G4int Af)
static int engine(pchar, pchar, double &, pchar &, const dic_type &)
static const double poissonTables[51 *((95-10)/5+1)]