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