| 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. )
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. )
110 ; 110 ;
111 return r; 111 return r;
112 } 112 }
113 113
114 #ifndef INCLXX_IN_GEANT4_MODE <<
115 G4double gauss(G4double sigma) { 114 G4double gauss(G4double sigma) {
116 // generate a Gaussian random number wit 115 // generate a Gaussian random number with standard deviation sigma
117 // uses the flat() and flat0() methods 116 // uses the flat() and flat0() methods
118 /* static G4ThreadLocal G4bool generated << 117 static G4ThreadLocal G4bool generated = false;
119 static G4ThreadLocal G4double u, v; 118 static G4ThreadLocal G4double u, v;
120 119
121 if( !generated ) 120 if( !generated )
122 { 121 {
123 u = shoot0(); 122 u = shoot0();
124 v = Math::twoPi*shoot(); 123 v = Math::twoPi*shoot();
125 generated = true; 124 generated = true;
126 return sigma*std::sqrt(-2*std::log(u)) 125 return sigma*std::sqrt(-2*std::log(u))*std::cos(v);
127 } 126 }
128 else 127 else
129 { 128 {
130 generated = false; 129 generated = false;
131 return sigma*std::sqrt(-2*std::log(u)) 130 return sigma*std::sqrt(-2*std::log(u))*std::sin(v);
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 } 131 }
>> 132 }
160 133
161 #endif <<
162 ThreeVector normVector(G4double norm) { 134 ThreeVector normVector(G4double norm) {
163 135
164 const G4double ctheta = (1.-2.*shoot()); 136 const G4double ctheta = (1.-2.*shoot());
165 const G4double stheta = std::sqrt(1.-cth 137 const G4double stheta = std::sqrt(1.-ctheta*ctheta);
166 const G4double phi = Math::twoPi*shoot() 138 const G4double phi = Math::twoPi*shoot();
167 return ThreeVector( 139 return ThreeVector(
168 norm * stheta * std:: 140 norm * stheta * std::cos(phi),
169 norm * stheta * std:: 141 norm * stheta * std::sin(phi),
170 norm * ctheta); 142 norm * ctheta);
171 143
172 } 144 }
173 145
174 ThreeVector sphereVector(G4double rmax) { 146 ThreeVector sphereVector(G4double rmax) {
175 return normVector( rmax*Math::pow13(shoo 147 return normVector( rmax*Math::pow13(shoot0()) );
176 } 148 }
177 149
178 ThreeVector gaussVector(G4double sigma) { 150 ThreeVector gaussVector(G4double sigma) {
179 const G4double sigmax = sigma * Math::on 151 const G4double sigmax = sigma * Math::oneOverSqrtThree;
180 return ThreeVector(gauss(sigmax), gauss( 152 return ThreeVector(gauss(sigmax), gauss(sigmax), gauss(sigmax));
181 } 153 }
182 154
183 std::pair<G4double,G4double> correlatedGau 155 std::pair<G4double,G4double> correlatedGaussian(const G4double corrCoeff, const G4double x0, const G4double sigma) {
184 // assert(corrCoeff<=1. && corrCoeff>=-1.); 156 // assert(corrCoeff<=1. && corrCoeff>=-1.);
185 G4double factor = 1.-corrCoeff*corrCoeff 157 G4double factor = 1.-corrCoeff*corrCoeff;
186 if(factor<=0.) 158 if(factor<=0.)
187 factor=0.; 159 factor=0.;
188 #ifndef INCLXX_IN_GEANT4_MODE <<
189 const G4double x = gauss(sigma) + x0; 160 const G4double x = gauss(sigma) + x0;
190 const G4double y = corrCoeff * x + gauss 161 const G4double y = corrCoeff * x + gauss(sigma*std::sqrt(factor)) + x0;
191 #else << 162 return std::make_pair(x, y);
192 const G4double x = gaussWithMemory(sigm <<
193 const G4double y = corrCoeff * x + gaus <<
194 <<
195 #endif <<
196 return std::make_pair(x, y); <<
197 } 163 }
198 164
199 std::pair<G4double,G4double> correlatedUni 165 std::pair<G4double,G4double> correlatedUniform(const G4double corrCoeff) {
200 std::pair<G4double,G4double> gaussians = 166 std::pair<G4double,G4double> gaussians = correlatedGaussian(corrCoeff);
201 return std::make_pair(Math::gaussianCDF( 167 return std::make_pair(Math::gaussianCDF(gaussians.first), Math::gaussianCDF(gaussians.second));
202 } 168 }
203 169
204 void deleteGenerator() { 170 void deleteGenerator() {
205 delete theGenerator; 171 delete theGenerator;
206 theGenerator = NULL; 172 theGenerator = NULL;
207 delete savedSeeds; 173 delete savedSeeds;
208 savedSeeds = NULL; 174 savedSeeds = NULL;
209 delete theAdapter; 175 delete theAdapter;
210 theAdapter = NULL; 176 theAdapter = NULL;
211 } 177 }
212 178
213 G4bool isInitialized() { 179 G4bool isInitialized() {
214 if(theGenerator == 0) return false; 180 if(theGenerator == 0) return false;
215 return true; 181 return true;
216 } 182 }
217 183
218 #ifdef INCL_COUNT_RND_CALLS 184 #ifdef INCL_COUNT_RND_CALLS
219 /// \brief Return the number of calls to t 185 /// \brief Return the number of calls to the RNG
220 unsigned long long getNumberOfCalls() { 186 unsigned long long getNumberOfCalls() {
221 return nCalls; 187 return nCalls;
222 } 188 }
223 #endif 189 #endif
224 190
225 void saveSeeds() { 191 void saveSeeds() {
226 if(!savedSeeds) 192 if(!savedSeeds)
227 savedSeeds = new SeedVector; 193 savedSeeds = new SeedVector;
228 194
229 (*savedSeeds) = theGenerator->getSeeds() 195 (*savedSeeds) = theGenerator->getSeeds();
230 } 196 }
231 197
232 SeedVector getSavedSeeds() { 198 SeedVector getSavedSeeds() {
233 if(!savedSeeds) 199 if(!savedSeeds)
234 savedSeeds = new SeedVector; 200 savedSeeds = new SeedVector;
235 201
236 return *savedSeeds; 202 return *savedSeeds;
237 } 203 }
238 204
239 void initialize(Config const * const 205 void initialize(Config const * const
240 #ifndef INCLXX_IN_GEANT4_MODE 206 #ifndef INCLXX_IN_GEANT4_MODE
241 theConfig 207 theConfig
242 #endif 208 #endif
243 ) { 209 ) {
244 #ifdef INCLXX_IN_GEANT4_MODE 210 #ifdef INCLXX_IN_GEANT4_MODE
245 Random::setGenerator(new Geant4RandomGen 211 Random::setGenerator(new Geant4RandomGenerator());
246 #else // INCLXX_IN_GEANT4_MODE 212 #else // INCLXX_IN_GEANT4_MODE
247 RNGType rng = theConfig->getRNGType(); 213 RNGType rng = theConfig->getRNGType();
248 if(rng == RanecuType) 214 if(rng == RanecuType)
249 setGenerator(new Ranecu(theConfig->get 215 setGenerator(new Ranecu(theConfig->getRandomSeeds()));
250 else if(rng == Ranecu3Type) 216 else if(rng == Ranecu3Type)
251 setGenerator(new Ranecu3(theConfig->ge 217 setGenerator(new Ranecu3(theConfig->getRandomSeeds()));
252 else 218 else
253 setGenerator(NULL); 219 setGenerator(NULL);
254 #endif // INCLXX_IN_GEANT4_MODE 220 #endif // INCLXX_IN_GEANT4_MODE
>> 221 }
>> 222
>> 223 G4int Adapter::operator()(const G4int n) const {
>> 224 return shootInteger(n);
255 } 225 }
256 226
257 Adapter const &getAdapter() { 227 Adapter const &getAdapter() {
258 return *theAdapter; 228 return *theAdapter;
259 } 229 }
260 230
261 } 231 }
262 232
263 } 233 }
264 234