| 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 26
27 /* 27 /*
28 * Implements calculation of the Fermi density 28 * Implements calculation of the Fermi density effect as per the method
29 * described in: 29 * described in:
30 * 30 *
31 * R. M. Sternheimer, M. J. Berger, and S. M 31 * R. M. Sternheimer, M. J. Berger, and S. M. Seltzer. Density
32 * effect for the ionization loss of charged 32 * effect for the ionization loss of charged particles in various sub-
33 * stances. Atom. Data Nucl. Data Tabl., 30: 33 * stances. Atom. Data Nucl. Data Tabl., 30:261, 1984.
34 * 34 *
35 * Which (among other Sternheimer references) 35 * Which (among other Sternheimer references) builds on:
36 * 36 *
37 * R. M. Sternheimer. The density effect for 37 * R. M. Sternheimer. The density effect for ionization loss in
38 * materials. Phys. Rev., 88:851859, 1952. 38 * materials. Phys. Rev., 88:851859, 1952.
39 * 39 *
40 * The returned values of delta are directly f 40 * The returned values of delta are directly from the Sternheimer calculation,
41 * and not Sternheimer's popular three-part ap 41 * and not Sternheimer's popular three-part approximate parameterization
42 * introduced in the same paper. 42 * introduced in the same paper.
43 * 43 *
44 * Author: Matthew Strait <straitm@umn.edu> 20 44 * Author: Matthew Strait <straitm@umn.edu> 2019
45 */ 45 */
46 46
47 #include "G4DensityEffectCalculator.hh" 47 #include "G4DensityEffectCalculator.hh"
48 <<
49 #include "G4AtomicShells.hh" 48 #include "G4AtomicShells.hh"
50 #include "G4NistManager.hh" 49 #include "G4NistManager.hh"
51 #include "G4Pow.hh" 50 #include "G4Pow.hh"
52 51
53 static G4Pow* gpow = G4Pow::GetInstance(); << 52 static G4Pow * gpow = G4Pow::GetInstance();
54 53
55 const G4int maxWarnings = 20; 54 const G4int maxWarnings = 20;
56 55
57 G4DensityEffectCalculator::G4DensityEffectCalc 56 G4DensityEffectCalculator::G4DensityEffectCalculator(const G4Material* mat, G4int n)
58 : fMaterial(mat), nlev(n) << 57 : fMaterial(mat), fVerbose(0), fWarnings(0), nlev(n)
59 { 58 {
60 fVerbose = std::max(fVerbose, G4NistManager: 59 fVerbose = std::max(fVerbose, G4NistManager::Instance()->GetVerbose());
61 60
62 sternf = new G4double[nlev]; << 61 sternf = new G4double [nlev];
63 levE = new G4double[nlev]; << 62 levE = new G4double [nlev];
64 sternl = new G4double[nlev]; << 63 sternl = new G4double [nlev];
65 sternEbar = new G4double[nlev]; << 64 sternEbar = new G4double [nlev];
66 for (G4int i = 0; i < nlev; ++i) { << 65 for(G4int i=0; i<nlev; ++i) {
67 sternf[i] = 0.0; << 66 sternf[i] = 0.0;
68 levE[i] = 0.0; << 67 levE[i] = 0.0;
69 sternl[i] = 0.0; << 68 sternl[i] = 0.0;
70 sternEbar[i] = 0.0; 69 sternEbar[i] = 0.0;
71 } 70 }
72 71
73 fConductivity = sternx = 0.0; 72 fConductivity = sternx = 0.0;
74 G4bool conductor = (fMaterial->GetFreeElectr 73 G4bool conductor = (fMaterial->GetFreeElectronDensity() > 0.0);
75 74
76 G4int sh = 0; 75 G4int sh = 0;
77 G4double sum = 0.; 76 G4double sum = 0.;
78 const G4double tot = fMaterial->GetTotNbOfAt 77 const G4double tot = fMaterial->GetTotNbOfAtomsPerVolume();
79 for (size_t j = 0; j < fMaterial->GetNumberO << 78 for(size_t j = 0; j < fMaterial->GetNumberOfElements(); ++j) {
80 // The last subshell is considered to cont 79 // The last subshell is considered to contain the conduction
81 // electrons. Sternheimer 1984 says "the l 80 // electrons. Sternheimer 1984 says "the lowest chemical valance of
82 // the element" is used to set the number 81 // the element" is used to set the number of conduction electrons.
83 // I'm not sure if that means the highest 82 // I'm not sure if that means the highest subshell or the whole
84 // shell, but in any case, he also says th 83 // shell, but in any case, he also says that the choice is arbitrary
85 // and offers a possible alternative. This 84 // and offers a possible alternative. This is one of the sources of
86 // uncertainty in the model. 85 // uncertainty in the model.
87 const G4double frac = fMaterial->GetVecNbO << 86 const G4double frac = fMaterial->GetVecNbOfAtomsPerVolume()[j]/tot;
88 const G4int Z = fMaterial->GetElement((G4i << 87 const G4int Z = fMaterial->GetElement(j)->GetZasInt();
89 const G4int nshell = G4AtomicShells::GetNu 88 const G4int nshell = G4AtomicShells::GetNumberOfShells(Z);
90 for (G4int i = 0; i < nshell; ++i) { << 89 for(G4int i = 0; i < nshell; ++i) {
91 // For conductors, put *all* top shell e 90 // For conductors, put *all* top shell electrons into the conduction
92 // band, regardless of element. 91 // band, regardless of element.
93 const G4double xx = frac * G4AtomicShell << 92 const G4double xx = frac*G4AtomicShells::GetNumberOfElectrons(Z, i);
94 if (i < nshell - 1 || ! conductor) { << 93 if(i < nshell-1 || !conductor) {
95 sternf[sh] += xx; << 94 sternf[sh] += xx;
96 } << 95 } else {
97 else { <<
98 fConductivity += xx; 96 fConductivity += xx;
99 } 97 }
100 levE[sh] = G4AtomicShells::GetBindingEne << 98 levE[sh] = G4AtomicShells::GetBindingEnergy(Z, i)/CLHEP::eV;
101 ++sh; 99 ++sh;
102 } 100 }
103 } 101 }
104 for (G4int i = 0; i < nlev; ++i) { << 102 for(G4int i=0; i<nlev; ++i) {
105 sum += sternf[i]; 103 sum += sternf[i];
106 } 104 }
107 sum += fConductivity; 105 sum += fConductivity;
108 106
109 const G4double invsum = (sum > 0.0) ? 1. / s << 107 const G4double invsum = (sum > 0.0) ? 1./sum : 0.0;
110 for (G4int i = 0; i < nlev; ++i) { << 108 for(G4int i=0; i<nlev; ++i) {
111 sternf[i] *= invsum; 109 sternf[i] *= invsum;
112 } 110 }
113 fConductivity *= invsum; << 111 fConductivity *= invsum;
114 plasmaE = fMaterial->GetIonisation()->GetPla << 112 plasmaE = fMaterial->GetIonisation()->GetPlasmaEnergy()/CLHEP::eV;
115 meanexcite = fMaterial->GetIonisation()->Get << 113 meanexcite = fMaterial->GetIonisation()->GetMeanExcitationEnergy()/CLHEP::eV;
116 } 114 }
117 115
118 G4DensityEffectCalculator::~G4DensityEffectCal 116 G4DensityEffectCalculator::~G4DensityEffectCalculator()
119 { 117 {
120 delete[] sternf; << 118 delete [] sternf;
121 delete[] levE; << 119 delete [] levE;
122 delete[] sternl; << 120 delete [] sternl;
123 delete[] sternEbar; << 121 delete [] sternEbar;
124 } 122 }
125 123
126 G4double G4DensityEffectCalculator::ComputeDen 124 G4double G4DensityEffectCalculator::ComputeDensityCorrection(G4double x)
127 { 125 {
128 if (fVerbose > 1) { << 126 if(fVerbose > 1) {
129 G4cout << "G4DensityEffectCalculator::Comp << 127 G4cout << "G4DensityEffectCalculator::ComputeDensityCorrection for "
130 << ", x= " << x << G4endl; << 128 << fMaterial->GetName() << ", x= " << x << G4endl;
131 } 129 }
132 const G4double approx = fMaterial->GetIonisa 130 const G4double approx = fMaterial->GetIonisation()->GetDensityCorrection(x);
133 const G4double exact = FermiDeltaCalculation << 131 const G4double exact = FermiDeltaCalculation(x);
134 132
135 if (fVerbose > 1) { << 133 if(fVerbose > 1) {
136 G4cout << " Delta: computed= " << exact << 134 G4cout << " Delta: computed= " << exact
>> 135 << ", parametrized= " << approx << G4endl;
137 } 136 }
138 if (approx >= 0. && exact < 0.) { << 137 if(approx >= 0. && exact < 0.) {
139 if (fVerbose > 0) { << 138 if(fVerbose > 0) {
140 ++fWarnings; 139 ++fWarnings;
141 if (fWarnings < maxWarnings) { << 140 if(fWarnings < maxWarnings) {
142 G4ExceptionDescription ed; << 141 G4ExceptionDescription ed;
143 ed << "Sternheimer fit failed for " << << 142 ed << "Sternheimer fit failed for " << fMaterial->GetName()
144 << ": Delta exact= " << exact << ", << 143 << ", x = " << x << ": Delta exact= "
145 G4Exception("G4DensityEffectCalculator << 144 << exact << ", approx= " << approx;
>> 145 G4Exception("G4DensityEffectCalculator::DensityCorrection", "mat008",
>> 146 JustWarning, ed);
146 } 147 }
147 } 148 }
148 return approx; 149 return approx;
149 } 150 }
150 // Fall back to approx if exact and approx a 151 // Fall back to approx if exact and approx are very different, under the
151 // assumption that this means the exact calc 152 // assumption that this means the exact calculation has gone haywire
152 // somehow, with the exception of the case w 153 // somehow, with the exception of the case where approx is negative. I
153 // have seen this clearly-wrong result occur 154 // have seen this clearly-wrong result occur for substances with extremely
154 // low density (1e-25 g/cc). 155 // low density (1e-25 g/cc).
155 if (approx >= 0. && std::abs(exact - approx) << 156 if(approx >= 0. && std::abs(exact - approx) > 1.) {
156 if (fVerbose > 0) { << 157 if(fVerbose > 0) {
157 ++fWarnings; 158 ++fWarnings;
158 if (fWarnings < maxWarnings) { << 159 if(fWarnings < maxWarnings) {
159 G4ExceptionDescription ed; << 160 G4ExceptionDescription ed;
160 ed << "Sternheimer exact= " << exact < << 161 ed << "Sternheimer exact= " << exact << " and approx= "
161 << " are too different for " << fMa << 162 << approx << " are too different for "
162 G4Exception("G4DensityEffectCalculator << 163 << fMaterial->GetName() << ", x = " << x;
>> 164 G4Exception("G4DensityEffectCalculator::DensityCorrection", "mat008",
>> 165 JustWarning, ed);
163 } 166 }
164 } 167 }
165 return approx; 168 return approx;
166 } 169 }
167 return exact; 170 return exact;
168 } 171 }
169 172
170 G4double G4DensityEffectCalculator::FermiDelta 173 G4double G4DensityEffectCalculator::FermiDeltaCalculation(G4double x)
171 { 174 {
172 // Above beta*gamma of 10^10, the exact trea 175 // Above beta*gamma of 10^10, the exact treatment is within machine
173 // precision of the limiting case, for ordin 176 // precision of the limiting case, for ordinary solids, at least. The
174 // convergence goes up as the density goes d 177 // convergence goes up as the density goes down, but even in a pretty
175 // hard vacuum it converges by 10^20. Also, 178 // hard vacuum it converges by 10^20. Also, it's hard to imagine how
176 // this energy is relevant (x = 20 -> 10^19 179 // this energy is relevant (x = 20 -> 10^19 GeV for muons). So this
177 // is mostly not here for physical reasons, 180 // is mostly not here for physical reasons, but rather to avoid ugly
178 // discontinuities in the return value. 181 // discontinuities in the return value.
179 if (x > 20.) { << 182 if(x > 20.) { return -1.; }
180 return -1.; <<
181 } <<
182 183
183 sternx = x; 184 sternx = x;
184 G4double sternrho = Newton(1.5, true); 185 G4double sternrho = Newton(1.5, true);
185 186
186 // Negative values, and values much larger t 187 // Negative values, and values much larger than unity are non-physical.
187 // Values between zero and one are also susp 188 // Values between zero and one are also suspect, but not as clearly wrong.
188 if (sternrho <= 0. || sternrho > 100.) { << 189 if(sternrho <= 0. || sternrho > 100.) {
189 if (fVerbose > 0) { << 190 if(fVerbose > 0) {
190 ++fWarnings; 191 ++fWarnings;
191 if (fWarnings < maxWarnings) { << 192 if(fWarnings < maxWarnings) {
192 G4ExceptionDescription ed; << 193 G4ExceptionDescription ed;
193 ed << "Sternheimer computation failed << 194 ed << "Sternheimer computation failed for " << fMaterial->GetName()
194 << ":\n" << 195 << ", x = " << x << ":\n"
195 << "Could not solve for Sternheimer << 196 << "Could not solve for Sternheimer rho. Probably you have a \n"
196 << "mean ionization energy which is << 197 << "mean ionization energy which is incompatible with your\n"
197 << "distribution of energy levels, << 198 << "distribution of energy levels, or an unusually dense material.\n"
198 << "Number of levels: " << nlev << << 199 << "Number of levels: " << nlev
199 << " Plasma energy(eV): " << plasma << 200 << " Mean ionization energy(eV): " << meanexcite
200 for (G4int i = 0; i < nlev; ++i) { << 201 << " Plasma energy(eV): " << plasmaE << "\n";
201 ed << "Level " << i << ": strength " << 202 for(G4int i = 0; i < nlev; ++i) {
202 } << 203 ed << "Level " << i << ": strength " << sternf[i]
203 G4Exception("G4DensityEffectCalculator << 204 << ": energy(eV)= " << levE[i] << "\n";
>> 205 }
>> 206 G4Exception("G4DensityEffectCalculator::SetupFermiDeltaCalc", "mat008",
>> 207 JustWarning, ed);
204 } 208 }
205 } 209 }
206 return -1.; 210 return -1.;
207 } 211 }
208 212
209 // Calculate the Sternheimer adjusted energy 213 // Calculate the Sternheimer adjusted energy levels and parameters l_i given
210 // the Sternheimer parameter rho. 214 // the Sternheimer parameter rho.
211 for (G4int i = 0; i < nlev; ++i) { << 215 for(G4int i=0; i<nlev; ++i) {
212 sternEbar[i] = levE[i] * (sternrho / plasm << 216 sternEbar[i] = levE[i] * (sternrho/plasmaE);
213 sternl[i] = std::sqrt(gpow->powN(sternEbar << 217 sternl[i] = std::sqrt(gpow->powN(sternEbar[i], 2) + (2./3.)*sternf[i]);
214 } 218 }
215 // The derivative of the function we are sol 219 // The derivative of the function we are solving for is strictly
216 // negative for positive (physical) values, 220 // negative for positive (physical) values, so if the value at
217 // zero is less than zero, it has no solutio 221 // zero is less than zero, it has no solution, and there is no
218 // density effect in the Sternheimer "exact" 222 // density effect in the Sternheimer "exact" treatment (which is
219 // still an approximation). 223 // still an approximation).
220 // 224 //
221 // For conductors, this test is not needed, 225 // For conductors, this test is not needed, because Ell(L) contains
222 // the term fConductivity/(L*L), so the valu 226 // the term fConductivity/(L*L), so the value at L=0 is always
223 // positive infinity. In the code we don't r 227 // positive infinity. In the code we don't return inf, though, but
224 // rather set that term to zero, which means 228 // rather set that term to zero, which means that if this test were
225 // used, it would give the wrong result for 229 // used, it would give the wrong result for some materials.
226 if (fConductivity == 0 && Ell(0) <= 0) { << 230 if(fConductivity == 0 && Ell(0) <= 0) return 0;
227 return 0; <<
228 } <<
229 231
230 // Attempt to find the root from 40 starting 232 // Attempt to find the root from 40 starting points evenly distributed
231 // in log space. Trying a single starting p 233 // in log space. Trying a single starting point is not sufficient for
232 // convergence in most cases. 234 // convergence in most cases.
233 for (G4int startLi = -10; startLi < 30; ++st << 235 for(G4int startLi = -10; startLi < 30; ++startLi){
234 const G4double sternL = Newton(gpow->powN( 236 const G4double sternL = Newton(gpow->powN(2, startLi), false);
235 if (sternL != -1.) { << 237 if(sternL != -1.) {
236 return DeltaOnceSolved(sternL); 238 return DeltaOnceSolved(sternL);
237 } 239 }
238 } 240 }
239 return -1.; // Signal the caller to use the << 241 return -1.; // Signal the caller to use the Sternheimer approximation,
240 // because we have been unable << 242 // because we have been unable to solve the exact form.
241 } 243 }
242 244
243 /* Newton's method for finding roots. Adapted 245 /* Newton's method for finding roots. Adapted from G4PolynominalSolver, but
244 * without the assumption that the input is a 246 * without the assumption that the input is a polynomial. Also, here we
245 * always expect the roots to be positive, so 247 * always expect the roots to be positive, so return -1 as an error value. */
246 G4double G4DensityEffectCalculator::Newton(G4d 248 G4double G4DensityEffectCalculator::Newton(G4double start, G4bool first)
247 { 249 {
248 const G4int maxIter = 100; 250 const G4int maxIter = 100;
249 G4int nbad = 0, ngood = 0; 251 G4int nbad = 0, ngood = 0;
250 252
251 G4double lambda(start), value(0.), dvalue(0. 253 G4double lambda(start), value(0.), dvalue(0.);
252 254
253 if (fVerbose > 2) { << 255 if(fVerbose > 2) {
254 G4cout << "G4DensityEffectCalculator::Newt << 256 G4cout << "G4DensityEffectCalculator::Newton: strat= " << start
>> 257 << " type: " << first << G4endl;
255 } 258 }
256 while (true) { << 259 while(true) {
257 if (first) { << 260 if(first) {
258 value = FRho(lambda); 261 value = FRho(lambda);
259 dvalue = DFRho(lambda); 262 dvalue = DFRho(lambda);
260 } << 263 } else {
261 else { <<
262 value = Ell(lambda); 264 value = Ell(lambda);
263 dvalue = DEll(lambda); 265 dvalue = DEll(lambda);
264 } 266 }
265 if (dvalue == 0.0) { << 267 if(dvalue == 0.0) { break; }
266 break; << 268 const G4double del = value/dvalue;
267 } <<
268 const G4double del = value / dvalue; <<
269 lambda -= del; 269 lambda -= del;
270 270
271 const G4double eps = std::abs(del / lambda << 271 const G4double eps = std::abs(del/lambda);
272 if (eps <= 1.e-12) { << 272 if(eps <= 1.e-12) {
273 ++ngood; 273 ++ngood;
274 if (ngood == 2) { << 274 if(ngood == 2) {
275 if (fVerbose > 2) { << 275 if(fVerbose > 2) {
276 G4cout << " Converged with result= << 276 G4cout << " Converged with result= " << lambda << G4endl;
277 } << 277 }
278 return lambda; << 278 return lambda;
279 } 279 }
280 } << 280 } else {
281 else { <<
282 ++nbad; 281 ++nbad;
283 } 282 }
284 if (nbad > maxIter || std::isnan(value) || << 283 if(nbad > maxIter || std::isnan(value) || std::isinf(value)) { break; }
285 break; <<
286 } <<
287 } 284 }
288 if (fVerbose > 2) { << 285 if(fVerbose > 2) {
289 G4cout << " Failed to converge last value << 286 G4cout << " Failed to converge last value= " << value
290 << " lambda= " << lambda << G4endl; << 287 << " dvalue= " << dvalue << " lambda= " << lambda << G4endl;
291 } 288 }
292 return -1.; 289 return -1.;
293 } 290 }
294 291
295 /* Return the derivative of the equation used 292 /* Return the derivative of the equation used
296 * to solve for the Sternheimer parameter rho. 293 * to solve for the Sternheimer parameter rho. */
297 G4double G4DensityEffectCalculator::DFRho(G4do 294 G4double G4DensityEffectCalculator::DFRho(G4double rho)
298 { 295 {
299 G4double ans = 0.0; 296 G4double ans = 0.0;
300 for (G4int i = 0; i < nlev; ++i) { << 297 for(G4int i = 0; i < nlev; ++i) {
301 if (sternf[i] > 0.) { << 298 if(sternf[i] > 0.) {
302 ans += sternf[i] * gpow->powN(levE[i], 2 299 ans += sternf[i] * gpow->powN(levE[i], 2) * rho /
303 (gpow->powN(levE[i] * rho, 2) + 2 << 300 (gpow->powN(levE[i] * rho, 2)
>> 301 + 2./3. * sternf[i] * gpow->powN(plasmaE, 2));
304 } 302 }
305 } 303 }
306 return ans; 304 return ans;
307 } 305 }
308 306
309 /* Return the functional value for the equatio 307 /* Return the functional value for the equation used
310 * to solve for the Sternheimer parameter rho. 308 * to solve for the Sternheimer parameter rho. */
311 G4double G4DensityEffectCalculator::FRho(G4dou 309 G4double G4DensityEffectCalculator::FRho(G4double rho)
312 { 310 {
313 G4double ans = 0.0; 311 G4double ans = 0.0;
314 for (G4int i = 0; i < nlev; ++i) { << 312 for(G4int i = 0; i<nlev; ++i) {
315 if (sternf[i] > 0.) { << 313 if(sternf[i] > 0.) {
316 ans += sternf[i] * << 314 ans += sternf[i] * G4Log(gpow->powN(levE[i]*rho, 2) +
317 G4Log(gpow->powN(levE[i] * rho, 2 << 315 2./3. * sternf[i]*gpow->powN(plasmaE, 2));
318 } << 316 }
319 } 317 }
320 ans *= 0.5; // pulled out of loop for effic << 318 ans *= 0.5; // pulled out of loop for efficiency
321 319
322 if (fConductivity > 0.) { << 320 if(fConductivity > 0.) {
323 ans += fConductivity * G4Log(plasmaE * std 321 ans += fConductivity * G4Log(plasmaE * std::sqrt(fConductivity));
324 } 322 }
325 ans -= G4Log(meanexcite); 323 ans -= G4Log(meanexcite);
326 return ans; 324 return ans;
327 } 325 }
328 326
329 /* Return the derivative for the equation used 327 /* Return the derivative for the equation used to
330 * solve for the Sternheimer parameter l, call 328 * solve for the Sternheimer parameter l, called 'L' here. */
331 G4double G4DensityEffectCalculator::DEll(G4dou 329 G4double G4DensityEffectCalculator::DEll(G4double L)
332 { 330 {
333 G4double ans = 0.; 331 G4double ans = 0.;
334 for (G4int i = 0; i < nlev; ++i) { << 332 for(G4int i=0; i<nlev; ++i) {
335 if (sternf[i] > 0 && (sternEbar[i] > 0. || << 333 if(sternf[i] > 0 && (sternEbar[i] > 0. || L != 0.)) {
336 const G4double y = gpow->powN(sternEbar[ 334 const G4double y = gpow->powN(sternEbar[i], 2);
337 ans += sternf[i] / gpow->powN(y + L * L, << 335 ans += sternf[i]/gpow->powN(y + L*L, 2);
338 } 336 }
339 } 337 }
340 ans += fConductivity / gpow->powN(L * L, 2); << 338 ans += fConductivity/gpow->powN(L*L, 2);
341 ans *= (-2 * L); // pulled out of the loop << 339 ans *= (-2*L); // pulled out of the loop for efficiency
342 return ans; 340 return ans;
343 } 341 }
344 342
345 /* Return the functional value for the equatio 343 /* Return the functional value for the equation used to
346 * solve for the Sternheimer parameter l, call 344 * solve for the Sternheimer parameter l, called 'L' here. */
347 G4double G4DensityEffectCalculator::Ell(G4doub 345 G4double G4DensityEffectCalculator::Ell(G4double L)
348 { 346 {
349 G4double ans = 0.; 347 G4double ans = 0.;
350 for (G4int i = 0; i < nlev; ++i) { << 348 for(G4int i=0; i<nlev; ++i) {
351 if (sternf[i] > 0. && (sternEbar[i] > 0. | << 349 if(sternf[i] > 0. && (sternEbar[i] > 0. || L != 0.)) {
352 ans += sternf[i] / (gpow->powN(sternEbar << 350 ans += sternf[i]/(gpow->powN(sternEbar[i], 2) + L*L);
353 } 351 }
354 } 352 }
355 if (fConductivity > 0. && L != 0.) { << 353 if(fConductivity > 0. && L != 0.) {
356 ans += fConductivity / (L * L); << 354 ans += fConductivity/(L*L);
357 } << 355 }
358 ans -= gpow->powZ(10, -2 * sternx); 356 ans -= gpow->powZ(10, -2 * sternx);
359 return ans; 357 return ans;
360 } 358 }
361 359
362 /** 360 /**
363 * Given the Sternheimer parameter l (called ' 361 * Given the Sternheimer parameter l (called 'sternL' here), and that
364 * the l_i and adjusted energies have been fou 362 * the l_i and adjusted energies have been found with SetupFermiDeltaCalc(),
365 * return the value of delta. Helper function 363 * return the value of delta. Helper function for DoFermiDeltaCalc().
366 */ 364 */
367 G4double G4DensityEffectCalculator::DeltaOnceS 365 G4double G4DensityEffectCalculator::DeltaOnceSolved(G4double sternL)
368 { 366 {
369 G4double ans = 0.; 367 G4double ans = 0.;
370 for (G4int i = 0; i < nlev; ++i) { << 368 for(G4int i=0; i<nlev; ++i) {
371 if (sternf[i] > 0.) { << 369 if(sternf[i] > 0.) {
372 ans += sternf[i] * << 370 ans += sternf[i] * G4Log((gpow->powN(sternl[i], 2)
373 G4Log((gpow->powN(sternl[i], 2) + << 371 + gpow->powN(sternL, 2))/gpow->powN(sternl[i], 2));
374 } 372 }
375 } 373 }
376 // sternl for the conduction electrons is sq 374 // sternl for the conduction electrons is sqrt(fConductivity), with
377 // no factor of 2./3 as with the other level 375 // no factor of 2./3 as with the other levels.
378 if (fConductivity > 0) { << 376 if(fConductivity > 0) {
379 ans += fConductivity * G4Log((fConductivit << 377 ans += fConductivity * G4Log((fConductivity
>> 378 + gpow->powN(sternL, 2))/fConductivity);
380 } 379 }
381 ans -= gpow->powN(sternL, 2) / (1 + gpow->po << 380 ans -= gpow->powN(sternL, 2)/(1 + gpow->powZ(10, 2 * sternx));
382 return ans; 381 return ans;
383 } 382 }
384 383