| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 // Geant4 class G4NuclearRadii
28 //
29 // Author V.Ivanchenko 27.05.2019
30 //
31 //
32
33 #include "G4NuclearRadii.hh"
34 #include "G4Pow.hh"
35 #include "G4PhysicalConstants.hh"
36 #include "G4ParticleDefinition.hh"
37 #include "G4NucleiProperties.hh"
38 #include "G4IsotopeList.hh"
39 #include "G4Log.hh"
40 #include "G4Exp.hh"
41
42 G4Pow* G4NuclearRadii::fG4pow = G4Pow::GetInst
43
44 namespace
45 {
46 const G4double llog10 = G4Log(10.);
47 const G4double fAlpha = 0.5*CLHEP::fine_stru
48 const G4double fInvep = 1.0/CLHEP::eplus;
49 }
50
51 G4double G4NuclearRadii::ExplicitRadius(G4int
52 {
53 G4double R = 0.0;
54 // Special rms radii for light nucleii
55 if(Z <= 4) {
56 if(A == 1) { R = 0.895*CLHE
57 else if(A == 2) { R = 2.13*CLHEP
58 else if(Z == 1 && A == 3) { R = 1.80*CLHEP
59 else if(Z == 2 && A == 3) { R = 1.96*CLHEP
60 else if(Z == 2 && A == 4) { R = 1.68*CLHEP
61 else if(Z == 3) { R = 2.40*CLHEP
62 else if(Z == 4) { R = 2.51*CLHEP
63 }
64 return R;
65 }
66
67 G4double G4NuclearRadii::Radius(G4int Z, G4int
68 {
69 G4double R = ExplicitRadius(Z, A);
70 if(0.0 == R) {
71 if (A <= 50) {
72 G4double y = 1.1;
73 if( A <= 15) { y = 1.26; }
74 else if( A <= 20) { y = 1.19; }
75 else if( A <= 30) { y = 1.12; }
76 G4double x = fG4pow->Z13(A);
77 R = y*(x - 1./x);
78 } else {
79 R = fG4pow->powZ(A, 0.27);
80 }
81 R *= CLHEP::fermi;
82 }
83 return R;
84 }
85
86 G4double G4NuclearRadii::RadiusRMS(G4int Z, G4
87 {
88 G4double R = ExplicitRadius(Z, A);
89 if(0.0 == R) {
90 R = 1.24*fG4pow->powZ(A, 0.28)*CLHEP::ferm
91 }
92 return R;
93 }
94
95 G4double G4NuclearRadii::RadiusNNGG(G4int Z, G
96 {
97 G4double R = ExplicitRadius(Z, A);
98 if(0.0 == R) {
99 if(A > 20) {
100 R = 1.08*fG4pow->Z13(A)*(0.85 + 0.15*G4E
101 } else {
102 R = 1.08*fG4pow->Z13(A)*(1.0 + 0.3*G4Exp
103 }
104 R *= CLHEP::fermi;
105 }
106 return R;
107 }
108
109 G4double G4NuclearRadii::RadiusECS(G4int Z, G4
110 {
111 G4double R=0.;
112 const G4double c[3]={0.77329745, 1.38206072,
113 const G4double c1=c[0];
114 const G4double c2=c[1];
115 const G4double c3=c[2];
116
117 // Special rms radii for light nuclei
118 if (A <= 30) {
119 G4double vn = 0.5*A + fG4pow->powN(0.028*A
120 G4double dev = vn - (A-Z);
121 R = c1*fG4pow->Z13(A) + c2/fG4pow->Z13(A)
122 } else if (A<=50){
123 G4double y = 1.1;
124 G4double x = fG4pow->Z13(A);
125 R = y*(x - 1./x);
126 }
127 return R*CLHEP::fermi;
128 }
129
130 G4double G4NuclearRadii::RadiusHNGG(G4int A)
131 {
132 G4double R = CLHEP::fermi;
133 if(A > 20) {
134 R *= 1.08*fG4pow->Z13(A)*(0.8 + 0.2*G4Exp(
135 } else {
136 R *= 1.08*fG4pow->Z13(A)*(1.0 + 0.1*G4Exp(
137 }
138 return R;
139 }
140
141 G4double G4NuclearRadii::RadiusKNGG(G4int A)
142 {
143 return 1.3*CLHEP::fermi*fG4pow->Z13(A);
144 }
145
146 G4double G4NuclearRadii::RadiusND(G4int A)
147 {
148 G4double R = CLHEP::fermi;
149 if (1 == A) {
150 R *= 0.895;
151 } else {
152 R *= fG4pow->Z13(A);
153 if (A <= 3.) { R *= 0.8; }
154 else { R *= 1.7; }
155 }
156 return R;
157 }
158
159 G4double G4NuclearRadii::RadiusCB(G4int Z, G4i
160 {
161 G4double R = ExplicitRadius(Z, A);
162 if(0.0 == R) {
163 G4int z = std::min(Z, 92);
164 R = r0[z]*fG4pow->Z13(A)*CLHEP::fermi;
165 }
166 return R;
167 }
168
169 G4double G4NuclearRadii::ParticleRadius(const
170 {
171 G4double R = CLHEP::fermi;
172 G4int pdg = std::abs(p->GetPDGEncoding());
173 if(pdg == 2112 || pdg == 2212) { R *= 0.89
174 else if(pdg == 211) { R *= 0.663; }
175 else if(pdg == 321) { R *= 0.340; }
176 else { R *= 0.5; }
177 return R;
178 }
179
180 G4double G4NuclearRadii::CoulombFactor(
181 const G4ParticleDefinition* thePartic
182 const G4ParticleDefinition* nucleon,
183 G4double ekin)
184 {
185 G4double tR = 0.895*CLHEP::fermi;
186 G4double pR = ParticleRadius(theParticle);
187
188 G4double pZ = theParticle->GetPDGCharge()*fI
189 G4double tZ = nucleon->GetPDGCharge()*fInvep
190
191 G4double pM = theParticle->GetPDGMass();
192 G4double tM = nucleon->GetPDGMass();
193
194 G4double pElab = ekin + pM;
195 G4double totTcm = std::sqrt(pM*pM + tM*tM +
196
197 G4double bC = fAlpha*pZ*tZ/(pR + tR);
198 return (totTcm > bC) ? 1. - bC/totTcm : 0.0;
199 }
200
201 G4double G4NuclearRadii::CoulombFactor(
202 G4int Z, G4int A,
203 const G4ParticleDefinition* theParticle,
204 G4double ekin)
205 {
206 G4double tR = RadiusCB(Z, A);
207 G4double pR = ParticleRadius(theParticle);
208
209 G4double pZ = theParticle->GetPDGCharge()*fI
210
211 G4double pM = theParticle->GetPDGMass();
212 G4double tM = G4NucleiProperties::GetNuclear
213
214 G4double pElab = ekin + pM;
215 G4double totTcm = std::sqrt(pM*pM + tM*tM +
216
217 G4double bC = fAlpha*pZ*Z/(pR + tR);
218 return (totTcm > bC) ? 1. - bC/totTcm : 0.0;
219 }
220
221 G4double
222 G4NuclearRadii::NeutronInelasticShape(G4int Z,
223 {
224 G4double A = (Z < 100) ? aeff[Z] : aeff[100]
225 G4double elog = G4Log(ekin/CLHEP::GeV)/llog1
226 G4double p3 = 0.6 + 13./A - 0.0005*A;
227 G4double p4 = 7.2449 - 0.018242*A;
228 G4double p5 = 1.36 + 1.8/A + 0.0005*A;
229 G4double p6 = 1. + 200./A + 0.02*A;
230 G4double p7 = 3.0 - (A - 70.)*(A - 200.)/110
231
232 G4double firstexp = G4Exp(-p4*(elog + p5));
233 G4double secondexp = G4Exp(-p6*(elog + p7));
234
235 return (1. + p3*firstexp/(1. + firstexp))/(1
236 }
237
238 G4double
239 G4NuclearRadii::ProtonInelasticShape(G4int Z,
240 {
241 G4double A = (Z < 100) ? aeff[Z] : aeff[100]
242 G4double elog = G4Log(ekin/CLHEP::GeV)/llog1
243 G4double ff1 = 5.6 - 0.016*A; // slope of t
244 G4double ff2 = 1.37 + 1.37/A; // start of t
245 G4double ff3 = 0.8 + 18./A - 0.002*A; //
246 G4double res = (1.0 + ff3*(1.0 - (1.0/(1+G4E
247 ff1 = 8. - 8./A - 0.008*A; // slope of the r
248 ff2 = 2.34 - 5.4/A - 0.0028*A; // start of t
249 res /= (1.0 + G4Exp(-ff1*(elog + ff2)));
250 return res;
251 }
252
253 const G4double G4NuclearRadii::r0[] = {
254 1.2,
255 1.3, 1.3, 1.3, 1.3,1.17,1.54,1.65,1.71, 1.7,1
256 1.7,1.57,1.53, 1.4, 1.3,1.30,1.44, 1.4, 1.4,
257 1.4, 1.4,1.46, 1.4, 1.4,1.46,1.55, 1.5,1.38,1
258 1.4, 1.4, 1.4,1.46, 1.4, 1.4, 1.4, 1.4, 1.4,1
259 1.4, 1.4, 1.4, 1.4, 1.4, 1.4,1.45,1.48, 1.4,1
260 1.46, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4,
261 1.4, 1.4, 1.4, 1.3, 1.3, 1.3, 1.3, 1.3, 1.3,
262 1.3, 1.3, 1.3, 1.3, 1.3, 1.3, 1.3, 1.3,1.33,1
263 1.3,1.32,1.34, 1.3, 1.3, 1.3, 1.3, 1.3, 1.3,
264 1.3, 1.3};
265