ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4InuclSpecialFunctions.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4InuclSpecialFunctions.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 //
27 // 20100114 M. Kelsey -- Remove G4CascadeMomentum, use G4LorentzVector directly
28 // 20100914 M. Kelsey -- Migrate to integer A and Z. Discard pointless
29 // verbosity.
30 // 20120608 M. Kelsey -- Fix variable-name "shadowing" compiler warnings.
31 // 20130308 M. Kelsey -- New function to compute INUCL-style random value
32 // 20130314 M. Kelsey -- Restore null initializer and if-block for _TLS_.
33 // 20130924 M. Kelsey -- Use G4Log, G4Exp, G4Pow for CPU speedup
34 // 20150619 M. Kelsey -- Define G4cbrt(int) to use G4Pow::Z13, rearrange
35 // FermiEnergy() to use G4Pow::Z23.
36 // 20150622 M. Kelsey -- Use G4AutoDelete for _TLS_ buffers.
37 
38 #include <cmath>
39 
40 #include "G4InuclSpecialFunctions.hh"
41 #include "G4AutoDelete.hh"
42 #include "G4Exp.hh"
43 #include "G4Log.hh"
44 #include "G4LorentzVector.hh"
45 #include "G4PhysicalConstants.hh"
46 #include "G4Pow.hh"
47 #include "G4ThreeVector.hh"
48 #include "Randomize.hh"
49 
50 
51 // Compute power series in random value, with powers-of-Ekin coeffciences
52 
53 G4double
54 G4InuclSpecialFunctions::randomInuclPowers(G4double ekin,
55  const G4double (&coeff)[4][4]) {
56  G4Pow* theG4Pow = G4Pow::GetInstance();
57 
58  G4double S = G4UniformRand(); // Random fraction for expansion
59 
60  G4double C, V;
61  G4double PQ=0., PR=0.;
62  for (G4int i=0; i<4; i++) {
63  V = 0.0;
64  for (G4int k=0; k<4; k++) {
65  C = coeff[i][k];
66  V += C * theG4Pow->powN(ekin, k);
67  }
68 
69  PQ += V;
70  PR += V * theG4Pow->powN(S, i);
71  }
72 
73  return std::sqrt(S) * (PR + (1-PQ)*(S*S*S*S));
74 }
75 
76 
77 
78 G4double G4InuclSpecialFunctions::getAL(G4int A) {
79  return 0.76 + 2.2 / G4cbrt(A);
80 }
81 
82 G4double G4InuclSpecialFunctions::csNN(G4double e) {
83  G4double snn;
84 
85  if (e < 40.0) {
86  snn = -1174.8 / (e * e) + 3088.5 / e + 5.3107;
87  } else {
88  snn = 93074.0 / (e * e) - 11.148 / e + 22.429;
89  }
90 
91  return snn;
92 }
93 
94 G4double G4InuclSpecialFunctions::csPN(G4double e) {
95  G4double spn;
96 
97  if (e < 40.0) {
98  spn = -5057.4 / (e * e) + 9069.2 / e + 6.9466;
99  } else {
100  spn = 239380.0 / (e * e) + 1802.0 / e + 27.147;
101  }
102 
103  return spn;
104 }
105 
106 // calculates the nuclei Fermi energy for 0 - neutron and 1 - proton
107 
108 G4double G4InuclSpecialFunctions::FermiEnergy(G4int A, G4int Z, G4int ntype) {
109  G4Pow* g4pow = G4Pow::GetInstance();
110  const G4double C = 55.4 / g4pow->Z23(A);
111  G4double arg = (ntype==0) ? g4pow->Z23(A-Z) : g4pow->Z23(Z);
112 
113  return C * arg;
114 }
115 
116 G4double G4InuclSpecialFunctions::G4cbrt(G4double x) {
117  return x==0 ? 0. : (x<0?-1.:1.)*G4Exp(G4Log(std::fabs(x))/3.);
118 }
119 
120 G4double G4InuclSpecialFunctions::G4cbrt(G4int n) {
121  return n==0 ? 0. : (n<0?-1.:1.)*G4Pow::GetInstance()->Z13(std::abs(n));
122 }
123 
124 G4double G4InuclSpecialFunctions::inuclRndm() {
125  return G4UniformRand();
126 }
127 
128 G4double G4InuclSpecialFunctions::randomGauss(G4double sigma) {
129  const G4double eps = 1.0e-6;
130  G4double r1 = inuclRndm();
131  r1 = r1 > eps ? r1 : eps;
132  G4double r2 = inuclRndm();
133  r2 = r2 > eps ? r2 : eps;
134  r2 = r2 < 1.0 - eps ? r2 : 1.0 - eps;
135 
136  return sigma * std::sin(twopi * r1) * std::sqrt(-2.0 * G4Log(r2));
137 }
138 
139 G4double G4InuclSpecialFunctions::randomPHI() {
140  return twopi * inuclRndm();
141 }
142 
143 std::pair<G4double, G4double> G4InuclSpecialFunctions::randomCOS_SIN() {
144  G4double CT = 1.0 - 2.0 * inuclRndm();
145 
146  return std::pair<G4double, G4double>(CT, std::sqrt(1.0 - CT*CT));
147 }
148 
149 G4LorentzVector
150 G4InuclSpecialFunctions::generateWithFixedTheta(G4double ct, G4double p,
151  G4double mass) {
152  G4double phi = randomPHI();
153  G4double pt = p * std::sqrt(std::fabs(1.0 - ct * ct));
154 
155  // Buffers to avoid memory thrashing
156  static G4ThreadLocal G4ThreeVector *pvec_G4MT_TLS_ = 0;
157  if (!pvec_G4MT_TLS_) {
158  pvec_G4MT_TLS_ = new G4ThreeVector;
159  G4AutoDelete::Register(pvec_G4MT_TLS_);
160  }
161  G4ThreeVector &pvec = *pvec_G4MT_TLS_;
162 
163  static G4ThreadLocal G4LorentzVector *momr_G4MT_TLS_ = 0;
164  if (!momr_G4MT_TLS_) {
165  momr_G4MT_TLS_ = new G4LorentzVector;
166  G4AutoDelete::Register(momr_G4MT_TLS_);
167  }
168  G4LorentzVector &momr = *momr_G4MT_TLS_;
169 
170  pvec.set(pt*std::cos(phi), pt*std::sin(phi), p*ct);
171  momr.setVectM(pvec, mass);
172 
173  return momr;
174 }
175 
176 G4LorentzVector
177 G4InuclSpecialFunctions::generateWithRandomAngles(G4double p, G4double mass) {
178  std::pair<G4double, G4double> COS_SIN = randomCOS_SIN();
179  G4double phi = randomPHI();
180  G4double pt = p * COS_SIN.second;
181 
182  // Buffers to avoid memory thrashing
183  static G4ThreadLocal G4ThreeVector *pvec_G4MT_TLS_ = 0;
184  if (!pvec_G4MT_TLS_) {
185  pvec_G4MT_TLS_ = new G4ThreeVector;
186  G4AutoDelete::Register(pvec_G4MT_TLS_);
187  }
188  G4ThreeVector &pvec = *pvec_G4MT_TLS_;
189 
190  static G4ThreadLocal G4LorentzVector *momr_G4MT_TLS_ = 0;
191  if (!momr_G4MT_TLS_) {
192  momr_G4MT_TLS_ = new G4LorentzVector;
193  G4AutoDelete::Register(momr_G4MT_TLS_);
194  }
195  G4LorentzVector &momr = *momr_G4MT_TLS_;
196 
197  pvec.set(pt*std::cos(phi), pt*std::sin(phi), p*COS_SIN.first);
198  momr.setVectM(pvec, mass);
199 
200  return momr;
201 }