CLHEP VERSION Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // $Id: RandPoisson.cc,v 1.7 2010/06/16 17:24:53 garren Exp $ 00002 // -*- C++ -*- 00003 // 00004 // ----------------------------------------------------------------------- 00005 // HEP Random 00006 // --- RandPoisson --- 00007 // class implementation file 00008 // ----------------------------------------------------------------------- 00009 // This file is part of Geant4 (simulation toolkit for HEP). 00010 00011 // ======================================================================= 00012 // Gabriele Cosmo - Created: 5th September 1995 00013 // - Added not static Shoot() method: 17th May 1996 00014 // - Algorithm now operates on doubles: 31st Oct 1996 00015 // - Added methods to shoot arrays: 28th July 1997 00016 // - Added check in case xm=-1: 4th February 1998 00017 // J.Marraffino - Added default mean as attribute and 00018 // operator() with mean: 16th Feb 1998 00019 // Gabriele Cosmo - Relocated static data from HepRandom: 5th Jan 1999 00020 // M Fischler - put and get to/from streams 12/15/04 00021 // M Fischler - put/get to/from streams uses pairs of ulongs when 00022 // + storing doubles avoid problems with precision 00023 // 4/14/05 00024 // Mark Fischler - Repaired BUG - when mean > 2 billion, was returning 00025 // mean instead of the proper value. 01/13/06 00026 // ======================================================================= 00027 00028 #include "CLHEP/Random/defs.h" 00029 #include "CLHEP/Random/RandPoisson.h" 00030 #include "CLHEP/Units/PhysicalConstants.h" 00031 #include "CLHEP/Random/DoubConv.hh" 00032 #include <cmath> // for std::floor() 00033 00034 namespace CLHEP { 00035 00036 std::string RandPoisson::name() const {return "RandPoisson";} 00037 HepRandomEngine & RandPoisson::engine() {return *localEngine;} 00038 00039 // Initialisation of static data 00040 double RandPoisson::status_st[3] = {0., 0., 0.}; 00041 double RandPoisson::oldm_st = -1.0; 00042 const double RandPoisson::meanMax_st = 2.0E9; 00043 00044 RandPoisson::~RandPoisson() { 00045 } 00046 00047 double RandPoisson::operator()() { 00048 return double(fire( defaultMean )); 00049 } 00050 00051 double RandPoisson::operator()( double mean ) { 00052 return double(fire( mean )); 00053 } 00054 00055 double gammln(double xx) { 00056 00057 // Returns the value ln(Gamma(xx) for xx > 0. Full accuracy is obtained for 00058 // xx > 1. For 0 < xx < 1. the reflection formula (6.1.4) can be used first. 00059 // (Adapted from Numerical Recipes in C) 00060 00061 static double cof[6] = {76.18009172947146,-86.50532032941677, 00062 24.01409824083091, -1.231739572450155, 00063 0.1208650973866179e-2, -0.5395239384953e-5}; 00064 int j; 00065 double x = xx - 1.0; 00066 double tmp = x + 5.5; 00067 tmp -= (x + 0.5) * std::log(tmp); 00068 double ser = 1.000000000190015; 00069 00070 for ( j = 0; j <= 5; j++ ) { 00071 x += 1.0; 00072 ser += cof[j]/x; 00073 } 00074 return -tmp + std::log(2.5066282746310005*ser); 00075 } 00076 00077 static 00078 double normal (HepRandomEngine* eptr) // mf 1/13/06 00079 { 00080 double r; 00081 double v1,v2,fac; 00082 do { 00083 v1 = 2.0 * eptr->flat() - 1.0; 00084 v2 = 2.0 * eptr->flat() - 1.0; 00085 r = v1*v1 + v2*v2; 00086 } while ( r > 1.0 ); 00087 00088 fac = std::sqrt(-2.0*std::log(r)/r); 00089 return v2*fac; 00090 } 00091 00092 long RandPoisson::shoot(double xm) { 00093 00094 // Returns as a floating-point number an integer value that is a random 00095 // deviation drawn from a Poisson distribution of mean xm, using flat() 00096 // as a source of uniform random numbers. 00097 // (Adapted from Numerical Recipes in C) 00098 00099 double em, t, y; 00100 double sq, alxm, g1; 00101 double om = getOldMean(); 00102 HepRandomEngine* anEngine = HepRandom::getTheEngine(); 00103 00104 double* status = getPStatus(); 00105 sq = status[0]; 00106 alxm = status[1]; 00107 g1 = status[2]; 00108 00109 if( xm == -1 ) return 0; 00110 if( xm < 12.0 ) { 00111 if( xm != om ) { 00112 setOldMean(xm); 00113 g1 = std::exp(-xm); 00114 } 00115 em = -1; 00116 t = 1.0; 00117 do { 00118 em += 1.0; 00119 t *= anEngine->flat(); 00120 } while( t > g1 ); 00121 } 00122 else if ( xm < getMaxMean() ) { 00123 if ( xm != om ) { 00124 setOldMean(xm); 00125 sq = std::sqrt(2.0*xm); 00126 alxm = std::log(xm); 00127 g1 = xm*alxm - gammln(xm + 1.0); 00128 } 00129 do { 00130 do { 00131 y = std::tan(CLHEP::pi*anEngine->flat()); 00132 em = sq*y + xm; 00133 } while( em < 0.0 ); 00134 em = std::floor(em); 00135 t = 0.9*(1.0 + y*y)* std::exp(em*alxm - gammln(em + 1.0) - g1); 00136 } while( anEngine->flat() > t ); 00137 } 00138 else { 00139 em = xm + std::sqrt(xm) * normal (anEngine); // mf 1/13/06 00140 if ( static_cast<long>(em) < 0 ) 00141 em = static_cast<long>(xm) >= 0 ? xm : getMaxMean(); 00142 } 00143 setPStatus(sq,alxm,g1); 00144 return long(em); 00145 } 00146 00147 void RandPoisson::shootArray(const int size, long* vect, double m1) 00148 { 00149 for( long* v = vect; v != vect + size; ++v ) 00150 *v = shoot(m1); 00151 } 00152 00153 long RandPoisson::shoot(HepRandomEngine* anEngine, double xm) { 00154 00155 // Returns as a floating-point number an integer value that is a random 00156 // deviation drawn from a Poisson distribution of mean xm, using flat() 00157 // of a given Random Engine as a source of uniform random numbers. 00158 // (Adapted from Numerical Recipes in C) 00159 00160 double em, t, y; 00161 double sq, alxm, g1; 00162 double om = getOldMean(); 00163 00164 double* status = getPStatus(); 00165 sq = status[0]; 00166 alxm = status[1]; 00167 g1 = status[2]; 00168 00169 if( xm == -1 ) return 0; 00170 if( xm < 12.0 ) { 00171 if( xm != om ) { 00172 setOldMean(xm); 00173 g1 = std::exp(-xm); 00174 } 00175 em = -1; 00176 t = 1.0; 00177 do { 00178 em += 1.0; 00179 t *= anEngine->flat(); 00180 } while( t > g1 ); 00181 } 00182 else if ( xm < getMaxMean() ) { 00183 if ( xm != om ) { 00184 setOldMean(xm); 00185 sq = std::sqrt(2.0*xm); 00186 alxm = std::log(xm); 00187 g1 = xm*alxm - gammln(xm + 1.0); 00188 } 00189 do { 00190 do { 00191 y = std::tan(CLHEP::pi*anEngine->flat()); 00192 em = sq*y + xm; 00193 } while( em < 0.0 ); 00194 em = std::floor(em); 00195 t = 0.9*(1.0 + y*y)* std::exp(em*alxm - gammln(em + 1.0) - g1); 00196 } while( anEngine->flat() > t ); 00197 } 00198 else { 00199 em = xm + std::sqrt(xm) * normal (anEngine); // mf 1/13/06 00200 if ( static_cast<long>(em) < 0 ) 00201 em = static_cast<long>(xm) >= 0 ? xm : getMaxMean(); 00202 } 00203 setPStatus(sq,alxm,g1); 00204 return long(em); 00205 } 00206 00207 void RandPoisson::shootArray(HepRandomEngine* anEngine, const int size, 00208 long* vect, double m1) 00209 { 00210 for( long* v = vect; v != vect + size; ++v ) 00211 *v = shoot(anEngine,m1); 00212 } 00213 00214 long RandPoisson::fire() { 00215 return long(fire( defaultMean )); 00216 } 00217 00218 long RandPoisson::fire(double xm) { 00219 00220 // Returns as a floating-point number an integer value that is a random 00221 // deviation drawn from a Poisson distribution of mean xm, using flat() 00222 // as a source of uniform random numbers. 00223 // (Adapted from Numerical Recipes in C) 00224 00225 double em, t, y; 00226 double sq, alxm, g1; 00227 00228 sq = status[0]; 00229 alxm = status[1]; 00230 g1 = status[2]; 00231 00232 if( xm == -1 ) return 0; 00233 if( xm < 12.0 ) { 00234 if( xm != oldm ) { 00235 oldm = xm; 00236 g1 = std::exp(-xm); 00237 } 00238 em = -1; 00239 t = 1.0; 00240 do { 00241 em += 1.0; 00242 t *= localEngine->flat(); 00243 } while( t > g1 ); 00244 } 00245 else if ( xm < meanMax ) { 00246 if ( xm != oldm ) { 00247 oldm = xm; 00248 sq = std::sqrt(2.0*xm); 00249 alxm = std::log(xm); 00250 g1 = xm*alxm - gammln(xm + 1.0); 00251 } 00252 do { 00253 do { 00254 y = std::tan(CLHEP::pi*localEngine->flat()); 00255 em = sq*y + xm; 00256 } while( em < 0.0 ); 00257 em = std::floor(em); 00258 t = 0.9*(1.0 + y*y)* std::exp(em*alxm - gammln(em + 1.0) - g1); 00259 } while( localEngine->flat() > t ); 00260 } 00261 else { 00262 em = xm + std::sqrt(xm) * normal (localEngine.get()); // mf 1/13/06 00263 if ( static_cast<long>(em) < 0 ) 00264 em = static_cast<long>(xm) >= 0 ? xm : getMaxMean(); 00265 } 00266 status[0] = sq; status[1] = alxm; status[2] = g1; 00267 return long(em); 00268 } 00269 00270 void RandPoisson::fireArray(const int size, long* vect ) 00271 { 00272 for( long* v = vect; v != vect + size; ++v ) 00273 *v = fire( defaultMean ); 00274 } 00275 00276 void RandPoisson::fireArray(const int size, long* vect, double m1) 00277 { 00278 for( long* v = vect; v != vect + size; ++v ) 00279 *v = fire( m1 ); 00280 } 00281 00282 std::ostream & RandPoisson::put ( std::ostream & os ) const { 00283 int pr=os.precision(20); 00284 std::vector<unsigned long> t(2); 00285 os << " " << name() << "\n"; 00286 os << "Uvec" << "\n"; 00287 t = DoubConv::dto2longs(meanMax); 00288 os << meanMax << " " << t[0] << " " << t[1] << "\n"; 00289 t = DoubConv::dto2longs(defaultMean); 00290 os << defaultMean << " " << t[0] << " " << t[1] << "\n"; 00291 t = DoubConv::dto2longs(status[0]); 00292 os << status[0] << " " << t[0] << " " << t[1] << "\n"; 00293 t = DoubConv::dto2longs(status[1]); 00294 os << status[1] << " " << t[0] << " " << t[1] << "\n"; 00295 t = DoubConv::dto2longs(status[2]); 00296 os << status[2] << " " << t[0] << " " << t[1] << "\n"; 00297 t = DoubConv::dto2longs(oldm); 00298 os << oldm << " " << t[0] << " " << t[1] << "\n"; 00299 os.precision(pr); 00300 return os; 00301 #ifdef REMOVED 00302 int pr=os.precision(20); 00303 os << " " << name() << "\n"; 00304 os << meanMax << " " << defaultMean << "\n"; 00305 os << status[0] << " " << status[1] << " " << status[2] << "\n"; 00306 os.precision(pr); 00307 return os; 00308 #endif 00309 } 00310 00311 std::istream & RandPoisson::get ( std::istream & is ) { 00312 std::string inName; 00313 is >> inName; 00314 if (inName != name()) { 00315 is.clear(std::ios::badbit | is.rdstate()); 00316 std::cerr << "Mismatch when expecting to read state of a " 00317 << name() << " distribution\n" 00318 << "Name found was " << inName 00319 << "\nistream is left in the badbit state\n"; 00320 return is; 00321 } 00322 if (possibleKeywordInput(is, "Uvec", meanMax)) { 00323 std::vector<unsigned long> t(2); 00324 is >> meanMax >> t[0] >> t[1]; meanMax = DoubConv::longs2double(t); 00325 is >> defaultMean >> t[0] >> t[1]; defaultMean = DoubConv::longs2double(t); 00326 is >> status[0] >> t[0] >> t[1]; status[0] = DoubConv::longs2double(t); 00327 is >> status[1] >> t[0] >> t[1]; status[1] = DoubConv::longs2double(t); 00328 is >> status[2] >> t[0] >> t[1]; status[2] = DoubConv::longs2double(t); 00329 is >> oldm >> t[0] >> t[1]; oldm = DoubConv::longs2double(t); 00330 return is; 00331 } 00332 // is >> meanMax encompassed by possibleKeywordInput 00333 is >> defaultMean >> status[0] >> status[1] >> status[2]; 00334 return is; 00335 } 00336 00337 } // namespace CLHEP 00338