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

Last change on this file since 1199 was 1196, checked in by garnier, 15 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.