28 #include "CLHEP/Random/RandPoisson.h"
29 #include "CLHEP/Units/PhysicalConstants.h"
30 #include "CLHEP/Random/DoubConv.h"
39 CLHEP_THREAD_LOCAL
double RandPoisson::status_st[3] = {0., 0., 0.};
40 CLHEP_THREAD_LOCAL
double RandPoisson::oldm_st = -1.0;
41 const double RandPoisson::meanMax_st = 2.0E9;
43 RandPoisson::~RandPoisson() {
46 double RandPoisson::operator()() {
47 return double(fire( defaultMean ));
50 double RandPoisson::operator()(
double mean ) {
51 return double(fire( mean ));
60 static const double cof[6] = {76.18009172947146,-86.50532032941677,
61 24.01409824083091, -1.231739572450155,
62 0.1208650973866179e-2, -0.5395239384953e-5};
66 tmp -= (x + 0.5) * std::log(tmp);
67 double ser = 1.000000000190015;
69 for ( j = 0; j <= 5; j++ ) {
73 return -tmp + std::log(2.5066282746310005*ser);
77 double normal (HepRandomEngine* eptr)
82 v1 = 2.0 * eptr->flat() - 1.0;
83 v2 = 2.0 * eptr->flat() - 1.0;
87 fac = std::sqrt(-2.0*std::log(r)/r);
100 double om = getOldMean();
101 HepRandomEngine* anEngine = HepRandom::getTheEngine();
103 double* status = getPStatus();
108 if( xm == -1 )
return 0;
118 t *= anEngine->flat();
121 else if ( xm < getMaxMean() ) {
124 sq = std::sqrt(2.0*xm);
126 g1 = xm*alxm -
gammln(xm + 1.0);
130 y = std::tan(
CLHEP::pi*anEngine->flat());
134 t = 0.9*(1.0 + y*y)* std::exp(em*alxm -
gammln(em + 1.0) - g1);
135 }
while( anEngine->flat() > t );
138 em = xm + std::sqrt(xm) *
normal (anEngine);
139 if ( static_cast<long>(em) < 0 )
140 em =
static_cast<long>(xm) >= 0 ? xm : getMaxMean();
142 setPStatus(sq,alxm,g1);
146 void RandPoisson::shootArray(
const int size,
long* vect,
double m1)
148 for(
long* v = vect; v != vect + size; ++v )
161 double om = getOldMean();
163 double* status = getPStatus();
168 if( xm == -1 )
return 0;
178 t *= anEngine->flat();
181 else if ( xm < getMaxMean() ) {
184 sq = std::sqrt(2.0*xm);
186 g1 = xm*alxm -
gammln(xm + 1.0);
190 y = std::tan(
CLHEP::pi*anEngine->flat());
194 t = 0.9*(1.0 + y*y)* std::exp(em*alxm -
gammln(em + 1.0) - g1);
195 }
while( anEngine->flat() > t );
198 em = xm + std::sqrt(xm) *
normal (anEngine);
199 if ( static_cast<long>(em) < 0 )
200 em =
static_cast<long>(xm) >= 0 ? xm : getMaxMean();
202 setPStatus(sq,alxm,g1);
206 void RandPoisson::shootArray(HepRandomEngine* anEngine,
const int size,
207 long* vect,
double m1)
209 for(
long* v = vect; v != vect + size; ++v )
210 *v =
shoot(anEngine,m1);
213 long RandPoisson::fire() {
214 return long(fire( defaultMean ));
217 long RandPoisson::fire(
double xm) {
231 if( xm == -1 )
return 0;
241 t *= localEngine->flat();
244 else if ( xm < meanMax ) {
247 sq = std::sqrt(2.0*xm);
249 g1 = xm*alxm -
gammln(xm + 1.0);
253 y = std::tan(
CLHEP::pi*localEngine->flat());
257 t = 0.9*(1.0 + y*y)* std::exp(em*alxm -
gammln(em + 1.0) - g1);
258 }
while( localEngine->flat() > t );
261 em = xm + std::sqrt(xm) *
normal (localEngine.get());
262 if ( static_cast<long>(em) < 0 )
263 em = static_cast<long>(xm) >= 0 ? xm : getMaxMean();
265 status[0] = sq; status[1] = alxm; status[2] = g1;
269 void RandPoisson::fireArray(
const int size,
long* vect )
271 for(
long* v = vect; v != vect + size; ++v )
272 *v = fire( defaultMean );
275 void RandPoisson::fireArray(
const int size,
long* vect,
double m1)
277 for(
long* v = vect; v != vect + size; ++v )
281 std::ostream & RandPoisson::put ( std::ostream & os )
const {
282 int pr=os.precision(20);
283 std::vector<unsigned long> t(2);
284 os <<
" " <<
name() <<
"\n";
285 os <<
"Uvec" <<
"\n";
286 t = DoubConv::dto2longs(meanMax);
287 os << meanMax <<
" " << t[0] <<
" " << t[1] <<
"\n";
288 t = DoubConv::dto2longs(defaultMean);
289 os << defaultMean <<
" " << t[0] <<
" " << t[1] <<
"\n";
290 t = DoubConv::dto2longs(status[0]);
291 os << status[0] <<
" " << t[0] <<
" " << t[1] <<
"\n";
292 t = DoubConv::dto2longs(status[1]);
293 os << status[1] <<
" " << t[0] <<
" " << t[1] <<
"\n";
294 t = DoubConv::dto2longs(status[2]);
295 os << status[2] <<
" " << t[0] <<
" " << t[1] <<
"\n";
296 t = DoubConv::dto2longs(oldm);
297 os << oldm <<
" " << t[0] <<
" " << t[1] <<
"\n";
302 std::istream & RandPoisson::get ( std::istream & is ) {
305 if (inName !=
name()) {
306 is.clear(std::ios::badbit | is.rdstate());
307 std::cerr <<
"Mismatch when expecting to read state of a "
308 <<
name() <<
" distribution\n"
309 <<
"Name found was " << inName
310 <<
"\nistream is left in the badbit state\n";
313 if (possibleKeywordInput(is,
"Uvec", meanMax)) {
314 std::vector<unsigned long> t(2);
315 is >> meanMax >> t[0] >> t[1]; meanMax = DoubConv::longs2double(t);
316 is >> defaultMean >> t[0] >> t[1]; defaultMean = DoubConv::longs2double(t);
317 is >> status[0] >> t[0] >> t[1]; status[0] = DoubConv::longs2double(t);
318 is >> status[1] >> t[0] >> t[1]; status[1] = DoubConv::longs2double(t);
319 is >> status[2] >> t[0] >> t[1]; status[2] = DoubConv::longs2double(t);
320 is >> oldm >> t[0] >> t[1]; oldm = DoubConv::longs2double(t);
324 is >> defaultMean >> status[0] >> status[1] >> status[2];
ThreeVector shoot(const G4int Ap, const G4int Af)
static int engine(pchar, pchar, double &, pchar &, const dic_type &)
static double normal(HepRandomEngine *eptr)
const G4double x[NPOINTSGL]
static const G4double fac