Geant4 Cross Reference |
1 // -*- C++ -*- 1 // -*- C++ -*- 2 // 2 // 3 // ------------------------------------------- 3 // ----------------------------------------------------------------------- 4 // HEP Random 4 // HEP Random 5 // --- RandPoisson --- 5 // --- RandPoisson --- 6 // class implementation f 6 // class implementation file 7 // ------------------------------------------- 7 // ----------------------------------------------------------------------- 8 // This file is part of Geant4 (simulation too 8 // This file is part of Geant4 (simulation toolkit for HEP). 9 9 10 // =========================================== 10 // ======================================================================= 11 // Gabriele Cosmo - Created: 5th September 199 11 // Gabriele Cosmo - Created: 5th September 1995 12 // - Added not static Shoot() m 12 // - Added not static Shoot() method: 17th May 1996 13 // - Algorithm now operates on 13 // - Algorithm now operates on doubles: 31st Oct 1996 14 // - Added methods to shoot arr 14 // - Added methods to shoot arrays: 28th July 1997 15 // - Added check in case xm=-1: 15 // - Added check in case xm=-1: 4th February 1998 16 // J.Marraffino - Added default mean as attr 16 // J.Marraffino - Added default mean as attribute and 17 // operator() with mean: 16th 17 // operator() with mean: 16th Feb 1998 18 // Gabriele Cosmo - Relocated static data from 18 // Gabriele Cosmo - Relocated static data from HepRandom: 5th Jan 1999 19 // M Fischler - put and get to/from stream 19 // M Fischler - put and get to/from streams 12/15/04 20 // M Fischler - put/get to/from streams 20 // M Fischler - put/get to/from streams uses pairs of ulongs when 21 // + storing doubles avoid problems with 21 // + storing doubles avoid problems with precision 22 // 4/14/05 22 // 4/14/05 23 // Mark Fischler - Repaired BUG - when mean > 23 // Mark Fischler - Repaired BUG - when mean > 2 billion, was returning 24 // mean instead of the proper 24 // mean instead of the proper value. 01/13/06 25 // =========================================== 25 // ======================================================================= 26 26 27 #include "CLHEP/Random/RandPoisson.h" 27 #include "CLHEP/Random/RandPoisson.h" 28 #include "CLHEP/Units/PhysicalConstants.h" 28 #include "CLHEP/Units/PhysicalConstants.h" 29 #include "CLHEP/Random/DoubConv.h" 29 #include "CLHEP/Random/DoubConv.h" 30 #include <cmath> // for std::floor() 30 #include <cmath> // for std::floor() 31 #include <iostream> << 32 #include <string> << 33 #include <vector> << 34 31 35 namespace CLHEP { 32 namespace CLHEP { 36 33 37 std::string RandPoisson::name() const {return 34 std::string RandPoisson::name() const {return "RandPoisson";} 38 HepRandomEngine & RandPoisson::engine() {retur 35 HepRandomEngine & RandPoisson::engine() {return *localEngine;} 39 36 40 // Initialisation of static data 37 // Initialisation of static data 41 CLHEP_THREAD_LOCAL double RandPoisson::status_ 38 CLHEP_THREAD_LOCAL double RandPoisson::status_st[3] = {0., 0., 0.}; 42 CLHEP_THREAD_LOCAL double RandPoisson::oldm_st 39 CLHEP_THREAD_LOCAL double RandPoisson::oldm_st = -1.0; 43 const double RandPoisson::meanMax_st = 2.0E9; 40 const double RandPoisson::meanMax_st = 2.0E9; 44 41 45 RandPoisson::~RandPoisson() { 42 RandPoisson::~RandPoisson() { 46 } 43 } 47 44 48 double RandPoisson::operator()() { 45 double RandPoisson::operator()() { 49 return double(fire( defaultMean )); 46 return double(fire( defaultMean )); 50 } 47 } 51 48 52 double RandPoisson::operator()( double mean ) 49 double RandPoisson::operator()( double mean ) { 53 return double(fire( mean )); 50 return double(fire( mean )); 54 } 51 } 55 52 56 double gammln(double xx) { 53 double gammln(double xx) { 57 54 58 // Returns the value ln(Gamma(xx) for xx > 0. 55 // Returns the value ln(Gamma(xx) for xx > 0. Full accuracy is obtained for 59 // xx > 1. For 0 < xx < 1. the reflection form 56 // xx > 1. For 0 < xx < 1. the reflection formula (6.1.4) can be used first. 60 // (Adapted from Numerical Recipes in C) 57 // (Adapted from Numerical Recipes in C) 61 58 62 static const double cof[6] = {76.18009172947 59 static const double cof[6] = {76.18009172947146,-86.50532032941677, 63 24.01409824083091 60 24.01409824083091, -1.231739572450155, 64 0.120865097386617 61 0.1208650973866179e-2, -0.5395239384953e-5}; 65 int j; 62 int j; 66 double x = xx - 1.0; 63 double x = xx - 1.0; 67 double tmp = x + 5.5; 64 double tmp = x + 5.5; 68 tmp -= (x + 0.5) * std::log(tmp); 65 tmp -= (x + 0.5) * std::log(tmp); 69 double ser = 1.000000000190015; 66 double ser = 1.000000000190015; 70 67 71 for ( j = 0; j <= 5; j++ ) { 68 for ( j = 0; j <= 5; j++ ) { 72 x += 1.0; 69 x += 1.0; 73 ser += cof[j]/x; 70 ser += cof[j]/x; 74 } 71 } 75 return -tmp + std::log(2.5066282746310005*se 72 return -tmp + std::log(2.5066282746310005*ser); 76 } 73 } 77 74 78 static 75 static 79 double normal (HepRandomEngine* eptr) // m 76 double normal (HepRandomEngine* eptr) // mf 1/13/06 80 { 77 { 81 double r; 78 double r; 82 double v1,v2,fac; 79 double v1,v2,fac; 83 do { 80 do { 84 v1 = 2.0 * eptr->flat() - 1.0; 81 v1 = 2.0 * eptr->flat() - 1.0; 85 v2 = 2.0 * eptr->flat() - 1.0; 82 v2 = 2.0 * eptr->flat() - 1.0; 86 r = v1*v1 + v2*v2; 83 r = v1*v1 + v2*v2; 87 } while ( r > 1.0 ); 84 } while ( r > 1.0 ); 88 85 89 fac = std::sqrt(-2.0*std::log(r)/r); 86 fac = std::sqrt(-2.0*std::log(r)/r); 90 return v2*fac; 87 return v2*fac; 91 } 88 } 92 89 93 long RandPoisson::shoot(double xm) { 90 long RandPoisson::shoot(double xm) { 94 91 95 // Returns as a floating-point number an integ 92 // Returns as a floating-point number an integer value that is a random 96 // deviation drawn from a Poisson distribution 93 // deviation drawn from a Poisson distribution of mean xm, using flat() 97 // as a source of uniform random numbers. 94 // as a source of uniform random numbers. 98 // (Adapted from Numerical Recipes in C) 95 // (Adapted from Numerical Recipes in C) 99 96 100 double em, t, y; 97 double em, t, y; 101 double sq, alxm, g1; 98 double sq, alxm, g1; 102 double om = getOldMean(); 99 double om = getOldMean(); 103 HepRandomEngine* anEngine = HepRandom::getTh 100 HepRandomEngine* anEngine = HepRandom::getTheEngine(); 104 101 105 double* pstatus = getPStatus(); << 102 double* status = getPStatus(); 106 sq = pstatus[0]; << 103 sq = status[0]; 107 alxm = pstatus[1]; << 104 alxm = status[1]; 108 g1 = pstatus[2]; << 105 g1 = status[2]; 109 106 110 if( xm == -1 ) return 0; 107 if( xm == -1 ) return 0; 111 if( xm < 12.0 ) { 108 if( xm < 12.0 ) { 112 if( xm != om ) { 109 if( xm != om ) { 113 setOldMean(xm); 110 setOldMean(xm); 114 g1 = std::exp(-xm); 111 g1 = std::exp(-xm); 115 } 112 } 116 em = -1; 113 em = -1; 117 t = 1.0; 114 t = 1.0; 118 do { 115 do { 119 em += 1.0; 116 em += 1.0; 120 t *= anEngine->flat(); 117 t *= anEngine->flat(); 121 } while( t > g1 ); 118 } while( t > g1 ); 122 } 119 } 123 else if ( xm < getMaxMean() ) { 120 else if ( xm < getMaxMean() ) { 124 if ( xm != om ) { 121 if ( xm != om ) { 125 setOldMean(xm); 122 setOldMean(xm); 126 sq = std::sqrt(2.0*xm); 123 sq = std::sqrt(2.0*xm); 127 alxm = std::log(xm); 124 alxm = std::log(xm); 128 g1 = xm*alxm - gammln(xm + 1.0); 125 g1 = xm*alxm - gammln(xm + 1.0); 129 } 126 } 130 do { 127 do { 131 do { 128 do { 132 y = std::tan(CLHEP::pi*anEngine->flat()); 129 y = std::tan(CLHEP::pi*anEngine->flat()); 133 em = sq*y + xm; 130 em = sq*y + xm; 134 } while( em < 0.0 ); 131 } while( em < 0.0 ); 135 em = std::floor(em); 132 em = std::floor(em); 136 t = 0.9*(1.0 + y*y)* std::exp(em*alxm - 133 t = 0.9*(1.0 + y*y)* std::exp(em*alxm - gammln(em + 1.0) - g1); 137 } while( anEngine->flat() > t ); 134 } while( anEngine->flat() > t ); 138 } 135 } 139 else { 136 else { 140 em = xm + std::sqrt(xm) * normal (anEngine 137 em = xm + std::sqrt(xm) * normal (anEngine); // mf 1/13/06 141 if ( static_cast<long>(em) < 0 ) 138 if ( static_cast<long>(em) < 0 ) 142 em = static_cast<long>(xm) >= 0 ? xm : g 139 em = static_cast<long>(xm) >= 0 ? xm : getMaxMean(); 143 } 140 } 144 setPStatus(sq,alxm,g1); 141 setPStatus(sq,alxm,g1); 145 return long(em); 142 return long(em); 146 } 143 } 147 144 148 void RandPoisson::shootArray(const int size, l 145 void RandPoisson::shootArray(const int size, long* vect, double m1) 149 { 146 { 150 for( long* v = vect; v != vect + size; ++v ) 147 for( long* v = vect; v != vect + size; ++v ) 151 *v = shoot(m1); 148 *v = shoot(m1); 152 } 149 } 153 150 154 long RandPoisson::shoot(HepRandomEngine* anEng 151 long RandPoisson::shoot(HepRandomEngine* anEngine, double xm) { 155 152 156 // Returns as a floating-point number an integ 153 // Returns as a floating-point number an integer value that is a random 157 // deviation drawn from a Poisson distribution 154 // deviation drawn from a Poisson distribution of mean xm, using flat() 158 // of a given Random Engine as a source of uni 155 // of a given Random Engine as a source of uniform random numbers. 159 // (Adapted from Numerical Recipes in C) 156 // (Adapted from Numerical Recipes in C) 160 157 161 double em, t, y; 158 double em, t, y; 162 double sq, alxm, g1; 159 double sq, alxm, g1; 163 double om = getOldMean(); 160 double om = getOldMean(); 164 161 165 double* pstatus = getPStatus(); << 162 double* status = getPStatus(); 166 sq = pstatus[0]; << 163 sq = status[0]; 167 alxm = pstatus[1]; << 164 alxm = status[1]; 168 g1 = pstatus[2]; << 165 g1 = status[2]; 169 166 170 if( xm == -1 ) return 0; 167 if( xm == -1 ) return 0; 171 if( xm < 12.0 ) { 168 if( xm < 12.0 ) { 172 if( xm != om ) { 169 if( xm != om ) { 173 setOldMean(xm); 170 setOldMean(xm); 174 g1 = std::exp(-xm); 171 g1 = std::exp(-xm); 175 } 172 } 176 em = -1; 173 em = -1; 177 t = 1.0; 174 t = 1.0; 178 do { 175 do { 179 em += 1.0; 176 em += 1.0; 180 t *= anEngine->flat(); 177 t *= anEngine->flat(); 181 } while( t > g1 ); 178 } while( t > g1 ); 182 } 179 } 183 else if ( xm < getMaxMean() ) { 180 else if ( xm < getMaxMean() ) { 184 if ( xm != om ) { 181 if ( xm != om ) { 185 setOldMean(xm); 182 setOldMean(xm); 186 sq = std::sqrt(2.0*xm); 183 sq = std::sqrt(2.0*xm); 187 alxm = std::log(xm); 184 alxm = std::log(xm); 188 g1 = xm*alxm - gammln(xm + 1.0); 185 g1 = xm*alxm - gammln(xm + 1.0); 189 } 186 } 190 do { 187 do { 191 do { 188 do { 192 y = std::tan(CLHEP::pi*anEngine->flat()); 189 y = std::tan(CLHEP::pi*anEngine->flat()); 193 em = sq*y + xm; 190 em = sq*y + xm; 194 } while( em < 0.0 ); 191 } while( em < 0.0 ); 195 em = std::floor(em); 192 em = std::floor(em); 196 t = 0.9*(1.0 + y*y)* std::exp(em*alxm - 193 t = 0.9*(1.0 + y*y)* std::exp(em*alxm - gammln(em + 1.0) - g1); 197 } while( anEngine->flat() > t ); 194 } while( anEngine->flat() > t ); 198 } 195 } 199 else { 196 else { 200 em = xm + std::sqrt(xm) * normal (anEngine 197 em = xm + std::sqrt(xm) * normal (anEngine); // mf 1/13/06 201 if ( static_cast<long>(em) < 0 ) 198 if ( static_cast<long>(em) < 0 ) 202 em = static_cast<long>(xm) >= 0 ? xm : g 199 em = static_cast<long>(xm) >= 0 ? xm : getMaxMean(); 203 } 200 } 204 setPStatus(sq,alxm,g1); 201 setPStatus(sq,alxm,g1); 205 return long(em); 202 return long(em); 206 } 203 } 207 204 208 void RandPoisson::shootArray(HepRandomEngine* 205 void RandPoisson::shootArray(HepRandomEngine* anEngine, const int size, 209 long* vect, doubl 206 long* vect, double m1) 210 { 207 { 211 for( long* v = vect; v != vect + size; ++v ) 208 for( long* v = vect; v != vect + size; ++v ) 212 *v = shoot(anEngine,m1); 209 *v = shoot(anEngine,m1); 213 } 210 } 214 211 215 long RandPoisson::fire() { 212 long RandPoisson::fire() { 216 return long(fire( defaultMean )); 213 return long(fire( defaultMean )); 217 } 214 } 218 215 219 long RandPoisson::fire(double xm) { 216 long RandPoisson::fire(double xm) { 220 217 221 // Returns as a floating-point number an integ 218 // Returns as a floating-point number an integer value that is a random 222 // deviation drawn from a Poisson distribution 219 // deviation drawn from a Poisson distribution of mean xm, using flat() 223 // as a source of uniform random numbers. 220 // as a source of uniform random numbers. 224 // (Adapted from Numerical Recipes in C) 221 // (Adapted from Numerical Recipes in C) 225 222 226 double em, t, y; 223 double em, t, y; 227 double sq, alxm, g1; 224 double sq, alxm, g1; 228 225 229 sq = status[0]; 226 sq = status[0]; 230 alxm = status[1]; 227 alxm = status[1]; 231 g1 = status[2]; 228 g1 = status[2]; 232 229 233 if( xm == -1 ) return 0; 230 if( xm == -1 ) return 0; 234 if( xm < 12.0 ) { 231 if( xm < 12.0 ) { 235 if( xm != oldm ) { 232 if( xm != oldm ) { 236 oldm = xm; 233 oldm = xm; 237 g1 = std::exp(-xm); 234 g1 = std::exp(-xm); 238 } 235 } 239 em = -1; 236 em = -1; 240 t = 1.0; 237 t = 1.0; 241 do { 238 do { 242 em += 1.0; 239 em += 1.0; 243 t *= localEngine->flat(); 240 t *= localEngine->flat(); 244 } while( t > g1 ); 241 } while( t > g1 ); 245 } 242 } 246 else if ( xm < meanMax ) { 243 else if ( xm < meanMax ) { 247 if ( xm != oldm ) { 244 if ( xm != oldm ) { 248 oldm = xm; 245 oldm = xm; 249 sq = std::sqrt(2.0*xm); 246 sq = std::sqrt(2.0*xm); 250 alxm = std::log(xm); 247 alxm = std::log(xm); 251 g1 = xm*alxm - gammln(xm + 1.0); 248 g1 = xm*alxm - gammln(xm + 1.0); 252 } 249 } 253 do { 250 do { 254 do { 251 do { 255 y = std::tan(CLHEP::pi*localEngine->flat()); 252 y = std::tan(CLHEP::pi*localEngine->flat()); 256 em = sq*y + xm; 253 em = sq*y + xm; 257 } while( em < 0.0 ); 254 } while( em < 0.0 ); 258 em = std::floor(em); 255 em = std::floor(em); 259 t = 0.9*(1.0 + y*y)* std::exp(em*alxm - 256 t = 0.9*(1.0 + y*y)* std::exp(em*alxm - gammln(em + 1.0) - g1); 260 } while( localEngine->flat() > t ); 257 } while( localEngine->flat() > t ); 261 } 258 } 262 else { 259 else { 263 em = xm + std::sqrt(xm) * normal (localEng 260 em = xm + std::sqrt(xm) * normal (localEngine.get()); // mf 1/13/06 264 if ( static_cast<long>(em) < 0 ) 261 if ( static_cast<long>(em) < 0 ) 265 em = static_cast<long>(xm) >= 0 ? xm : g 262 em = static_cast<long>(xm) >= 0 ? xm : getMaxMean(); 266 } 263 } 267 status[0] = sq; status[1] = alxm; status[2] 264 status[0] = sq; status[1] = alxm; status[2] = g1; 268 return long(em); 265 return long(em); 269 } 266 } 270 267 271 void RandPoisson::fireArray(const int size, lo 268 void RandPoisson::fireArray(const int size, long* vect ) 272 { 269 { 273 for( long* v = vect; v != vect + size; ++v ) 270 for( long* v = vect; v != vect + size; ++v ) 274 *v = fire( defaultMean ); 271 *v = fire( defaultMean ); 275 } 272 } 276 273 277 void RandPoisson::fireArray(const int size, lo 274 void RandPoisson::fireArray(const int size, long* vect, double m1) 278 { 275 { 279 for( long* v = vect; v != vect + size; ++v ) 276 for( long* v = vect; v != vect + size; ++v ) 280 *v = fire( m1 ); 277 *v = fire( m1 ); 281 } 278 } 282 279 283 std::ostream & RandPoisson::put ( std::ostream 280 std::ostream & RandPoisson::put ( std::ostream & os ) const { 284 long pr=os.precision(20); << 281 int pr=os.precision(20); 285 std::vector<unsigned long> t(2); 282 std::vector<unsigned long> t(2); 286 os << " " << name() << "\n"; 283 os << " " << name() << "\n"; 287 os << "Uvec" << "\n"; 284 os << "Uvec" << "\n"; 288 t = DoubConv::dto2longs(meanMax); 285 t = DoubConv::dto2longs(meanMax); 289 os << meanMax << " " << t[0] << " " << t[1] 286 os << meanMax << " " << t[0] << " " << t[1] << "\n"; 290 t = DoubConv::dto2longs(defaultMean); 287 t = DoubConv::dto2longs(defaultMean); 291 os << defaultMean << " " << t[0] << " " << t 288 os << defaultMean << " " << t[0] << " " << t[1] << "\n"; 292 t = DoubConv::dto2longs(status[0]); 289 t = DoubConv::dto2longs(status[0]); 293 os << status[0] << " " << t[0] << " " << t[1 290 os << status[0] << " " << t[0] << " " << t[1] << "\n"; 294 t = DoubConv::dto2longs(status[1]); 291 t = DoubConv::dto2longs(status[1]); 295 os << status[1] << " " << t[0] << " " << t[1 292 os << status[1] << " " << t[0] << " " << t[1] << "\n"; 296 t = DoubConv::dto2longs(status[2]); 293 t = DoubConv::dto2longs(status[2]); 297 os << status[2] << " " << t[0] << " " << t[1 294 os << status[2] << " " << t[0] << " " << t[1] << "\n"; 298 t = DoubConv::dto2longs(oldm); 295 t = DoubConv::dto2longs(oldm); 299 os << oldm << " " << t[0] << " " << t[1] << 296 os << oldm << " " << t[0] << " " << t[1] << "\n"; 300 os.precision(pr); 297 os.precision(pr); 301 return os; 298 return os; 302 } 299 } 303 300 304 std::istream & RandPoisson::get ( std::istream 301 std::istream & RandPoisson::get ( std::istream & is ) { 305 std::string inName; 302 std::string inName; 306 is >> inName; 303 is >> inName; 307 if (inName != name()) { 304 if (inName != name()) { 308 is.clear(std::ios::badbit | is.rdstate()); 305 is.clear(std::ios::badbit | is.rdstate()); 309 std::cerr << "Mismatch when expecting to r 306 std::cerr << "Mismatch when expecting to read state of a " 310 << name() << " distribution\n" 307 << name() << " distribution\n" 311 << "Name found was " << inName 308 << "Name found was " << inName 312 << "\nistream is left in the badbit st 309 << "\nistream is left in the badbit state\n"; 313 return is; 310 return is; 314 } 311 } 315 if (possibleKeywordInput(is, "Uvec", meanMax 312 if (possibleKeywordInput(is, "Uvec", meanMax)) { 316 std::vector<unsigned long> t(2); 313 std::vector<unsigned long> t(2); 317 is >> meanMax >> t[0] >> t[1]; meanMax 314 is >> meanMax >> t[0] >> t[1]; meanMax = DoubConv::longs2double(t); 318 is >> defaultMean >> t[0] >> t[1]; default 315 is >> defaultMean >> t[0] >> t[1]; defaultMean = DoubConv::longs2double(t); 319 is >> status[0] >> t[0] >> t[1]; status[ 316 is >> status[0] >> t[0] >> t[1]; status[0] = DoubConv::longs2double(t); 320 is >> status[1] >> t[0] >> t[1]; status[ 317 is >> status[1] >> t[0] >> t[1]; status[1] = DoubConv::longs2double(t); 321 is >> status[2] >> t[0] >> t[1]; status[ 318 is >> status[2] >> t[0] >> t[1]; status[2] = DoubConv::longs2double(t); 322 is >> oldm >> t[0] >> t[1]; oldm 319 is >> oldm >> t[0] >> t[1]; oldm = DoubConv::longs2double(t); 323 return is; 320 return is; 324 } 321 } 325 // is >> meanMax encompassed by possibleKeyw 322 // is >> meanMax encompassed by possibleKeywordInput 326 is >> defaultMean >> status[0] >> status[1] 323 is >> defaultMean >> status[0] >> status[1] >> status[2]; 327 return is; 324 return is; 328 } 325 } 329 326 330 } // namespace CLHEP 327 } // namespace CLHEP 331 328 332 329