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