Geant4: G4NucleiProperties.cc Source File

00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 //
00027 // $Id$
00028 //
00029 // 
00030 // ------------------------------------------------------------
00031 //      GEANT 4 class header file 
00032 //
00033 // ------------------------------------------------------------
00034 //
00035 // Hadronic Process: Nuclear De-excitations
00036 // by V. Lara (Oct 1998)
00037 // Migrate into particles category by H.Kurashige (17 Nov. 98)
00038 // Added Shell-Pairing corrections to the Cameron mass 
00039 // excess formula by V.Lara (9 May 99)
00040 // 090331 Migrate to AME03 by Koi, Tatsumi 
00041 
00042 #include "G4NucleiProperties.hh"
00043 #include "G4PhysicalConstants.hh"
00044 #include "G4SystemOfUnits.hh"
00045 
00046 G4double G4NucleiProperties::mass_proton = -1.;
00047 G4double G4NucleiProperties::mass_neutron = -1.;
00048 G4double G4NucleiProperties::mass_deuteron = -1.;
00049 G4double G4NucleiProperties::mass_triton = -1.;
00050 G4double G4NucleiProperties::mass_alpha = -1.;
00051 G4double G4NucleiProperties::mass_He3 = -1.;
00052 
00053 G4double G4NucleiProperties::GetNuclearMass(const G4double A, const G4double Z)
00054 {
00055   G4double mass=0.0;
00056 
00057   if (std::fabs(A - G4int(A)) > 1.e-10) {
00058     mass = NuclearMass(A,Z);
00059  
00060   } else {
00061     // use mass table
00062     G4int iZ = G4int(Z);
00063     G4int iA = G4int(A);
00064     mass =GetNuclearMass(iA,iZ);
00065   }
00066   
00067    return mass;
00068 }
00069 
00070 
00071 G4double G4NucleiProperties::GetNuclearMass(const G4int A, const G4int Z)
00072 {
00073   if (mass_proton  <= 0.0 ) {
00074     const G4ParticleDefinition * nucleus = 0;
00075     nucleus = G4ParticleTable::GetParticleTable()->FindParticle("proton"); // proton 
00076     if (nucleus!=0) mass_proton = nucleus->GetPDGMass();
00077     nucleus = G4ParticleTable::GetParticleTable()->FindParticle("neutron"); // neutron 
00078     if (nucleus!=0) mass_neutron = nucleus->GetPDGMass();
00079     nucleus = G4ParticleTable::GetParticleTable()->FindParticle("deuteron"); // deuteron 
00080     if (nucleus!=0) mass_deuteron = nucleus->GetPDGMass();
00081     nucleus = G4ParticleTable::GetParticleTable()->FindParticle("triton"); // triton 
00082     if (nucleus!=0) mass_triton = nucleus->GetPDGMass();
00083     nucleus = G4ParticleTable::GetParticleTable()->FindParticle("alpha"); // alpha 
00084     if (nucleus!=0) mass_alpha = nucleus->GetPDGMass();
00085     nucleus = G4ParticleTable::GetParticleTable()->FindParticle("He3"); // He3 
00086     if (nucleus!=0) mass_He3 = nucleus->GetPDGMass();
00087 
00088   }
00089 
00090   if (A < 1 || Z < 0 || Z > A) {
00091 #ifdef G4VERBOSE
00092     if (G4ParticleTable::GetParticleTable()->GetVerboseLevel()>0) {
00093       G4cerr << "G4NucleiProperties::GetNuclearMass: Wrong values for A = " << A 
00094              << " and Z = " << Z << G4endl;
00095     }
00096 #endif    
00097     return 0.0;
00098   }
00099   
00100   G4double mass= -1.;
00101   if ( (Z<=2) ) {
00102     // light nuclei
00103     if ( (Z==1)&&(A==1) ) {
00104       mass = mass_proton;
00105     } else if ( (Z==0)&&(A==1) ) {
00106       mass = mass_neutron;
00107     } else if ( (Z==1)&&(A==2) ) {
00108       mass = mass_deuteron;
00109     } else if ( (Z==1)&&(A==3) ) {
00110       mass = mass_triton;
00111     } else if ( (Z==2)&&(A==4) ) {
00112       mass = mass_alpha;
00113     } else if ( (Z==2)&&(A==3) ) {
00114       mass = mass_He3;
00115     }
00116   }
00117   
00118   if (mass < 0.) {
00119     if (G4NucleiPropertiesTableAME03::IsInTable(Z,A)) {
00120       // AME 03 table
00121       mass = G4NucleiPropertiesTableAME03::GetNuclearMass(Z,A);
00122     } else if (G4NucleiPropertiesTheoreticalTable::IsInTable(Z,A)){
00123       // Theoretical table
00124       mass = G4NucleiPropertiesTheoreticalTable::GetNuclearMass(Z,A);
00125     } else {
00126       mass = NuclearMass(G4double(A),G4double(Z));
00127     }
00128   }
00129 
00130   if (mass < 0.) mass = 0.0;
00131   return mass;
00132 }
00133 
00134 G4bool G4NucleiProperties::IsInStableTable(const G4double A, const G4double Z)
00135 {
00136   G4int iA = G4int(A);
00137   G4int iZ = G4int(Z);
00138   return IsInStableTable(iA, iZ);
00139 }
00140 
00141 G4bool G4NucleiProperties::IsInStableTable(const G4int A, const int Z)
00142 {
00143   if (A < 1 || Z < 0 || Z > A) {
00144 #ifdef G4VERBOSE
00145     if (G4ParticleTable::GetParticleTable()->GetVerboseLevel()>0) {
00146       G4cerr << "G4NucleiProperties::IsInStableTable: Wrong values for A = " 
00147              << A << " and Z = " << Z << G4endl;        
00148     }
00149 #endif 
00150     return false;
00151   } 
00152 
00153   return G4NucleiPropertiesTableAME03::IsInTable(Z,A);
00154 
00155 }
00156 
00157 G4double G4NucleiProperties::GetMassExcess(const G4double A, const G4double Z)
00158 {
00159   G4int iA = G4int(A);
00160   G4int iZ = G4int(Z);
00161   return GetMassExcess(iA,iZ);
00162 }
00163 
00164 G4double G4NucleiProperties::GetMassExcess(const G4int A, const G4int Z)
00165 {
00166   if (A < 1 || Z < 0 || Z > A) {
00167 #ifdef G4VERBOSE
00168     if (G4ParticleTable::GetParticleTable()->GetVerboseLevel()>0) {
00169       G4cerr << "G4NucleiProperties::GetMassExccess: Wrong values for A = " 
00170              << A << " and Z = " << Z << G4endl;
00171     }
00172 #endif    
00173     return 0.0;
00174     
00175   } else {
00176 
00177     if (G4NucleiPropertiesTableAME03::IsInTable(Z,A)){
00178       return G4NucleiPropertiesTableAME03::GetMassExcess(Z,A);
00179     } else if (G4NucleiPropertiesTheoreticalTable::IsInTable(Z,A)){
00180       return G4NucleiPropertiesTheoreticalTable::GetMassExcess(Z,A);
00181     } else {
00182       return MassExcess(A,Z);
00183     }
00184   }
00185 
00186 }
00187 
00188 
00189 G4double G4NucleiProperties::GetAtomicMass(const G4double A, const G4double Z)
00190 {
00191   if (A < 1 || Z < 0 || Z > A) {
00192 #ifdef G4VERBOSE
00193     if (G4ParticleTable::GetParticleTable()->GetVerboseLevel()>0) {
00194       G4cerr << "G4NucleiProperties::GetAtomicMass: Wrong values for A = " 
00195              << A << " and Z = " << Z << G4endl;        
00196     }
00197 #endif 
00198     return 0.0;
00199 
00200   } else if (std::fabs(A - G4int(A)) > 1.e-10) {
00201     return AtomicMass(A,Z);
00202 
00203   } else {
00204     G4int iA = G4int(A);
00205     G4int iZ = G4int(Z);
00206     if (G4NucleiPropertiesTableAME03::IsInTable(iZ,iA)) {
00207       return G4NucleiPropertiesTableAME03::GetAtomicMass(iZ,iA);
00208     } else if (G4NucleiPropertiesTheoreticalTable::IsInTable(iZ,iA)){
00209       return G4NucleiPropertiesTheoreticalTable::GetAtomicMass(iZ,iA);
00210     } else {
00211       return AtomicMass(A,Z);
00212     }
00213   }
00214 }
00215 
00216 G4double G4NucleiProperties::GetBindingEnergy(const G4double A, const G4double Z)
00217 {
00218   G4int iA = G4int(A);
00219   G4int iZ = G4int(Z);
00220   return GetBindingEnergy(iA,iZ);
00221 }
00222 
00223 G4double G4NucleiProperties::GetBindingEnergy(const G4int A, const G4int Z)
00224 {
00225   if (A < 1 || Z < 0 || Z > A) {
00226 #ifdef G4VERBOSE
00227     if (G4ParticleTable::GetParticleTable()->GetVerboseLevel()>0) {
00228       G4cerr << "G4NucleiProperties::GetMassExccess: Wrong values for A = " 
00229              << A << " and Z = " << Z << G4endl;
00230     }
00231 #endif
00232     return 0.0;
00233 
00234   } else {
00235     if (G4NucleiPropertiesTableAME03::IsInTable(Z,A)) {
00236       return G4NucleiPropertiesTableAME03::GetBindingEnergy(Z,A);
00237     } else if (G4NucleiPropertiesTheoreticalTable::IsInTable(Z,A)) {
00238       return G4NucleiPropertiesTheoreticalTable::GetBindingEnergy(Z,A);
00239     }else {
00240       return BindingEnergy(A,Z);
00241     }
00242 
00243   }
00244 }
00245 
00246 
00247 G4double G4NucleiProperties::MassExcess(G4double A, G4double Z) 
00248 {
00249   return GetAtomicMass(A,Z) - A*amu_c2;
00250 }
00251 
00252 G4double  G4NucleiProperties::AtomicMass(G4double A, G4double Z)
00253 {
00254   const G4double hydrogen_mass_excess = G4NucleiPropertiesTableAME03::GetMassExcess(1,1);
00255   const G4double neutron_mass_excess =  G4NucleiPropertiesTableAME03::GetMassExcess(0,1);
00256 
00257   G4double mass =
00258       (A-Z)*neutron_mass_excess + Z*hydrogen_mass_excess - BindingEnergy(A,Z) + A*amu_c2;
00259 
00260   return mass;
00261 }
00262 
00263 G4double  G4NucleiProperties::NuclearMass(G4double A, G4double Z)
00264 {
00265   if (A < 1 || Z < 0 || Z > A) {
00266 #ifdef G4VERBOSE
00267     if (G4ParticleTable::GetParticleTable()->GetVerboseLevel()>0) {
00268       G4cerr << "G4NucleiProperties::NuclearMass: Wrong values for A = " 
00269              << A << " and Z = " << Z << G4endl;
00270     }
00271 #endif 
00272     return 0.0;
00273   }
00274 
00275   G4double mass = AtomicMass(A,Z);
00276   // atomic mass is converted to nuclear mass according formula in  AME03
00277   mass -= Z*electron_mass_c2;
00278   mass += ( 14.4381*std::pow ( Z , 2.39 ) + 1.55468*1e-6*std::pow ( Z , 5.35 ) )*eV;      
00279 
00280   return mass;
00281 }
00282 
00283 G4double  G4NucleiProperties::BindingEnergy(G4double A, G4double Z)
00284 { 
00285   //
00286   // Weitzsaecker's Mass formula
00287   //
00288   G4int Npairing = G4int(A-Z)%2;                  // pairing
00289   G4int Zpairing = G4int(Z)%2;
00290   G4double binding =
00291       - 15.67*A                           // nuclear volume
00292       + 17.23*std::pow(A,2./3.)                // surface energy
00293       + 93.15*((A/2.-Z)*(A/2.-Z))/A       // asymmetry
00294       + 0.6984523*Z*Z*std::pow(A,-1./3.);      // coulomb
00295   if( Npairing == Zpairing ) binding += (Npairing+Zpairing-1) * 12.0 / std::sqrt(A);  // pairing
00296 
00297   return -binding*MeV;
00298 }
00299