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

Last change on this file since 1335 was 1315, checked in by garnier, 15 years ago

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