Geant4-11
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4DensityEffectCalculator.cc
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26
27/*
28 * Implements calculation of the Fermi density effect as per the method
29 * described in:
30 *
31 * R. M. Sternheimer, M. J. Berger, and S. M. Seltzer. Density
32 * effect for the ionization loss of charged particles in various sub-
33 * stances. Atom. Data Nucl. Data Tabl., 30:261, 1984.
34 *
35 * Which (among other Sternheimer references) builds on:
36 *
37 * R. M. Sternheimer. The density effect for ionization loss in
38 * materials. Phys. Rev., 88:851­859, 1952.
39 *
40 * The returned values of delta are directly from the Sternheimer calculation,
41 * and not Sternheimer's popular three-part approximate parameterization
42 * introduced in the same paper.
43 *
44 * Author: Matthew Strait <straitm@umn.edu> 2019
45 */
46
48#include "G4AtomicShells.hh"
49#include "G4NistManager.hh"
50#include "G4Pow.hh"
51
53
54const G4int maxWarnings = 20;
55
57 : fMaterial(mat), fVerbose(0), fWarnings(0), nlev(n)
58{
60
61 sternf = new G4double [nlev];
62 levE = new G4double [nlev];
63 sternl = new G4double [nlev];
64 sternEbar = new G4double [nlev];
65 for(G4int i=0; i<nlev; ++i) {
66 sternf[i] = 0.0;
67 levE[i] = 0.0;
68 sternl[i] = 0.0;
69 sternEbar[i] = 0.0;
70 }
71
72 fConductivity = sternx = 0.0;
73 G4bool conductor = (fMaterial->GetFreeElectronDensity() > 0.0);
74
75 G4int sh = 0;
76 G4double sum = 0.;
78 for(size_t j = 0; j < fMaterial->GetNumberOfElements(); ++j) {
79 // The last subshell is considered to contain the conduction
80 // electrons. Sternheimer 1984 says "the lowest chemical valance of
81 // the element" is used to set the number of conduction electrons.
82 // I'm not sure if that means the highest subshell or the whole
83 // shell, but in any case, he also says that the choice is arbitrary
84 // and offers a possible alternative. This is one of the sources of
85 // uncertainty in the model.
86 const G4double frac = fMaterial->GetVecNbOfAtomsPerVolume()[j]/tot;
87 const G4int Z = fMaterial->GetElement(j)->GetZasInt();
89 for(G4int i = 0; i < nshell; ++i) {
90 // For conductors, put *all* top shell electrons into the conduction
91 // band, regardless of element.
93 if(i < nshell-1 || !conductor) {
94 sternf[sh] += xx;
95 } else {
96 fConductivity += xx;
97 }
99 ++sh;
100 }
101 }
102 for(G4int i=0; i<nlev; ++i) {
103 sum += sternf[i];
104 }
105 sum += fConductivity;
106
107 const G4double invsum = (sum > 0.0) ? 1./sum : 0.0;
108 for(G4int i=0; i<nlev; ++i) {
109 sternf[i] *= invsum;
110 }
111 fConductivity *= invsum;
114}
115
117{
118 delete [] sternf;
119 delete [] levE;
120 delete [] sternl;
121 delete [] sternEbar;
122}
123
125{
126 if(fVerbose > 1) {
127 G4cout << "G4DensityEffectCalculator::ComputeDensityCorrection for "
128 << fMaterial->GetName() << ", x= " << x << G4endl;
129 }
131 const G4double exact = FermiDeltaCalculation(x);
132
133 if(fVerbose > 1) {
134 G4cout << " Delta: computed= " << exact
135 << ", parametrized= " << approx << G4endl;
136 }
137 if(approx >= 0. && exact < 0.) {
138 if(fVerbose > 0) {
139 ++fWarnings;
140 if(fWarnings < maxWarnings) {
142 ed << "Sternheimer fit failed for " << fMaterial->GetName()
143 << ", x = " << x << ": Delta exact= "
144 << exact << ", approx= " << approx;
145 G4Exception("G4DensityEffectCalculator::DensityCorrection", "mat008",
146 JustWarning, ed);
147 }
148 }
149 return approx;
150 }
151 // Fall back to approx if exact and approx are very different, under the
152 // assumption that this means the exact calculation has gone haywire
153 // somehow, with the exception of the case where approx is negative. I
154 // have seen this clearly-wrong result occur for substances with extremely
155 // low density (1e-25 g/cc).
156 if(approx >= 0. && std::abs(exact - approx) > 1.) {
157 if(fVerbose > 0) {
158 ++fWarnings;
159 if(fWarnings < maxWarnings) {
161 ed << "Sternheimer exact= " << exact << " and approx= "
162 << approx << " are too different for "
163 << fMaterial->GetName() << ", x = " << x;
164 G4Exception("G4DensityEffectCalculator::DensityCorrection", "mat008",
165 JustWarning, ed);
166 }
167 }
168 return approx;
169 }
170 return exact;
171}
172
174{
175 // Above beta*gamma of 10^10, the exact treatment is within machine
176 // precision of the limiting case, for ordinary solids, at least. The
177 // convergence goes up as the density goes down, but even in a pretty
178 // hard vacuum it converges by 10^20. Also, it's hard to imagine how
179 // this energy is relevant (x = 20 -> 10^19 GeV for muons). So this
180 // is mostly not here for physical reasons, but rather to avoid ugly
181 // discontinuities in the return value.
182 if(x > 20.) { return -1.; }
183
184 sternx = x;
185 G4double sternrho = Newton(1.5, true);
186
187 // Negative values, and values much larger than unity are non-physical.
188 // Values between zero and one are also suspect, but not as clearly wrong.
189 if(sternrho <= 0. || sternrho > 100.) {
190 if(fVerbose > 0) {
191 ++fWarnings;
192 if(fWarnings < maxWarnings) {
194 ed << "Sternheimer computation failed for " << fMaterial->GetName()
195 << ", x = " << x << ":\n"
196 << "Could not solve for Sternheimer rho. Probably you have a \n"
197 << "mean ionization energy which is incompatible with your\n"
198 << "distribution of energy levels, or an unusually dense material.\n"
199 << "Number of levels: " << nlev
200 << " Mean ionization energy(eV): " << meanexcite
201 << " Plasma energy(eV): " << plasmaE << "\n";
202 for(G4int i = 0; i < nlev; ++i) {
203 ed << "Level " << i << ": strength " << sternf[i]
204 << ": energy(eV)= " << levE[i] << "\n";
205 }
206 G4Exception("G4DensityEffectCalculator::SetupFermiDeltaCalc", "mat008",
207 JustWarning, ed);
208 }
209 }
210 return -1.;
211 }
212
213 // Calculate the Sternheimer adjusted energy levels and parameters l_i given
214 // the Sternheimer parameter rho.
215 for(G4int i=0; i<nlev; ++i) {
216 sternEbar[i] = levE[i] * (sternrho/plasmaE);
217 sternl[i] = std::sqrt(gpow->powN(sternEbar[i], 2) + (2./3.)*sternf[i]);
218 }
219 // The derivative of the function we are solving for is strictly
220 // negative for positive (physical) values, so if the value at
221 // zero is less than zero, it has no solution, and there is no
222 // density effect in the Sternheimer "exact" treatment (which is
223 // still an approximation).
224 //
225 // For conductors, this test is not needed, because Ell(L) contains
226 // the term fConductivity/(L*L), so the value at L=0 is always
227 // positive infinity. In the code we don't return inf, though, but
228 // rather set that term to zero, which means that if this test were
229 // used, it would give the wrong result for some materials.
230 if(fConductivity == 0 && Ell(0) <= 0) return 0;
231
232 // Attempt to find the root from 40 starting points evenly distributed
233 // in log space. Trying a single starting point is not sufficient for
234 // convergence in most cases.
235 for(G4int startLi = -10; startLi < 30; ++startLi){
236 const G4double sternL = Newton(gpow->powN(2, startLi), false);
237 if(sternL != -1.) {
238 return DeltaOnceSolved(sternL);
239 }
240 }
241 return -1.; // Signal the caller to use the Sternheimer approximation,
242 // because we have been unable to solve the exact form.
243}
244
245/* Newton's method for finding roots. Adapted from G4PolynominalSolver, but
246 * without the assumption that the input is a polynomial. Also, here we
247 * always expect the roots to be positive, so return -1 as an error value. */
249{
250 const G4int maxIter = 100;
251 G4int nbad = 0, ngood = 0;
252
253 G4double lambda(start), value(0.), dvalue(0.);
254
255 if(fVerbose > 2) {
256 G4cout << "G4DensityEffectCalculator::Newton: strat= " << start
257 << " type: " << first << G4endl;
258 }
259 while(true) {
260 if(first) {
261 value = FRho(lambda);
262 dvalue = DFRho(lambda);
263 } else {
264 value = Ell(lambda);
265 dvalue = DEll(lambda);
266 }
267 if(dvalue == 0.0) { break; }
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= " << lambda << G4endl;
277 }
278 return lambda;
279 }
280 } else {
281 ++nbad;
282 }
283 if(nbad > maxIter || std::isnan(value) || std::isinf(value)) { break; }
284 }
285 if(fVerbose > 2) {
286 G4cout << " Failed to converge last value= " << value
287 << " dvalue= " << dvalue << " lambda= " << lambda << G4endl;
288 }
289 return -1.;
290}
291
292/* Return the derivative of the equation used
293 * to solve for the Sternheimer parameter rho. */
295{
296 G4double ans = 0.0;
297 for(G4int i = 0; i < nlev; ++i) {
298 if(sternf[i] > 0.) {
299 ans += sternf[i] * gpow->powN(levE[i], 2) * rho /
300 (gpow->powN(levE[i] * rho, 2)
301 + 2./3. * sternf[i] * gpow->powN(plasmaE, 2));
302 }
303 }
304 return ans;
305}
306
307/* Return the functional value for the equation used
308 * to solve for the Sternheimer parameter rho. */
310{
311 G4double ans = 0.0;
312 for(G4int i = 0; i<nlev; ++i) {
313 if(sternf[i] > 0.) {
314 ans += sternf[i] * G4Log(gpow->powN(levE[i]*rho, 2) +
315 2./3. * sternf[i]*gpow->powN(plasmaE, 2));
316 }
317 }
318 ans *= 0.5; // pulled out of loop for efficiency
319
320 if(fConductivity > 0.) {
321 ans += fConductivity * G4Log(plasmaE * std::sqrt(fConductivity));
322 }
323 ans -= G4Log(meanexcite);
324 return ans;
325}
326
327/* Return the derivative for the equation used to
328 * solve for the Sternheimer parameter l, called 'L' here. */
330{
331 G4double ans = 0.;
332 for(G4int i=0; i<nlev; ++i) {
333 if(sternf[i] > 0 && (sternEbar[i] > 0. || L != 0.)) {
334 const G4double y = gpow->powN(sternEbar[i], 2);
335 ans += sternf[i]/gpow->powN(y + L*L, 2);
336 }
337 }
338 ans += fConductivity/gpow->powN(L*L, 2);
339 ans *= (-2*L); // pulled out of the loop for efficiency
340 return ans;
341}
342
343/* Return the functional value for the equation used to
344 * solve for the Sternheimer parameter l, called 'L' here. */
346{
347 G4double ans = 0.;
348 for(G4int i=0; i<nlev; ++i) {
349 if(sternf[i] > 0. && (sternEbar[i] > 0. || L != 0.)) {
350 ans += sternf[i]/(gpow->powN(sternEbar[i], 2) + L*L);
351 }
352 }
353 if(fConductivity > 0. && L != 0.) {
354 ans += fConductivity/(L*L);
355 }
356 ans -= gpow->powZ(10, -2 * sternx);
357 return ans;
358}
359
366{
367 G4double ans = 0.;
368 for(G4int i=0; i<nlev; ++i) {
369 if(sternf[i] > 0.) {
370 ans += sternf[i] * G4Log((gpow->powN(sternl[i], 2)
371 + gpow->powN(sternL, 2))/gpow->powN(sternl[i], 2));
372 }
373 }
374 // sternl for the conduction electrons is sqrt(fConductivity), with
375 // no factor of 2./3 as with the other levels.
376 if(fConductivity > 0) {
378 + gpow->powN(sternL, 2))/fConductivity);
379 }
380 ans -= gpow->powN(sternL, 2)/(1 + gpow->powZ(10, 2 * sternx));
381 return ans;
382}
const G4int maxWarnings
static G4Pow * gpow
static const G4double eps
@ JustWarning
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:35
std::ostringstream G4ExceptionDescription
Definition: G4Exception.hh:40
G4double G4Log(G4double x)
Definition: G4Log.hh:226
static constexpr double L
Definition: G4SIunits.hh:104
double G4double
Definition: G4Types.hh:83
bool G4bool
Definition: G4Types.hh:86
int G4int
Definition: G4Types.hh:85
const G4int Z[17]
#define G4endl
Definition: G4ios.hh:57
G4GLOB_DLL std::ostream G4cout
static G4int GetNumberOfElectrons(G4int Z, G4int SubshellNb)
static G4double GetBindingEnergy(G4int Z, G4int SubshellNb)
static G4int GetNumberOfShells(G4int Z)
G4DensityEffectCalculator(const G4Material *, G4int)
G4double Newton(G4double x0, G4bool first)
G4double ComputeDensityCorrection(G4double x)
G4double FermiDeltaCalculation(G4double x)
G4int GetZasInt() const
Definition: G4Element.hh:132
G4double GetDensityCorrection(G4double x)
G4double GetMeanExcitationEnergy() const
G4double GetPlasmaEnergy() const
G4double GetTotNbOfAtomsPerVolume() const
Definition: G4Material.hh:205
const G4Element * GetElement(G4int iel) const
Definition: G4Material.hh:198
G4IonisParamMat * GetIonisation() const
Definition: G4Material.hh:222
G4double GetFreeElectronDensity() const
Definition: G4Material.hh:175
size_t GetNumberOfElements() const
Definition: G4Material.hh:182
const G4double * GetVecNbOfAtomsPerVolume() const
Definition: G4Material.hh:202
const G4String & GetName() const
Definition: G4Material.hh:173
static G4NistManager * Instance()
Definition: G4Pow.hh:49
static G4Pow * GetInstance()
Definition: G4Pow.cc:41
G4double powZ(G4int Z, G4double y) const
Definition: G4Pow.hh:225
G4double powN(G4double x, G4int n) const
Definition: G4Pow.cc:166
static constexpr double eV
T max(const T t1, const T t2)
brief Return the largest of the two arguments