| Geant4 Cross Reference |
1 #include "CLHEP/Random/RandGaussZiggurat.h" 1
2 #include "CLHEP/Units/PhysicalConstants.h"
3 #include <iostream>
4 #include <cmath> // for std::log()
5
6 namespace CLHEP {
7
8 CLHEP_THREAD_LOCAL unsigned long RandGaussZigg
9 CLHEP_THREAD_LOCAL float RandGaussZiggurat::wn
10 CLHEP_THREAD_LOCAL bool RandGaussZiggurat::zig
11
12 HepRandomEngine & RandGaussZiggurat::engine()
13
14 RandGaussZiggurat::~RandGaussZiggurat() {
15 }
16
17 std::string RandGaussZiggurat::name() const
18 {
19 return "RandGaussZiggurat";
20 }
21
22 bool RandGaussZiggurat::ziggurat_init()
23 {
24 const double rzm1 = 2147483648.0, rzm2 = 429
25 double dn=3.442619855899,tn=dn,vn=9.91256303
26 double de=7.697117470131487, te=de, ve=3.949
27 int i;
28
29 /* Set up tables for RNOR */
30 q=vn/std::exp(-.5*dn*dn);
31 kn[0]=(unsigned long)((dn/q)*rzm1);
32 kn[1]=0;
33
34 wn[0]=q/rzm1;
35 wn[127]=dn/rzm1;
36
37 fn[0]=1.;
38 fn[127]=std::exp(-.5*dn*dn);
39
40 for(i=126;i>=1;i--) {
41 dn=std::sqrt(-2.*std::log(vn/dn+std::exp(-
42 kn[i+1]=(unsigned long)((dn/tn)*rzm1);
43 tn=dn;
44 fn[i]=std::exp(-.5*dn*dn);
45 wn[i]=dn/rzm1;
46 }
47
48 /* Set up tables for REXP */
49 q = ve/std::exp(-de);
50 ke[0]=(unsigned long)((de/q)*rzm2);
51 ke[1]=0;
52
53 we[0]=q/rzm2;
54 we[255]=de/rzm2;
55
56 fe[0]=1.;
57 fe[255]=std::exp(-de);
58
59 for(i=254;i>=1;i--) {
60 de=-std::log(ve/de+std::exp(-de));
61 ke[i+1]= (unsigned long)((de/te)*rzm2);
62 te=de;
63 fe[i]=std::exp(-de);
64 we[i]=de/rzm2;
65 }
66 ziggurat_is_init=true;
67
68 //std::cout<<"Done RandGaussZiggurat::ziggur
69
70 return true;
71 }
72
73 float RandGaussZiggurat::ziggurat_nfix(long hz
74 {
75 if(!ziggurat_is_init) ziggurat_init();
76 const float r = 3.442620f; /* The start
77 float x, y;
78 unsigned long iz=hz&127;
79 for(;;)
80 {
81 x=hz*wn[iz]; /* iz==0, handles the ba
82 if(iz==0) {
83 do {
84 /* change to (1.0 - UNI) as argument t
85 while the original UNI generates (0
86 x=-std::log(1.0 - ziggurat_UNI(anEngin
87 y=-std::log(1.0 - ziggurat_UNI(anEngin
88 } while(y+y<x*x);
89 return (hz>0)? r+x : -r-x;
90 }
91 /* iz>0, handle the wedges of other strips
92 if( fn[iz]+(1.0 - ziggurat_UNI(anEngine))*
93
94 /* initiate, try to exit for(;;) for loop*
95 hz=(signed)ziggurat_SHR3(anEngine);
96 iz=hz&127;
97 if((unsigned long)std::abs(hz)<kn[iz]) ret
98 }
99 }
100
101 double RandGaussZiggurat::operator()() {
102 return ziggurat_RNOR(localEngine.get()) * de
103 }
104
105 double RandGaussZiggurat::operator()( double m
106 return ziggurat_RNOR(localEngine.get()) * st
107 }
108
109 void RandGaussZiggurat::shootArray( const int
110 {
111 for (int i=0; i<size; ++i) {
112 vect[i] = shoot(mean,stdDev);
113 }
114 }
115
116 void RandGaussZiggurat::shootArray( const int
117 {
118 for (int i=0; i<size; ++i) {
119 vect[i] = shoot(mean,stdDev);
120 }
121 }
122
123 void RandGaussZiggurat::shootArray( HepRandomE
124 {
125 for (int i=0; i<size; ++i) {
126 vect[i] = shoot(anEngine,mean,stdDev);
127 }
128 }
129
130 void RandGaussZiggurat::shootArray( HepRandomE
131 {
132 for (int i=0; i<size; ++i) {
133 vect[i] = shoot(anEngine,mean,stdDev);
134 }
135 }
136
137 void RandGaussZiggurat::fireArray( const int s
138 {
139 for (int i=0; i<size; ++i) {
140 vect[i] = fire( defaultMean, defaultStdDe
141 }
142 }
143
144 void RandGaussZiggurat::fireArray( const int s
145 {
146 for (int i=0; i<size; ++i) {
147 vect[i] = fire( defaultMean, defaultStdDe
148 }
149 }
150
151 void RandGaussZiggurat::fireArray( const int s
152 {
153 for (int i=0; i<size; ++i) {
154 vect[i] = fire( mean, stdDev );
155 }
156 }
157
158 void RandGaussZiggurat::fireArray( const int s
159 {
160 for (int i=0; i<size; ++i) {
161 vect[i] = fire( mean, stdDev );
162 }
163 }
164
165 std::ostream & RandGaussZiggurat::put ( std::o
166 int pr=os.precision(20);
167 os << " " << name() << "\n";
168 RandGauss::put(os);
169 os.precision(pr);
170 return os;
171 }
172
173 std::istream & RandGaussZiggurat::get ( std::i
174 std::string inName;
175 is >> inName;
176 if (inName != name()) {
177 is.clear(std::ios::badbit | is.rdstate());
178 std::cerr << "Mismatch when expecting to r
179 << name() << " distribution\n"
180 << "Name found was " << inName
181 << "\nistream is left in the badbit st
182 return is;
183 }
184 RandGauss::get(is);
185 return is;
186 }
187
188 } // namespace CLHEP
189