| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // INCL++ intra-nuclear cascade model 26 // INCL++ intra-nuclear cascade model
27 // Alain Boudard, CEA-Saclay, France 27 // Alain Boudard, CEA-Saclay, France
28 // Joseph Cugnon, University of Liege, Belgium 28 // Joseph Cugnon, University of Liege, Belgium
29 // Jean-Christophe David, CEA-Saclay, France 29 // Jean-Christophe David, CEA-Saclay, France
30 // Pekka Kaitaniemi, CEA-Saclay, France, and H 30 // Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland
31 // Sylvie Leray, CEA-Saclay, France 31 // Sylvie Leray, CEA-Saclay, France
32 // Davide Mancusi, CEA-Saclay, France 32 // Davide Mancusi, CEA-Saclay, France
33 // 33 //
34 #define INCLXX_IN_GEANT4_MODE 1 34 #define INCLXX_IN_GEANT4_MODE 1
35 35
36 #include "globals.hh" 36 #include "globals.hh"
37 37
38 /* 38 /*
39 * G4INCLRandom.cc 39 * G4INCLRandom.cc
40 * 40 *
41 * \date 7 June 2009 41 * \date 7 June 2009
42 * \author Pekka Kaitaniemi 42 * \author Pekka Kaitaniemi
43 */ 43 */
44 44
45 #include "G4INCLRandom.hh" 45 #include "G4INCLRandom.hh"
46 #include "G4INCLGlobals.hh" 46 #include "G4INCLGlobals.hh"
47 // #include <cassert> 47 // #include <cassert>
48 48
49 #include "G4INCLRanecu.hh" 49 #include "G4INCLRanecu.hh"
50 #include "G4INCLRanecu3.hh" 50 #include "G4INCLRanecu3.hh"
51 #include "G4INCLGeant4Random.hh" 51 #include "G4INCLGeant4Random.hh"
52 #include "G4INCLLogger.hh" 52 #include "G4INCLLogger.hh"
53 53
54 namespace G4INCL { 54 namespace G4INCL {
55 55
56 namespace Random { 56 namespace Random {
57 57
58 namespace { 58 namespace {
59 59
60 G4ThreadLocal IRandomGenerator* theGener 60 G4ThreadLocal IRandomGenerator* theGenerator = NULL;
61 61
62 #ifdef INCL_COUNT_RND_CALLS 62 #ifdef INCL_COUNT_RND_CALLS
63 G4ThreadLocal unsigned long long nCalls; 63 G4ThreadLocal unsigned long long nCalls;
64 #endif 64 #endif
65 65
66 G4ThreadLocal SeedVector *savedSeeds = N 66 G4ThreadLocal SeedVector *savedSeeds = NULL;
67 67
68 G4ThreadLocal Adapter *theAdapter = NULL 68 G4ThreadLocal Adapter *theAdapter = NULL;
69 69
70 } 70 }
71 71
72 void setGenerator(G4INCL::IRandomGenerator 72 void setGenerator(G4INCL::IRandomGenerator *aGenerator) {
73 if(isInitialized()) { 73 if(isInitialized()) {
74 INCL_ERROR("INCL random number generat 74 INCL_ERROR("INCL random number generator already initialized." << '\n');
75 } else { 75 } else {
76 #ifdef INCL_COUNT_RND_CALLS 76 #ifdef INCL_COUNT_RND_CALLS
77 nCalls = 0; 77 nCalls = 0;
78 #endif 78 #endif
79 theGenerator = aGenerator; 79 theGenerator = aGenerator;
80 } 80 }
81 if(!theAdapter) 81 if(!theAdapter)
82 theAdapter = new Adapter(); 82 theAdapter = new Adapter();
83 } 83 }
84 84
85 void setSeeds(const SeedVector &sv) { 85 void setSeeds(const SeedVector &sv) {
86 theGenerator->setSeeds(sv); 86 theGenerator->setSeeds(sv);
87 } 87 }
88 88
89 SeedVector getSeeds() { 89 SeedVector getSeeds() {
90 return theGenerator->getSeeds(); 90 return theGenerator->getSeeds();
91 } 91 }
92 92
93 G4double shoot() { 93 G4double shoot() {
94 #ifdef INCL_COUNT_RND_CALLS 94 #ifdef INCL_COUNT_RND_CALLS
95 nCalls++; 95 nCalls++;
96 #endif 96 #endif
97 return theGenerator->flat(); 97 return theGenerator->flat();
98 } 98 }
99 99
100 G4double shoot0() { 100 G4double shoot0() {
101 G4double r; 101 G4double r;
102 while( (r=shoot()) <= 0. ) /* Loop check 102 while( (r=shoot()) <= 0. ) /* Loop checking, 10.07.2015, D.Mancusi */
103 ; 103 ;
104 return r; 104 return r;
105 } 105 }
106 106
107 G4double shoot1() { 107 G4double shoot1() {
108 G4double r; 108 G4double r;
109 while( (r=shoot()) >= 1. ) /* Loop check 109 while( (r=shoot()) >= 1. ) /* Loop checking, 10.07.2015, D.Mancusi */
110 ; 110 ;
111 return r; 111 return r;
112 } 112 }
113 113
114 #ifndef INCLXX_IN_GEANT4_MODE 114 #ifndef INCLXX_IN_GEANT4_MODE
115 G4double gauss(G4double sigma) { 115 G4double gauss(G4double sigma) {
116 // generate a Gaussian random number wit 116 // generate a Gaussian random number with standard deviation sigma
117 // uses the flat() and flat0() methods 117 // uses the flat() and flat0() methods
118 /* static G4ThreadLocal G4bool generated 118 /* static G4ThreadLocal G4bool generated = false;
119 static G4ThreadLocal G4double u, v; 119 static G4ThreadLocal G4double u, v;
120 120
121 if( !generated ) 121 if( !generated )
122 { 122 {
123 u = shoot0(); 123 u = shoot0();
124 v = Math::twoPi*shoot(); 124 v = Math::twoPi*shoot();
125 generated = true; 125 generated = true;
126 return sigma*std::sqrt(-2*std::log(u)) 126 return sigma*std::sqrt(-2*std::log(u))*std::cos(v);
127 } 127 }
128 else 128 else
129 { 129 {
130 generated = false; 130 generated = false;
131 return sigma*std::sqrt(-2*std::log(u)) 131 return sigma*std::sqrt(-2*std::log(u))*std::sin(v);
132 }*/ 132 }*/
133 133
134 return sigma*std::sqrt(-2*std::log(shoot 134 return sigma*std::sqrt(-2*std::log(shoot0()))*std::cos(Math::twoPi*shoot());
135 } 135 }
136 #else 136 #else
137 G4double gauss(G4double sigma) { 137 G4double gauss(G4double sigma) {
138 return G4RandGauss::shoot(0.,sigma); 138 return G4RandGauss::shoot(0.,sigma);
139 } 139 }
140 140
141 G4double gaussWithMemory(G4double sigma) 141 G4double gaussWithMemory(G4double sigma) {
142 // generate a Gaussian random number 142 // generate a Gaussian random number with standard deviation sigma
143 // uses the flat() and flat0() metho 143 // uses the flat() and flat0() methods
144 static G4ThreadLocal G4bool generate 144 static G4ThreadLocal G4bool generated = false;
145 static G4ThreadLocal G4double u, v; 145 static G4ThreadLocal G4double u, v;
146 146
147 if( !generated ) 147 if( !generated )
148 { 148 {
149 u = shoot0(); 149 u = shoot0();
150 v = Math::twoPi*shoot(); 150 v = Math::twoPi*shoot();
151 generated = true; 151 generated = true;
152 return sigma*std::sqrt(-2*std::l 152 return sigma*std::sqrt(-2*std::log(u))*std::cos(v);
153 } 153 }
154 else 154 else
155 { 155 {
156 generated = false; 156 generated = false;
157 return sigma*std::sqrt(-2*std::l 157 return sigma*std::sqrt(-2*std::log(u))*std::sin(v);
158 } 158 }
159 } 159 }
160 160
161 #endif 161 #endif
162 ThreeVector normVector(G4double norm) { 162 ThreeVector normVector(G4double norm) {
163 163
164 const G4double ctheta = (1.-2.*shoot()); 164 const G4double ctheta = (1.-2.*shoot());
165 const G4double stheta = std::sqrt(1.-cth 165 const G4double stheta = std::sqrt(1.-ctheta*ctheta);
166 const G4double phi = Math::twoPi*shoot() 166 const G4double phi = Math::twoPi*shoot();
167 return ThreeVector( 167 return ThreeVector(
168 norm * stheta * std:: 168 norm * stheta * std::cos(phi),
169 norm * stheta * std:: 169 norm * stheta * std::sin(phi),
170 norm * ctheta); 170 norm * ctheta);
171 171
172 } 172 }
173 173
174 ThreeVector sphereVector(G4double rmax) { 174 ThreeVector sphereVector(G4double rmax) {
175 return normVector( rmax*Math::pow13(shoo 175 return normVector( rmax*Math::pow13(shoot0()) );
176 } 176 }
177 177
178 ThreeVector gaussVector(G4double sigma) { 178 ThreeVector gaussVector(G4double sigma) {
179 const G4double sigmax = sigma * Math::on 179 const G4double sigmax = sigma * Math::oneOverSqrtThree;
180 return ThreeVector(gauss(sigmax), gauss( 180 return ThreeVector(gauss(sigmax), gauss(sigmax), gauss(sigmax));
181 } 181 }
182 182
183 std::pair<G4double,G4double> correlatedGau 183 std::pair<G4double,G4double> correlatedGaussian(const G4double corrCoeff, const G4double x0, const G4double sigma) {
184 // assert(corrCoeff<=1. && corrCoeff>=-1.); 184 // assert(corrCoeff<=1. && corrCoeff>=-1.);
185 G4double factor = 1.-corrCoeff*corrCoeff 185 G4double factor = 1.-corrCoeff*corrCoeff;
186 if(factor<=0.) 186 if(factor<=0.)
187 factor=0.; 187 factor=0.;
188 #ifndef INCLXX_IN_GEANT4_MODE 188 #ifndef INCLXX_IN_GEANT4_MODE
189 const G4double x = gauss(sigma) + x0; 189 const G4double x = gauss(sigma) + x0;
190 const G4double y = corrCoeff * x + gauss 190 const G4double y = corrCoeff * x + gauss(sigma*std::sqrt(factor)) + x0;
191 #else 191 #else
192 const G4double x = gaussWithMemory(sigm 192 const G4double x = gaussWithMemory(sigma) + x0;
193 const G4double y = corrCoeff * x + gaus 193 const G4double y = corrCoeff * x + gaussWithMemory(sigma*std::sqrt(factor)) + x0;
194 194
195 #endif 195 #endif
196 return std::make_pair(x, y); 196 return std::make_pair(x, y);
197 } 197 }
198 198
199 std::pair<G4double,G4double> correlatedUni 199 std::pair<G4double,G4double> correlatedUniform(const G4double corrCoeff) {
200 std::pair<G4double,G4double> gaussians = 200 std::pair<G4double,G4double> gaussians = correlatedGaussian(corrCoeff);
201 return std::make_pair(Math::gaussianCDF( 201 return std::make_pair(Math::gaussianCDF(gaussians.first), Math::gaussianCDF(gaussians.second));
202 } 202 }
203 203
204 void deleteGenerator() { 204 void deleteGenerator() {
205 delete theGenerator; 205 delete theGenerator;
206 theGenerator = NULL; 206 theGenerator = NULL;
207 delete savedSeeds; 207 delete savedSeeds;
208 savedSeeds = NULL; 208 savedSeeds = NULL;
209 delete theAdapter; 209 delete theAdapter;
210 theAdapter = NULL; 210 theAdapter = NULL;
211 } 211 }
212 212
213 G4bool isInitialized() { 213 G4bool isInitialized() {
214 if(theGenerator == 0) return false; 214 if(theGenerator == 0) return false;
215 return true; 215 return true;
216 } 216 }
217 217
218 #ifdef INCL_COUNT_RND_CALLS 218 #ifdef INCL_COUNT_RND_CALLS
219 /// \brief Return the number of calls to t 219 /// \brief Return the number of calls to the RNG
220 unsigned long long getNumberOfCalls() { 220 unsigned long long getNumberOfCalls() {
221 return nCalls; 221 return nCalls;
222 } 222 }
223 #endif 223 #endif
224 224
225 void saveSeeds() { 225 void saveSeeds() {
226 if(!savedSeeds) 226 if(!savedSeeds)
227 savedSeeds = new SeedVector; 227 savedSeeds = new SeedVector;
228 228
229 (*savedSeeds) = theGenerator->getSeeds() 229 (*savedSeeds) = theGenerator->getSeeds();
230 } 230 }
231 231
232 SeedVector getSavedSeeds() { 232 SeedVector getSavedSeeds() {
233 if(!savedSeeds) 233 if(!savedSeeds)
234 savedSeeds = new SeedVector; 234 savedSeeds = new SeedVector;
235 235
236 return *savedSeeds; 236 return *savedSeeds;
237 } 237 }
238 238
239 void initialize(Config const * const 239 void initialize(Config const * const
240 #ifndef INCLXX_IN_GEANT4_MODE 240 #ifndef INCLXX_IN_GEANT4_MODE
241 theConfig 241 theConfig
242 #endif 242 #endif
243 ) { 243 ) {
244 #ifdef INCLXX_IN_GEANT4_MODE 244 #ifdef INCLXX_IN_GEANT4_MODE
245 Random::setGenerator(new Geant4RandomGen 245 Random::setGenerator(new Geant4RandomGenerator());
246 #else // INCLXX_IN_GEANT4_MODE 246 #else // INCLXX_IN_GEANT4_MODE
247 RNGType rng = theConfig->getRNGType(); 247 RNGType rng = theConfig->getRNGType();
248 if(rng == RanecuType) 248 if(rng == RanecuType)
249 setGenerator(new Ranecu(theConfig->get 249 setGenerator(new Ranecu(theConfig->getRandomSeeds()));
250 else if(rng == Ranecu3Type) 250 else if(rng == Ranecu3Type)
251 setGenerator(new Ranecu3(theConfig->ge 251 setGenerator(new Ranecu3(theConfig->getRandomSeeds()));
252 else 252 else
253 setGenerator(NULL); 253 setGenerator(NULL);
254 #endif // INCLXX_IN_GEANT4_MODE 254 #endif // INCLXX_IN_GEANT4_MODE
255 } 255 }
256 256
257 Adapter const &getAdapter() { 257 Adapter const &getAdapter() {
258 return *theAdapter; 258 return *theAdapter;
259 } 259 }
260 260
261 } 261 }
262 262
263 } 263 }
264 264