| 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 // Pekka Kaitaniemi, CEA and Helsinki Institute of Physics
28 // Joseph Cugnon, University of Liege, Belgium << 28 // Davide Mancusi, CEA
29 // Jean-Christophe David, CEA-Saclay, France << 29 // Alain Boudard, CEA
30 // Pekka Kaitaniemi, CEA-Saclay, France, and H << 30 // Sylvie Leray, CEA
31 // Sylvie Leray, CEA-Saclay, France << 31 // Joseph Cugnon, University of Liege
32 // Davide Mancusi, CEA-Saclay, France <<
33 // 32 //
34 #define INCLXX_IN_GEANT4_MODE 1 33 #define INCLXX_IN_GEANT4_MODE 1
35 34
36 #include "globals.hh" 35 #include "globals.hh"
37 36
38 /* 37 /*
39 * G4INCLRandom.cc 38 * G4INCLRandom.cc
40 * 39 *
41 * \date 7 June 2009 40 * \date 7 June 2009
42 * \author Pekka Kaitaniemi 41 * \author Pekka Kaitaniemi
43 */ 42 */
44 43
45 #include "G4INCLRandom.hh" 44 #include "G4INCLRandom.hh"
46 #include "G4INCLGlobals.hh" 45 #include "G4INCLGlobals.hh"
47 // #include <cassert> <<
48 <<
49 #include "G4INCLRanecu.hh" <<
50 #include "G4INCLRanecu3.hh" <<
51 #include "G4INCLGeant4Random.hh" <<
52 #include "G4INCLLogger.hh" <<
53 46
54 namespace G4INCL { 47 namespace G4INCL {
55 48
56 namespace Random { << 49 G4INCL::IRandomGenerator* Random::theGenerator;
57 <<
58 namespace { <<
59 <<
60 G4ThreadLocal IRandomGenerator* theGener <<
61 <<
62 #ifdef INCL_COUNT_RND_CALLS <<
63 G4ThreadLocal unsigned long long nCalls; <<
64 #endif <<
65 <<
66 G4ThreadLocal SeedVector *savedSeeds = N <<
67 <<
68 G4ThreadLocal Adapter *theAdapter = NULL <<
69 <<
70 } <<
71 <<
72 void setGenerator(G4INCL::IRandomGenerator <<
73 if(isInitialized()) { <<
74 INCL_ERROR("INCL random number generat <<
75 } else { <<
76 #ifdef INCL_COUNT_RND_CALLS <<
77 nCalls = 0; <<
78 #endif <<
79 theGenerator = aGenerator; <<
80 } <<
81 if(!theAdapter) <<
82 theAdapter = new Adapter(); <<
83 } <<
84 <<
85 void setSeeds(const SeedVector &sv) { <<
86 theGenerator->setSeeds(sv); <<
87 } <<
88 <<
89 SeedVector getSeeds() { <<
90 return theGenerator->getSeeds(); <<
91 } <<
92 <<
93 G4double shoot() { <<
94 #ifdef INCL_COUNT_RND_CALLS <<
95 nCalls++; <<
96 #endif <<
97 return theGenerator->flat(); <<
98 } <<
99 <<
100 G4double shoot0() { <<
101 G4double r; <<
102 while( (r=shoot()) <= 0. ) /* Loop check <<
103 ; <<
104 return r; <<
105 } <<
106 <<
107 G4double shoot1() { <<
108 G4double r; <<
109 while( (r=shoot()) >= 1. ) /* Loop check <<
110 ; <<
111 return r; <<
112 } <<
113 <<
114 #ifndef INCLXX_IN_GEANT4_MODE <<
115 G4double gauss(G4double sigma) { <<
116 // generate a Gaussian random number wit <<
117 // uses the flat() and flat0() methods <<
118 /* static G4ThreadLocal G4bool generated <<
119 static G4ThreadLocal G4double u, v; <<
120 <<
121 if( !generated ) <<
122 { <<
123 u = shoot0(); <<
124 v = Math::twoPi*shoot(); <<
125 generated = true; <<
126 return sigma*std::sqrt(-2*std::log(u)) <<
127 } <<
128 else <<
129 { <<
130 generated = false; <<
131 return sigma*std::sqrt(-2*std::log(u)) <<
132 }*/ <<
133 <<
134 return sigma*std::sqrt(-2*std::log(shoot <<
135 } <<
136 #else <<
137 G4double gauss(G4double sigma) { <<
138 return G4RandGauss::shoot(0.,sigma); <<
139 } <<
140 <<
141 G4double gaussWithMemory(G4double sigma) <<
142 // generate a Gaussian random number <<
143 // uses the flat() and flat0() metho <<
144 static G4ThreadLocal G4bool generate <<
145 static G4ThreadLocal G4double u, v; <<
146 <<
147 if( !generated ) <<
148 { <<
149 u = shoot0(); <<
150 v = Math::twoPi*shoot(); <<
151 generated = true; <<
152 return sigma*std::sqrt(-2*std::l <<
153 } <<
154 else <<
155 { <<
156 generated = false; <<
157 return sigma*std::sqrt(-2*std::l <<
158 } <<
159 } <<
160 <<
161 #endif <<
162 ThreeVector normVector(G4double norm) { <<
163 <<
164 const G4double ctheta = (1.-2.*shoot()); <<
165 const G4double stheta = std::sqrt(1.-cth <<
166 const G4double phi = Math::twoPi*shoot() <<
167 return ThreeVector( <<
168 norm * stheta * std:: <<
169 norm * stheta * std:: <<
170 norm * ctheta); <<
171 <<
172 } <<
173 50
174 ThreeVector sphereVector(G4double rmax) { << 51 G4double Random::gauss(G4double sigma) {
175 return normVector( rmax*Math::pow13(shoo << 52 // generate a Gaussian random number with standard deviation sigma
176 } << 53 // uses the flat() and flat0() methods
177 << 54 static G4bool generated = false;
178 ThreeVector gaussVector(G4double sigma) { << 55 static G4double u, v;
179 const G4double sigmax = sigma * Math::on << 56
180 return ThreeVector(gauss(sigmax), gauss( << 57 if( !generated )
181 } << 58 {
182 << 59 u = shoot0();
183 std::pair<G4double,G4double> correlatedGau << 60 v = Math::twoPi*shoot();
184 // assert(corrCoeff<=1. && corrCoeff>=-1.); << 61 generated = true;
185 G4double factor = 1.-corrCoeff*corrCoeff << 62 return sigma*std::sqrt(-2*std::log(u))*std::cos(v);
186 if(factor<=0.) << 63 }
187 factor=0.; << 64 else
188 #ifndef INCLXX_IN_GEANT4_MODE << 65 {
189 const G4double x = gauss(sigma) + x0; << 66 generated = false;
190 const G4double y = corrCoeff * x + gauss << 67 return sigma*std::sqrt(-2*std::log(u))*std::sin(v);
191 #else << 68 }
192 const G4double x = gaussWithMemory(sigm << 69 }
193 const G4double y = corrCoeff * x + gaus << 70
194 << 71 ThreeVector Random::normVector(G4double norm) {
195 #endif << 72
196 return std::make_pair(x, y); << 73 const G4double ctheta = (1.-2.*shoot());
197 } << 74 const G4double stheta = std::sqrt(1.-ctheta*ctheta);
198 << 75 const G4double phi = Math::twoPi*shoot();
199 std::pair<G4double,G4double> correlatedUni << 76 return ThreeVector(
200 std::pair<G4double,G4double> gaussians = << 77 norm * stheta * std::cos(phi),
201 return std::make_pair(Math::gaussianCDF( << 78 norm * stheta * std::sin(phi),
202 } << 79 norm * ctheta);
203 <<
204 void deleteGenerator() { <<
205 delete theGenerator; <<
206 theGenerator = NULL; <<
207 delete savedSeeds; <<
208 savedSeeds = NULL; <<
209 delete theAdapter; <<
210 theAdapter = NULL; <<
211 } <<
212 <<
213 G4bool isInitialized() { <<
214 if(theGenerator == 0) return false; <<
215 return true; <<
216 } <<
217 <<
218 #ifdef INCL_COUNT_RND_CALLS <<
219 /// \brief Return the number of calls to t <<
220 unsigned long long getNumberOfCalls() { <<
221 return nCalls; <<
222 } <<
223 #endif <<
224 <<
225 void saveSeeds() { <<
226 if(!savedSeeds) <<
227 savedSeeds = new SeedVector; <<
228 <<
229 (*savedSeeds) = theGenerator->getSeeds() <<
230 } <<
231 <<
232 SeedVector getSavedSeeds() { <<
233 if(!savedSeeds) <<
234 savedSeeds = new SeedVector; <<
235 <<
236 return *savedSeeds; <<
237 } <<
238 <<
239 void initialize(Config const * const <<
240 #ifndef INCLXX_IN_GEANT4_MODE <<
241 theConfig <<
242 #endif <<
243 ) { <<
244 #ifdef INCLXX_IN_GEANT4_MODE <<
245 Random::setGenerator(new Geant4RandomGen <<
246 #else // INCLXX_IN_GEANT4_MODE <<
247 RNGType rng = theConfig->getRNGType(); <<
248 if(rng == RanecuType) <<
249 setGenerator(new Ranecu(theConfig->get <<
250 else if(rng == Ranecu3Type) <<
251 setGenerator(new Ranecu3(theConfig->ge <<
252 else <<
253 setGenerator(NULL); <<
254 #endif // INCLXX_IN_GEANT4_MODE <<
255 } <<
256 <<
257 Adapter const &getAdapter() { <<
258 return *theAdapter; <<
259 } <<
260 80
261 } 81 }
262 82
263 } 83 }
264 84