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