| 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 // G4NucleiProperties class implementation
27 //
28 // Author: V.Lara, October 1998
29 // History:
30 // - 17.11.1998, H.Kurashige - Migrated into particles category
31 // - 31.03.2009, T.Koi - Migrated to AME03
32 // --------------------------------------------------------------------
33
34 #include "G4NucleiProperties.hh"
35
36 #include "G4NucleiPropertiesTableAME12.hh"
37 #include "G4NucleiPropertiesTheoreticalTable.hh"
38 #include "G4ParticleTable.hh"
39 #include "G4PhysicalConstants.hh"
40 #include "G4SystemOfUnits.hh"
41
42 G4ThreadLocal G4double G4NucleiProperties::mass_proton = -1.;
43 G4ThreadLocal G4double G4NucleiProperties::mass_neutron = -1.;
44 G4ThreadLocal G4double G4NucleiProperties::mass_deuteron = -1.;
45 G4ThreadLocal G4double G4NucleiProperties::mass_triton = -1.;
46 G4ThreadLocal G4double G4NucleiProperties::mass_alpha = -1.;
47 G4ThreadLocal G4double G4NucleiProperties::mass_He3 = -1.;
48
49 G4double G4NucleiProperties::GetNuclearMass(const G4double A, const G4double Z)
50 {
51 G4double mass = 0.0;
52
53 if (std::fabs(A - G4int(A)) > 1.e-10) {
54 mass = NuclearMass(A, Z);
55 }
56 else {
57 // use mass table
58 auto iZ = G4int(Z);
59 auto iA = G4int(A);
60 mass = GetNuclearMass(iA, iZ);
61 }
62 return mass;
63 }
64
65 G4double G4NucleiProperties::GetNuclearMass(const G4int A, const G4int Z)
66 {
67 if (mass_proton <= 0.0) {
68 const G4ParticleDefinition* nucleus = nullptr;
69 nucleus = G4ParticleTable::GetParticleTable()->FindParticle("neutron");
70 if (nucleus != nullptr) mass_neutron = nucleus->GetPDGMass();
71
72 nucleus = G4ParticleTable::GetParticleTable()->FindParticle("deuteron");
73 if (nucleus != nullptr) mass_deuteron = nucleus->GetPDGMass();
74
75 nucleus = G4ParticleTable::GetParticleTable()->FindParticle("triton");
76 if (nucleus != nullptr) mass_triton = nucleus->GetPDGMass();
77
78 nucleus = G4ParticleTable::GetParticleTable()->FindParticle("alpha");
79 if (nucleus != nullptr) mass_alpha = nucleus->GetPDGMass();
80
81 nucleus = G4ParticleTable::GetParticleTable()->FindParticle("He3");
82 if (nucleus != nullptr) mass_He3 = nucleus->GetPDGMass();
83
84 nucleus = G4ParticleTable::GetParticleTable()->FindParticle("proton");
85 if (nucleus != nullptr) mass_proton = nucleus->GetPDGMass();
86 }
87
88 if (A < 1 || Z < 0 || Z > A) {
89 #ifdef G4VERBOSE
90 if (G4ParticleTable::GetParticleTable()->GetVerboseLevel() > 0) {
91 G4cout << "G4NucleiProperties::GetNuclearMass: Wrong values for A = " << A << " and Z = " << Z
92 << G4endl;
93 }
94 #endif
95 return 0.0;
96 }
97
98 G4double mass = -1.;
99 if ((Z <= 2)) {
100 // light nuclei
101 if ((Z == 1) && (A == 1)) {
102 mass = mass_proton;
103 }
104 else if ((Z == 0) && (A == 1)) {
105 mass = mass_neutron;
106 }
107 else if ((Z == 1) && (A == 2)) {
108 mass = mass_deuteron;
109 }
110 else if ((Z == 1) && (A == 3)) {
111 mass = mass_triton;
112 }
113 else if ((Z == 2) && (A == 4)) {
114 mass = mass_alpha;
115 }
116 else if ((Z == 2) && (A == 3)) {
117 mass = mass_He3;
118 }
119 }
120
121 if (mass < 0.) {
122 if (G4NucleiPropertiesTableAME12::IsInTable(Z, A)) {
123 // AME table
124 mass = G4NucleiPropertiesTableAME12::GetNuclearMass(Z, A);
125 }
126 else if (G4NucleiPropertiesTheoreticalTable::IsInTable(Z, A)) {
127 // Theoretical table
128 mass = G4NucleiPropertiesTheoreticalTable::GetNuclearMass(Z, A);
129 }
130 else if (Z == A) {
131 mass = A * mass_proton;
132 }
133 else if (0 == Z) {
134 mass = A * mass_neutron;
135 }
136 else {
137 mass = NuclearMass(G4double(A), G4double(Z));
138 }
139 }
140
141 if (mass < 0.) mass = 0.0;
142 return mass;
143 }
144
145 G4bool G4NucleiProperties::IsInStableTable(const G4double A, const G4double Z)
146 {
147 auto iA = G4int(A);
148 auto iZ = G4int(Z);
149 return IsInStableTable(iA, iZ);
150 }
151
152 G4bool G4NucleiProperties::IsInStableTable(const G4int A, const G4int Z)
153 {
154 if (A < 1 || Z < 0 || Z > A) {
155 #ifdef G4VERBOSE
156 if (G4ParticleTable::GetParticleTable()->GetVerboseLevel() > 0) {
157 G4cout << "G4NucleiProperties::IsInStableTable: Wrong values for A = " << A
158 << " and Z = " << Z << G4endl;
159 }
160 #endif
161 return false;
162 }
163
164 return G4NucleiPropertiesTableAME12::IsInTable(Z, A);
165 }
166
167 G4double G4NucleiProperties::GetMassExcess(const G4double A, const G4double Z)
168 {
169 auto iA = G4int(A);
170 auto iZ = G4int(Z);
171 return GetMassExcess(iA, iZ);
172 }
173
174 G4double G4NucleiProperties::GetMassExcess(const G4int A, const G4int Z)
175 {
176 if (A < 1 || Z < 0 || Z > A) {
177 #ifdef G4VERBOSE
178 if (G4ParticleTable::GetParticleTable()->GetVerboseLevel() > 0) {
179 G4cout << "G4NucleiProperties::GetMassExccess: Wrong values for A = " << A << " and Z = " << Z
180 << G4endl;
181 }
182 #endif
183 return 0.0;
184 }
185
186 if (G4NucleiPropertiesTableAME12::IsInTable(Z, A)) {
187 // AME table
188 return G4NucleiPropertiesTableAME12::GetMassExcess(Z, A);
189 }
190 if (G4NucleiPropertiesTheoreticalTable::IsInTable(Z, A)) {
191 return G4NucleiPropertiesTheoreticalTable::GetMassExcess(Z, A);
192 }
193 return MassExcess(A, Z);
194 }
195
196 G4double G4NucleiProperties::GetAtomicMass(const G4double A, const G4double Z)
197 {
198 if (A < 1 || Z < 0 || Z > A) {
199 #ifdef G4VERBOSE
200 if (G4ParticleTable::GetParticleTable()->GetVerboseLevel() > 0) {
201 G4cout << "G4NucleiProperties::GetAtomicMass: Wrong values for A = " << A << " and Z = " << Z
202 << G4endl;
203 }
204 #endif
205 return 0.0;
206 }
207 if (std::fabs(A - G4int(A)) > 1.e-10) {
208 return AtomicMass(A, Z);
209 }
210
211 auto iA = G4int(A);
212 auto iZ = G4int(Z);
213 if (G4NucleiPropertiesTableAME12::IsInTable(Z, A)) {
214 return G4NucleiPropertiesTableAME12::GetAtomicMass(Z, A);
215 }
216 if (G4NucleiPropertiesTheoreticalTable::IsInTable(iZ, iA)) {
217 return G4NucleiPropertiesTheoreticalTable::GetAtomicMass(iZ, iA);
218 }
219 return AtomicMass(A, Z);
220 }
221
222 G4double G4NucleiProperties::GetBindingEnergy(const G4double A, const G4double Z)
223 {
224 auto iA = G4int(A);
225 auto iZ = G4int(Z);
226 return GetBindingEnergy(iA, iZ);
227 }
228
229 G4double G4NucleiProperties::GetBindingEnergy(const G4int A, const G4int Z)
230 {
231 if (A < 1 || Z < 0 || Z > A) {
232 #ifdef G4VERBOSE
233 if (G4ParticleTable::GetParticleTable()->GetVerboseLevel() > 0) {
234 G4cout << "G4NucleiProperties::GetMassExccess: Wrong values for A = " << A << " and Z = " << Z
235 << G4endl;
236 }
237 #endif
238 return 0.0;
239 }
240
241 if (G4NucleiPropertiesTableAME12::IsInTable(Z, A)) {
242 return G4NucleiPropertiesTableAME12::GetBindingEnergy(Z, A);
243 }
244 if (G4NucleiPropertiesTheoreticalTable::IsInTable(Z, A)) {
245 return G4NucleiPropertiesTheoreticalTable::GetBindingEnergy(Z, A);
246 }
247 return BindingEnergy(A, Z);
248 }
249
250 G4double G4NucleiProperties::MassExcess(G4double A, G4double Z)
251 {
252 return GetAtomicMass(A, Z) - A * amu_c2;
253 }
254
255 G4double G4NucleiProperties::AtomicMass(G4double A, G4double Z)
256 {
257 G4double hydrogen_mass_excess;
258 G4double neutron_mass_excess;
259 hydrogen_mass_excess = G4NucleiPropertiesTableAME12::GetMassExcess(1, 1);
260 neutron_mass_excess = G4NucleiPropertiesTableAME12::GetMassExcess(0, 1);
261 G4double mass =
262 (A - Z) * neutron_mass_excess + Z * hydrogen_mass_excess - BindingEnergy(A, Z) + A * amu_c2;
263 return mass;
264 }
265
266 G4double G4NucleiProperties::NuclearMass(G4double A, G4double Z)
267 {
268 if (A < 1 || Z < 0 || Z > A) {
269 #ifdef G4VERBOSE
270 if (G4ParticleTable::GetParticleTable()->GetVerboseLevel() > 0) {
271 G4cout << "G4NucleiProperties::NuclearMass: Wrong values for A = " << A << " and Z = " << Z
272 << G4endl;
273 }
274 #endif
275 return 0.0;
276 }
277
278 G4double mass = AtomicMass(A, Z);
279
280 // atomic mass is converted to nuclear mass according to
281 // formula in AME03 and 12
282 //
283 mass -= Z * electron_mass_c2;
284 mass += (14.4381 * std::pow(Z, 2.39) + 1.55468 * 1e-6 * std::pow(Z, 5.35)) * eV;
285
286 return mass;
287 }
288
289 G4double G4NucleiProperties::BindingEnergy(G4double A, G4double Z)
290 {
291 //
292 // Weitzsaecker's Mass formula
293 //
294 G4int Npairing = G4int(A - Z) % 2; // pairing
295 G4int Zpairing = G4int(Z) % 2;
296 G4double binding = -15.67 * A // nuclear volume
297 + 17.23 * std::pow(A, 2. / 3.) // surface energy
298 + 93.15 * ((A / 2. - Z) * (A / 2. - Z)) / A // asymmetry
299 + 0.6984523 * Z * Z * std::pow(A, -1. / 3.); // coulomb
300 if (Npairing == Zpairing) {
301 binding += (Npairing + Zpairing - 1) * 12.0 / std::sqrt(A); // pairing
302 }
303
304 return -binding * MeV;
305 }
306