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