| 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.18 2008/11/06 13:17:36 kurasige Exp $
29 // History: << 28 // GEANT4 tag $Name: geant4-09-02 $
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)
>> 41
33 42
34 #include "G4NucleiProperties.hh" 43 #include "G4NucleiProperties.hh"
35 44
36 #include "G4NucleiPropertiesTableAME12.hh" << 45 G4bool G4NucleiProperties::isIntialized = false;
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 46
49 G4double G4NucleiProperties::GetNuclearMass(co << 47 G4double G4NucleiProperties::mass_proton = -1.;
>> 48 G4double G4NucleiProperties::mass_neutron = -1.;
>> 49 G4double G4NucleiProperties::mass_deuteron = -1.;
>> 50 G4double G4NucleiProperties::mass_triton = -1.;
>> 51 G4double G4NucleiProperties::mass_alpha = -1.;
>> 52 G4double G4NucleiProperties::mass_He3 = -1.;
>> 53 G4double G4NucleiProperties::electronMass[MaxZ];
>> 54
>> 55 G4double G4NucleiProperties::AtomicMass(G4double A, G4double Z)
50 { 56 {
51 G4double mass = 0.0; << 57 const G4double hydrogen_mass_excess = G4NucleiPropertiesTable::GetMassExcess(1,1);
>> 58 const G4double neutron_mass_excess = G4NucleiPropertiesTable::GetMassExcess(0,1);
>> 59
>> 60 G4double mass =
>> 61 (A-Z)*neutron_mass_excess + Z*hydrogen_mass_excess - BindingEnergy(A,Z) + A*amu_c2;
52 62
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 return mass;
63 } 64 }
64 65
65 G4double G4NucleiProperties::GetNuclearMass(co << 66 G4double G4NucleiProperties::BindingEnergy(G4double A, G4double Z)
66 { << 67 {
67 if (mass_proton <= 0.0) { << 68 //
68 const G4ParticleDefinition* nucleus = null << 69 // Weitzsaecker's Mass formula
69 nucleus = G4ParticleTable::GetParticleTabl << 70 //
70 if (nucleus != nullptr) mass_neutron = nuc << 71 G4int Npairing = G4int(A-Z)%2; // pairing
71 << 72 G4int Zpairing = G4int(Z)%2;
72 nucleus = G4ParticleTable::GetParticleTabl << 73 G4double binding =
73 if (nucleus != nullptr) mass_deuteron = nu << 74 - 15.67*A // nuclear volume
74 << 75 + 17.23*std::pow(A,2./3.) // surface energy
75 nucleus = G4ParticleTable::GetParticleTabl << 76 + 93.15*((A/2.-Z)*(A/2.-Z))/A // asymmetry
76 if (nucleus != nullptr) mass_triton = nucl << 77 + 0.6984523*Z*Z*std::pow(A,-1./3.); // coulomb
>> 78 if( Npairing == Zpairing ) binding += (Npairing+Zpairing-1) * 12.0 / std::sqrt(A); // pairing
77 79
78 nucleus = G4ParticleTable::GetParticleTabl << 80 return -binding*MeV;
79 if (nucleus != nullptr) mass_alpha = nucle << 81 }
80 82
81 nucleus = G4ParticleTable::GetParticleTabl <<
82 if (nucleus != nullptr) mass_He3 = nucleus <<
83 83
84 nucleus = G4ParticleTable::GetParticleTabl << 84 G4double G4NucleiProperties::GetNuclearMass(const G4double A, const G4double Z)
85 if (nucleus != nullptr) mass_proton = nucl << 85 {
>> 86 if (!isIntialized) {
>> 87 isIntialized = true;
>> 88 G4ParticleDefinition * nucleus = 0;
>> 89 nucleus = G4ParticleTable::GetParticleTable()->FindParticle("proton"); // proton
>> 90 if (nucleus!=0) mass_proton = nucleus->GetPDGMass();
>> 91 nucleus = G4ParticleTable::GetParticleTable()->FindParticle("neutron"); // neutron
>> 92 if (nucleus!=0) mass_neutron = nucleus->GetPDGMass();
>> 93 nucleus = G4ParticleTable::GetParticleTable()->FindParticle("deuteron"); // deuteron
>> 94 if (nucleus!=0) mass_deuteron = nucleus->GetPDGMass();
>> 95 nucleus = G4ParticleTable::GetParticleTable()->FindParticle("triton"); // triton
>> 96 if (nucleus!=0) mass_triton = nucleus->GetPDGMass();
>> 97 nucleus = G4ParticleTable::GetParticleTable()->FindParticle("alpha"); // alpha
>> 98 if (nucleus!=0) mass_alpha = nucleus->GetPDGMass();
>> 99 nucleus = G4ParticleTable::GetParticleTable()->FindParticle("He3"); // He3
>> 100 if (nucleus!=0) mass_He3 = nucleus->GetPDGMass();
>> 101
>> 102 for (int iz=1; iz<MaxZ; iz+=1){
>> 103 electronMass[iz] = iz*electron_mass_c2
>> 104 - 1.433e-5*MeV*std::pow(G4double(iz),2.39); ;
>> 105 }
86 } 106 }
87 107
88 if (A < 1 || Z < 0 || Z > A) { 108 if (A < 1 || Z < 0 || Z > A) {
89 #ifdef G4VERBOSE 109 #ifdef G4VERBOSE
90 if (G4ParticleTable::GetParticleTable()->G << 110 if (G4ParticleTable::GetParticleTable()->GetVerboseLevel()>0) {
91 G4cout << "G4NucleiProperties::GetNuclea << 111 G4cout << "G4NucleiProperties::GetNuclearMass: Wrong values for A = " << A
92 << G4endl; << 112 << " and Z = " << Z << G4endl;
93 } 113 }
94 #endif << 114 #endif
95 return 0.0; 115 return 0.0;
96 } << 116 } else {
97 << 117 G4double mass= -1.;
98 G4double mass = -1.; << 118 if ( (Z<=2) ) {
99 if ((Z <= 2)) { << 119 if ( (Z==1)&&(A==1) ) {
100 // light nuclei << 120 mass = mass_proton;
101 if ((Z == 1) && (A == 1)) { << 121 } else if ( (Z==0)&&(A==1) ) {
102 mass = mass_proton; << 122 mass = mass_neutron;
103 } << 123 } else if ( (Z==1)&&(A==2) ) {
104 else if ((Z == 0) && (A == 1)) { << 124 mass = mass_deuteron;
105 mass = mass_neutron; << 125 } else if ( (Z==1)&&(A==3) ) {
106 } << 126 mass = mass_triton;
107 else if ((Z == 1) && (A == 2)) { << 127 } else if ( (Z==2)&&(A==4) ) {
108 mass = mass_deuteron; << 128 mass = mass_alpha;
109 } << 129 } else if ( (Z==2)&&(A==3) ) {
110 else if ((Z == 1) && (A == 3)) { << 130 mass = mass_He3;
111 mass = mass_triton; << 131 }
112 } << 132 }
113 else if ((Z == 2) && (A == 4)) { << 133 if (mass < 0.) {
114 mass = mass_alpha; << 134 if (Z >= MaxZ) {
115 } << 135 mass = GetAtomicMass(A,Z) - Z*electron_mass_c2 + 1.433e-5*MeV*std::pow(Z,2.39);
116 else if ((Z == 2) && (A == 3)) { << 136 } else {
117 mass = mass_He3; << 137 mass = GetAtomicMass(A,Z) - electronMass[G4int(Z)];
118 } << 138 }
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 if (mass < 0.) mass = 0.0;
>> 141 return mass;
139 } 142 }
140 <<
141 if (mass < 0.) mass = 0.0; <<
142 return mass; <<
143 } 143 }
144 144
145 G4bool G4NucleiProperties::IsInStableTable(con 145 G4bool G4NucleiProperties::IsInStableTable(const G4double A, const G4double Z)
146 { 146 {
147 auto iA = G4int(A); << 147 if (Z < 0 || Z > A) {
148 auto iZ = G4int(Z); <<
149 return IsInStableTable(iA, iZ); <<
150 } <<
151 <<
152 G4bool G4NucleiProperties::IsInStableTable(con <<
153 { <<
154 if (A < 1 || Z < 0 || Z > A) { <<
155 #ifdef G4VERBOSE 148 #ifdef G4VERBOSE
156 if (G4ParticleTable::GetParticleTable()->G << 149 if (G4ParticleTable::GetParticleTable()->GetVerboseLevel()>0) {
157 G4cout << "G4NucleiProperties::IsInStabl << 150 G4cout << "G4NucleiProperties::IsInStableTable: Wrong values for A = "
158 << " and Z = " << Z << G4endl; << 151 << A << " and Z = " << Z << G4endl;
159 } 152 }
160 #endif << 153 #endif
161 return false; 154 return false;
162 } <<
163 155
164 return G4NucleiPropertiesTableAME12::IsInTab << 156 } else {
>> 157 G4int iA = G4int(A);
>> 158 G4int iZ = G4int(Z);
>> 159 return G4NucleiPropertiesTable::IsInTable(iZ,iA);
>> 160 }
165 } 161 }
166 162
167 G4double G4NucleiProperties::GetMassExcess(con 163 G4double G4NucleiProperties::GetMassExcess(const G4double A, const G4double Z)
168 { 164 {
169 auto iA = G4int(A); << 165 G4int iA = G4int(A);
170 auto iZ = G4int(Z); << 166 G4int iZ = G4int(Z);
171 return GetMassExcess(iA, iZ); << 167 return GetMassExcess(iA,iZ);
172 } 168 }
173 169
174 G4double G4NucleiProperties::GetMassExcess(con 170 G4double G4NucleiProperties::GetMassExcess(const G4int A, const G4int Z)
175 { 171 {
176 if (A < 1 || Z < 0 || Z > A) { 172 if (A < 1 || Z < 0 || Z > A) {
177 #ifdef G4VERBOSE 173 #ifdef G4VERBOSE
178 if (G4ParticleTable::GetParticleTable()->G << 174 if (G4ParticleTable::GetParticleTable()->GetVerboseLevel()>0) {
179 G4cout << "G4NucleiProperties::GetMassEx << 175 G4cout << "G4NucleiProperties::GetMassExccess: Wrong values for A = "
180 << G4endl; << 176 << A << " and Z = " << Z << G4endl;
181 } 177 }
182 #endif << 178 #endif
183 return 0.0; 179 return 0.0;
184 } << 180
>> 181 } else {
185 182
186 if (G4NucleiPropertiesTableAME12::IsInTable( << 183 if (G4NucleiPropertiesTable::IsInTable(Z,A)){
187 // AME table << 184 return G4NucleiPropertiesTable::GetMassExcess(Z,A);
188 return G4NucleiPropertiesTableAME12::GetMa << 185 } else if (G4NucleiPropertiesTheoreticalTable::IsInTable(Z,A)){
189 } << 186 return G4NucleiPropertiesTheoreticalTable::GetMassExcess(Z,A);
190 if (G4NucleiPropertiesTheoreticalTable::IsIn << 187 } else {
191 return G4NucleiPropertiesTheoreticalTable: << 188 return MassExcess(A,Z);
>> 189 }
192 } 190 }
193 return MassExcess(A, Z); << 191
194 } 192 }
195 193
>> 194
196 G4double G4NucleiProperties::GetAtomicMass(con 195 G4double G4NucleiProperties::GetAtomicMass(const G4double A, const G4double Z)
197 { 196 {
198 if (A < 1 || Z < 0 || Z > A) { << 197 if (Z < 0 || Z > A) {
199 #ifdef G4VERBOSE 198 #ifdef G4VERBOSE
200 if (G4ParticleTable::GetParticleTable()->G << 199 if (G4ParticleTable::GetParticleTable()->GetVerboseLevel()>0) {
201 G4cout << "G4NucleiProperties::GetAtomic << 200 G4cout << "G4NucleiProperties::GetAtomicMass: Wrong values for A = "
202 << G4endl; << 201 << A << " and Z = " << Z << G4endl;
203 } 202 }
204 #endif << 203 #endif
205 return 0.0; 204 return 0.0;
206 } <<
207 if (std::fabs(A - G4int(A)) > 1.e-10) { <<
208 return AtomicMass(A, Z); <<
209 } <<
210 205
211 auto iA = G4int(A); << 206 } else if (std::abs(A - G4int(A)) > 1.e-10) {
212 auto iZ = G4int(Z); << 207 return AtomicMass(A,Z);
213 if (G4NucleiPropertiesTableAME12::IsInTable( << 208
214 return G4NucleiPropertiesTableAME12::GetAt << 209 } else {
215 } << 210 G4int iA = G4int(A);
216 if (G4NucleiPropertiesTheoreticalTable::IsIn << 211 G4int iZ = G4int(Z);
217 return G4NucleiPropertiesTheoreticalTable: << 212 if (G4NucleiPropertiesTable::IsInTable(iZ,iA)) {
>> 213 return G4NucleiPropertiesTable::GetAtomicMass(iZ,iA);
>> 214 } else if (G4NucleiPropertiesTheoreticalTable::IsInTable(iZ,iA)){
>> 215 return G4NucleiPropertiesTheoreticalTable::GetAtomicMass(iZ,iA);
>> 216 } else {
>> 217 return AtomicMass(A,Z);
>> 218 }
218 } 219 }
219 return AtomicMass(A, Z); <<
220 } 220 }
221 221
222 G4double G4NucleiProperties::GetBindingEnergy( 222 G4double G4NucleiProperties::GetBindingEnergy(const G4double A, const G4double Z)
223 { 223 {
224 auto iA = G4int(A); << 224 G4int iA = G4int(A);
225 auto iZ = G4int(Z); << 225 G4int iZ = G4int(Z);
226 return GetBindingEnergy(iA, iZ); << 226 return GetBindingEnergy(iA,iZ);
227 } 227 }
228 228
229 G4double G4NucleiProperties::GetBindingEnergy( 229 G4double G4NucleiProperties::GetBindingEnergy(const G4int A, const G4int Z)
230 { 230 {
231 if (A < 1 || Z < 0 || Z > A) { 231 if (A < 1 || Z < 0 || Z > A) {
232 #ifdef G4VERBOSE 232 #ifdef G4VERBOSE
233 if (G4ParticleTable::GetParticleTable()->G << 233 if (G4ParticleTable::GetParticleTable()->GetVerboseLevel()>0) {
234 G4cout << "G4NucleiProperties::GetMassEx << 234 G4cout << "G4NucleiProperties::GetMassExccess: Wrong values for A = "
235 << G4endl; << 235 << A << " and Z = " << Z << G4endl;
236 } 236 }
237 #endif 237 #endif
238 return 0.0; 238 return 0.0;
239 } <<
240 <<
241 if (G4NucleiPropertiesTableAME12::IsInTable( <<
242 return G4NucleiPropertiesTableAME12::GetBi <<
243 } <<
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 } <<
265 239
266 G4double G4NucleiProperties::NuclearMass(G4dou << 240 } else {
267 { << 241 if (G4NucleiPropertiesTable::IsInTable(Z,A)) {
268 if (A < 1 || Z < 0 || Z > A) { << 242 return G4NucleiPropertiesTable::GetBindingEnergy(Z,A);
269 #ifdef G4VERBOSE << 243 } else if (G4NucleiPropertiesTheoreticalTable::IsInTable(Z,A)) {
270 if (G4ParticleTable::GetParticleTable()->G << 244 return G4NucleiPropertiesTheoreticalTable::GetBindingEnergy(Z,A);
271 G4cout << "G4NucleiProperties::NuclearMa << 245 }else {
272 << G4endl; << 246 return BindingEnergy(A,Z);
273 } 247 }
274 #endif <<
275 return 0.0; <<
276 } <<
277 <<
278 G4double mass = AtomicMass(A, Z); <<
279 248
280 // atomic mass is converted to nuclear mass << 249 }
281 // formula in AME03 and 12 <<
282 // <<
283 mass -= Z * electron_mass_c2; <<
284 mass += (14.4381 * std::pow(Z, 2.39) + 1.554 <<
285 <<
286 return mass; <<
287 } 250 }
288 251
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 252
304 return -binding * MeV; << 253 G4double G4NucleiProperties::MassExcess(G4double A, G4double Z)
>> 254 {
>> 255 return GetAtomicMass(A,Z) - A*amu_c2;
305 } 256 }
>> 257
306 258