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

Last change on this file since 1340 was 1340, checked in by garnier, 14 years ago

update ti head
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.23 2010/08/10 15:47:42 kurasige Exp $
28// GEANT4 tag $Name: particles-V09-03-15 $
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    const 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.