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

Last change on this file since 859 was 850, checked in by garnier, 17 years ago

geant4.8.2 beta
File size: 12.0 KB
RevLine 
[824]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.13 2007/09/14 07:04:09 kurasige Exp $
[850]28// GEANT4 tag $Name: HEAD $
[824]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
42
43#include "G4NucleiProperties.hh"
44
45
46
47G4double G4NucleiProperties::AtomicMass(G4double A, G4double Z)
48{
49 const G4double hydrogen_mass_excess = G4NucleiPropertiesTable::GetMassExcess(1,1);
50 const G4double neutron_mass_excess = G4NucleiPropertiesTable::GetMassExcess(0,1);
51
52 G4double mass =
53 (A-Z)*neutron_mass_excess + Z*hydrogen_mass_excess - BindingEnergy(A,Z) + A*amu_c2;
54
55 return mass;
56}
57
58G4double G4NucleiProperties::BindingEnergy(G4double A, G4double Z)
59{
60 //
61 // Weitzsaecker's Mass formula
62 //
63 G4int Npairing = G4int(A-Z)%2; // pairing
64 G4int Zpairing = G4int(Z)%2;
65 G4double binding =
66 - 15.67*A // nuclear volume
67 + 17.23*std::pow(A,2./3.) // surface energy
68 + 93.15*((A/2.-Z)*(A/2.-Z))/A // asymmetry
69 + 0.6984523*Z*Z*std::pow(A,-1./3.); // coulomb
70 if( Npairing == Zpairing ) binding += (Npairing+Zpairing-1) * 12.0 / std::sqrt(A); // pairing
71
72 return -binding*MeV;
73}
74
75
76G4double G4NucleiProperties::GetNuclearMass(const G4double A, const G4double Z)
77{
78 if (A < 1 || Z < 0 || Z > A) {
79#ifdef G4VERBOSE
80 if (G4ParticleTable::GetParticleTable()->GetVerboseLevel()>0) {
81 G4cout << "G4NucleiProperties::GetNuclearMass: Wrong values for A = " << A
82 << " and Z = " << Z << G4endl;
83 }
84#endif
85 return 0.0;
86 } else {
87 G4ParticleDefinition * nucleus = 0;
88 if ( (Z<=2) ) {
89 if ( (Z==1)&&(A==1) ) {
90 nucleus = G4ParticleTable::GetParticleTable()->FindParticle("proton"); // proton
91 } else if ( (Z==0)&&(A==1) ) {
92 nucleus = G4ParticleTable::GetParticleTable()->FindParticle("neutron"); // neutron
93 } else if ( (Z==1)&&(A==2) ) {
94 nucleus = G4ParticleTable::GetParticleTable()->FindParticle("deuteron"); // deuteron
95 } else if ( (Z==1)&&(A==3) ) {
96 nucleus = G4ParticleTable::GetParticleTable()->FindParticle("triton"); // triton
97 } else if ( (Z==2)&&(A==4) ) {
98 nucleus = G4ParticleTable::GetParticleTable()->FindParticle("alpha"); // alpha
99 } else if ( (Z==2)&&(A==3) ) {
100 nucleus = G4ParticleTable::GetParticleTable()->FindParticle("He3"); // He3
101 }
102 }
103 if (nucleus!=0) {
104 return nucleus->GetPDGMass();
105 }else {
106 return GetAtomicMass(A,Z) - Z*electron_mass_c2 + 1.433e-5*MeV*std::pow(Z,2.39);
107 }
108 }
109}
110
111
112// G4double G4NucleiProperties::CameronMassExcess(const G4int A, const G4int Z)
113// {
114// const G4double alpha = -17.0354*MeV;
115// const G4double beta = -31.4506*MeV;
116// const G4double phi = 44.2355*MeV;
117// const G4double gamma = 25.8357*MeV;
118//
119// const G4double A13 = std::pow(G4double(A),1.0/3.0);
120// const G4double A23 = A13*A13;
121// const G4double A43 = A23*A23;
122// const G4double Z43 = std::pow(G4double(Z),4.0/3.0);
123// G4double D = (G4double(A) - 2.0*G4double(Z))/G4double(A);
124// D *= D; // D = std::pow((A-2Z)/A,2)
125//
126//
127// // Surface term
128// G4double SurfaceEnergy = (gamma - phi*D)*(1.0 - 0.62025/A23)*(1.0 - 0.62025/A23)*A23;
129//
130// // Coulomb term
131// G4double CoulombEnergy = 0.779*MeV*(G4double(Z*(Z-1))/A13)*(1.0-1.5849/A23+1.2273/A+1.5772/A43);
132//
133// // Exchange term
134// G4double ExchangeEnergy = -0.4323*MeV*(Z43/A13)*(1.0-0.57811/A13-0.14518/A23+0.49597/A);
135//
136// // Volume term
137// G4double VolumeEnergy = G4double(A)*(alpha-beta*D);
138//
139// // Shell+Pairing corrections for protons
140// G4double SPcorrectionZ;
141// if (Z <= TableSize) SPcorrectionZ = SPZTable[Z-1]*MeV;
142// else SPcorrectionZ = 0.0;
143//
144// // Shell+Pairing corrections for protons
145// G4int N = A - Z;
146// G4double SPcorrectionN;
147// if (N <= TableSize) SPcorrectionN = SPNTable[N-1]*MeV;
148// else SPcorrectionN = 0.0;
149//
150//
151// // Mass Excess
152// // First two terms give the mass excess of the neutrons and protons in the nucleus
153// // (see Cameron, Canadian Journal of Physics, 35, 1957 page 1022)
154// G4double MassExcess = 8.367*MeV*A - 0.783*MeV*Z +
155// SurfaceEnergy + CoulombEnergy + ExchangeEnergy + VolumeEnergy +
156// SPcorrectionZ + SPcorrectionN;
157//
158// return MassExcess;
159// }
160
161
162// S(Z)+P(Z) from Tab. 1 from A.G.W. Cameron, Canad. J. Phys., 35(1957)1021
163// or Delta M(Z) from Tab. 97 of book [1]
164// const G4double G4NucleiProperties::SPZTable[TableSize] = {
165// 20.80, 15.80, 21.00, 16.80, 19.80, 16.50, 18.80, 16.50, 18.50, 17.20, // 1 - 10
166// 18.26, 15.05, 16.01, 12.04, 13.27, 11.09, 12.17, 10.26, 11.04, 8.41, // 11 - 20
167// 9.79, 7.36, 8.15, 5.63, 5.88, 3.17, 3.32, .82, 1.83, .97, // 21 - 30
168// 2.33, 1.27, 2.92, 1.61, 2.91, 1.35, 2.40, .89, 1.74, .36, // 31
169// 0.95, -0.65, -0.04, -1.73, -0.96, -2.87, -2.05, -4.05, -3.40, -5.72, // 41
170// -3.75, -4.13, -2.42, -2.85, -1.01, -1.33, 0.54, -0.02, 1.74, 0.75, // 51
171// 2.24, 1.00, 1.98, 0.79, 1.54, 0.39, 1.08, 0.00, 0.78, -0.35, // 61
172// 0.58, -0.55, 0.59, -0.61, 0.59, -0.35, 0.32, -0.96, -0.52, -2.08, // 71
173// -2.46, -3.64, -1.55, -0.96, 0.97, 0.88, 2.37, 1.75, 2.72, 1.90, // 81
174// 2.55, 1.46, 1.93, 0.86, 1.17, 0.08, 0.39, -0.76, -0.39, -1.51, // 91 - 100
175// -1.17, -2.36, -1.95, -3.06, -2.62, -3.55, -2.95, -3.75, -3.07, -3.79, // 101 - 110
176// -3.06, -3.77, -3.05, -3.78, -3.12, -3.90, -3.35, -4.24, -3.86, -4.92, // 111 - 120
177// -5.06, -6.77, -7.41, -9.18,-10.16,-11.12, -9.76, -9.23, -7.96, -7.65, // 121 - 130
178// // --------- from this point there are not tabulated values -----------------------
179// 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, // 131 - 140
180// 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, // 141 - 150
181// 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, // 151
182// 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, // 161
183// 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, // 171
184// 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, // 181
185// 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 // 191 - 200
186// };
187
188// S(N)+P(N) from Tab. 1 from A.G.W. Cameron, Canad. J. Phys., 35(1957)1021
189// or Delta M(N) from Tab. 97 of book [1]
190// const G4double G4NucleiProperties::SPNTable[TableSize] = {
191// -8.40,-12.90, -8.00, 11.90, -9.20,-12.50,-10.80,-13.60,-11.20,-12.20, // 1 - 10
192// -12.81,-15.40,-13.07,-15.80,-13.81,-14.98,-12.63,-13.76,-11.37,-12.38, // 11 - 20
193// -9.23, -9.65, -7.64, -9.17, -8.05, -9.72, -8.87,-10.76, -8.64, -8.89, // 21 - 30
194// -6.60, -7.13, -4.77, -5.33, -3.06, -3.79, -1.72, -2.79, -0.93, -2.19, // 31
195// -0.52, -1.90, -0.45, -2.20, -1.22, -3.07, -2.42, -4.37, -3.94, -6.08, // 41
196// -4.49, -4.50, -3.14, -2.93, -1.04, -1.36, 0.69, 0.21, 2.11, 1.33, // 51
197// 3.29, 2.46, 4.30, 3.32, 4.79, 3.62, 4.97, 3.64, 4.63, 3.07, // 61
198// 4.06, 2.49, 3.30, 1.46, 2.06, 0.51, 0.74, -1.18, -1.26, -3.54, // 71
199// -3.97, -5.26, -4.18, -3.71, -2.10, -1.70, -0.08, -0.18, 0.94, 0.27, // 81
200// 1.13, 0.08, 0.91, -0.31, 0.49, -0.78, 0.08, -1.15, -0.23, -1.41, // 91 - 100
201// -0.42, -1.55, -0.55, -1.66, -0.66, -1.73, -0.75, -1.74, -0.78, -1.69, // 101 - 110
202// -0.78, -1.60, -0.75, -1.46, -0.67, -1.26, -0.51, -1.04, -0.53, -1.84, // 111 - 120
203// -2.42, -4.52, -4.76, -6.33, -6.76, -7.81, -5.80, -5.37, -3.63, -3.35, // 121 - 130
204// -1.75, -1.88, -0.61, -0.90, 0.09, -0.32, 0.55, -0.13, 0.70, -0.06, // 131 - 140
205// 0.49, -0.20, 0.40, -0.22, 0.36, -0.09, 0.58, 0.12, 0.75, 0.15, // 141 - 150
206// 0.70, 0.17, 1.11, 0.89, 1.85, 1.62, 2.54, 2.29, 3.20, 2.91, // 151
207// 3.84, 3.53, 4.48, 4.15, 5.12, 4.78, 5.75, 5.39, 6.31, 5.91, // 161
208// 6.87, 6.33, 7.13, 6.61, 7.30, 6.31, 6.27, 4.83, 4.49, 2.85, // 171
209// 2.32, 0.58, -0.11, -0.98, 0.81, 1.77, 3.37, 4.13, 5.60, 6.15, // 181
210// 7.29, 7.35, 7.95, 7.67, 8.16, 7.83, 8.31, 8.01, 8.53, 8.27 // 191 - 200
211// };
212
213
214
215G4double G4NucleiProperties::GetMassExcess(const G4int A, const G4int Z)
216{
217 if (A < 1 || Z < 0 || Z > A) {
218#ifdef G4VERBOSE
219 if (G4ParticleTable::GetParticleTable()->GetVerboseLevel()>0) {
220 G4cout << "G4NucleiProperties::GetMassExccess: Wrong values for A = "
221 << A << " and Z = " << Z << G4endl;
222 }
223#endif
224 return 0.0;
225
226 } else {
227
228 if (G4NucleiPropertiesTable::IsInTable(Z,A)){
229 return G4NucleiPropertiesTable::GetMassExcess(Z,A);
230 } else if (G4NucleiPropertiesTheoreticalTable::IsInTable(Z,A)){
231 return G4NucleiPropertiesTheoreticalTable::GetMassExcess(Z,A);
232 } else {
233 return MassExcess(A,Z);
234 }
235 }
236
237}
238
239
240G4double G4NucleiProperties::GetAtomicMass(const G4double A, const G4double Z)
241{
242 if (Z < 0 || Z > A) {
243#ifdef G4VERBOSE
244 if (G4ParticleTable::GetParticleTable()->GetVerboseLevel()>0) {
245 G4cout << "G4NucleiProperties::GetAtomicMass: Wrong values for A = "
246 << A << " and Z = " << Z << G4endl;
247 }
248#endif
249 return 0.0;
250
251 } else if (std::abs(A - G4int(A)) > 1.e-10) {
252 return AtomicMass(A,Z);
253
254 } else {
255 G4int iA = G4int(A);
256 G4int iZ = G4int(Z);
257 if (G4NucleiPropertiesTable::IsInTable(iZ,iA)) {
258 return G4NucleiPropertiesTable::GetAtomicMass(iZ,iA);
259 } else if (G4NucleiPropertiesTheoreticalTable::IsInTable(iZ,iA)){
260 return G4NucleiPropertiesTheoreticalTable::GetAtomicMass(iZ,iA);
261 } else {
262 return AtomicMass(A,Z);
263 }
264 }
265}
266
267G4double G4NucleiProperties::GetBindingEnergy(const G4int A, const G4int Z)
268{
269 if (A < 1 || Z < 0 || Z > A) {
270#ifdef G4VERBOSE
271 if (G4ParticleTable::GetParticleTable()->GetVerboseLevel()>0) {
272 G4cout << "G4NucleiProperties::GetMassExccess: Wrong values for A = "
273 << A << " and Z = " << Z << G4endl;
274 }
275#endif
276 return 0.0;
277
278 } else {
279 if (G4NucleiPropertiesTable::IsInTable(Z,A)) {
280 return G4NucleiPropertiesTable::GetBindingEnergy(Z,A);
281 } else if (G4NucleiPropertiesTheoreticalTable::IsInTable(Z,A)) {
282 return G4NucleiPropertiesTheoreticalTable::GetBindingEnergy(Z,A);
283 }else {
284 return BindingEnergy(A,Z);
285 }
286
287 }
288}
289
290
291G4double G4NucleiProperties::MassExcess(G4double A, G4double Z)
292{
293 return GetAtomicMass(A,Z) - A*amu_c2;
294}
295
Note: See TracBrowser for help on using the repository browser.