28 #include "CLHEP/Random/defs.h"
29 #include "CLHEP/Random/RandPoisson.h"
30 #include "CLHEP/Units/PhysicalConstants.h"
31 #include "CLHEP/Random/DoubConv.hh"
40 double RandPoisson::status_st[3] = {0., 0., 0.};
41 double RandPoisson::oldm_st = -1.0;
42 const double RandPoisson::meanMax_st = 2.0E9;
61 static double cof[6] = {76.18009172947146,-86.50532032941677,
62 24.01409824083091, -1.231739572450155,
63 0.1208650973866179e-2, -0.5395239384953e-5};
67 tmp -= (x + 0.5) * std::log(tmp);
68 double ser = 1.000000000190015;
70 for ( j = 0; j <= 5; j++ ) {
74 return -tmp + std::log(2.5066282746310005*ser);
78 double normal (HepRandomEngine* eptr)
83 v1 = 2.0 * eptr->flat() - 1.0;
84 v2 = 2.0 * eptr->flat() - 1.0;
88 fac = std::sqrt(-2.0*std::log(r)/r);
109 if( xm == -1 )
return 0;
119 t *= anEngine->
flat();
125 sq = std::sqrt(2.0*xm);
127 g1 = xm*alxm -
gammln(xm + 1.0);
131 y = std::tan(CLHEP::pi*anEngine->
flat());
135 t = 0.9*(1.0 + y*y)* std::exp(em*alxm -
gammln(em + 1.0) - g1);
136 }
while( anEngine->
flat() > t );
139 em = xm + std::sqrt(xm) * normal (anEngine);
140 if (
static_cast<long>(em) < 0 )
141 em =
static_cast<long>(xm) >= 0 ? xm :
getMaxMean();
149 for(
long* v = vect; v != vect + size; ++v )
169 if( xm == -1 )
return 0;
179 t *= anEngine->
flat();
185 sq = std::sqrt(2.0*xm);
187 g1 = xm*alxm -
gammln(xm + 1.0);
191 y = std::tan(CLHEP::pi*anEngine->
flat());
195 t = 0.9*(1.0 + y*y)* std::exp(em*alxm -
gammln(em + 1.0) - g1);
196 }
while( anEngine->
flat() > t );
199 em = xm + std::sqrt(xm) * normal (anEngine);
200 if (
static_cast<long>(em) < 0 )
201 em =
static_cast<long>(xm) >= 0 ? xm :
getMaxMean();
208 long* vect,
double m1)
210 for(
long* v = vect; v != vect + size; ++v )
211 *v =
shoot(anEngine,m1);
232 if( xm == -1 )
return 0;
242 t *= localEngine->flat();
248 sq = std::sqrt(2.0*xm);
250 g1 = xm*alxm -
gammln(xm + 1.0);
254 y = std::tan(CLHEP::pi*localEngine->flat());
258 t = 0.9*(1.0 + y*y)* std::exp(em*alxm -
gammln(em + 1.0) - g1);
259 }
while( localEngine->flat() > t );
262 em = xm + std::sqrt(xm) * normal (localEngine.get());
263 if (
static_cast<long>(em) < 0 )
264 em =
static_cast<long>(xm) >= 0 ? xm :
getMaxMean();
266 status[0] = sq; status[1] = alxm; status[2] = g1;
272 for(
long* v = vect; v != vect + size; ++v )
278 for(
long* v = vect; v != vect + size; ++v )
283 int pr=os.precision(20);
284 std::vector<unsigned long> t(2);
285 os <<
" " <<
name() <<
"\n";
286 os <<
"Uvec" <<
"\n";
288 os <<
meanMax <<
" " << t[0] <<
" " << t[1] <<
"\n";
290 os <<
defaultMean <<
" " << t[0] <<
" " << t[1] <<
"\n";
292 os << status[0] <<
" " << t[0] <<
" " << t[1] <<
"\n";
294 os << status[1] <<
" " << t[0] <<
" " << t[1] <<
"\n";
296 os << status[2] <<
" " << t[0] <<
" " << t[1] <<
"\n";
298 os << oldm <<
" " << t[0] <<
" " << t[1] <<
"\n";
302 int pr=os.precision(20);
303 os <<
" " <<
name() <<
"\n";
305 os << status[0] <<
" " << status[1] <<
" " << status[2] <<
"\n";
314 if (inName !=
name()) {
315 is.clear(std::ios::badbit | is.rdstate());
316 std::cerr <<
"Mismatch when expecting to read state of a "
317 <<
name() <<
" distribution\n"
318 <<
"Name found was " << inName
319 <<
"\nistream is left in the badbit state\n";
323 std::vector<unsigned long> t(2);
333 is >>
defaultMean >> status[0] >> status[1] >> status[2];
static double longs2double(const std::vector< unsigned long > &v)
static std::vector< unsigned long > dto2longs(double d)
static HepRandomEngine * getTheEngine()
std::ostream & put(std::ostream &os) const
static double getOldMean()
static double * getPStatus()
static double getMaxMean()
static long shoot(double m=1.0)
static void setOldMean(double val)
static void setPStatus(double sq, double alxm, double g1)
static void shootArray(const int size, long *vect, double m=1.0)
void fireArray(const int size, long *vect)
std::istream & get(std::istream &is)
HepRandomEngine & engine()
bool possibleKeywordInput(IS &is, const std::string &key, T &t)