CLHEP 2.0.4.7 Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // $Id: RandPoisson.cc,v 1.5.4.3 2006/01/17 19:28:46 fischler 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 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 if ( deleteEngine ) delete localEngine; 00046 } 00047 00048 RandPoisson::RandPoisson(const RandPoisson& right) 00049 : meanMax(right.meanMax), defaultMean(right.defaultMean) 00050 {;} 00051 00052 double RandPoisson::operator()() { 00053 return double(fire( defaultMean )); 00054 } 00055 00056 double RandPoisson::operator()( double mean ) { 00057 return double(fire( mean )); 00058 } 00059 00060 double gammln(double xx) { 00061 00062 // Returns the value ln(Gamma(xx) for xx > 0. Full accuracy is obtained for 00063 // xx > 1. For 0 < xx < 1. the reflection formula (6.1.4) can be used first. 00064 // (Adapted from Numerical Recipes in C) 00065 00066 static double cof[6] = {76.18009172947146,-86.50532032941677, 00067 24.01409824083091, -1.231739572450155, 00068 0.1208650973866179e-2, -0.5395239384953e-5}; 00069 int j; 00070 double x = xx - 1.0; 00071 double tmp = x + 5.5; 00072 tmp -= (x + 0.5) * log(tmp); 00073 double ser = 1.000000000190015; 00074 00075 for ( j = 0; j <= 5; j++ ) { 00076 x += 1.0; 00077 ser += cof[j]/x; 00078 } 00079 return -tmp + log(2.5066282746310005*ser); 00080 } 00081 00082 static 00083 double normal (HepRandomEngine* eptr) // mf 1/13/06 00084 { 00085 double r; 00086 double v1,v2,fac; 00087 do { 00088 v1 = 2.0 * eptr->flat() - 1.0; 00089 v2 = 2.0 * eptr->flat() - 1.0; 00090 r = v1*v1 + v2*v2; 00091 } while ( r > 1.0 ); 00092 00093 fac = sqrt(-2.0*log(r)/r); 00094 return v2*fac; 00095 } 00096 00097 long RandPoisson::shoot(double xm) { 00098 00099 // Returns as a floating-point number an integer value that is a random 00100 // deviation drawn from a Poisson distribution of mean xm, using flat() 00101 // as a source of uniform random numbers. 00102 // (Adapted from Numerical Recipes in C) 00103 00104 double em, t, y; 00105 double sq, alxm, g; 00106 double om = getOldMean(); 00107 HepRandomEngine* anEngine = HepRandom::getTheEngine(); 00108 00109 double* status = getPStatus(); 00110 sq = status[0]; 00111 alxm = status[1]; 00112 g = status[2]; 00113 00114 if( xm == -1 ) return 0; 00115 if( xm < 12.0 ) { 00116 if( xm != om ) { 00117 setOldMean(xm); 00118 g = exp(-xm); 00119 } 00120 em = -1; 00121 t = 1.0; 00122 do { 00123 em += 1.0; 00124 t *= anEngine->flat(); 00125 } while( t > g ); 00126 } 00127 else if ( xm < getMaxMean() ) { 00128 if ( xm != om ) { 00129 setOldMean(xm); 00130 sq = sqrt(2.0*xm); 00131 alxm = log(xm); 00132 g = xm*alxm - gammln(xm + 1.0); 00133 } 00134 do { 00135 do { 00136 y = tan(CLHEP::pi*anEngine->flat()); 00137 em = sq*y + xm; 00138 } while( em < 0.0 ); 00139 em = floor(em); 00140 t = 0.9*(1.0 + y*y)* exp(em*alxm - gammln(em + 1.0) - g); 00141 } while( anEngine->flat() > t ); 00142 } 00143 else { 00144 em = xm + sqrt(xm) * normal (anEngine); // mf 1/13/06 00145 if ( static_cast<long>(em) < 0 ) 00146 em = static_cast<long>(xm) >= 0 ? xm : getMaxMean(); 00147 } 00148 setPStatus(sq,alxm,g); 00149 return long(em); 00150 } 00151 00152 void RandPoisson::shootArray(const int size, long* vect, double m) 00153 { 00154 int i; 00155 00156 for (i=0; i<size; ++i) 00157 vect[i] = shoot(m); 00158 } 00159 00160 long RandPoisson::shoot(HepRandomEngine* anEngine, double xm) { 00161 00162 // Returns as a floating-point number an integer value that is a random 00163 // deviation drawn from a Poisson distribution of mean xm, using flat() 00164 // of a given Random Engine as a source of uniform random numbers. 00165 // (Adapted from Numerical Recipes in C) 00166 00167 double em, t, y; 00168 double sq, alxm, g; 00169 double om = getOldMean(); 00170 00171 double* status = getPStatus(); 00172 sq = status[0]; 00173 alxm = status[1]; 00174 g = status[2]; 00175 00176 if( xm == -1 ) return 0; 00177 if( xm < 12.0 ) { 00178 if( xm != om ) { 00179 setOldMean(xm); 00180 g = exp(-xm); 00181 } 00182 em = -1; 00183 t = 1.0; 00184 do { 00185 em += 1.0; 00186 t *= anEngine->flat(); 00187 } while( t > g ); 00188 } 00189 else if ( xm < getMaxMean() ) { 00190 if ( xm != om ) { 00191 setOldMean(xm); 00192 sq = sqrt(2.0*xm); 00193 alxm = log(xm); 00194 g = xm*alxm - gammln(xm + 1.0); 00195 } 00196 do { 00197 do { 00198 y = tan(CLHEP::pi*anEngine->flat()); 00199 em = sq*y + xm; 00200 } while( em < 0.0 ); 00201 em = floor(em); 00202 t = 0.9*(1.0 + y*y)* exp(em*alxm - gammln(em + 1.0) - g); 00203 } while( anEngine->flat() > t ); 00204 } 00205 else { 00206 em = xm + sqrt(xm) * normal (anEngine); // mf 1/13/06 00207 if ( static_cast<long>(em) < 0 ) 00208 em = static_cast<long>(xm) >= 0 ? xm : getMaxMean(); 00209 } 00210 setPStatus(sq,alxm,g); 00211 return long(em); 00212 } 00213 00214 void RandPoisson::shootArray(HepRandomEngine* anEngine, const int size, 00215 long* vect, double m) 00216 { 00217 int i; 00218 00219 for (i=0; i<size; ++i) 00220 vect[i] = shoot(anEngine,m); 00221 } 00222 00223 long RandPoisson::fire() { 00224 return long(fire( defaultMean )); 00225 } 00226 00227 long RandPoisson::fire(double xm) { 00228 00229 // Returns as a floating-point number an integer value that is a random 00230 // deviation drawn from a Poisson distribution of mean xm, using flat() 00231 // as a source of uniform random numbers. 00232 // (Adapted from Numerical Recipes in C) 00233 00234 double em, t, y; 00235 double sq, alxm, g; 00236 00237 sq = status[0]; 00238 alxm = status[1]; 00239 g = status[2]; 00240 00241 if( xm == -1 ) return 0; 00242 if( xm < 12.0 ) { 00243 if( xm != oldm ) { 00244 oldm = xm; 00245 g = exp(-xm); 00246 } 00247 em = -1; 00248 t = 1.0; 00249 do { 00250 em += 1.0; 00251 t *= localEngine->flat(); 00252 } while( t > g ); 00253 } 00254 else if ( xm < meanMax ) { 00255 if ( xm != oldm ) { 00256 oldm = xm; 00257 sq = sqrt(2.0*xm); 00258 alxm = log(xm); 00259 g = xm*alxm - gammln(xm + 1.0); 00260 } 00261 do { 00262 do { 00263 y = tan(CLHEP::pi*localEngine->flat()); 00264 em = sq*y + xm; 00265 } while( em < 0.0 ); 00266 em = floor(em); 00267 t = 0.9*(1.0 + y*y)* exp(em*alxm - gammln(em + 1.0) - g); 00268 } while( localEngine->flat() > t ); 00269 } 00270 else { 00271 em = xm + sqrt(xm) * normal (localEngine); // mf 1/13/06 00272 if ( static_cast<long>(em) < 0 ) 00273 em = static_cast<long>(xm) >= 0 ? xm : getMaxMean(); 00274 } 00275 status[0] = sq; status[1] = alxm; status[2] = g; 00276 return long(em); 00277 } 00278 00279 void RandPoisson::fireArray(const int size, long* vect ) 00280 { 00281 int i; 00282 00283 for (i=0; i<size; ++i) 00284 vect[i] = fire( defaultMean ); 00285 } 00286 00287 void RandPoisson::fireArray(const int size, long* vect, double m) 00288 { 00289 int i; 00290 00291 for (i=0; i<size; ++i) 00292 vect[i] = fire( m ); 00293 } 00294 00295 std::ostream & RandPoisson::put ( std::ostream & os ) const { 00296 int pr=os.precision(20); 00297 std::vector<unsigned long> t(2); 00298 os << " " << name() << "\n"; 00299 os << "Uvec" << "\n"; 00300 t = DoubConv::dto2longs(meanMax); 00301 os << meanMax << " " << t[0] << " " << t[1] << "\n"; 00302 t = DoubConv::dto2longs(defaultMean); 00303 os << defaultMean << " " << t[0] << " " << t[1] << "\n"; 00304 t = DoubConv::dto2longs(status[0]); 00305 os << status[0] << " " << t[0] << " " << t[1] << "\n"; 00306 t = DoubConv::dto2longs(status[1]); 00307 os << status[1] << " " << t[0] << " " << t[1] << "\n"; 00308 t = DoubConv::dto2longs(status[2]); 00309 os << status[2] << " " << t[0] << " " << t[1] << "\n"; 00310 t = DoubConv::dto2longs(oldm); 00311 os << oldm << " " << t[0] << " " << t[1] << "\n"; 00312 os.precision(pr); 00313 return os; 00314 #ifdef REMOVED 00315 int pr=os.precision(20); 00316 os << " " << name() << "\n"; 00317 os << meanMax << " " << defaultMean << "\n"; 00318 os << status[0] << " " << status[1] << " " << status[2] << "\n"; 00319 os.precision(pr); 00320 return os; 00321 #endif 00322 } 00323 00324 std::istream & RandPoisson::get ( std::istream & is ) { 00325 std::string inName; 00326 is >> inName; 00327 if (inName != name()) { 00328 is.clear(std::ios::badbit | is.rdstate()); 00329 std::cerr << "Mismatch when expecting to read state of a " 00330 << name() << " distribution\n" 00331 << "Name found was " << inName 00332 << "\nistream is left in the badbit state\n"; 00333 return is; 00334 } 00335 if (possibleKeywordInput(is, "Uvec", meanMax)) { 00336 std::vector<unsigned long> t(2); 00337 is >> meanMax >> t[0] >> t[1]; meanMax = DoubConv::longs2double(t); 00338 is >> defaultMean >> t[0] >> t[1]; defaultMean = DoubConv::longs2double(t); 00339 is >> status[0] >> t[0] >> t[1]; status[0] = DoubConv::longs2double(t); 00340 is >> status[1] >> t[0] >> t[1]; status[1] = DoubConv::longs2double(t); 00341 is >> status[2] >> t[0] >> t[1]; status[2] = DoubConv::longs2double(t); 00342 is >> oldm >> t[0] >> t[1]; oldm = DoubConv::longs2double(t); 00343 return is; 00344 } 00345 // is >> meanMax encompassed by possibleKeywordInput 00346 is >> defaultMean >> status[0] >> status[1] >> status[2]; 00347 return is; 00348 } 00349 00350 } // namespace CLHEP 00351