| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // G4NucleiProperties class implementation 26 // G4NucleiProperties class implementation
27 // 27 //
28 // Author: V.Lara, October 1998 28 // Author: V.Lara, October 1998
29 // History: 29 // History:
30 // - 17.11.1998, H.Kurashige - Migrated into p 30 // - 17.11.1998, H.Kurashige - Migrated into particles category
31 // - 31.03.2009, T.Koi - Migrated to AME03 31 // - 31.03.2009, T.Koi - Migrated to AME03
32 // ------------------------------------------- 32 // --------------------------------------------------------------------
33 33
34 #include "G4NucleiProperties.hh" 34 #include "G4NucleiProperties.hh"
35 35
36 #include "G4NucleiPropertiesTableAME12.hh" 36 #include "G4NucleiPropertiesTableAME12.hh"
37 #include "G4NucleiPropertiesTheoreticalTable.h 37 #include "G4NucleiPropertiesTheoreticalTable.hh"
38 #include "G4ParticleTable.hh" 38 #include "G4ParticleTable.hh"
>> 39
39 #include "G4PhysicalConstants.hh" 40 #include "G4PhysicalConstants.hh"
40 #include "G4SystemOfUnits.hh" 41 #include "G4SystemOfUnits.hh"
41 42
42 G4ThreadLocal G4double G4NucleiProperties::mas 43 G4ThreadLocal G4double G4NucleiProperties::mass_proton = -1.;
43 G4ThreadLocal G4double G4NucleiProperties::mas 44 G4ThreadLocal G4double G4NucleiProperties::mass_neutron = -1.;
44 G4ThreadLocal G4double G4NucleiProperties::mas 45 G4ThreadLocal G4double G4NucleiProperties::mass_deuteron = -1.;
45 G4ThreadLocal G4double G4NucleiProperties::mas 46 G4ThreadLocal G4double G4NucleiProperties::mass_triton = -1.;
46 G4ThreadLocal G4double G4NucleiProperties::mas 47 G4ThreadLocal G4double G4NucleiProperties::mass_alpha = -1.;
47 G4ThreadLocal G4double G4NucleiProperties::mas 48 G4ThreadLocal G4double G4NucleiProperties::mass_He3 = -1.;
48 49
49 G4double G4NucleiProperties::GetNuclearMass(co 50 G4double G4NucleiProperties::GetNuclearMass(const G4double A, const G4double Z)
50 { 51 {
51 G4double mass = 0.0; << 52 G4double mass =0.0;
52 53
53 if (std::fabs(A - G4int(A)) > 1.e-10) { << 54 if (std::fabs(A - G4int(A)) > 1.e-10)
54 mass = NuclearMass(A, Z); << 55 {
>> 56 mass = NuclearMass(A,Z);
>> 57
55 } 58 }
56 else { << 59 else
>> 60 {
57 // use mass table 61 // use mass table
58 auto iZ = G4int(Z); << 62 G4int iZ = G4int(Z);
59 auto iA = G4int(A); << 63 G4int iA = G4int(A);
60 mass = GetNuclearMass(iA, iZ); << 64 mass =GetNuclearMass(iA,iZ);
61 } 65 }
62 return mass; 66 return mass;
63 } 67 }
64 68
>> 69
65 G4double G4NucleiProperties::GetNuclearMass(co 70 G4double G4NucleiProperties::GetNuclearMass(const G4int A, const G4int Z)
66 { 71 {
67 if (mass_proton <= 0.0) { << 72 if (mass_proton <= 0.0 )
68 const G4ParticleDefinition* nucleus = null << 73 {
>> 74 const G4ParticleDefinition * nucleus = nullptr;
69 nucleus = G4ParticleTable::GetParticleTabl 75 nucleus = G4ParticleTable::GetParticleTable()->FindParticle("neutron");
70 if (nucleus != nullptr) mass_neutron = nuc << 76 if (nucleus!=nullptr) mass_neutron = nucleus->GetPDGMass();
71 77
72 nucleus = G4ParticleTable::GetParticleTabl 78 nucleus = G4ParticleTable::GetParticleTable()->FindParticle("deuteron");
73 if (nucleus != nullptr) mass_deuteron = nu << 79 if (nucleus!=nullptr) mass_deuteron = nucleus->GetPDGMass();
74 80
75 nucleus = G4ParticleTable::GetParticleTabl 81 nucleus = G4ParticleTable::GetParticleTable()->FindParticle("triton");
76 if (nucleus != nullptr) mass_triton = nucl << 82 if (nucleus!=nullptr) mass_triton = nucleus->GetPDGMass();
77 83
78 nucleus = G4ParticleTable::GetParticleTabl 84 nucleus = G4ParticleTable::GetParticleTable()->FindParticle("alpha");
79 if (nucleus != nullptr) mass_alpha = nucle << 85 if (nucleus!=nullptr) mass_alpha = nucleus->GetPDGMass();
80 86
81 nucleus = G4ParticleTable::GetParticleTabl 87 nucleus = G4ParticleTable::GetParticleTable()->FindParticle("He3");
82 if (nucleus != nullptr) mass_He3 = nucleus << 88 if (nucleus!=nullptr) mass_He3 = nucleus->GetPDGMass();
83 89
84 nucleus = G4ParticleTable::GetParticleTabl 90 nucleus = G4ParticleTable::GetParticleTable()->FindParticle("proton");
85 if (nucleus != nullptr) mass_proton = nucl << 91 if (nucleus!=nullptr) mass_proton = nucleus->GetPDGMass();
86 } 92 }
87 93
88 if (A < 1 || Z < 0 || Z > A) { << 94 if (A < 1 || Z < 0 || Z > A)
>> 95 {
89 #ifdef G4VERBOSE 96 #ifdef G4VERBOSE
90 if (G4ParticleTable::GetParticleTable()->G << 97 if (G4ParticleTable::GetParticleTable()->GetVerboseLevel()>0)
91 G4cout << "G4NucleiProperties::GetNuclea << 98 {
92 << G4endl; << 99 G4cout << "G4NucleiProperties::GetNuclearMass: Wrong values for A = "
>> 100 << A << " and Z = " << Z << G4endl;
93 } 101 }
94 #endif << 102 #endif
95 return 0.0; 103 return 0.0;
96 } 104 }
97 << 105
98 G4double mass = -1.; << 106 G4double mass= -1.;
99 if ((Z <= 2)) { << 107 if ( (Z<=2) )
>> 108 {
100 // light nuclei 109 // light nuclei
101 if ((Z == 1) && (A == 1)) { << 110 if ( (Z==1)&&(A==1) ) {
102 mass = mass_proton; 111 mass = mass_proton;
103 } << 112 } else if ( (Z==0)&&(A==1) ) {
104 else if ((Z == 0) && (A == 1)) { <<
105 mass = mass_neutron; 113 mass = mass_neutron;
106 } << 114 } else if ( (Z==1)&&(A==2) ) {
107 else if ((Z == 1) && (A == 2)) { <<
108 mass = mass_deuteron; 115 mass = mass_deuteron;
109 } << 116 } else if ( (Z==1)&&(A==3) ) {
110 else if ((Z == 1) && (A == 3)) { <<
111 mass = mass_triton; 117 mass = mass_triton;
112 } << 118 } else if ( (Z==2)&&(A==4) ) {
113 else if ((Z == 2) && (A == 4)) { <<
114 mass = mass_alpha; 119 mass = mass_alpha;
115 } << 120 } else if ( (Z==2)&&(A==3) ) {
116 else if ((Z == 2) && (A == 3)) { <<
117 mass = mass_He3; 121 mass = mass_He3;
118 } 122 }
119 } 123 }
120 << 124
121 if (mass < 0.) { << 125 if (mass < 0.)
122 if (G4NucleiPropertiesTableAME12::IsInTabl << 126 {
>> 127 if ( G4NucleiPropertiesTableAME12::IsInTable(Z,A) ) {
123 // AME table 128 // AME table
124 mass = G4NucleiPropertiesTableAME12::Get << 129 mass = G4NucleiPropertiesTableAME12::GetNuclearMass(Z,A);
125 } << 130 } else if (G4NucleiPropertiesTheoreticalTable::IsInTable(Z,A)){
126 else if (G4NucleiPropertiesTheoreticalTabl <<
127 // Theoretical table 131 // Theoretical table
128 mass = G4NucleiPropertiesTheoreticalTabl << 132 mass = G4NucleiPropertiesTheoreticalTable::GetNuclearMass(Z,A);
129 } << 133 } else if ( Z == A ) {
130 else if (Z == A) { << 134 mass = A*mass_proton;
131 mass = A * mass_proton; << 135 } else if( 0 == Z ) {
132 } << 136 mass = A*mass_neutron;
133 else if (0 == Z) { << 137 } else {
134 mass = A * mass_neutron; << 138 mass = NuclearMass(G4double(A),G4double(Z));
135 } <<
136 else { <<
137 mass = NuclearMass(G4double(A), G4double <<
138 } 139 }
139 } 140 }
140 141
141 if (mass < 0.) mass = 0.0; 142 if (mass < 0.) mass = 0.0;
142 return mass; 143 return mass;
143 } 144 }
144 145
145 G4bool G4NucleiProperties::IsInStableTable(con 146 G4bool G4NucleiProperties::IsInStableTable(const G4double A, const G4double Z)
146 { 147 {
147 auto iA = G4int(A); << 148 G4int iA = G4int(A);
148 auto iZ = G4int(Z); << 149 G4int iZ = G4int(Z);
149 return IsInStableTable(iA, iZ); 150 return IsInStableTable(iA, iZ);
150 } 151 }
151 152
152 G4bool G4NucleiProperties::IsInStableTable(con 153 G4bool G4NucleiProperties::IsInStableTable(const G4int A, const G4int Z)
153 { 154 {
154 if (A < 1 || Z < 0 || Z > A) { << 155 if (A < 1 || Z < 0 || Z > A)
>> 156 {
155 #ifdef G4VERBOSE 157 #ifdef G4VERBOSE
156 if (G4ParticleTable::GetParticleTable()->G << 158 if (G4ParticleTable::GetParticleTable()->GetVerboseLevel()>0)
157 G4cout << "G4NucleiProperties::IsInStabl << 159 {
158 << " and Z = " << Z << G4endl; << 160 G4cout << "G4NucleiProperties::IsInStableTable: Wrong values for A = "
>> 161 << A << " and Z = " << Z << G4endl;
159 } 162 }
160 #endif << 163 #endif
161 return false; 164 return false;
162 } << 165 }
163 166
164 return G4NucleiPropertiesTableAME12::IsInTab << 167 return G4NucleiPropertiesTableAME12::IsInTable(Z,A);
165 } 168 }
166 169
167 G4double G4NucleiProperties::GetMassExcess(con 170 G4double G4NucleiProperties::GetMassExcess(const G4double A, const G4double Z)
168 { 171 {
169 auto iA = G4int(A); << 172 G4int iA = G4int(A);
170 auto iZ = G4int(Z); << 173 G4int iZ = G4int(Z);
171 return GetMassExcess(iA, iZ); << 174 return GetMassExcess(iA,iZ);
172 } 175 }
173 176
174 G4double G4NucleiProperties::GetMassExcess(con 177 G4double G4NucleiProperties::GetMassExcess(const G4int A, const G4int Z)
175 { 178 {
176 if (A < 1 || Z < 0 || Z > A) { << 179 if (A < 1 || Z < 0 || Z > A)
>> 180 {
177 #ifdef G4VERBOSE 181 #ifdef G4VERBOSE
178 if (G4ParticleTable::GetParticleTable()->G << 182 if (G4ParticleTable::GetParticleTable()->GetVerboseLevel()>0)
179 G4cout << "G4NucleiProperties::GetMassEx << 183 {
180 << G4endl; << 184 G4cout << "G4NucleiProperties::GetMassExccess: Wrong values for A = "
>> 185 << A << " and Z = " << Z << G4endl;
181 } 186 }
182 #endif << 187 #endif
183 return 0.0; 188 return 0.0;
>> 189
184 } 190 }
>> 191 else
>> 192 {
185 193
186 if (G4NucleiPropertiesTableAME12::IsInTable( << 194 if ( G4NucleiPropertiesTableAME12::IsInTable(Z,A) ){
187 // AME table << 195 // AME table
188 return G4NucleiPropertiesTableAME12::GetMa << 196 return G4NucleiPropertiesTableAME12::GetMassExcess(Z,A);
189 } << 197 } else if (G4NucleiPropertiesTheoreticalTable::IsInTable(Z,A)){
190 if (G4NucleiPropertiesTheoreticalTable::IsIn << 198 return G4NucleiPropertiesTheoreticalTable::GetMassExcess(Z,A);
191 return G4NucleiPropertiesTheoreticalTable: << 199 } else {
>> 200 return MassExcess(A,Z);
>> 201 }
192 } 202 }
193 return MassExcess(A, Z); <<
194 } 203 }
195 204
>> 205
196 G4double G4NucleiProperties::GetAtomicMass(con 206 G4double G4NucleiProperties::GetAtomicMass(const G4double A, const G4double Z)
197 { 207 {
198 if (A < 1 || Z < 0 || Z > A) { << 208 if (A < 1 || Z < 0 || Z > A)
>> 209 {
199 #ifdef G4VERBOSE 210 #ifdef G4VERBOSE
200 if (G4ParticleTable::GetParticleTable()->G << 211 if (G4ParticleTable::GetParticleTable()->GetVerboseLevel()>0)
201 G4cout << "G4NucleiProperties::GetAtomic << 212 {
202 << G4endl; << 213 G4cout << "G4NucleiProperties::GetAtomicMass: Wrong values for A = "
>> 214 << A << " and Z = " << Z << G4endl;
203 } 215 }
204 #endif << 216 #endif
205 return 0.0; 217 return 0.0;
206 } <<
207 if (std::fabs(A - G4int(A)) > 1.e-10) { <<
208 return AtomicMass(A, Z); <<
209 } <<
210 218
211 auto iA = G4int(A); <<
212 auto iZ = G4int(Z); <<
213 if (G4NucleiPropertiesTableAME12::IsInTable( <<
214 return G4NucleiPropertiesTableAME12::GetAt <<
215 } 219 }
216 if (G4NucleiPropertiesTheoreticalTable::IsIn << 220 else if (std::fabs(A - G4int(A)) > 1.e-10)
217 return G4NucleiPropertiesTheoreticalTable: << 221 {
>> 222 return AtomicMass(A,Z);
>> 223 }
>> 224 else
>> 225 {
>> 226 G4int iA = G4int(A);
>> 227 G4int iZ = G4int(Z);
>> 228 if ( G4NucleiPropertiesTableAME12::IsInTable(Z,A) ) {
>> 229 return G4NucleiPropertiesTableAME12::GetAtomicMass(Z,A);
>> 230 } else if (G4NucleiPropertiesTheoreticalTable::IsInTable(iZ,iA)){
>> 231 return G4NucleiPropertiesTheoreticalTable::GetAtomicMass(iZ,iA);
>> 232 } else {
>> 233 return AtomicMass(A,Z);
>> 234 }
218 } 235 }
219 return AtomicMass(A, Z); <<
220 } 236 }
221 237
222 G4double G4NucleiProperties::GetBindingEnergy( << 238 G4double
223 { << 239 G4NucleiProperties::GetBindingEnergy(const G4double A, const G4double Z)
224 auto iA = G4int(A); << 240 {
225 auto iZ = G4int(Z); << 241 G4int iA = G4int(A);
226 return GetBindingEnergy(iA, iZ); << 242 G4int iZ = G4int(Z);
>> 243 return GetBindingEnergy(iA,iZ);
227 } 244 }
228 245
229 G4double G4NucleiProperties::GetBindingEnergy( 246 G4double G4NucleiProperties::GetBindingEnergy(const G4int A, const G4int Z)
230 { 247 {
231 if (A < 1 || Z < 0 || Z > A) { << 248 if (A < 1 || Z < 0 || Z > A)
>> 249 {
232 #ifdef G4VERBOSE 250 #ifdef G4VERBOSE
233 if (G4ParticleTable::GetParticleTable()->G << 251 if (G4ParticleTable::GetParticleTable()->GetVerboseLevel()>0)
234 G4cout << "G4NucleiProperties::GetMassEx << 252 {
235 << G4endl; << 253 G4cout << "G4NucleiProperties::GetMassExccess: Wrong values for A = "
>> 254 << A << " and Z = " << Z << G4endl;
236 } 255 }
237 #endif 256 #endif
238 return 0.0; 257 return 0.0;
239 } <<
240 258
241 if (G4NucleiPropertiesTableAME12::IsInTable( <<
242 return G4NucleiPropertiesTableAME12::GetBi <<
243 } 259 }
244 if (G4NucleiPropertiesTheoreticalTable::IsIn << 260 else
245 return G4NucleiPropertiesTheoreticalTable: << 261 {
>> 262 if ( G4NucleiPropertiesTableAME12::IsInTable(Z,A) ) {
>> 263 return G4NucleiPropertiesTableAME12::GetBindingEnergy(Z,A);
>> 264 } else if (G4NucleiPropertiesTheoreticalTable::IsInTable(Z,A)) {
>> 265 return G4NucleiPropertiesTheoreticalTable::GetBindingEnergy(Z,A);
>> 266 }else {
>> 267 return BindingEnergy(A,Z);
>> 268 }
246 } 269 }
247 return BindingEnergy(A, Z); <<
248 } 270 }
249 271
250 G4double G4NucleiProperties::MassExcess(G4doub << 272 G4double G4NucleiProperties::MassExcess(G4double A, G4double Z)
251 { 273 {
252 return GetAtomicMass(A, Z) - A * amu_c2; << 274 return GetAtomicMass(A,Z) - A*amu_c2;
253 } 275 }
254 276
255 G4double G4NucleiProperties::AtomicMass(G4doub << 277 G4double G4NucleiProperties::AtomicMass(G4double A, G4double Z)
256 { 278 {
257 G4double hydrogen_mass_excess; 279 G4double hydrogen_mass_excess;
258 G4double neutron_mass_excess; << 280 G4double neutron_mass_excess;
259 hydrogen_mass_excess = G4NucleiPropertiesTab << 281 hydrogen_mass_excess = G4NucleiPropertiesTableAME12::GetMassExcess(1,1);
260 neutron_mass_excess = G4NucleiPropertiesTabl << 282 neutron_mass_excess = G4NucleiPropertiesTableAME12::GetMassExcess(0,1);
261 G4double mass = << 283 G4double mass = (A-Z)*neutron_mass_excess
262 (A - Z) * neutron_mass_excess + Z * hydrog << 284 + Z*hydrogen_mass_excess - BindingEnergy(A,Z) + A*amu_c2;
263 return mass; 285 return mass;
264 } 286 }
265 287
266 G4double G4NucleiProperties::NuclearMass(G4dou << 288 G4double G4NucleiProperties::NuclearMass(G4double A, G4double Z)
267 { 289 {
268 if (A < 1 || Z < 0 || Z > A) { << 290 if (A < 1 || Z < 0 || Z > A)
>> 291 {
269 #ifdef G4VERBOSE 292 #ifdef G4VERBOSE
270 if (G4ParticleTable::GetParticleTable()->G << 293 if (G4ParticleTable::GetParticleTable()->GetVerboseLevel()>0)
271 G4cout << "G4NucleiProperties::NuclearMa << 294 {
272 << G4endl; << 295 G4cout << "G4NucleiProperties::NuclearMass: Wrong values for A = "
>> 296 << A << " and Z = " << Z << G4endl;
273 } 297 }
274 #endif << 298 #endif
275 return 0.0; 299 return 0.0;
276 } 300 }
277 301
278 G4double mass = AtomicMass(A, Z); << 302 G4double mass = AtomicMass(A,Z);
279 303
280 // atomic mass is converted to nuclear mass 304 // atomic mass is converted to nuclear mass according to
281 // formula in AME03 and 12 305 // formula in AME03 and 12
282 // 306 //
283 mass -= Z * electron_mass_c2; << 307 mass -= Z*electron_mass_c2;
284 mass += (14.4381 * std::pow(Z, 2.39) + 1.554 << 308 mass += ( 14.4381*std::pow ( Z , 2.39 )
>> 309 + 1.55468*1e-6*std::pow ( Z , 5.35 ) )*eV;
285 310
286 return mass; 311 return mass;
287 } 312 }
288 313
289 G4double G4NucleiProperties::BindingEnergy(G4d 314 G4double G4NucleiProperties::BindingEnergy(G4double A, G4double Z)
290 { << 315 {
291 // 316 //
292 // Weitzsaecker's Mass formula 317 // Weitzsaecker's Mass formula
293 // 318 //
294 G4int Npairing = G4int(A - Z) % 2; // pairi << 319 G4int Npairing = G4int(A-Z)%2; // pairing
295 G4int Zpairing = G4int(Z) % 2; << 320 G4int Zpairing = G4int(Z)%2;
296 G4double binding = -15.67 * A // nuclear vo << 321 G4double binding =
297 + 17.23 * std::pow(A, 2. << 322 - 15.67*A // nuclear volume
298 + 93.15 * ((A / 2. - Z) * << 323 + 17.23*std::pow(A,2./3.) // surface energy
299 + 0.6984523 * Z * Z * std << 324 + 93.15*((A/2.-Z)*(A/2.-Z))/A // asymmetry
300 if (Npairing == Zpairing) { << 325 + 0.6984523*Z*Z*std::pow(A,-1./3.); // coulomb
301 binding += (Npairing + Zpairing - 1) * 12. << 326 if( Npairing == Zpairing )
>> 327 {
>> 328 binding += (Npairing+Zpairing-1) * 12.0 / std::sqrt(A); // pairing
302 } 329 }
303 330
304 return -binding * MeV; << 331 return -binding*MeV;
305 } 332 }
306 333