40 #include "CLHEP/Random/RandPoissonQ.h"
41 #include "CLHEP/Random/RandGaussQ.h"
42 #include "CLHEP/Random/DoubConv.h"
43 #include "CLHEP/Random/Stat.h"
44 #include "CLHEP/Utility/thread_local.h"
58 const double RandPoissonQ::FIRST_MU = 10;
59 const double RandPoissonQ::LAST_MU = 95;
61 const int RandPoissonQ::BELOW = 30;
62 const int RandPoissonQ::ENTRIES = 51;
64 const double RandPoissonQ::MAXIMUM_POISSON_DEVIATE = 2.0E9;
73 #include "CLHEP/Random/poissonTables.cdat"
80 RandPoissonQ::~RandPoissonQ() {
83 void RandPoissonQ::setupForDefaultMu() {
87 double sig2 = defaultMean * (.9998654 - .08346/defaultMean);
88 sigma = std::sqrt(sig2);
98 a1 = std::sqrt (1-2*
a2*
a2*sig2);
99 a0 = defaultMean + .5 - sig2 *
a2;
118 return shoot(getTheEngine(), xm);
121 double RandPoissonQ::operator()() {
122 return (
double) fire();
125 double RandPoissonQ::operator()(
double mean ) {
126 return (
double) fire(mean);
129 long RandPoissonQ::fire(
double mean) {
130 return shoot(getLocalEngine(), mean);
133 long RandPoissonQ::fire() {
134 if ( defaultMean < LAST_MU +
S ) {
135 return poissonDeviateSmall ( getLocalEngine(), defaultMean );
137 return poissonDeviateQuick ( getLocalEngine(),
a0,
a1,
a2, sigma );
147 static CLHEP_THREAD_LOCAL
double lastLargeMean = -1.;
149 static CLHEP_THREAD_LOCAL
double lastA0;
150 static CLHEP_THREAD_LOCAL
double lastA1;
151 static CLHEP_THREAD_LOCAL
double lastA2;
152 static CLHEP_THREAD_LOCAL
double lastSigma;
154 if ( mean < LAST_MU +
S ) {
155 return poissonDeviateSmall ( anEngine, 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;
167 return poissonDeviateQuick ( anEngine, lastA0, lastA1, lastA2, lastSigma );
172 void RandPoissonQ::shootArray(
const int size,
long* vect,
double m) {
173 for(
long* v = vect; v != vect + size; ++v )
180 void RandPoissonQ::fireArray(
const int size,
long* vect,
double m) {
181 for(
long* v = vect; v != vect + size; ++v )
185 void RandPoissonQ::fireArray(
const int size,
long* vect) {
186 for(
long* v = vect; v != vect + size; ++v )
187 *v = fire( defaultMean );
193 long RandPoissonQ::poissonDeviateQuick ( HepRandomEngine *e,
double mu ) {
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;
213 return poissonDeviateQuick ( e, sa0, sa1, sa2, sig );
217 long RandPoissonQ::poissonDeviateQuick ( HepRandomEngine *e,
218 double A0,
double A1,
double A2,
double sig) {
243 double p = A2*g*g + A1*g + A0;
244 if ( p < 0 )
return 0;
247 if ( p > MAXIMUM_POISSON_DEVIATE ) p = MAXIMUM_POISSON_DEVIATE;
255 long RandPoissonQ::poissonDeviateSmall (HepRandomEngine * e,
double mean) {
261 double rRemainder = 0;
266 if ( mean > LAST_MU +
S ) {
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. };
294 if ( mean < FIRST_MU ) {
297 double term = std::exp(-mean);
300 if ( r < (1 - 1.0E-9) ) {
306 const double* oneOverNptr = oneOverN;
310 term *= ( mean * (*oneOverNptr) );
325 term *= ( mean / N );
328 if (cdf == cdf0)
break;
339 int rowNumber = int((mean - FIRST_MU)/
S);
341 double mu = FIRST_MU + rowNumber*
S;
342 double deltaMu = mean - mu;
343 int Nmin = int(mu - BELOW);
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];
450 double cdf = cdfs[ENTRIES-1];
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.0E-10) ) {
518 const double* oneOverNptr = oneOverN;
522 term *= ( deltaMu * (*oneOverNptr) );
528 term *= ( deltaMu / N2 );
545 std::ostream & RandPoissonQ::put ( std::ostream & os )
const {
546 int pr=os.precision(20);
547 std::vector<unsigned long> t(2);
548 os <<
" " <<
name() <<
"\n";
549 os <<
"Uvec" <<
"\n";
550 t = DoubConv::dto2longs(
a0);
551 os <<
a0 <<
" " << t[0] <<
" " << t[1] <<
"\n";
552 t = DoubConv::dto2longs(
a1);
553 os <<
a1 <<
" " << t[0] <<
" " << t[1] <<
"\n";
554 t = DoubConv::dto2longs(
a2);
555 os <<
a2 <<
" " << t[0] <<
" " << t[1] <<
"\n";
556 t = DoubConv::dto2longs(sigma);
557 os << sigma <<
" " << t[0] <<
" " << t[1] <<
"\n";
558 RandPoisson::put(os);
562 int pr=os.precision(20);
563 os <<
" " <<
name() <<
"\n";
564 os <<
a0 <<
" " <<
a1 <<
" " <<
a2 <<
"\n";
566 RandPoisson::put(os);
572 std::istream & RandPoissonQ::get ( std::istream & is ) {
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";
583 if (possibleKeywordInput(is,
"Uvec",
a0)) {
584 std::vector<unsigned long> t(2);
585 is >>
a0 >> t[0] >> t[1];
a0 = DoubConv::longs2double(t);
586 is >>
a1 >> t[0] >> t[1];
a1 = DoubConv::longs2double(t);
587 is >> a2 >> t[0] >> t[1]; a2 = DoubConv::longs2double(t);
588 is >> sigma >> t[0] >> t[1]; sigma = DoubConv::longs2double(t);
589 RandPoisson::get(is);
593 is >>
a1 >> a2 >> sigma;
594 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)]