ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4INCLRandom.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4INCLRandom.cc
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 
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 
90  return theGenerator->getSeeds();
91  }
92 
94 #ifdef INCL_COUNT_RND_CALLS
95  nCalls++;
96 #endif
97  return theGenerator->flat();
98  }
99 
101  G4double r;
102  while( (r=shoot()) <= 0. ) /* Loop checking, 10.07.2015, D.Mancusi */
103  ;
104  return r;
105  }
106 
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
138  return G4RandGauss::shoot(0.,sigma);
139  }
140 
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
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 
175  return normVector( rmax*Math::pow13(shoot0()) );
176  }
177 
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 
205  delete theGenerator;
206  theGenerator = NULL;
207  delete savedSeeds;
208  savedSeeds = NULL;
209  delete theAdapter;
210  theAdapter = NULL;
211  }
212 
214  if(theGenerator == 0) return false;
215  return true;
216  }
217 
218 #ifdef INCL_COUNT_RND_CALLS
219 
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 
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
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 }