| Geant4 Cross Reference |
1 #include "CLHEP/Random/DoubConv.h" 1 #include "CLHEP/Random/DoubConv.h"
2 2
3 #include "CLHEP/Random/RandExpZiggurat.h" 3 #include "CLHEP/Random/RandExpZiggurat.h"
4 4
5 #include <cmath> // for std::log() 5 #include <cmath> // for std::log()
6 #include <iostream> 6 #include <iostream>
7 #include <vector> 7 #include <vector>
8 8
9 namespace CLHEP { 9 namespace CLHEP {
10 10
11 CLHEP_THREAD_LOCAL unsigned long RandExpZiggur 11 CLHEP_THREAD_LOCAL unsigned long RandExpZiggurat::kn[128], RandExpZiggurat::ke[256];
12 CLHEP_THREAD_LOCAL float RandExpZiggurat::wn[1 12 CLHEP_THREAD_LOCAL float RandExpZiggurat::wn[128],RandExpZiggurat::fn[128],RandExpZiggurat::we[256],RandExpZiggurat::fe[256];
13 CLHEP_THREAD_LOCAL bool RandExpZiggurat::ziggu 13 CLHEP_THREAD_LOCAL bool RandExpZiggurat::ziggurat_is_init = false;
14 14
15 std::string RandExpZiggurat::name() const {ret 15 std::string RandExpZiggurat::name() const {return "RandExpZiggurat";}
16 16
17 HepRandomEngine & RandExpZiggurat::engine() {r 17 HepRandomEngine & RandExpZiggurat::engine() {return *localEngine;}
18 18
19 RandExpZiggurat::~RandExpZiggurat() { 19 RandExpZiggurat::~RandExpZiggurat() {
20 } 20 }
21 21
22 RandExpZiggurat::RandExpZiggurat(const RandExp 22 RandExpZiggurat::RandExpZiggurat(const RandExpZiggurat& right) : HepRandom(right),defaultMean(right.defaultMean)
23 { 23 {
24 } 24 }
25 25
26 double RandExpZiggurat::operator()() 26 double RandExpZiggurat::operator()()
27 { 27 {
28 return fire( defaultMean ); 28 return fire( defaultMean );
29 } 29 }
30 30
31 void RandExpZiggurat::shootArray( const int si 31 void RandExpZiggurat::shootArray( const int size, float* vect, float mean )
32 { 32 {
33 for (int i=0; i<size; ++i) vect[i] = shoot( 33 for (int i=0; i<size; ++i) vect[i] = shoot(mean);
34 } 34 }
35 35
36 void RandExpZiggurat::shootArray( const int si 36 void RandExpZiggurat::shootArray( const int size, double* vect, double mean )
37 { 37 {
38 for (int i=0; i<size; ++i) vect[i] = shoot( 38 for (int i=0; i<size; ++i) vect[i] = shoot(mean);
39 } 39 }
40 40
41 void RandExpZiggurat::shootArray(HepRandomEngi 41 void RandExpZiggurat::shootArray(HepRandomEngine* anEngine, const int size, float* vect, float mean )
42 { 42 {
43 for (int i=0; i<size; ++i) vect[i] = shoot( 43 for (int i=0; i<size; ++i) vect[i] = shoot(anEngine, mean);
44 } 44 }
45 45
46 void RandExpZiggurat::shootArray(HepRandomEngi 46 void RandExpZiggurat::shootArray(HepRandomEngine* anEngine, const int size, double* vect, double mean )
47 { 47 {
48 for (int i=0; i<size; ++i) vect[i] = shoot( 48 for (int i=0; i<size; ++i) vect[i] = shoot(anEngine, mean);
49 } 49 }
50 50
51 void RandExpZiggurat::fireArray( const int siz 51 void RandExpZiggurat::fireArray( const int size, float* vect)
52 { 52 {
53 for (int i=0; i<size; ++i) vect[i] = fire( 53 for (int i=0; i<size; ++i) vect[i] = fire( defaultMean );
54 } 54 }
55 55
56 void RandExpZiggurat::fireArray( const int siz 56 void RandExpZiggurat::fireArray( const int size, double* vect)
57 { 57 {
58 for (int i=0; i<size; ++i) vect[i] = fire( 58 for (int i=0; i<size; ++i) vect[i] = fire( defaultMean );
59 } 59 }
60 60
61 void RandExpZiggurat::fireArray( const int siz 61 void RandExpZiggurat::fireArray( const int size, float* vect, float mean )
62 { 62 {
63 for (int i=0; i<size; ++i) vect[i] = fire( 63 for (int i=0; i<size; ++i) vect[i] = fire( mean );
64 } 64 }
65 65
66 void RandExpZiggurat::fireArray( const int siz 66 void RandExpZiggurat::fireArray( const int size, double* vect, double mean )
67 { 67 {
68 for (int i=0; i<size; ++i) vect[i] = fire( 68 for (int i=0; i<size; ++i) vect[i] = fire( mean );
69 } 69 }
70 70
71 std::ostream & RandExpZiggurat::put ( std::ost 71 std::ostream & RandExpZiggurat::put ( std::ostream & os ) const {
72 int pr=os.precision(20); 72 int pr=os.precision(20);
73 std::vector<unsigned long> t(2); 73 std::vector<unsigned long> t(2);
74 os << " " << name() << "\n"; 74 os << " " << name() << "\n";
75 os << "Uvec" << "\n"; 75 os << "Uvec" << "\n";
76 t = DoubConv::dto2longs(defaultMean); 76 t = DoubConv::dto2longs(defaultMean);
77 os << defaultMean << " " << t[0] << " " << t 77 os << defaultMean << " " << t[0] << " " << t[1] << "\n";
78 os.precision(pr); 78 os.precision(pr);
79 return os; 79 return os;
80 #ifdef REMOVED 80 #ifdef REMOVED
81 int pr=os.precision(20); 81 int pr=os.precision(20);
82 os << " " << name() << "\n"; 82 os << " " << name() << "\n";
83 os << defaultMean << "\n"; 83 os << defaultMean << "\n";
84 os.precision(pr); 84 os.precision(pr);
85 return os; 85 return os;
86 #endif 86 #endif
87 } 87 }
88 88
89 std::istream & RandExpZiggurat::get ( std::ist 89 std::istream & RandExpZiggurat::get ( std::istream & is ) {
90 std::string inName; 90 std::string inName;
91 is >> inName; 91 is >> inName;
92 if (inName != name()) { 92 if (inName != name()) {
93 is.clear(std::ios::badbit | is.rdstate()); 93 is.clear(std::ios::badbit | is.rdstate());
94 std::cerr << "Mismatch when expecting to r 94 std::cerr << "Mismatch when expecting to read state of a "
95 << name() << " distribution\n" 95 << name() << " distribution\n"
96 << "Name found was " << inName 96 << "Name found was " << inName
97 << "\nistream is left in the badbit st 97 << "\nistream is left in the badbit state\n";
98 return is; 98 return is;
99 } 99 }
100 if (possibleKeywordInput(is, "Uvec", default 100 if (possibleKeywordInput(is, "Uvec", defaultMean)) {
101 std::vector<unsigned long> t(2); 101 std::vector<unsigned long> t(2);
102 is >> defaultMean >> t[0] >> t[1]; default 102 is >> defaultMean >> t[0] >> t[1]; defaultMean = DoubConv::longs2double(t);
103 return is; 103 return is;
104 } 104 }
105 // is >> defaultMean encompassed by possible 105 // is >> defaultMean encompassed by possibleKeywordInput
106 return is; 106 return is;
107 } 107 }
108 108
109 109
110 float RandExpZiggurat::ziggurat_efix(unsigned 110 float RandExpZiggurat::ziggurat_efix(unsigned long jz,HepRandomEngine* anEngine)
111 { 111 {
112 if(!ziggurat_is_init) ziggurat_init(); 112 if(!ziggurat_is_init) ziggurat_init();
113 113
114 unsigned long iz=jz&255; 114 unsigned long iz=jz&255;
115 115
116 float x; 116 float x;
117 for(;;) 117 for(;;)
118 { 118 {
119 if(iz==0) return (7.69711-std::log(ziggura 119 if(iz==0) return (7.69711-std::log(ziggurat_UNI(anEngine))); /* iz==0 */
120 x=jz*we[iz]; 120 x=jz*we[iz];
121 if( fe[iz]+ziggurat_UNI(anEngine)*(fe[iz-1 121 if( fe[iz]+ziggurat_UNI(anEngine)*(fe[iz-1]-fe[iz]) < std::exp(-x) ) return (x);
122 122
123 /* initiate, try to exit for(;;) loop */ 123 /* initiate, try to exit for(;;) loop */
124 jz=ziggurat_SHR3(anEngine); 124 jz=ziggurat_SHR3(anEngine);
125 iz=(jz&255); 125 iz=(jz&255);
126 if(jz<ke[iz]) return (jz*we[iz]); 126 if(jz<ke[iz]) return (jz*we[iz]);
127 } 127 }
128 } 128 }
129 129
130 bool RandExpZiggurat::ziggurat_init() 130 bool RandExpZiggurat::ziggurat_init()
131 { 131 {
132 const double rzm1 = 2147483648.0, rzm2 = 429 132 const double rzm1 = 2147483648.0, rzm2 = 4294967296.;
133 double dn=3.442619855899,tn=dn,vn=9.91256303 133 double dn=3.442619855899,tn=dn,vn=9.91256303526217e-3, q;
134 double de=7.697117470131487, te=de, ve=3.949 134 double de=7.697117470131487, te=de, ve=3.949659822581572e-3;
135 int i; 135 int i;
136 136
137 /* Set up tables for RNOR */ 137 /* Set up tables for RNOR */
138 q=vn/std::exp(-.5*dn*dn); 138 q=vn/std::exp(-.5*dn*dn);
139 kn[0]=(unsigned long)((dn/q)*rzm1); 139 kn[0]=(unsigned long)((dn/q)*rzm1);
140 kn[1]=0; 140 kn[1]=0;
141 141
142 wn[0]=q/rzm1; 142 wn[0]=q/rzm1;
143 wn[127]=dn/rzm1; 143 wn[127]=dn/rzm1;
144 144
145 fn[0]=1.; 145 fn[0]=1.;
146 fn[127]=std::exp(-.5*dn*dn); 146 fn[127]=std::exp(-.5*dn*dn);
147 147
148 for(i=126;i>=1;i--) { 148 for(i=126;i>=1;i--) {
149 dn=std::sqrt(-2.*std::log(vn/dn+std::exp(- 149 dn=std::sqrt(-2.*std::log(vn/dn+std::exp(-.5*dn*dn)));
150 kn[i+1]=(unsigned long)((dn/tn)*rzm1); 150 kn[i+1]=(unsigned long)((dn/tn)*rzm1);
151 tn=dn; 151 tn=dn;
152 fn[i]=std::exp(-.5*dn*dn); 152 fn[i]=std::exp(-.5*dn*dn);
153 wn[i]=dn/rzm1; 153 wn[i]=dn/rzm1;
154 } 154 }
155 155
156 /* Set up tables for REXP */ 156 /* Set up tables for REXP */
157 q = ve/std::exp(-de); 157 q = ve/std::exp(-de);
158 ke[0]=(unsigned long)((de/q)*rzm2); 158 ke[0]=(unsigned long)((de/q)*rzm2);
159 ke[1]=0; 159 ke[1]=0;
160 160
161 we[0]=q/rzm2; 161 we[0]=q/rzm2;
162 we[255]=de/rzm2; 162 we[255]=de/rzm2;
163 163
164 fe[0]=1.; 164 fe[0]=1.;
165 fe[255]=std::exp(-de); 165 fe[255]=std::exp(-de);
166 166
167 for(i=254;i>=1;i--) { 167 for(i=254;i>=1;i--) {
168 de=-std::log(ve/de+std::exp(-de)); 168 de=-std::log(ve/de+std::exp(-de));
169 ke[i+1]= (unsigned long)((de/te)*rzm2); 169 ke[i+1]= (unsigned long)((de/te)*rzm2);
170 te=de; 170 te=de;
171 fe[i]=std::exp(-de); 171 fe[i]=std::exp(-de);
172 we[i]=de/rzm2; 172 we[i]=de/rzm2;
173 } 173 }
174 ziggurat_is_init=true; 174 ziggurat_is_init=true;
175 return true; 175 return true;
176 } 176 }
177 177
178 } // namespace CLHEP 178 } // namespace CLHEP
179 179