| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // INCL++ intra-nuclear cascade model
27 // Alain Boudard, CEA-Saclay, France
28 // Joseph Cugnon, University of Liege, Belgium
29 // Jean-Christophe David, CEA-Saclay, France
30 // Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland
31 // Sylvie Leray, CEA-Saclay, France
32 // Davide Mancusi, CEA-Saclay, France
33 //
34 #define INCLXX_IN_GEANT4_MODE 1
35
36 #include "globals.hh"
37
38 /*
39 * G4INCLRandom.cc
40 *
41 * \date 7 June 2009
42 * \author Pekka Kaitaniemi
43 */
44
45 #include "G4INCLRandom.hh"
46 #include "G4INCLGlobals.hh"
47 // #include <cassert>
48
49 #include "G4INCLRanecu.hh"
50 #include "G4INCLRanecu3.hh"
51 #include "G4INCLGeant4Random.hh"
52 #include "G4INCLLogger.hh"
53
54 namespace G4INCL {
55
56 namespace Random {
57
58 namespace {
59
60 G4ThreadLocal IRandomGenerator* theGenerator = NULL;
61
62 #ifdef INCL_COUNT_RND_CALLS
63 G4ThreadLocal unsigned long long nCalls;
64 #endif
65
66 G4ThreadLocal SeedVector *savedSeeds = NULL;
67
68 G4ThreadLocal Adapter *theAdapter = NULL;
69
70 }
71
72 void setGenerator(G4INCL::IRandomGenerator *aGenerator) {
73 if(isInitialized()) {
74 INCL_ERROR("INCL random number generator already initialized." << '\n');
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 checking, 10.07.2015, D.Mancusi */
103 ;
104 return r;
105 }
106
107 G4double shoot1() {
108 G4double r;
109 while( (r=shoot()) >= 1. ) /* Loop checking, 10.07.2015, D.Mancusi */
110 ;
111 return r;
112 }
113
114 #ifndef INCLXX_IN_GEANT4_MODE
115 G4double gauss(G4double sigma) {
116 // generate a Gaussian random number with standard deviation sigma
117 // uses the flat() and flat0() methods
118 /* static G4ThreadLocal G4bool generated = false;
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))*std::cos(v);
127 }
128 else
129 {
130 generated = false;
131 return sigma*std::sqrt(-2*std::log(u))*std::sin(v);
132 }*/
133
134 return sigma*std::sqrt(-2*std::log(shoot0()))*std::cos(Math::twoPi*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 with standard deviation sigma
143 // uses the flat() and flat0() methods
144 static G4ThreadLocal G4bool generated = false;
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::log(u))*std::cos(v);
153 }
154 else
155 {
156 generated = false;
157 return sigma*std::sqrt(-2*std::log(u))*std::sin(v);
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.-ctheta*ctheta);
166 const G4double phi = Math::twoPi*shoot();
167 return ThreeVector(
168 norm * stheta * std::cos(phi),
169 norm * stheta * std::sin(phi),
170 norm * ctheta);
171
172 }
173
174 ThreeVector sphereVector(G4double rmax) {
175 return normVector( rmax*Math::pow13(shoot0()) );
176 }
177
178 ThreeVector gaussVector(G4double sigma) {
179 const G4double sigmax = sigma * Math::oneOverSqrtThree;
180 return ThreeVector(gauss(sigmax), gauss(sigmax), gauss(sigmax));
181 }
182
183 std::pair<G4double,G4double> correlatedGaussian(const G4double corrCoeff, const G4double x0, const G4double sigma) {
184 // assert(corrCoeff<=1. && corrCoeff>=-1.);
185 G4double factor = 1.-corrCoeff*corrCoeff;
186 if(factor<=0.)
187 factor=0.;
188 #ifndef INCLXX_IN_GEANT4_MODE
189 const G4double x = gauss(sigma) + x0;
190 const G4double y = corrCoeff * x + gauss(sigma*std::sqrt(factor)) + x0;
191 #else
192 const G4double x = gaussWithMemory(sigma) + x0;
193 const G4double y = corrCoeff * x + gaussWithMemory(sigma*std::sqrt(factor)) + x0;
194
195 #endif
196 return std::make_pair(x, y);
197 }
198
199 std::pair<G4double,G4double> correlatedUniform(const G4double corrCoeff) {
200 std::pair<G4double,G4double> gaussians = correlatedGaussian(corrCoeff);
201 return std::make_pair(Math::gaussianCDF(gaussians.first), Math::gaussianCDF(gaussians.second));
202 }
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 the RNG
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 Geant4RandomGenerator());
246 #else // INCLXX_IN_GEANT4_MODE
247 RNGType rng = theConfig->getRNGType();
248 if(rng == RanecuType)
249 setGenerator(new Ranecu(theConfig->getRandomSeeds()));
250 else if(rng == Ranecu3Type)
251 setGenerator(new Ranecu3(theConfig->getRandomSeeds()));
252 else
253 setGenerator(NULL);
254 #endif // INCLXX_IN_GEANT4_MODE
255 }
256
257 Adapter const &getAdapter() {
258 return *theAdapter;
259 }
260
261 }
262
263 }
264