Garfield++ v1r0
A toolkit for the detailed simulation of particle detectors based on ionisation measurement in gases and semiconductors
Loading...
Searching...
No Matches
pois.cpp
Go to the documentation of this file.
1#include <cmath>
2
7
8/*
9// to avoid one else file put it here
10#ifdef PRINT_RANLUX
11unsigned long num_ranlux = 0;
12#endif
13
14PoisState pois_state;
15*/
16
17long pois(double AMU, int &IERROR)
18 //C
19 //C POISSON GENERATOR
20 //C CODED FROM LOS ALAMOS REPORT LA-5061-MS
21 //C PROB(N)=EXP(-AMU)*AMU**N/FACT(N)
22 //C WHERE FACT(N) STANDS FOR FACTORIAL OF N
23 //C ON RETURN IERROR.EQ.0 NORMALLY
24 //C IERROR.EQ.1 IF AMU.LE.0.
25 //C
26 {
27 double AMUOL = -1.;
28 double AMAX = 100.;
29 double EXPMA = 0.;
30 double PIR = 0;
31 double RAN;
32 //static double SECOND_RAN;
33 //static int s_inited_SECOND_RAN = 0;
34 long N = 0;
35 IERROR = 0;
36 if (AMU > AMAX) goto m500;
37 if (AMU == AMUOL) goto m200;
38 if (AMU > 0.0) goto m100;
39 //C MEAN SHOULD BE POSITIVE
40 IERROR = 1;
41 return 0;
42//C SAVE EXPONENTIAL FOR FURTHER IDENTICAL REQUESTS
43m100:
44 IERROR = 0;
45 AMUOL = AMU;
46 EXPMA = exp(-AMU);
47m200:
48 PIR = 1.;
49 N = -1;
50m300:
51 N = N + 1;
52 //{ // for debug
53 // double x = SRANLUX();
54 // mcout<<"pois: x = "<<x<<'\n';
55 // PIR=PIR * x;
56 //}
57 PIR = PIR * SRANLUX(); // working variant
58 if (PIR > EXPMA) goto m300;
59 return N;
60//C NORMAL APPROXIMATION FOR AMU.GT.AMAX
61m500:
62 /*
63 //Iprint2n(mcout,
64 // pois_state.s_inited_SECOND_RAN,
65 // pois_state.SECOND_RAN);
66 if(pois_state.s_inited_SECOND_RAN == 1)
67 {
68 pois_state.s_inited_SECOND_RAN = 0;
69 N = long(pois_state.SECOND_RAN*sqrt(AMU) + AMU + .5);
70 return N;
71 }
72 else
73 {
74 rnorm_double(SRANLUX(), SRANLUX(), RAN, pois_state.SECOND_RAN);
75 N = long(RAN*sqrt(AMU) + AMU + .5);
76 pois_state.s_inited_SECOND_RAN = 1;
77 return N;
78 }
79 */
80 RAN = rnorm_improved();
81 N = long(RAN * sqrt(AMU) + AMU + .5);
82 return N;
83
84}
85
86/*
87long pois (float AMU,int &IERROR)
88//C
89//C POISSON GENERATOR
90//C CODED FROM LOS ALAMOS REPORT LA-5061-MS
91//C PROB(N)=EXP(-AMU)*AMU**N/FACT(N)
92//C WHERE FACT(N) STANDS FOR FACTORIAL OF N
93//C ON RETURN IERROR.EQ.0 NORMALLY
94//C IERROR.EQ.1 IF AMU.LE.0.
95//C
96{
97 float AMUOL=-1.;
98 float AMAX=100.;
99 float EXPMA;
100 float PIR=0;
101 float RAN,DUMMY;
102 long N=0;
103 IERROR=0;
104 if(AMU>AMAX) goto m500;
105 if(AMU==AMUOL) goto m200;
106 if(AMU>0.) goto m100 ;
107//C MEAN SHOULD BE POSITIVE
108 IERROR=1;
109 return 0;
110//C SAVE EXPONENTIAL FOR FURTHER IDENTICAL REQUESTS
111 m100: IERROR=0;
112 AMUOL=AMU;
113 EXPMA=exp(-AMU);
114m200: PIR=1.;
115 N=-1;
116m300: N=N+1;
117 //{ // for debug
118 // double x = SRANLUX();
119 // mcout<<"pois: x = "<<x<<'\n';
120 // PIR=PIR * x;
121 //}
122 PIR=PIR * SRANLUX(); // working variant
123 if(PIR>EXPMA) goto m300;
124 return N;
125//C NORMAL APPROXIMATION FOR AMU.GT.AMAX
126m500:
127 rnorm(SRANLUX(), SRANLUX(), RAN,DUMMY);
128 N = long(RAN*sqrt(AMU) + AMU + .5);
129 return N;
130}
131*/
DoubleAc sqrt(const DoubleAc &f)
Definition: DoubleAc.cpp:313
DoubleAc exp(const DoubleAc &f)
Definition: DoubleAc.cpp:376
long pois(double AMU, int &IERROR)
Definition: pois.cpp:17
ffloat SRANLUX(void)
Definition: ranluxint.h:262
double rnorm_improved(void)
Definition: rnorm.cpp:15