source: trunk/source/particles/management/src/G4NucleiProperties.cc@ 1260

Last change on this file since 1260 was 1196, checked in by garnier, 16 years ago

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