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 | // $Id: G4NucleiProperties.cc,v 1.13 2007/09/14 07:04:09 kurasige Exp $ |
---|
28 | // GEANT4 tag $Name: $ |
---|
29 | // |
---|
30 | // |
---|
31 | // ------------------------------------------------------------ |
---|
32 | // GEANT 4 class header file |
---|
33 | // |
---|
34 | // ------------------------------------------------------------ |
---|
35 | // |
---|
36 | // Hadronic Process: Nuclear De-excitations |
---|
37 | // by V. Lara (Oct 1998) |
---|
38 | // Migrate into particles category by H.Kurashige (17 Nov. 98) |
---|
39 | // Added Shell-Pairing corrections to the Cameron mass |
---|
40 | // excess formula by V.Lara (9 May 99) |
---|
41 | |
---|
42 | |
---|
43 | #include "G4NucleiProperties.hh" |
---|
44 | |
---|
45 | |
---|
46 | |
---|
47 | G4double G4NucleiProperties::AtomicMass(G4double A, G4double Z) |
---|
48 | { |
---|
49 | const G4double hydrogen_mass_excess = G4NucleiPropertiesTable::GetMassExcess(1,1); |
---|
50 | const G4double neutron_mass_excess = G4NucleiPropertiesTable::GetMassExcess(0,1); |
---|
51 | |
---|
52 | G4double mass = |
---|
53 | (A-Z)*neutron_mass_excess + Z*hydrogen_mass_excess - BindingEnergy(A,Z) + A*amu_c2; |
---|
54 | |
---|
55 | return mass; |
---|
56 | } |
---|
57 | |
---|
58 | G4double G4NucleiProperties::BindingEnergy(G4double A, G4double Z) |
---|
59 | { |
---|
60 | // |
---|
61 | // Weitzsaecker's Mass formula |
---|
62 | // |
---|
63 | G4int Npairing = G4int(A-Z)%2; // pairing |
---|
64 | G4int Zpairing = G4int(Z)%2; |
---|
65 | G4double binding = |
---|
66 | - 15.67*A // nuclear volume |
---|
67 | + 17.23*std::pow(A,2./3.) // surface energy |
---|
68 | + 93.15*((A/2.-Z)*(A/2.-Z))/A // asymmetry |
---|
69 | + 0.6984523*Z*Z*std::pow(A,-1./3.); // coulomb |
---|
70 | if( Npairing == Zpairing ) binding += (Npairing+Zpairing-1) * 12.0 / std::sqrt(A); // pairing |
---|
71 | |
---|
72 | return -binding*MeV; |
---|
73 | } |
---|
74 | |
---|
75 | |
---|
76 | G4double G4NucleiProperties::GetNuclearMass(const G4double A, const G4double Z) |
---|
77 | { |
---|
78 | if (A < 1 || Z < 0 || Z > A) { |
---|
79 | #ifdef G4VERBOSE |
---|
80 | if (G4ParticleTable::GetParticleTable()->GetVerboseLevel()>0) { |
---|
81 | G4cout << "G4NucleiProperties::GetNuclearMass: Wrong values for A = " << A |
---|
82 | << " and Z = " << Z << G4endl; |
---|
83 | } |
---|
84 | #endif |
---|
85 | return 0.0; |
---|
86 | } else { |
---|
87 | G4ParticleDefinition * nucleus = 0; |
---|
88 | if ( (Z<=2) ) { |
---|
89 | if ( (Z==1)&&(A==1) ) { |
---|
90 | nucleus = G4ParticleTable::GetParticleTable()->FindParticle("proton"); // proton |
---|
91 | } else if ( (Z==0)&&(A==1) ) { |
---|
92 | nucleus = G4ParticleTable::GetParticleTable()->FindParticle("neutron"); // neutron |
---|
93 | } else if ( (Z==1)&&(A==2) ) { |
---|
94 | nucleus = G4ParticleTable::GetParticleTable()->FindParticle("deuteron"); // deuteron |
---|
95 | } else if ( (Z==1)&&(A==3) ) { |
---|
96 | nucleus = G4ParticleTable::GetParticleTable()->FindParticle("triton"); // triton |
---|
97 | } else if ( (Z==2)&&(A==4) ) { |
---|
98 | nucleus = G4ParticleTable::GetParticleTable()->FindParticle("alpha"); // alpha |
---|
99 | } else if ( (Z==2)&&(A==3) ) { |
---|
100 | nucleus = G4ParticleTable::GetParticleTable()->FindParticle("He3"); // He3 |
---|
101 | } |
---|
102 | } |
---|
103 | if (nucleus!=0) { |
---|
104 | return nucleus->GetPDGMass(); |
---|
105 | }else { |
---|
106 | return GetAtomicMass(A,Z) - Z*electron_mass_c2 + 1.433e-5*MeV*std::pow(Z,2.39); |
---|
107 | } |
---|
108 | } |
---|
109 | } |
---|
110 | |
---|
111 | |
---|
112 | // G4double G4NucleiProperties::CameronMassExcess(const G4int A, const G4int Z) |
---|
113 | // { |
---|
114 | // const G4double alpha = -17.0354*MeV; |
---|
115 | // const G4double beta = -31.4506*MeV; |
---|
116 | // const G4double phi = 44.2355*MeV; |
---|
117 | // const G4double gamma = 25.8357*MeV; |
---|
118 | // |
---|
119 | // const G4double A13 = std::pow(G4double(A),1.0/3.0); |
---|
120 | // const G4double A23 = A13*A13; |
---|
121 | // const G4double A43 = A23*A23; |
---|
122 | // const G4double Z43 = std::pow(G4double(Z),4.0/3.0); |
---|
123 | // G4double D = (G4double(A) - 2.0*G4double(Z))/G4double(A); |
---|
124 | // D *= D; // D = std::pow((A-2Z)/A,2) |
---|
125 | // |
---|
126 | // |
---|
127 | // // Surface term |
---|
128 | // G4double SurfaceEnergy = (gamma - phi*D)*(1.0 - 0.62025/A23)*(1.0 - 0.62025/A23)*A23; |
---|
129 | // |
---|
130 | // // Coulomb term |
---|
131 | // G4double CoulombEnergy = 0.779*MeV*(G4double(Z*(Z-1))/A13)*(1.0-1.5849/A23+1.2273/A+1.5772/A43); |
---|
132 | // |
---|
133 | // // Exchange term |
---|
134 | // G4double ExchangeEnergy = -0.4323*MeV*(Z43/A13)*(1.0-0.57811/A13-0.14518/A23+0.49597/A); |
---|
135 | // |
---|
136 | // // Volume term |
---|
137 | // G4double VolumeEnergy = G4double(A)*(alpha-beta*D); |
---|
138 | // |
---|
139 | // // Shell+Pairing corrections for protons |
---|
140 | // G4double SPcorrectionZ; |
---|
141 | // if (Z <= TableSize) SPcorrectionZ = SPZTable[Z-1]*MeV; |
---|
142 | // else SPcorrectionZ = 0.0; |
---|
143 | // |
---|
144 | // // Shell+Pairing corrections for protons |
---|
145 | // G4int N = A - Z; |
---|
146 | // G4double SPcorrectionN; |
---|
147 | // if (N <= TableSize) SPcorrectionN = SPNTable[N-1]*MeV; |
---|
148 | // else SPcorrectionN = 0.0; |
---|
149 | // |
---|
150 | // |
---|
151 | // // Mass Excess |
---|
152 | // // First two terms give the mass excess of the neutrons and protons in the nucleus |
---|
153 | // // (see Cameron, Canadian Journal of Physics, 35, 1957 page 1022) |
---|
154 | // G4double MassExcess = 8.367*MeV*A - 0.783*MeV*Z + |
---|
155 | // SurfaceEnergy + CoulombEnergy + ExchangeEnergy + VolumeEnergy + |
---|
156 | // SPcorrectionZ + SPcorrectionN; |
---|
157 | // |
---|
158 | // return MassExcess; |
---|
159 | // } |
---|
160 | |
---|
161 | |
---|
162 | // S(Z)+P(Z) from Tab. 1 from A.G.W. Cameron, Canad. J. Phys., 35(1957)1021 |
---|
163 | // or Delta M(Z) from Tab. 97 of book [1] |
---|
164 | // const G4double G4NucleiProperties::SPZTable[TableSize] = { |
---|
165 | // 20.80, 15.80, 21.00, 16.80, 19.80, 16.50, 18.80, 16.50, 18.50, 17.20, // 1 - 10 |
---|
166 | // 18.26, 15.05, 16.01, 12.04, 13.27, 11.09, 12.17, 10.26, 11.04, 8.41, // 11 - 20 |
---|
167 | // 9.79, 7.36, 8.15, 5.63, 5.88, 3.17, 3.32, .82, 1.83, .97, // 21 - 30 |
---|
168 | // 2.33, 1.27, 2.92, 1.61, 2.91, 1.35, 2.40, .89, 1.74, .36, // 31 |
---|
169 | // 0.95, -0.65, -0.04, -1.73, -0.96, -2.87, -2.05, -4.05, -3.40, -5.72, // 41 |
---|
170 | // -3.75, -4.13, -2.42, -2.85, -1.01, -1.33, 0.54, -0.02, 1.74, 0.75, // 51 |
---|
171 | // 2.24, 1.00, 1.98, 0.79, 1.54, 0.39, 1.08, 0.00, 0.78, -0.35, // 61 |
---|
172 | // 0.58, -0.55, 0.59, -0.61, 0.59, -0.35, 0.32, -0.96, -0.52, -2.08, // 71 |
---|
173 | // -2.46, -3.64, -1.55, -0.96, 0.97, 0.88, 2.37, 1.75, 2.72, 1.90, // 81 |
---|
174 | // 2.55, 1.46, 1.93, 0.86, 1.17, 0.08, 0.39, -0.76, -0.39, -1.51, // 91 - 100 |
---|
175 | // -1.17, -2.36, -1.95, -3.06, -2.62, -3.55, -2.95, -3.75, -3.07, -3.79, // 101 - 110 |
---|
176 | // -3.06, -3.77, -3.05, -3.78, -3.12, -3.90, -3.35, -4.24, -3.86, -4.92, // 111 - 120 |
---|
177 | // -5.06, -6.77, -7.41, -9.18,-10.16,-11.12, -9.76, -9.23, -7.96, -7.65, // 121 - 130 |
---|
178 | // // --------- from this point there are not tabulated values ----------------------- |
---|
179 | // 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, // 131 - 140 |
---|
180 | // 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, // 141 - 150 |
---|
181 | // 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, // 151 |
---|
182 | // 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, // 161 |
---|
183 | // 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, // 171 |
---|
184 | // 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, // 181 |
---|
185 | // 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 // 191 - 200 |
---|
186 | // }; |
---|
187 | |
---|
188 | // S(N)+P(N) from Tab. 1 from A.G.W. Cameron, Canad. J. Phys., 35(1957)1021 |
---|
189 | // or Delta M(N) from Tab. 97 of book [1] |
---|
190 | // const G4double G4NucleiProperties::SPNTable[TableSize] = { |
---|
191 | // -8.40,-12.90, -8.00, 11.90, -9.20,-12.50,-10.80,-13.60,-11.20,-12.20, // 1 - 10 |
---|
192 | // -12.81,-15.40,-13.07,-15.80,-13.81,-14.98,-12.63,-13.76,-11.37,-12.38, // 11 - 20 |
---|
193 | // -9.23, -9.65, -7.64, -9.17, -8.05, -9.72, -8.87,-10.76, -8.64, -8.89, // 21 - 30 |
---|
194 | // -6.60, -7.13, -4.77, -5.33, -3.06, -3.79, -1.72, -2.79, -0.93, -2.19, // 31 |
---|
195 | // -0.52, -1.90, -0.45, -2.20, -1.22, -3.07, -2.42, -4.37, -3.94, -6.08, // 41 |
---|
196 | // -4.49, -4.50, -3.14, -2.93, -1.04, -1.36, 0.69, 0.21, 2.11, 1.33, // 51 |
---|
197 | // 3.29, 2.46, 4.30, 3.32, 4.79, 3.62, 4.97, 3.64, 4.63, 3.07, // 61 |
---|
198 | // 4.06, 2.49, 3.30, 1.46, 2.06, 0.51, 0.74, -1.18, -1.26, -3.54, // 71 |
---|
199 | // -3.97, -5.26, -4.18, -3.71, -2.10, -1.70, -0.08, -0.18, 0.94, 0.27, // 81 |
---|
200 | // 1.13, 0.08, 0.91, -0.31, 0.49, -0.78, 0.08, -1.15, -0.23, -1.41, // 91 - 100 |
---|
201 | // -0.42, -1.55, -0.55, -1.66, -0.66, -1.73, -0.75, -1.74, -0.78, -1.69, // 101 - 110 |
---|
202 | // -0.78, -1.60, -0.75, -1.46, -0.67, -1.26, -0.51, -1.04, -0.53, -1.84, // 111 - 120 |
---|
203 | // -2.42, -4.52, -4.76, -6.33, -6.76, -7.81, -5.80, -5.37, -3.63, -3.35, // 121 - 130 |
---|
204 | // -1.75, -1.88, -0.61, -0.90, 0.09, -0.32, 0.55, -0.13, 0.70, -0.06, // 131 - 140 |
---|
205 | // 0.49, -0.20, 0.40, -0.22, 0.36, -0.09, 0.58, 0.12, 0.75, 0.15, // 141 - 150 |
---|
206 | // 0.70, 0.17, 1.11, 0.89, 1.85, 1.62, 2.54, 2.29, 3.20, 2.91, // 151 |
---|
207 | // 3.84, 3.53, 4.48, 4.15, 5.12, 4.78, 5.75, 5.39, 6.31, 5.91, // 161 |
---|
208 | // 6.87, 6.33, 7.13, 6.61, 7.30, 6.31, 6.27, 4.83, 4.49, 2.85, // 171 |
---|
209 | // 2.32, 0.58, -0.11, -0.98, 0.81, 1.77, 3.37, 4.13, 5.60, 6.15, // 181 |
---|
210 | // 7.29, 7.35, 7.95, 7.67, 8.16, 7.83, 8.31, 8.01, 8.53, 8.27 // 191 - 200 |
---|
211 | // }; |
---|
212 | |
---|
213 | |
---|
214 | |
---|
215 | G4double G4NucleiProperties::GetMassExcess(const G4int A, const G4int Z) |
---|
216 | { |
---|
217 | if (A < 1 || Z < 0 || Z > A) { |
---|
218 | #ifdef G4VERBOSE |
---|
219 | if (G4ParticleTable::GetParticleTable()->GetVerboseLevel()>0) { |
---|
220 | G4cout << "G4NucleiProperties::GetMassExccess: Wrong values for A = " |
---|
221 | << A << " and Z = " << Z << G4endl; |
---|
222 | } |
---|
223 | #endif |
---|
224 | return 0.0; |
---|
225 | |
---|
226 | } else { |
---|
227 | |
---|
228 | if (G4NucleiPropertiesTable::IsInTable(Z,A)){ |
---|
229 | return G4NucleiPropertiesTable::GetMassExcess(Z,A); |
---|
230 | } else if (G4NucleiPropertiesTheoreticalTable::IsInTable(Z,A)){ |
---|
231 | return G4NucleiPropertiesTheoreticalTable::GetMassExcess(Z,A); |
---|
232 | } else { |
---|
233 | return MassExcess(A,Z); |
---|
234 | } |
---|
235 | } |
---|
236 | |
---|
237 | } |
---|
238 | |
---|
239 | |
---|
240 | G4double G4NucleiProperties::GetAtomicMass(const G4double A, const G4double Z) |
---|
241 | { |
---|
242 | if (Z < 0 || Z > A) { |
---|
243 | #ifdef G4VERBOSE |
---|
244 | if (G4ParticleTable::GetParticleTable()->GetVerboseLevel()>0) { |
---|
245 | G4cout << "G4NucleiProperties::GetAtomicMass: Wrong values for A = " |
---|
246 | << A << " and Z = " << Z << G4endl; |
---|
247 | } |
---|
248 | #endif |
---|
249 | return 0.0; |
---|
250 | |
---|
251 | } else if (std::abs(A - G4int(A)) > 1.e-10) { |
---|
252 | return AtomicMass(A,Z); |
---|
253 | |
---|
254 | } else { |
---|
255 | G4int iA = G4int(A); |
---|
256 | G4int iZ = G4int(Z); |
---|
257 | if (G4NucleiPropertiesTable::IsInTable(iZ,iA)) { |
---|
258 | return G4NucleiPropertiesTable::GetAtomicMass(iZ,iA); |
---|
259 | } else if (G4NucleiPropertiesTheoreticalTable::IsInTable(iZ,iA)){ |
---|
260 | return G4NucleiPropertiesTheoreticalTable::GetAtomicMass(iZ,iA); |
---|
261 | } else { |
---|
262 | return AtomicMass(A,Z); |
---|
263 | } |
---|
264 | } |
---|
265 | } |
---|
266 | |
---|
267 | G4double G4NucleiProperties::GetBindingEnergy(const G4int A, const G4int Z) |
---|
268 | { |
---|
269 | if (A < 1 || Z < 0 || Z > A) { |
---|
270 | #ifdef G4VERBOSE |
---|
271 | if (G4ParticleTable::GetParticleTable()->GetVerboseLevel()>0) { |
---|
272 | G4cout << "G4NucleiProperties::GetMassExccess: Wrong values for A = " |
---|
273 | << A << " and Z = " << Z << G4endl; |
---|
274 | } |
---|
275 | #endif |
---|
276 | return 0.0; |
---|
277 | |
---|
278 | } else { |
---|
279 | if (G4NucleiPropertiesTable::IsInTable(Z,A)) { |
---|
280 | return G4NucleiPropertiesTable::GetBindingEnergy(Z,A); |
---|
281 | } else if (G4NucleiPropertiesTheoreticalTable::IsInTable(Z,A)) { |
---|
282 | return G4NucleiPropertiesTheoreticalTable::GetBindingEnergy(Z,A); |
---|
283 | }else { |
---|
284 | return BindingEnergy(A,Z); |
---|
285 | } |
---|
286 | |
---|
287 | } |
---|
288 | } |
---|
289 | |
---|
290 | |
---|
291 | G4double G4NucleiProperties::MassExcess(G4double A, G4double Z) |
---|
292 | { |
---|
293 | return GetAtomicMass(A,Z) - A*amu_c2; |
---|
294 | } |
---|
295 | |
---|