Geant4 Cross Reference |
>> 1 // $Id:$ 1 // -*- C++ -*- 2 // -*- C++ -*- 2 // 3 // 3 // ------------------------------------------- 4 // ----------------------------------------------------------------------- 4 // HEP Random 5 // HEP Random 5 // --- RandPoissonQ -- 6 // --- RandPoissonQ --- 6 // class implementation f 7 // class implementation file 7 // ------------------------------------------- 8 // ----------------------------------------------------------------------- 8 9 9 // =========================================== 10 // ======================================================================= 10 // M. Fischler - Implemented new, much fast 11 // M. Fischler - Implemented new, much faster table-driven algorithm 11 // applicable for mu < 100 12 // applicable for mu < 100 12 // - Implemented "quick()" methods, shich 13 // - Implemented "quick()" methods, shich are the same as the 13 // new methods for mu < 100 and are a s 14 // new methods for mu < 100 and are a skew-corrected gaussian 14 // approximation for large mu. 15 // approximation for large mu. 15 // M. Fischler - Removed mean=100 from the 16 // M. Fischler - Removed mean=100 from the table-driven set, since it 16 // uses a value just off the end of the 17 // uses a value just off the end of the table. (April 2004) 17 // M Fischler - put and get to/from stream 18 // M Fischler - put and get to/from streams 12/15/04 18 // M Fischler - Utilize RandGaussQ rather 19 // M Fischler - Utilize RandGaussQ rather than RandGauss, as clearly 19 // intended by the inclusion of RandGau 20 // intended by the inclusion of RandGaussQ.h. Using RandGauss 20 // introduces a subtle trap in that the 21 // introduces a subtle trap in that the state of RandPoissonQ 21 // can never be properly captured witho 22 // can never be properly captured without also saveing the 22 // state of RandGauss! RandGaussQ is, 23 // state of RandGauss! RandGaussQ is, on the other hand, 23 // stateless except for the engine used 24 // stateless except for the engine used. 24 // M Fischler - Modified use of wrong engine << 25 // M Fisculer - Modified use of wrong engine when shoot (anEngine, mean) 25 // is called. This flaw was preventing 26 // is called. This flaw was preventing any hope of proper 26 // saving and restoring in the instance 27 // saving and restoring in the instance cases. 27 // M Fischler - fireArray using defaultMea 28 // M Fischler - fireArray using defaultMean 2/10/05 28 // M Fischler - put/get to/from streams 29 // M Fischler - put/get to/from streams uses pairs of ulongs when 29 // + storing doubles avoid problems with 30 // + storing doubles avoid problems with precision 30 // 4/14/05 31 // 4/14/05 31 // M Fischler - Modified use of shoot (mean) << 32 // M Fisculer - Modified use of shoot (mean) instead of 32 // shoot(getLocalEngine(), mean) when f 33 // shoot(getLocalEngine(), mean) when fire(mean) is called. 33 // This flaw was causing bad "cross-tal 34 // This flaw was causing bad "cross-talk" between modules 34 // in CMS, where one used its own engin 35 // in CMS, where one used its own engine, and the other 35 // used the static generator. 10/18/07 36 // used the static generator. 10/18/07 36 // 37 // 37 // =========================================== 38 // ======================================================================= 38 39 39 #include "CLHEP/Random/RandPoissonQ.h" 40 #include "CLHEP/Random/RandPoissonQ.h" 40 #include "CLHEP/Random/RandGaussQ.h" 41 #include "CLHEP/Random/RandGaussQ.h" 41 #include "CLHEP/Random/DoubConv.h" 42 #include "CLHEP/Random/DoubConv.h" 42 #include "CLHEP/Random/Stat.h" 43 #include "CLHEP/Random/Stat.h" 43 #include "CLHEP/Utility/thread_local.h" 44 #include "CLHEP/Utility/thread_local.h" 44 #include <cmath> // for std::pow() 45 #include <cmath> // for std::pow() 45 #include <iostream> << 46 #include <string> << 47 #include <vector> << 48 46 49 namespace CLHEP { 47 namespace CLHEP { 50 48 51 std::string RandPoissonQ::name() const {return 49 std::string RandPoissonQ::name() const {return "RandPoissonQ";} 52 HepRandomEngine & RandPoissonQ::engine() {retu 50 HepRandomEngine & RandPoissonQ::engine() {return RandPoisson::engine();} 53 51 54 // Initialization of static data: Note that t 52 // Initialization of static data: Note that this is all const static data, 55 // so that saveEngineStatus properly saves all 53 // so that saveEngineStatus properly saves all needed information. 56 54 57 // The following MUST MATCH the correspondin 55 // The following MUST MATCH the corresponding values used (in 58 // poissonTables.cc) when poissonTables.cdat 56 // poissonTables.cc) when poissonTables.cdat was created. 59 57 60 const double RandPoissonQ::FIRST_MU = 10;// lo 58 const double RandPoissonQ::FIRST_MU = 10;// lowest mu value in table 61 const double RandPoissonQ::LAST_MU = 95;// hi 59 const double RandPoissonQ::LAST_MU = 95;// highest mu value 62 const double RandPoissonQ::S = 5; // Sp 60 const double RandPoissonQ::S = 5; // Spacing between mu values 63 const int RandPoissonQ::BELOW = 30; // St 61 const int RandPoissonQ::BELOW = 30; // Starting point for N is at mu - BELOW 64 const int RandPoissonQ::ENTRIES = 51; // Nu 62 const int RandPoissonQ::ENTRIES = 51; // Number of entries in each mu row 65 63 66 const double RandPoissonQ::MAXIMUM_POISSON_DEV 64 const double RandPoissonQ::MAXIMUM_POISSON_DEVIATE = 2.0E9; 67 // Careful -- this is NOT the maximum number 65 // Careful -- this is NOT the maximum number that can be held in 68 // a long. It actually should be some large 66 // a long. It actually should be some large number of sigma below 69 // that. 67 // that. 70 68 71 // Here comes the big (9K bytes) table, kept 69 // Here comes the big (9K bytes) table, kept in a file of 72 // ENTRIES * (FIRST_MU - LAST_MU + 1)/S doub 70 // ENTRIES * (FIRST_MU - LAST_MU + 1)/S doubles 73 71 74 static const double poissonTables [ 51 * ( (95 72 static const double poissonTables [ 51 * ( (95-10)/5 + 1 ) ] = { 75 #include "CLHEP/Random/poissonTables.cdat" 73 #include "CLHEP/Random/poissonTables.cdat" 76 }; 74 }; 77 75 78 // 76 // 79 // Constructors and destructors: 77 // Constructors and destructors: 80 // 78 // 81 79 82 RandPoissonQ::~RandPoissonQ() { 80 RandPoissonQ::~RandPoissonQ() { 83 } 81 } 84 82 85 void RandPoissonQ::setupForDefaultMu() { 83 void RandPoissonQ::setupForDefaultMu() { 86 84 87 // The following are useful for quick approx 85 // The following are useful for quick approximation, for large mu 88 86 89 double sig2 = defaultMean * (.9998654 - .083 87 double sig2 = defaultMean * (.9998654 - .08346/defaultMean); 90 sigma = std::sqrt(sig2); 88 sigma = std::sqrt(sig2); 91 // sigma for the Guassian which approximates 89 // sigma for the Guassian which approximates the Poisson -- naively 92 // sqrt (defaultMean). 90 // sqrt (defaultMean). 93 // 91 // 94 // The multiplier corrects for fact that dis 92 // The multiplier corrects for fact that discretization of the form 95 // [gaussian+.5] increases the second moment 93 // [gaussian+.5] increases the second moment by a small amount. 96 94 97 double t = 1./(sig2); 95 double t = 1./(sig2); 98 96 99 a2 = t/6 + t*t/324; 97 a2 = t/6 + t*t/324; 100 a1 = std::sqrt (1-2*a2*a2*sig2); 98 a1 = std::sqrt (1-2*a2*a2*sig2); 101 a0 = defaultMean + .5 - sig2 * a2; 99 a0 = defaultMean + .5 - sig2 * a2; 102 100 103 // The formula will be a0 + a1*x + a2*x*x wh 101 // The formula will be a0 + a1*x + a2*x*x where x has 2nd moment of sigma. 104 // The coeffeicients are chosen to match the 102 // The coeffeicients are chosen to match the first THREE moments of the 105 // true Poisson distribution. 103 // true Poisson distribution. 106 // 104 // 107 // Actually, if the correction for discretiz 105 // Actually, if the correction for discretization were not needed, then 108 // a2 could be taken one order higher by add 106 // a2 could be taken one order higher by adding t*t*t/5832. However, 109 // the discretization correction is not perf 107 // the discretization correction is not perfect, leading to inaccuracy 110 // on the order to 1/mu**2, so adding a thir 108 // on the order to 1/mu**2, so adding a third term is overkill. 111 109 112 } // setupForDefaultMu() 110 } // setupForDefaultMu() 113 111 114 112 115 // 113 // 116 // fire, quick, operator(), and shoot methods: 114 // fire, quick, operator(), and shoot methods: 117 // 115 // 118 116 119 long RandPoissonQ::shoot(double xm) { 117 long RandPoissonQ::shoot(double xm) { 120 return shoot(getTheEngine(), xm); 118 return shoot(getTheEngine(), xm); 121 } 119 } 122 120 123 double RandPoissonQ::operator()() { 121 double RandPoissonQ::operator()() { 124 return (double) fire(); 122 return (double) fire(); 125 } 123 } 126 124 127 double RandPoissonQ::operator()( double mean ) 125 double RandPoissonQ::operator()( double mean ) { 128 return (double) fire(mean); 126 return (double) fire(mean); 129 } 127 } 130 128 131 long RandPoissonQ::fire(double mean) { 129 long RandPoissonQ::fire(double mean) { 132 return shoot(getLocalEngine(), mean); 130 return shoot(getLocalEngine(), mean); 133 } 131 } 134 132 135 long RandPoissonQ::fire() { 133 long RandPoissonQ::fire() { 136 if ( defaultMean < LAST_MU + S ) { 134 if ( defaultMean < LAST_MU + S ) { 137 return poissonDeviateSmall ( getLocalEngin 135 return poissonDeviateSmall ( getLocalEngine(), defaultMean ); 138 } else { 136 } else { 139 return poissonDeviateQuick ( getLocalEngin 137 return poissonDeviateQuick ( getLocalEngine(), a0, a1, a2, sigma ); 140 } 138 } 141 } // fire() 139 } // fire() 142 140 143 long RandPoissonQ::shoot(HepRandomEngine* anEn 141 long RandPoissonQ::shoot(HepRandomEngine* anEngine, double mean) { 144 142 145 // The following variables, static to this m 143 // The following variables, static to this method, apply to the 146 // last time a large mean was supplied; they 144 // last time a large mean was supplied; they obviate certain calculations 147 // if consecutive calls use the same mean. 145 // if consecutive calls use the same mean. 148 146 149 static CLHEP_THREAD_LOCAL double lastLargeMe 147 static CLHEP_THREAD_LOCAL double lastLargeMean = -1.; // Mean from previous shoot 150 // requiring poissonDeviateQuick() 148 // requiring poissonDeviateQuick() 151 static CLHEP_THREAD_LOCAL double lastA0; 149 static CLHEP_THREAD_LOCAL double lastA0; 152 static CLHEP_THREAD_LOCAL double lastA1; 150 static CLHEP_THREAD_LOCAL double lastA1; 153 static CLHEP_THREAD_LOCAL double lastA2; 151 static CLHEP_THREAD_LOCAL double lastA2; 154 static CLHEP_THREAD_LOCAL double lastSigma; 152 static CLHEP_THREAD_LOCAL double lastSigma; 155 153 156 if ( mean < LAST_MU + S ) { 154 if ( mean < LAST_MU + S ) { 157 return poissonDeviateSmall ( anEngine, mea 155 return poissonDeviateSmall ( anEngine, mean ); 158 } else { 156 } else { 159 if ( mean != lastLargeMean ) { 157 if ( mean != lastLargeMean ) { 160 // Compute the coefficients defining the 158 // Compute the coefficients defining the quadratic transformation from a 161 // Gaussian to a Poisson for this mean. 159 // Gaussian to a Poisson for this mean. Also save these for next time. 162 double sig2 = mean * (.9998654 - .08346/ 160 double sig2 = mean * (.9998654 - .08346/mean); 163 lastSigma = std::sqrt(sig2); 161 lastSigma = std::sqrt(sig2); 164 double t = 1./sig2; 162 double t = 1./sig2; 165 lastA2 = t*(1./6.) + t*t*(1./324.); 163 lastA2 = t*(1./6.) + t*t*(1./324.); 166 lastA1 = std::sqrt (1-2*lastA2*lastA2*si 164 lastA1 = std::sqrt (1-2*lastA2*lastA2*sig2); 167 lastA0 = mean + .5 - sig2 * lastA2; 165 lastA0 = mean + .5 - sig2 * lastA2; 168 } 166 } 169 return poissonDeviateQuick ( anEngine, las 167 return poissonDeviateQuick ( anEngine, lastA0, lastA1, lastA2, lastSigma ); 170 } 168 } 171 169 172 } // shoot (anEngine, mean) 170 } // shoot (anEngine, mean) 173 171 174 void RandPoissonQ::shootArray(const int size, 172 void RandPoissonQ::shootArray(const int size, long* vect, double m) { 175 for( long* v = vect; v != vect + size; ++v ) 173 for( long* v = vect; v != vect + size; ++v ) 176 *v = shoot(m); 174 *v = shoot(m); 177 // Note: We could test for m > 100, and i 175 // Note: We could test for m > 100, and if it is, precompute a0, a1, a2, 178 // and sigma and call the appropriate for 176 // and sigma and call the appropriate form of poissonDeviateQuick. 179 // But since those are cached anyway, not 177 // But since those are cached anyway, not much time would be saved. 180 } 178 } 181 179 182 void RandPoissonQ::shootArray(HepRandomEngine* << 183 long* vect, doubl << 184 for( long* v = vect; v != vect + size; ++v ) << 185 *v = shoot(anEngine,m1); << 186 } << 187 << 188 void RandPoissonQ::fireArray(const int size, l 180 void RandPoissonQ::fireArray(const int size, long* vect, double m) { 189 for( long* v = vect; v != vect + size; ++v ) 181 for( long* v = vect; v != vect + size; ++v ) 190 *v = fire( m ); 182 *v = fire( m ); 191 } 183 } 192 184 193 void RandPoissonQ::fireArray(const int size, l 185 void RandPoissonQ::fireArray(const int size, long* vect) { 194 for( long* v = vect; v != vect + size; ++v ) 186 for( long* v = vect; v != vect + size; ++v ) 195 *v = fire( defaultMean ); 187 *v = fire( defaultMean ); 196 } 188 } 197 189 198 190 199 // Quick Poisson deviate algorithm used by qui 191 // Quick Poisson deviate algorithm used by quick for large mu: 200 192 201 long RandPoissonQ::poissonDeviateQuick ( HepRa 193 long RandPoissonQ::poissonDeviateQuick ( HepRandomEngine *e, double mu ) { 202 194 203 // Compute the coefficients defining the qua 195 // Compute the coefficients defining the quadratic transformation from a 204 // Gaussian to a Poisson: 196 // Gaussian to a Poisson: 205 197 206 double sig2 = mu * (.9998654 - .08346/mu); 198 double sig2 = mu * (.9998654 - .08346/mu); 207 double sig = std::sqrt(sig2); 199 double sig = std::sqrt(sig2); 208 // The multiplier corrects for fact that dis 200 // The multiplier corrects for fact that discretization of the form 209 // [gaussian+.5] increases the second moment 201 // [gaussian+.5] increases the second moment by a small amount. 210 202 211 double t = 1./sig2; 203 double t = 1./sig2; 212 204 213 double sa2 = t*(1./6.) + t*t*(1./324.); 205 double sa2 = t*(1./6.) + t*t*(1./324.); 214 double sa1 = std::sqrt (1-2*sa2*sa2*sig2); 206 double sa1 = std::sqrt (1-2*sa2*sa2*sig2); 215 double sa0 = mu + .5 - sig2 * sa2; 207 double sa0 = mu + .5 - sig2 * sa2; 216 208 217 // The formula will be sa0 + sa1*x + sa2*x*x 209 // The formula will be sa0 + sa1*x + sa2*x*x where x has sigma of sq. 218 // The coeffeicients are chosen to match the 210 // The coeffeicients are chosen to match the first THREE moments of the 219 // true Poisson distribution. 211 // true Poisson distribution. 220 212 221 return poissonDeviateQuick ( e, sa0, sa1, sa 213 return poissonDeviateQuick ( e, sa0, sa1, sa2, sig ); 222 } 214 } 223 215 224 216 225 long RandPoissonQ::poissonDeviateQuick ( HepRa 217 long RandPoissonQ::poissonDeviateQuick ( HepRandomEngine *e, 226 double A0, double A1, double A2, double si 218 double A0, double A1, double A2, double sig) { 227 // 219 // 228 // Quick Poisson deviate algorithm used by qui 220 // Quick Poisson deviate algorithm used by quick for large mu: 229 // 221 // 230 // The principle: For very large mu, a poisso 222 // The principle: For very large mu, a poisson distribution can be approximated 231 // by a gaussian: return the integer part of m 223 // by a gaussian: return the integer part of mu + .5 + g where g is a unit 232 // normal. However, this yelds a miserable ap 224 // normal. However, this yelds a miserable approximation at values as 233 // "large" as 100. The primary problem is tha 225 // "large" as 100. The primary problem is that the poisson distribution is 234 // supposed to have a skew of 1/mu**2, and the 226 // supposed to have a skew of 1/mu**2, and the zero skew of the Guassian 235 // leads to errors of order as big as 1/mu**2. 227 // leads to errors of order as big as 1/mu**2. 236 // 228 // 237 // We substitute for the gaussian a quadratic 229 // We substitute for the gaussian a quadratic function of that gaussian random. 238 // The expression looks very nearly like mu + 230 // The expression looks very nearly like mu + .5 - 1/6 + g + g**2/(6*mu). 239 // The small positive quadratic term causes th 231 // The small positive quadratic term causes the resulting variate to have 240 // a positive skew; the -1/6 constant term is 232 // a positive skew; the -1/6 constant term is there to correct for this bias 241 // in the mean. By adjusting these two and th 233 // in the mean. By adjusting these two and the linear term, we can match the 242 // first three moments to high accuracy in 1/m 234 // first three moments to high accuracy in 1/mu. 243 // 235 // 244 // The sigma used is not precisely sqrt(mu) si 236 // The sigma used is not precisely sqrt(mu) since a rounded-off Gaussian 245 // has a second moment which is slightly large 237 // has a second moment which is slightly larger than that of the Gaussian. 246 // To compensate, sig is multiplied by a facto 238 // To compensate, sig is multiplied by a factor which is slightly less than 1. 247 239 248 // double g = RandGauss::shootQuick( e ); 240 // double g = RandGauss::shootQuick( e ); // TEMPORARY MOD: 249 double g = RandGaussQ::shoot( e ); // Unit 241 double g = RandGaussQ::shoot( e ); // Unit normal 250 g *= sig; 242 g *= sig; 251 double p = A2*g*g + A1*g + A0; 243 double p = A2*g*g + A1*g + A0; 252 if ( p < 0 ) return 0; // Shouldn't ever po 244 if ( p < 0 ) return 0; // Shouldn't ever possibly happen since 253 // mean should not be less than 100, b 245 // mean should not be less than 100, but 254 // we check due to paranoia. 246 // we check due to paranoia. 255 if ( p > MAXIMUM_POISSON_DEVIATE ) p = MAXIM 247 if ( p > MAXIMUM_POISSON_DEVIATE ) p = MAXIMUM_POISSON_DEVIATE; 256 248 257 return long(p); 249 return long(p); 258 250 259 } // poissonDeviateQuick () 251 } // poissonDeviateQuick () 260 252 261 253 262 254 263 long RandPoissonQ::poissonDeviateSmall (HepRan 255 long RandPoissonQ::poissonDeviateSmall (HepRandomEngine * e, double mean) { 264 long N1; 256 long N1; 265 long N2; 257 long N2; 266 // The following are for later use to form a 258 // The following are for later use to form a secondary random s: 267 double rRange; // This will hold the in 259 double rRange; // This will hold the interval between cdf for the 268 // computed N1 and cdf for N1+1. 260 // computed N1 and cdf for N1+1. 269 double rRemainder = 0; // This will hold the 261 double rRemainder = 0; // This will hold the length into that interval. 270 262 271 // Coming in, mean should not be more than L 263 // Coming in, mean should not be more than LAST_MU + S. However, we will 272 // be paranoid and test for this: 264 // be paranoid and test for this: 273 265 274 if ( mean > LAST_MU + S ) { 266 if ( mean > LAST_MU + S ) { 275 return RandPoisson::shoot(e, mean); 267 return RandPoisson::shoot(e, mean); 276 } 268 } 277 269 278 if (mean <= 0) { 270 if (mean <= 0) { 279 return 0; // Perhaps we ought to balk 271 return 0; // Perhaps we ought to balk harder here! 280 } 272 } 281 273 282 // >>> 1 <<< 274 // >>> 1 <<< 283 // Generate the first random, which we alwa 275 // Generate the first random, which we always will need. 284 276 285 double r = e->flat(); 277 double r = e->flat(); 286 278 287 // >>> 2 <<< 279 // >>> 2 <<< 288 // For small mean, below the start of the t 280 // For small mean, below the start of the tables, 289 // do the series for cdf directly. 281 // do the series for cdf directly. 290 282 291 // In this case, since we know the series wi 283 // In this case, since we know the series will terminate relatively quickly, 292 // almost alwaye use a precomputed 1/N array 284 // almost alwaye use a precomputed 1/N array without fear of overrunning it. 293 285 294 static const double oneOverN[50] = 286 static const double oneOverN[50] = 295 { 0, 1., 1/2., 1/3., 1/4., 1/5., 287 { 0, 1., 1/2., 1/3., 1/4., 1/5., 1/6., 1/7., 1/8., 1/9., 296 1/10., 1/11., 1/12., 1/13., 1/14., 1/15., 288 1/10., 1/11., 1/12., 1/13., 1/14., 1/15., 1/16., 1/17., 1/18., 1/19., 297 1/20., 1/21., 1/22., 1/23., 1/24., 1/25., 289 1/20., 1/21., 1/22., 1/23., 1/24., 1/25., 1/26., 1/27., 1/28., 1/29., 298 1/30., 1/31., 1/32., 1/33., 1/34., 1/35., 290 1/30., 1/31., 1/32., 1/33., 1/34., 1/35., 1/36., 1/37., 1/38., 1/39., 299 1/40., 1/41., 1/42., 1/43., 1/44., 1/45., 291 1/40., 1/41., 1/42., 1/43., 1/44., 1/45., 1/46., 1/47., 1/48., 1/49. }; 300 292 301 293 302 if ( mean < FIRST_MU ) { 294 if ( mean < FIRST_MU ) { 303 295 304 long N = 0; 296 long N = 0; 305 double term = std::exp(-mean); 297 double term = std::exp(-mean); 306 double cdf = term; 298 double cdf = term; 307 299 308 if ( r < (1 - 1.0E-9) ) { 300 if ( r < (1 - 1.0E-9) ) { 309 // 301 // 310 // **** This is a normal path: **** 302 // **** This is a normal path: **** 311 // 303 // 312 // Except when r is very close to 1, it 304 // Except when r is very close to 1, it is certain that we will exceed r 313 // before the 30-th term in the series, 305 // before the 30-th term in the series, so a simple while loop is OK. 314 const double* oneOverNptr = oneOverN; 306 const double* oneOverNptr = oneOverN; 315 while( cdf <= r ) { 307 while( cdf <= r ) { 316 N++ ; 308 N++ ; 317 oneOverNptr++; 309 oneOverNptr++; 318 term *= ( mean * (*oneOverNptr) ); 310 term *= ( mean * (*oneOverNptr) ); 319 cdf += term; 311 cdf += term; 320 } 312 } 321 return N; 313 return N; 322 // 314 // 323 // **** **** 315 // **** **** 324 // 316 // 325 } else { // r is almost 1... 317 } else { // r is almost 1... 326 // For r very near to 1 we would have to 318 // For r very near to 1 we would have to check that we don't fall 327 // off the end of the table of 1/N. Sin 319 // off the end of the table of 1/N. Since this is very rare, we just 328 // ignore the table and do the identical 320 // ignore the table and do the identical while loop, using explicit 329 // division. 321 // division. 330 double cdf0; 322 double cdf0; 331 while ( cdf <= r ) { 323 while ( cdf <= r ) { 332 N++ ; 324 N++ ; 333 term *= ( mean / N ); 325 term *= ( mean / N ); 334 cdf0 = cdf; 326 cdf0 = cdf; 335 cdf += term; 327 cdf += term; 336 if (cdf == cdf0) break; // Can't happen, but 328 if (cdf == cdf0) break; // Can't happen, but just in case... 337 } 329 } 338 return N; 330 return N; 339 } // end of if ( r compared to (1 - 1.0E-9 331 } // end of if ( r compared to (1 - 1.0E-9) ) 340 332 341 } // End of the code for mean < FIRST_MU 333 } // End of the code for mean < FIRST_MU 342 334 343 // >>> 3 <<< 335 // >>> 3 <<< 344 // Find the row of the tables corresponding 336 // Find the row of the tables corresponding to the highest tabulated mu 345 // which is no greater than our actual mean 337 // which is no greater than our actual mean. 346 338 347 int rowNumber = int((mean - FIRST_MU)/S); 339 int rowNumber = int((mean - FIRST_MU)/S); 348 const double * cdfs = &poissonTables [rowNum 340 const double * cdfs = &poissonTables [rowNumber*ENTRIES]; 349 double mu = FIRST_MU + rowNumber*S; 341 double mu = FIRST_MU + rowNumber*S; 350 double deltaMu = mean - mu; 342 double deltaMu = mean - mu; 351 int Nmin = int(mu - BELOW); 343 int Nmin = int(mu - BELOW); 352 if (Nmin < 1) Nmin = 1; 344 if (Nmin < 1) Nmin = 1; 353 int Nmax = Nmin + (ENTRIES - 1); 345 int Nmax = Nmin + (ENTRIES - 1); 354 346 355 347 356 // >>> 4 <<< 348 // >>> 4 <<< 357 // If r is less that the smallest entry in 349 // If r is less that the smallest entry in the row, then 358 // generate the deviate directly from the s 350 // generate the deviate directly from the series. 359 351 360 if ( r < cdfs[0] ) { 352 if ( r < cdfs[0] ) { 361 353 362 // In this case, we are tempted to use the 354 // In this case, we are tempted to use the actual mean, and not 363 // generate a second deviate to account fo 355 // generate a second deviate to account for the leftover part mean - mu. 364 // That would be an error, generating a di 356 // That would be an error, generating a distribution with enough excess 365 // at Nmin + (mean-mu)/2 to be detectable 357 // at Nmin + (mean-mu)/2 to be detectable in 4,000,000 trials. 366 358 367 // Since this case is very rare (never mor 359 // Since this case is very rare (never more than .2% of the r values) 368 // and can happen where N will be large (u 360 // and can happen where N will be large (up to 65 for the mu=95 row) 369 // we use explicit division so as to avoid 361 // we use explicit division so as to avoid having to worry about running 370 // out of oneOverN table. 362 // out of oneOverN table. 371 363 372 long N = 0; 364 long N = 0; 373 double term = std::exp(-mu); 365 double term = std::exp(-mu); 374 double cdf = term; 366 double cdf = term; 375 double cdf0; 367 double cdf0; 376 368 377 while(cdf <= r) { 369 while(cdf <= r) { 378 N++ ; 370 N++ ; 379 term *= ( mu / N ); 371 term *= ( mu / N ); 380 cdf0 = cdf; 372 cdf0 = cdf; 381 cdf += term; 373 cdf += term; 382 if (cdf == cdf0) break; // Can't happen, but 374 if (cdf == cdf0) break; // Can't happen, but just in case... 383 } 375 } 384 N1 = N; 376 N1 = N; 385 // std::cout << r << " " << N << " "; 377 // std::cout << r << " " << N << " "; 386 // DBG_small = true; 378 // DBG_small = true; 387 rRange = 0; // In this case there is alw 379 rRange = 0; // In this case there is always a second r needed 388 380 389 } // end of small-r case 381 } // end of small-r case 390 382 391 383 392 // >>> 5 <<< 384 // >>> 5 <<< 393 // Assuming r lies within the scope of the 385 // Assuming r lies within the scope of the row for this mu, find the 394 // largest entry not greater than r. N1 is 386 // largest entry not greater than r. N1 is the N corresponding to the 395 // index a. 387 // index a. 396 388 397 else if ( r < cdfs[ENTRIES-1] ) { // r i 389 else if ( r < cdfs[ENTRIES-1] ) { // r is also >= cdfs[0] 398 390 399 // 391 // 400 // **** This is the normal code path **** 392 // **** This is the normal code path **** 401 // 393 // 402 394 403 int a = 0; // Highest val 395 int a = 0; // Highest value of index such that cdfs[a] 404 // is known NOT to be greater than r. 396 // is known NOT to be greater than r. 405 int b = ENTRIES - 1; // Lowest value of i 397 int b = ENTRIES - 1; // Lowest value of index such that cdfs[b] is 406 // known to exeed r. 398 // known to exeed r. 407 399 408 while (b != (a+1) ) { 400 while (b != (a+1) ) { 409 int c = (a+b+1)>>1; 401 int c = (a+b+1)>>1; 410 if (r > cdfs[c]) { 402 if (r > cdfs[c]) { 411 a = c; 403 a = c; 412 } else { 404 } else { 413 b = c; 405 b = c; 414 } 406 } 415 } 407 } 416 408 417 N1 = Nmin + a; 409 N1 = Nmin + a; 418 rRange = cdfs[a+1] - cdfs[a]; 410 rRange = cdfs[a+1] - cdfs[a]; 419 rRemainder = r - cdfs[a]; 411 rRemainder = r - cdfs[a]; 420 412 421 // 413 // 422 // **** **** 414 // **** **** 423 // 415 // 424 416 425 } // end of medium-r (normal) case 417 } // end of medium-r (normal) case 426 418 427 419 428 // >>> 6 <<< 420 // >>> 6 <<< 429 // If r exceeds the greatest entry in the t 421 // If r exceeds the greatest entry in the table for this mu, then start 430 // from that cdf, and use the series to com 422 // from that cdf, and use the series to compute from there until r is 431 // exceeded. 423 // exceeded. 432 424 433 else { // if ( r >= cdfs[ENTRIES-1] ) { 425 else { // if ( r >= cdfs[ENTRIES-1] ) { 434 426 435 // Here, division must be done explicitly, 427 // Here, division must be done explicitly, and we must also protect against 436 // roundoff preventing termination. 428 // roundoff preventing termination. 437 429 438 // 430 // 439 //+++ cdfs[ENTRIES-1] is exp(-mu) sum (mu**m 431 //+++ cdfs[ENTRIES-1] is exp(-mu) sum (mu**m/m! , m=0 to Nmax) 440 //+++ (where Nmax = mu - BELOW + ENTRIES - 1 432 //+++ (where Nmax = mu - BELOW + ENTRIES - 1) 441 //+++ cdfs[ENTRIES-1]-cdfs[ENTRIES-2] is exp 433 //+++ cdfs[ENTRIES-1]-cdfs[ENTRIES-2] is exp(-mu) mu**(Nmax)/(Nmax)! 442 //+++ If the sum up to k-1 <= r < sum up to 434 //+++ If the sum up to k-1 <= r < sum up to k, then N = k-1 443 //+++ Consider k = Nmax in the above stateme 435 //+++ Consider k = Nmax in the above statement: 444 //+++ If cdfs[ENTRIES-2] <= r < cdfs[ENTRIES 436 //+++ If cdfs[ENTRIES-2] <= r < cdfs[ENTRIES-1], N would be Nmax-1 445 //+++ But here r >= cdfs[ENTRIES-1] so N >= 437 //+++ But here r >= cdfs[ENTRIES-1] so N >= Nmax 446 // 438 // 447 439 448 // Erroneous: 440 // Erroneous: 449 //+++ cdfs[ENTRIES-1] is exp(-mu) sum (mu**m 441 //+++ cdfs[ENTRIES-1] is exp(-mu) sum (mu**m/m! , m=0 to Nmax-1) 450 //+++ cdfs[ENTRIES-1]-cdfs[ENTRIES-2] is exp 442 //+++ cdfs[ENTRIES-1]-cdfs[ENTRIES-2] is exp(-mu) mu**(Nmax-1)/(Nmax-1)! 451 //+++ If a sum up to k-1 <= r < sum up to k, 443 //+++ If a sum up to k-1 <= r < sum up to k, then N = k-1 452 //+++ So if cdfs[ENTRIES-1] were > r, N woul 444 //+++ So if cdfs[ENTRIES-1] were > r, N would be Nmax-1 (or less) 453 //+++ But here r >= cdfs[ENTRIES-1] so N >= 445 //+++ But here r >= cdfs[ENTRIES-1] so N >= Nmax 454 // 446 // 455 447 456 // std::cout << "r = " << r << " mu = " < 448 // std::cout << "r = " << r << " mu = " << mu << "\n"; 457 long N = Nmax -1; 449 long N = Nmax -1; 458 double cdf = cdfs[ENTRIES-1]; 450 double cdf = cdfs[ENTRIES-1]; 459 double term = cdf - cdfs[ENTRIES-2]; 451 double term = cdf - cdfs[ENTRIES-2]; 460 double cdf0; 452 double cdf0; 461 while(cdf <= r) { 453 while(cdf <= r) { 462 N++ ; 454 N++ ; 463 // std::cout << " N " << N << " term " < 455 // std::cout << " N " << N << " term " << 464 // term << " cdf " << cdf << "\n"; 456 // term << " cdf " << cdf << "\n"; 465 term *= ( mu / N ); 457 term *= ( mu / N ); 466 cdf0 = cdf; 458 cdf0 = cdf; 467 cdf += term; 459 cdf += term; 468 if (cdf == cdf0) break; // If term gets so s 460 if (cdf == cdf0) break; // If term gets so small cdf stops increasing, 469 // terminate using that value of N sin 461 // terminate using that value of N since we 470 // would never reach r. 462 // would never reach r. 471 } 463 } 472 N1 = N; 464 N1 = N; 473 rRange = 0; // We can't validly omit the 465 rRange = 0; // We can't validly omit the second true random 474 466 475 // N = Nmax -1; 467 // N = Nmax -1; 476 // cdf = cdfs[ENTRIES-1]; 468 // cdf = cdfs[ENTRIES-1]; 477 // term = cdf - cdfs[ENTRIES-2]; 469 // term = cdf - cdfs[ENTRIES-2]; 478 // for (int isxz=0; isxz < 100; isxz++) 470 // for (int isxz=0; isxz < 100; isxz++) { 479 // N++ ; 471 // N++ ; 480 // term *= ( mu / N ); 472 // term *= ( mu / N ); 481 // cdf0 = cdf; 473 // cdf0 = cdf; 482 // cdf += term; 474 // cdf += term; 483 // } 475 // } 484 // std::cout.precision(20); 476 // std::cout.precision(20); 485 // std::cout << "Final sum is " << cdf 477 // std::cout << "Final sum is " << cdf << "\n"; 486 478 487 } // end of large-r case 479 } // end of large-r case 488 480 489 481 490 482 491 // >>> 7 <<< 483 // >>> 7 <<< 492 // Form a second random, s, based on the po 484 // Form a second random, s, based on the position of r within the range 493 // of this table entry to the next entry. 485 // of this table entry to the next entry. 494 486 495 // However, if this range is very small, the 487 // However, if this range is very small, then we lose too many bits of 496 // randomness. In that situation, we genera 488 // randomness. In that situation, we generate a second random for s. 497 489 498 double s; 490 double s; 499 491 500 static const double MINRANGE = .01; // Sac 492 static const double MINRANGE = .01; // Sacrifice up to two digits of 501 // randomness when using r to 493 // randomness when using r to produce 502 // a second random s. Leads to up t 494 // a second random s. Leads to up to 503 // .09 extra randoms each time. 495 // .09 extra randoms each time. 504 496 505 if ( rRange > MINRANGE ) { 497 if ( rRange > MINRANGE ) { 506 // 498 // 507 // **** This path taken 90% of the time ** 499 // **** This path taken 90% of the time **** 508 // 500 // 509 s = rRemainder / rRange; 501 s = rRemainder / rRange; 510 } else { 502 } else { 511 s = e->flat(); // extra true random neede 503 s = e->flat(); // extra true random needed about one time in 10. 512 } 504 } 513 505 514 // >>> 8 <<< 506 // >>> 8 <<< 515 // Use the direct summation method to form 507 // Use the direct summation method to form a second poisson deviate N2 516 // from deltaMu and s. 508 // from deltaMu and s. 517 509 518 N2 = 0; 510 N2 = 0; 519 double term = std::exp(-deltaMu); 511 double term = std::exp(-deltaMu); 520 double cdf = term; 512 double cdf = term; 521 513 522 if ( s < (1 - 1.0E-10) ) { 514 if ( s < (1 - 1.0E-10) ) { 523 // 515 // 524 // This is the normal path: 516 // This is the normal path: 525 // 517 // 526 const double* oneOverNptr = oneOverN; 518 const double* oneOverNptr = oneOverN; 527 while( cdf <= s ) { 519 while( cdf <= s ) { 528 N2++ ; 520 N2++ ; 529 oneOverNptr++; 521 oneOverNptr++; 530 term *= ( deltaMu * (*oneOverNptr) ); 522 term *= ( deltaMu * (*oneOverNptr) ); 531 cdf += term; 523 cdf += term; 532 } 524 } 533 } else { // s is almost 1... 525 } else { // s is almost 1... 534 while( cdf <= s ) { 526 while( cdf <= s ) { 535 N2++ ; 527 N2++ ; 536 term *= ( deltaMu / N2 ); 528 term *= ( deltaMu / N2 ); 537 cdf += term; 529 cdf += term; 538 } 530 } 539 } // end of if ( s compared to (1 - 1.0E-10) 531 } // end of if ( s compared to (1 - 1.0E-10) ) 540 532 541 // >>> 9 <<< 533 // >>> 9 <<< 542 // The result is the sum of those two devia 534 // The result is the sum of those two deviates 543 535 544 // if (DBG_small) { 536 // if (DBG_small) { 545 // std::cout << N2 << " " << N1+N2 << 537 // std::cout << N2 << " " << N1+N2 << "\n"; 546 // DBG_small = false; 538 // DBG_small = false; 547 // } 539 // } 548 540 549 return N1 + N2; 541 return N1 + N2; 550 542 551 } // poissonDeviate() 543 } // poissonDeviate() 552 544 553 std::ostream & RandPoissonQ::put ( std::ostrea 545 std::ostream & RandPoissonQ::put ( std::ostream & os ) const { 554 long pr=os.precision(20); << 546 int pr=os.precision(20); 555 std::vector<unsigned long> t(2); 547 std::vector<unsigned long> t(2); 556 os << " " << name() << "\n"; 548 os << " " << name() << "\n"; 557 os << "Uvec" << "\n"; 549 os << "Uvec" << "\n"; 558 t = DoubConv::dto2longs(a0); 550 t = DoubConv::dto2longs(a0); 559 os << a0 << " " << t[0] << " " << t[1] << "\ 551 os << a0 << " " << t[0] << " " << t[1] << "\n"; 560 t = DoubConv::dto2longs(a1); 552 t = DoubConv::dto2longs(a1); 561 os << a1 << " " << t[0] << " " << t[1] << "\ 553 os << a1 << " " << t[0] << " " << t[1] << "\n"; 562 t = DoubConv::dto2longs(a2); 554 t = DoubConv::dto2longs(a2); 563 os << a2 << " " << t[0] << " " << t[1] << "\ 555 os << a2 << " " << t[0] << " " << t[1] << "\n"; 564 t = DoubConv::dto2longs(sigma); 556 t = DoubConv::dto2longs(sigma); 565 os << sigma << " " << t[0] << " " << t[1] << 557 os << sigma << " " << t[0] << " " << t[1] << "\n"; 566 RandPoisson::put(os); 558 RandPoisson::put(os); 567 os.precision(pr); 559 os.precision(pr); 568 return os; 560 return os; 569 #ifdef REMOVED 561 #ifdef REMOVED 570 int pr=os.precision(20); 562 int pr=os.precision(20); 571 os << " " << name() << "\n"; 563 os << " " << name() << "\n"; 572 os << a0 << " " << a1 << " " << a2 << "\n"; 564 os << a0 << " " << a1 << " " << a2 << "\n"; 573 os << sigma << "\n"; 565 os << sigma << "\n"; 574 RandPoisson::put(os); 566 RandPoisson::put(os); 575 os.precision(pr); 567 os.precision(pr); 576 return os; 568 return os; 577 #endif 569 #endif 578 } 570 } 579 571 580 std::istream & RandPoissonQ::get ( std::istrea 572 std::istream & RandPoissonQ::get ( std::istream & is ) { 581 std::string inName; 573 std::string inName; 582 is >> inName; 574 is >> inName; 583 if (inName != name()) { 575 if (inName != name()) { 584 is.clear(std::ios::badbit | is.rdstate()); 576 is.clear(std::ios::badbit | is.rdstate()); 585 std::cerr << "Mismatch when expecting to r 577 std::cerr << "Mismatch when expecting to read state of a " 586 << name() << " distribution\n" 578 << name() << " distribution\n" 587 << "Name found was " << inName 579 << "Name found was " << inName 588 << "\nistream is left in the badbit st 580 << "\nistream is left in the badbit state\n"; 589 return is; 581 return is; 590 } 582 } 591 if (possibleKeywordInput(is, "Uvec", a0)) { 583 if (possibleKeywordInput(is, "Uvec", a0)) { 592 std::vector<unsigned long> t(2); 584 std::vector<unsigned long> t(2); 593 is >> a0 >> t[0] >> t[1]; a0 = DoubConv::l 585 is >> a0 >> t[0] >> t[1]; a0 = DoubConv::longs2double(t); 594 is >> a1 >> t[0] >> t[1]; a1 = DoubConv::l 586 is >> a1 >> t[0] >> t[1]; a1 = DoubConv::longs2double(t); 595 is >> a2 >> t[0] >> t[1]; a2 = DoubConv::l 587 is >> a2 >> t[0] >> t[1]; a2 = DoubConv::longs2double(t); 596 is >> sigma >> t[0] >> t[1]; sigma = DoubC 588 is >> sigma >> t[0] >> t[1]; sigma = DoubConv::longs2double(t); 597 RandPoisson::get(is); 589 RandPoisson::get(is); 598 return is; 590 return is; 599 } 591 } 600 // is >> a0 encompassed by possibleKeywordIn 592 // is >> a0 encompassed by possibleKeywordInput 601 is >> a1 >> a2 >> sigma; 593 is >> a1 >> a2 >> sigma; 602 RandPoisson::get(is); 594 RandPoisson::get(is); 603 return is; 595 return is; 604 } 596 } 605 597 606 } // namespace CLHEP 598 } // namespace CLHEP 607 599 608 600