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

Last change on this file since 835 was 824, checked in by garnier, 16 years ago

import all except CVS
File size: 12.0 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.13 2007/09/14 07:04:09 kurasige Exp $
28// GEANT4 tag $Name:  $
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.