Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer << 3 // * DISCLAIMER * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers * 6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which * 7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: * 8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license * 9 // * include a list of copyright holders. << 10 // * 9 // * * 11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file << 14 // * use. * 16 // * for the full disclaimer and the limitatio << 17 // * 15 // * * 18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the * 19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. * 20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work * 21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this * 22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. * 23 // * acceptance of all terms of the Geant4 Sof << 24 // ******************************************* 21 // ******************************************************************** 25 // 22 // 26 // G4NucleiProperties class implementation << 27 // 23 // 28 // Author: V.Lara, October 1998 << 24 // $Id: G4NucleiProperties.cc,v 1.10 2002/05/30 02:13:10 kurasige Exp $ 29 // History: << 25 // GEANT4 tag $Name: geant4-05-00 $ 30 // - 17.11.1998, H.Kurashige - Migrated into p << 26 // 31 // - 31.03.2009, T.Koi - Migrated to AME03 << 27 // 32 // ------------------------------------------- << 28 // ------------------------------------------------------------ >> 29 // GEANT 4 class header file >> 30 // >> 31 // ------------------------------------------------------------ >> 32 // >> 33 // Hadronic Process: Nuclear De-excitations >> 34 // by V. Lara (Oct 1998) >> 35 // Migrate into particles category by H.Kurashige (17 Nov. 98) >> 36 // Added Shell-Pairing corrections to the Cameron mass >> 37 // excess formula by V.Lara (9 May 99) 33 38 34 #include "G4NucleiProperties.hh" << 35 39 36 #include "G4NucleiPropertiesTableAME12.hh" << 40 #include "G4NucleiProperties.hh" 37 #include "G4NucleiPropertiesTheoreticalTable.h << 38 #include "G4ParticleTable.hh" << 39 #include "G4PhysicalConstants.hh" << 40 #include "G4SystemOfUnits.hh" << 41 << 42 G4ThreadLocal G4double G4NucleiProperties::mas << 43 G4ThreadLocal G4double G4NucleiProperties::mas << 44 G4ThreadLocal G4double G4NucleiProperties::mas << 45 G4ThreadLocal G4double G4NucleiProperties::mas << 46 G4ThreadLocal G4double G4NucleiProperties::mas << 47 G4ThreadLocal G4double G4NucleiProperties::mas << 48 41 49 G4double G4NucleiProperties::GetNuclearMass(co << 50 { << 51 G4double mass = 0.0; << 52 42 53 if (std::fabs(A - G4int(A)) > 1.e-10) { << 54 mass = NuclearMass(A, Z); << 55 } << 56 else { << 57 // use mass table << 58 auto iZ = G4int(Z); << 59 auto iA = G4int(A); << 60 mass = GetNuclearMass(iA, iZ); << 61 } << 62 return mass; << 63 } << 64 43 65 G4double G4NucleiProperties::GetNuclearMass(co << 44 G4double G4NucleiProperties::AtomicMass(G4double A, G4double Z) 66 { 45 { 67 if (mass_proton <= 0.0) { << 46 const G4double hydrogen_mass_excess = G4NucleiPropertiesTable::GetMassExcess(1,1); 68 const G4ParticleDefinition* nucleus = null << 47 const G4double neutron_mass_excess = G4NucleiPropertiesTable::GetMassExcess(0,1); 69 nucleus = G4ParticleTable::GetParticleTabl << 70 if (nucleus != nullptr) mass_neutron = nuc << 71 << 72 nucleus = G4ParticleTable::GetParticleTabl << 73 if (nucleus != nullptr) mass_deuteron = nu << 74 48 75 nucleus = G4ParticleTable::GetParticleTabl << 49 G4double mass = 76 if (nucleus != nullptr) mass_triton = nucl << 50 (A-Z)*neutron_mass_excess + Z*hydrogen_mass_excess - BindingEnergy(A,Z) + A*amu_c2; 77 << 78 nucleus = G4ParticleTable::GetParticleTabl << 79 if (nucleus != nullptr) mass_alpha = nucle << 80 << 81 nucleus = G4ParticleTable::GetParticleTabl << 82 if (nucleus != nullptr) mass_He3 = nucleus << 83 << 84 nucleus = G4ParticleTable::GetParticleTabl << 85 if (nucleus != nullptr) mass_proton = nucl << 86 } << 87 << 88 if (A < 1 || Z < 0 || Z > A) { << 89 #ifdef G4VERBOSE << 90 if (G4ParticleTable::GetParticleTable()->G << 91 G4cout << "G4NucleiProperties::GetNuclea << 92 << G4endl; << 93 } << 94 #endif << 95 return 0.0; << 96 } << 97 << 98 G4double mass = -1.; << 99 if ((Z <= 2)) { << 100 // light nuclei << 101 if ((Z == 1) && (A == 1)) { << 102 mass = mass_proton; << 103 } << 104 else if ((Z == 0) && (A == 1)) { << 105 mass = mass_neutron; << 106 } << 107 else if ((Z == 1) && (A == 2)) { << 108 mass = mass_deuteron; << 109 } << 110 else if ((Z == 1) && (A == 3)) { << 111 mass = mass_triton; << 112 } << 113 else if ((Z == 2) && (A == 4)) { << 114 mass = mass_alpha; << 115 } << 116 else if ((Z == 2) && (A == 3)) { << 117 mass = mass_He3; << 118 } << 119 } << 120 << 121 if (mass < 0.) { << 122 if (G4NucleiPropertiesTableAME12::IsInTabl << 123 // AME table << 124 mass = G4NucleiPropertiesTableAME12::Get << 125 } << 126 else if (G4NucleiPropertiesTheoreticalTabl << 127 // Theoretical table << 128 mass = G4NucleiPropertiesTheoreticalTabl << 129 } << 130 else if (Z == A) { << 131 mass = A * mass_proton; << 132 } << 133 else if (0 == Z) { << 134 mass = A * mass_neutron; << 135 } << 136 else { << 137 mass = NuclearMass(G4double(A), G4double << 138 } << 139 } << 140 51 141 if (mass < 0.) mass = 0.0; << 142 return mass; 52 return mass; 143 } 53 } 144 54 145 G4bool G4NucleiProperties::IsInStableTable(con << 55 G4double G4NucleiProperties::BindingEnergy(G4double A, G4double Z) 146 { << 56 { 147 auto iA = G4int(A); << 57 // 148 auto iZ = G4int(Z); << 58 // Weitzsaecker's Mass formula 149 return IsInStableTable(iA, iZ); << 59 // >> 60 G4int Npairing = G4int(A-Z)%2; // pairing >> 61 G4int Zpairing = G4int(Z)%2; >> 62 G4double binding = >> 63 - 15.67*A // nuclear volume >> 64 + 17.23*pow(A,2./3.) // surface energy >> 65 + 93.15*((A/2.-Z)*(A/2.-Z))/A // asymmetry >> 66 + 0.6984523*Z*Z*pow(A,-1./3.); // coulomb >> 67 if( Npairing == Zpairing ) binding += (Npairing+Zpairing-1) * 12.0 / sqrt(A); // pairing >> 68 >> 69 return -binding*MeV; 150 } 70 } 151 71 152 G4bool G4NucleiProperties::IsInStableTable(con << 72 >> 73 G4double G4NucleiProperties::GetNuclearMass(const G4double A, const G4double Z) 153 { 74 { 154 if (A < 1 || Z < 0 || Z > A) { << 75 if (A < 1 || Z < 0 || Z > A) { 155 #ifdef G4VERBOSE << 76 G4cout << "G4NucleiProperties::GetNuclearMass: Wrong values for A = " << A 156 if (G4ParticleTable::GetParticleTable()->G << 77 << " and Z = " << Z << G4endl; 157 G4cout << "G4NucleiProperties::IsInStabl << 78 return 0.0; 158 << " and Z = " << Z << G4endl; << 79 } else { 159 } << 80 G4ParticleDefinition * nucleus = 0; 160 #endif << 81 if ( (Z<=2) ) { 161 return false; << 82 if ( (Z==1)&&(A==1) ) { 162 } << 83 nucleus = G4ParticleTable::GetParticleTable()->FindParticle("proton"); // proton >> 84 } else if ( (Z==0)&&(A==1) ) { >> 85 nucleus = G4ParticleTable::GetParticleTable()->FindParticle("neutron"); // neutron >> 86 } else if ( (Z==1)&&(A==2) ) { >> 87 nucleus = G4ParticleTable::GetParticleTable()->FindParticle("deuteron"); // deuteron >> 88 } else if ( (Z==1)&&(A==3) ) { >> 89 nucleus = G4ParticleTable::GetParticleTable()->FindParticle("triton"); // triton >> 90 } else if ( (Z==2)&&(A==4) ) { >> 91 nucleus = G4ParticleTable::GetParticleTable()->FindParticle("alpha"); // alpha >> 92 } else if ( (Z==2)&&(A==3) ) { >> 93 nucleus = G4ParticleTable::GetParticleTable()->FindParticle("He3"); // He3 >> 94 } >> 95 } >> 96 if (nucleus!=0) { >> 97 return nucleus->GetPDGMass(); >> 98 }else { >> 99 return GetAtomicMass(A,Z) - Z*electron_mass_c2 + 1.433e-5*MeV*pow(Z,2.39); >> 100 } >> 101 } >> 102 } >> 103 >> 104 >> 105 // G4double G4NucleiProperties::CameronMassExcess(const G4int A, const G4int Z) >> 106 // { >> 107 // const G4double alpha = -17.0354*MeV; >> 108 // const G4double beta = -31.4506*MeV; >> 109 // const G4double phi = 44.2355*MeV; >> 110 // const G4double gamma = 25.8357*MeV; >> 111 // >> 112 // const G4double A13 = pow(G4double(A),1.0/3.0); >> 113 // const G4double A23 = A13*A13; >> 114 // const G4double A43 = A23*A23; >> 115 // const G4double Z43 = pow(G4double(Z),4.0/3.0); >> 116 // G4double D = (G4double(A) - 2.0*G4double(Z))/G4double(A); >> 117 // D *= D; // D = pow((A-2Z)/A,2) >> 118 // >> 119 // >> 120 // // Surface term >> 121 // G4double SurfaceEnergy = (gamma - phi*D)*(1.0 - 0.62025/A23)*(1.0 - 0.62025/A23)*A23; >> 122 // >> 123 // // Coulomb term >> 124 // G4double CoulombEnergy = 0.779*MeV*(G4double(Z*(Z-1))/A13)*(1.0-1.5849/A23+1.2273/A+1.5772/A43); >> 125 // >> 126 // // Exchange term >> 127 // G4double ExchangeEnergy = -0.4323*MeV*(Z43/A13)*(1.0-0.57811/A13-0.14518/A23+0.49597/A); >> 128 // >> 129 // // Volume term >> 130 // G4double VolumeEnergy = G4double(A)*(alpha-beta*D); >> 131 // >> 132 // // Shell+Pairing corrections for protons >> 133 // G4double SPcorrectionZ; >> 134 // if (Z <= TableSize) SPcorrectionZ = SPZTable[Z-1]*MeV; >> 135 // else SPcorrectionZ = 0.0; >> 136 // >> 137 // // Shell+Pairing corrections for protons >> 138 // G4int N = A - Z; >> 139 // G4double SPcorrectionN; >> 140 // if (N <= TableSize) SPcorrectionN = SPNTable[N-1]*MeV; >> 141 // else SPcorrectionN = 0.0; >> 142 // >> 143 // >> 144 // // Mass Excess >> 145 // // First two terms give the mass excess of the neutrons and protons in the nucleus >> 146 // // (see Cameron, Canadian Journal of Physics, 35, 1957 page 1022) >> 147 // G4double MassExcess = 8.367*MeV*A - 0.783*MeV*Z + >> 148 // SurfaceEnergy + CoulombEnergy + ExchangeEnergy + VolumeEnergy + >> 149 // SPcorrectionZ + SPcorrectionN; >> 150 // >> 151 // return MassExcess; >> 152 // } >> 153 >> 154 >> 155 // S(Z)+P(Z) from Tab. 1 from A.G.W. Cameron, Canad. J. Phys., 35(1957)1021 >> 156 // or Delta M(Z) from Tab. 97 of book [1] >> 157 // const G4double G4NucleiProperties::SPZTable[TableSize] = { >> 158 // 20.80, 15.80, 21.00, 16.80, 19.80, 16.50, 18.80, 16.50, 18.50, 17.20, // 1 - 10 >> 159 // 18.26, 15.05, 16.01, 12.04, 13.27, 11.09, 12.17, 10.26, 11.04, 8.41, // 11 - 20 >> 160 // 9.79, 7.36, 8.15, 5.63, 5.88, 3.17, 3.32, .82, 1.83, .97, // 21 - 30 >> 161 // 2.33, 1.27, 2.92, 1.61, 2.91, 1.35, 2.40, .89, 1.74, .36, // 31 >> 162 // 0.95, -0.65, -0.04, -1.73, -0.96, -2.87, -2.05, -4.05, -3.40, -5.72, // 41 >> 163 // -3.75, -4.13, -2.42, -2.85, -1.01, -1.33, 0.54, -0.02, 1.74, 0.75, // 51 >> 164 // 2.24, 1.00, 1.98, 0.79, 1.54, 0.39, 1.08, 0.00, 0.78, -0.35, // 61 >> 165 // 0.58, -0.55, 0.59, -0.61, 0.59, -0.35, 0.32, -0.96, -0.52, -2.08, // 71 >> 166 // -2.46, -3.64, -1.55, -0.96, 0.97, 0.88, 2.37, 1.75, 2.72, 1.90, // 81 >> 167 // 2.55, 1.46, 1.93, 0.86, 1.17, 0.08, 0.39, -0.76, -0.39, -1.51, // 91 - 100 >> 168 // -1.17, -2.36, -1.95, -3.06, -2.62, -3.55, -2.95, -3.75, -3.07, -3.79, // 101 - 110 >> 169 // -3.06, -3.77, -3.05, -3.78, -3.12, -3.90, -3.35, -4.24, -3.86, -4.92, // 111 - 120 >> 170 // -5.06, -6.77, -7.41, -9.18,-10.16,-11.12, -9.76, -9.23, -7.96, -7.65, // 121 - 130 >> 171 // // --------- from this point there are not tabulated values ----------------------- >> 172 // 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, // 131 - 140 >> 173 // 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, // 141 - 150 >> 174 // 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, // 151 >> 175 // 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, // 161 >> 176 // 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, // 171 >> 177 // 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, // 181 >> 178 // 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 // 191 - 200 >> 179 // }; >> 180 >> 181 // S(N)+P(N) from Tab. 1 from A.G.W. Cameron, Canad. J. Phys., 35(1957)1021 >> 182 // or Delta M(N) from Tab. 97 of book [1] >> 183 // const G4double G4NucleiProperties::SPNTable[TableSize] = { >> 184 // -8.40,-12.90, -8.00, 11.90, -9.20,-12.50,-10.80,-13.60,-11.20,-12.20, // 1 - 10 >> 185 // -12.81,-15.40,-13.07,-15.80,-13.81,-14.98,-12.63,-13.76,-11.37,-12.38, // 11 - 20 >> 186 // -9.23, -9.65, -7.64, -9.17, -8.05, -9.72, -8.87,-10.76, -8.64, -8.89, // 21 - 30 >> 187 // -6.60, -7.13, -4.77, -5.33, -3.06, -3.79, -1.72, -2.79, -0.93, -2.19, // 31 >> 188 // -0.52, -1.90, -0.45, -2.20, -1.22, -3.07, -2.42, -4.37, -3.94, -6.08, // 41 >> 189 // -4.49, -4.50, -3.14, -2.93, -1.04, -1.36, 0.69, 0.21, 2.11, 1.33, // 51 >> 190 // 3.29, 2.46, 4.30, 3.32, 4.79, 3.62, 4.97, 3.64, 4.63, 3.07, // 61 >> 191 // 4.06, 2.49, 3.30, 1.46, 2.06, 0.51, 0.74, -1.18, -1.26, -3.54, // 71 >> 192 // -3.97, -5.26, -4.18, -3.71, -2.10, -1.70, -0.08, -0.18, 0.94, 0.27, // 81 >> 193 // 1.13, 0.08, 0.91, -0.31, 0.49, -0.78, 0.08, -1.15, -0.23, -1.41, // 91 - 100 >> 194 // -0.42, -1.55, -0.55, -1.66, -0.66, -1.73, -0.75, -1.74, -0.78, -1.69, // 101 - 110 >> 195 // -0.78, -1.60, -0.75, -1.46, -0.67, -1.26, -0.51, -1.04, -0.53, -1.84, // 111 - 120 >> 196 // -2.42, -4.52, -4.76, -6.33, -6.76, -7.81, -5.80, -5.37, -3.63, -3.35, // 121 - 130 >> 197 // -1.75, -1.88, -0.61, -0.90, 0.09, -0.32, 0.55, -0.13, 0.70, -0.06, // 131 - 140 >> 198 // 0.49, -0.20, 0.40, -0.22, 0.36, -0.09, 0.58, 0.12, 0.75, 0.15, // 141 - 150 >> 199 // 0.70, 0.17, 1.11, 0.89, 1.85, 1.62, 2.54, 2.29, 3.20, 2.91, // 151 >> 200 // 3.84, 3.53, 4.48, 4.15, 5.12, 4.78, 5.75, 5.39, 6.31, 5.91, // 161 >> 201 // 6.87, 6.33, 7.13, 6.61, 7.30, 6.31, 6.27, 4.83, 4.49, 2.85, // 171 >> 202 // 2.32, 0.58, -0.11, -0.98, 0.81, 1.77, 3.37, 4.13, 5.60, 6.15, // 181 >> 203 // 7.29, 7.35, 7.95, 7.67, 8.16, 7.83, 8.31, 8.01, 8.53, 8.27 // 191 - 200 >> 204 // }; 163 205 164 return G4NucleiPropertiesTableAME12::IsInTab << 165 } << 166 206 167 G4double G4NucleiProperties::GetMassExcess(con << 168 { << 169 auto iA = G4int(A); << 170 auto iZ = G4int(Z); << 171 return GetMassExcess(iA, iZ); << 172 } << 173 207 174 G4double G4NucleiProperties::GetMassExcess(con 208 G4double G4NucleiProperties::GetMassExcess(const G4int A, const G4int Z) 175 { 209 { 176 if (A < 1 || Z < 0 || Z > A) { 210 if (A < 1 || Z < 0 || Z > A) { 177 #ifdef G4VERBOSE << 211 G4cout << "G4NucleiProperties::GetMassExccess: Wrong values for A = " << A 178 if (G4ParticleTable::GetParticleTable()->G << 212 << " and Z = " << Z << G4endl; 179 G4cout << "G4NucleiProperties::GetMassEx << 180 << G4endl; << 181 } << 182 #endif << 183 return 0.0; 213 return 0.0; 184 } << 185 214 186 if (G4NucleiPropertiesTableAME12::IsInTable( << 215 } else { 187 // AME table << 188 return G4NucleiPropertiesTableAME12::GetMa << 189 } << 190 if (G4NucleiPropertiesTheoreticalTable::IsIn << 191 return G4NucleiPropertiesTheoreticalTable: << 192 } << 193 return MassExcess(A, Z); << 194 } << 195 216 196 G4double G4NucleiProperties::GetAtomicMass(con << 217 if (G4NucleiPropertiesTable::IsInTable(Z,A)){ 197 { << 218 return G4NucleiPropertiesTable::GetMassExcess(Z,A); 198 if (A < 1 || Z < 0 || Z > A) { << 219 } else if (G4NucleiPropertiesTheoreticalTable::IsInTable(Z,A)){ 199 #ifdef G4VERBOSE << 220 return G4NucleiPropertiesTheoreticalTable::GetMassExcess(Z,A); 200 if (G4ParticleTable::GetParticleTable()->G << 221 } else { 201 G4cout << "G4NucleiProperties::GetAtomic << 222 return MassExcess(A,Z); 202 << G4endl; << 203 } 223 } 204 #endif << 205 return 0.0; << 206 } << 207 if (std::fabs(A - G4int(A)) > 1.e-10) { << 208 return AtomicMass(A, Z); << 209 } 224 } 210 225 211 auto iA = G4int(A); << 212 auto iZ = G4int(Z); << 213 if (G4NucleiPropertiesTableAME12::IsInTable( << 214 return G4NucleiPropertiesTableAME12::GetAt << 215 } << 216 if (G4NucleiPropertiesTheoreticalTable::IsIn << 217 return G4NucleiPropertiesTheoreticalTable: << 218 } << 219 return AtomicMass(A, Z); << 220 } 226 } 221 227 222 G4double G4NucleiProperties::GetBindingEnergy( << 223 { << 224 auto iA = G4int(A); << 225 auto iZ = G4int(Z); << 226 return GetBindingEnergy(iA, iZ); << 227 } << 228 228 229 G4double G4NucleiProperties::GetBindingEnergy( << 229 G4double G4NucleiProperties::GetAtomicMass(const G4double A, const G4double Z) 230 { 230 { 231 if (A < 1 || Z < 0 || Z > A) { << 231 if (Z < 0 || Z > A) { 232 #ifdef G4VERBOSE << 232 G4cout << "G4NucleiProperties::GetAtomicMass: Wrong values for A = " << A 233 if (G4ParticleTable::GetParticleTable()->G << 233 << " and Z = " << Z << G4endl; return 0.0; 234 G4cout << "G4NucleiProperties::GetMassEx << 234 235 << G4endl; << 235 } else if (abs(A - G4int(A)) > 1.e-10) { >> 236 return AtomicMass(A,Z); >> 237 >> 238 } else { >> 239 G4int iA = G4int(A); >> 240 G4int iZ = G4int(Z); >> 241 if (G4NucleiPropertiesTable::IsInTable(iZ,iA)) { >> 242 return G4NucleiPropertiesTable::GetAtomicMass(iZ,iA); >> 243 } else if (G4NucleiPropertiesTheoreticalTable::IsInTable(iZ,iA)){ >> 244 return G4NucleiPropertiesTheoreticalTable::GetAtomicMass(iZ,iA); >> 245 } else { >> 246 return AtomicMass(A,Z); 236 } 247 } 237 #endif << 238 return 0.0; << 239 } << 240 << 241 if (G4NucleiPropertiesTableAME12::IsInTable( << 242 return G4NucleiPropertiesTableAME12::GetBi << 243 } 248 } 244 if (G4NucleiPropertiesTheoreticalTable::IsIn << 245 return G4NucleiPropertiesTheoreticalTable: << 246 } << 247 return BindingEnergy(A, Z); << 248 } << 249 << 250 G4double G4NucleiProperties::MassExcess(G4doub << 251 { << 252 return GetAtomicMass(A, Z) - A * amu_c2; << 253 } << 254 << 255 G4double G4NucleiProperties::AtomicMass(G4doub << 256 { << 257 G4double hydrogen_mass_excess; << 258 G4double neutron_mass_excess; << 259 hydrogen_mass_excess = G4NucleiPropertiesTab << 260 neutron_mass_excess = G4NucleiPropertiesTabl << 261 G4double mass = << 262 (A - Z) * neutron_mass_excess + Z * hydrog << 263 return mass; << 264 } 249 } 265 250 266 G4double G4NucleiProperties::NuclearMass(G4dou << 251 G4double G4NucleiProperties::GetBindingEnergy(const G4int A, const G4int Z) 267 { 252 { 268 if (A < 1 || Z < 0 || Z > A) { 253 if (A < 1 || Z < 0 || Z > A) { 269 #ifdef G4VERBOSE << 254 G4cout << "G4NucleiProperties::GetMassExccess: Wrong values for A = " << A 270 if (G4ParticleTable::GetParticleTable()->G << 255 << " and Z = " << Z << G4endl; 271 G4cout << "G4NucleiProperties::NuclearMa << 272 << G4endl; << 273 } << 274 #endif << 275 return 0.0; 256 return 0.0; 276 } << 277 257 278 G4double mass = AtomicMass(A, Z); << 258 } else { 279 << 259 if (G4NucleiPropertiesTable::IsInTable(Z,A)) { 280 // atomic mass is converted to nuclear mass << 260 return G4NucleiPropertiesTable::GetBindingEnergy(Z,A); 281 // formula in AME03 and 12 << 261 } else if (G4NucleiPropertiesTheoreticalTable::IsInTable(Z,A)) { 282 // << 262 return G4NucleiPropertiesTheoreticalTable::GetBindingEnergy(Z,A); 283 mass -= Z * electron_mass_c2; << 263 }else { 284 mass += (14.4381 * std::pow(Z, 2.39) + 1.554 << 264 return BindingEnergy(A,Z); >> 265 } 285 266 286 return mass; << 267 } 287 } 268 } 288 269 289 G4double G4NucleiProperties::BindingEnergy(G4d << 290 { << 291 // << 292 // Weitzsaecker's Mass formula << 293 // << 294 G4int Npairing = G4int(A - Z) % 2; // pairi << 295 G4int Zpairing = G4int(Z) % 2; << 296 G4double binding = -15.67 * A // nuclear vo << 297 + 17.23 * std::pow(A, 2. << 298 + 93.15 * ((A / 2. - Z) * << 299 + 0.6984523 * Z * Z * std << 300 if (Npairing == Zpairing) { << 301 binding += (Npairing + Zpairing - 1) * 12. << 302 } << 303 270 304 return -binding * MeV; << 271 G4double G4NucleiProperties::MassExcess(G4double A, G4double Z) >> 272 { >> 273 return GetAtomicMass(A,Z) - A*amu_c2; 305 } 274 } >> 275 306 276