Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4NucleiProperties.cc
Go to the documentation of this file.
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 99159 2016-09-07 08:11:50Z gcosmo $
28 //
29 //
30 // ------------------------------------------------------------
31 // GEANT 4 class header file
32 //
33 // ------------------------------------------------------------
34 //
35 // Hadronic Process: Nuclear De-excitations
36 // by V. Lara (Oct 1998)
37 // Migrate into particles category by H.Kurashige (17 Nov. 98)
38 // Added Shell-Pairing corrections to the Cameron mass
39 // excess formula by V.Lara (9 May 99)
40 // 090331 Migrate to AME03 by Koi, Tatsumi
41 
42 #include "G4NucleiProperties.hh"
43 
47 #include "G4ParticleTable.hh"
48 
49 #include "G4PhysicalConstants.hh"
50 #include "G4SystemOfUnits.hh"
51 
52 G4ThreadLocal G4double G4NucleiProperties::mass_proton = -1.;
53 G4ThreadLocal G4double G4NucleiProperties::mass_neutron = -1.;
54 G4ThreadLocal G4double G4NucleiProperties::mass_deuteron = -1.;
55 G4ThreadLocal G4double G4NucleiProperties::mass_triton = -1.;
56 G4ThreadLocal G4double G4NucleiProperties::mass_alpha = -1.;
57 G4ThreadLocal G4double G4NucleiProperties::mass_He3 = -1.;
58 #ifndef G4NucleiProperties_USE_OLD_AME_TABLE
59 G4bool G4NucleiProperties::use_old_evaluation = false;
60 #else
61 G4bool G4NucleiProperties::use_old_evaluation = true;
62 #endif
63 
65 {
66  G4double mass=0.0;
67 
68  if (std::fabs(A - G4int(A)) > 1.e-10) {
69  mass = NuclearMass(A,Z);
70 
71  } else {
72  // use mass table
73  G4int iZ = G4int(Z);
74  G4int iA = G4int(A);
75  mass =GetNuclearMass(iA,iZ);
76  }
77 
78  return mass;
79 }
80 
81 
83 {
84  if (mass_proton <= 0.0 ) {
85  const G4ParticleDefinition * nucleus = 0;
86  nucleus = G4ParticleTable::GetParticleTable()->FindParticle("proton"); // proton
87  if (nucleus!=0) mass_proton = nucleus->GetPDGMass();
88  nucleus = G4ParticleTable::GetParticleTable()->FindParticle("neutron"); // neutron
89  if (nucleus!=0) mass_neutron = nucleus->GetPDGMass();
90  nucleus = G4ParticleTable::GetParticleTable()->FindParticle("deuteron"); // deuteron
91  if (nucleus!=0) mass_deuteron = nucleus->GetPDGMass();
92  nucleus = G4ParticleTable::GetParticleTable()->FindParticle("triton"); // triton
93  if (nucleus!=0) mass_triton = nucleus->GetPDGMass();
94  nucleus = G4ParticleTable::GetParticleTable()->FindParticle("alpha"); // alpha
95  if (nucleus!=0) mass_alpha = nucleus->GetPDGMass();
96  nucleus = G4ParticleTable::GetParticleTable()->FindParticle("He3"); // He3
97  if (nucleus!=0) mass_He3 = nucleus->GetPDGMass();
98 
99  }
100 
101  if (A < 1 || Z < 0 || Z > A) {
102 #ifdef G4VERBOSE
103  if (G4ParticleTable::GetParticleTable()->GetVerboseLevel()>0) {
104  G4cout << "G4NucleiProperties::GetNuclearMass: Wrong values for A = " << A
105  << " and Z = " << Z << G4endl;
106  }
107 #endif
108  return 0.0;
109  }
110 
111  G4double mass= -1.;
112  if ( (Z<=2) ) {
113  // light nuclei
114  if ( (Z==1)&&(A==1) ) {
115  mass = mass_proton;
116  } else if ( (Z==0)&&(A==1) ) {
117  mass = mass_neutron;
118  } else if ( (Z==1)&&(A==2) ) {
119  mass = mass_deuteron;
120  } else if ( (Z==1)&&(A==3) ) {
121  mass = mass_triton;
122  } else if ( (Z==2)&&(A==4) ) {
123  mass = mass_alpha;
124  } else if ( (Z==2)&&(A==3) ) {
125  mass = mass_He3;
126  }
127  }
128 
129  if (mass < 0.) {
130  G4bool inAMETable = false;
131  if ( ! use_old_evaluation ) {
132  inAMETable = G4NucleiPropertiesTableAME12::IsInTable(Z,A);
133  } else {
134  inAMETable = G4NucleiPropertiesTableAME03::IsInTable(Z,A);
135  }
136  if ( inAMETable ) {
137  // AME table
138  if ( ! use_old_evaluation ) {
139  mass = G4NucleiPropertiesTableAME12::GetNuclearMass(Z,A);
140  } else {
141  mass = G4NucleiPropertiesTableAME03::GetNuclearMass(Z,A);
142  }
143  } else if (G4NucleiPropertiesTheoreticalTable::IsInTable(Z,A)){
144  // Theoretical table
145  mass = G4NucleiPropertiesTheoreticalTable::GetNuclearMass(Z,A);
146  } else {
147  mass = NuclearMass(G4double(A),G4double(Z));
148  }
149  }
150 
151  if (mass < 0.) mass = 0.0;
152  return mass;
153 }
154 
156 {
157  G4int iA = G4int(A);
158  G4int iZ = G4int(Z);
159  return IsInStableTable(iA, iZ);
160 }
161 
163 {
164  if (A < 1 || Z < 0 || Z > A) {
165 #ifdef G4VERBOSE
166  if (G4ParticleTable::GetParticleTable()->GetVerboseLevel()>0) {
167  G4cout << "G4NucleiProperties::IsInStableTable: Wrong values for A = "
168  << A << " and Z = " << Z << G4endl;
169  }
170 #endif
171  return false;
172  }
173 
174  if ( ! use_old_evaluation ) {
175  return G4NucleiPropertiesTableAME12::IsInTable(Z,A);
176  } else {
177  return G4NucleiPropertiesTableAME03::IsInTable(Z,A);
178  }
179 }
180 
182 {
183  G4int iA = G4int(A);
184  G4int iZ = G4int(Z);
185  return GetMassExcess(iA,iZ);
186 }
187 
189 {
190  if (A < 1 || Z < 0 || Z > A) {
191 #ifdef G4VERBOSE
192  if (G4ParticleTable::GetParticleTable()->GetVerboseLevel()>0) {
193  G4cout << "G4NucleiProperties::GetMassExccess: Wrong values for A = "
194  << A << " and Z = " << Z << G4endl;
195  }
196 #endif
197  return 0.0;
198 
199  } else {
200 
201  G4bool inAMETable = false;
202  if ( ! use_old_evaluation ) {
203  inAMETable = G4NucleiPropertiesTableAME12::IsInTable(Z,A);
204  } else {
205  inAMETable = G4NucleiPropertiesTableAME03::IsInTable(Z,A);
206  }
207  if (inAMETable){
208  // AME table
209  if ( ! use_old_evaluation ) {
210  return G4NucleiPropertiesTableAME12::GetMassExcess(Z,A);
211  } else {
212  return G4NucleiPropertiesTableAME03::GetMassExcess(Z,A);
213  }
214  } else if (G4NucleiPropertiesTheoreticalTable::IsInTable(Z,A)){
215  return G4NucleiPropertiesTheoreticalTable::GetMassExcess(Z,A);
216  } else {
217  return MassExcess(A,Z);
218  }
219  }
220 
221 }
222 
223 
224 G4double G4NucleiProperties::GetAtomicMass(const G4double A, const G4double Z)
225 {
226  if (A < 1 || Z < 0 || Z > A) {
227 #ifdef G4VERBOSE
228  if (G4ParticleTable::GetParticleTable()->GetVerboseLevel()>0) {
229  G4cout << "G4NucleiProperties::GetAtomicMass: Wrong values for A = "
230  << A << " and Z = " << Z << G4endl;
231  }
232 #endif
233  return 0.0;
234 
235  } else if (std::fabs(A - G4int(A)) > 1.e-10) {
236  return AtomicMass(A,Z);
237 
238  } else {
239  G4int iA = G4int(A);
240  G4int iZ = G4int(Z);
241  G4bool inAMETable = false;
242  if ( ! use_old_evaluation ) {
243  inAMETable = G4NucleiPropertiesTableAME12::IsInTable(Z,A);
244  } else {
245  inAMETable = G4NucleiPropertiesTableAME03::IsInTable(Z,A);
246  }
247  if (inAMETable) {
248  if ( ! use_old_evaluation ) {
249  return G4NucleiPropertiesTableAME12::GetAtomicMass(Z,A);
250  } else {
251  return G4NucleiPropertiesTableAME03::GetAtomicMass(Z,A);
252  }
253  } else if (G4NucleiPropertiesTheoreticalTable::IsInTable(iZ,iA)){
254  return G4NucleiPropertiesTheoreticalTable::GetAtomicMass(iZ,iA);
255  } else {
256  return AtomicMass(A,Z);
257  }
258  }
259 }
260 
262 {
263  G4int iA = G4int(A);
264  G4int iZ = G4int(Z);
265  return GetBindingEnergy(iA,iZ);
266 }
267 
269 {
270  if (A < 1 || Z < 0 || Z > A) {
271 #ifdef G4VERBOSE
272  if (G4ParticleTable::GetParticleTable()->GetVerboseLevel()>0) {
273  G4cout << "G4NucleiProperties::GetMassExccess: Wrong values for A = "
274  << A << " and Z = " << Z << G4endl;
275  }
276 #endif
277  return 0.0;
278 
279  } else {
280  G4bool inAMETable = false;
281  if ( ! use_old_evaluation ) {
282  inAMETable = G4NucleiPropertiesTableAME12::IsInTable(Z,A);
283  } else {
284  inAMETable = G4NucleiPropertiesTableAME03::IsInTable(Z,A);
285  }
286  if (inAMETable) {
287  if ( ! use_old_evaluation ) {
288  return G4NucleiPropertiesTableAME12::GetBindingEnergy(Z,A);
289  } else {
290  return G4NucleiPropertiesTableAME03::GetBindingEnergy(Z,A);
291  }
292  } else if (G4NucleiPropertiesTheoreticalTable::IsInTable(Z,A)) {
293  return G4NucleiPropertiesTheoreticalTable::GetBindingEnergy(Z,A);
294  }else {
295  return BindingEnergy(A,Z);
296  }
297 
298  }
299 }
300 
301 
302 G4double G4NucleiProperties::MassExcess(G4double A, G4double Z)
303 {
304  return GetAtomicMass(A,Z) - A*amu_c2;
305 }
306 
307 G4double G4NucleiProperties::AtomicMass(G4double A, G4double Z)
308 {
309  //const G4double hydrogen_mass_excess;
310  //const G4double neutron_mass_excess;
311  G4double hydrogen_mass_excess;
312  G4double neutron_mass_excess;
313  if ( ! use_old_evaluation ) {
314  hydrogen_mass_excess = G4NucleiPropertiesTableAME12::GetMassExcess(1,1);
315  neutron_mass_excess = G4NucleiPropertiesTableAME12::GetMassExcess(0,1);
316  } else {
317  hydrogen_mass_excess = G4NucleiPropertiesTableAME03::GetMassExcess(1,1);
318  neutron_mass_excess = G4NucleiPropertiesTableAME03::GetMassExcess(0,1);
319  }
320 
321  G4double mass =
322  (A-Z)*neutron_mass_excess + Z*hydrogen_mass_excess - BindingEnergy(A,Z) + A*amu_c2;
323 
324  return mass;
325 }
326 
327 G4double G4NucleiProperties::NuclearMass(G4double A, G4double Z)
328 {
329  if (A < 1 || Z < 0 || Z > A) {
330 #ifdef G4VERBOSE
331  if (G4ParticleTable::GetParticleTable()->GetVerboseLevel()>0) {
332  G4cout << "G4NucleiProperties::NuclearMass: Wrong values for A = "
333  << A << " and Z = " << Z << G4endl;
334  }
335 #endif
336  return 0.0;
337  }
338 
339  G4double mass = AtomicMass(A,Z);
340  // atomic mass is converted to nuclear mass according formula in AME03 and 12
341  mass -= Z*electron_mass_c2;
342  mass += ( 14.4381*std::pow ( Z , 2.39 ) + 1.55468*1e-6*std::pow ( Z , 5.35 ) )*eV;
343 
344  return mass;
345 }
346 
347 G4double G4NucleiProperties::BindingEnergy(G4double A, G4double Z)
348 {
349  //
350  // Weitzsaecker's Mass formula
351  //
352  G4int Npairing = G4int(A-Z)%2; // pairing
353  G4int Zpairing = G4int(Z)%2;
354  G4double binding =
355  - 15.67*A // nuclear volume
356  + 17.23*std::pow(A,2./3.) // surface energy
357  + 93.15*((A/2.-Z)*(A/2.-Z))/A // asymmetry
358  + 0.6984523*Z*Z*std::pow(A,-1./3.); // coulomb
359  if( Npairing == Zpairing ) binding += (Npairing+Zpairing-1) * 12.0 / std::sqrt(A); // pairing
360 
361  return -binding*MeV;
362 }
363 
365 {
366  use_old_evaluation = val;
367 }
static G4double GetNuclearMass(const G4double A, const G4double Z)
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
#define G4ThreadLocal
Definition: tls.hh:89
int G4int
Definition: G4Types.hh:78
static void UseOldAMETable(G4bool val=true)
static G4double GetMassExcess(const G4int A, const G4int Z)
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
bool G4bool
Definition: G4Types.hh:79
static bool IsInStableTable(const G4double A, const G4double Z)
static constexpr double eV
Definition: G4SIunits.hh:215
float electron_mass_c2
Definition: hepunit.py:274
G4double GetPDGMass() const
static G4ParticleTable * GetParticleTable()
static G4double GetBindingEnergy(const G4int A, const G4int Z)
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
double G4double
Definition: G4Types.hh:76
float amu_c2
Definition: hepunit.py:277