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