ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DensityEffectCalculator.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4DensityEffectCalculator.cc
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 
54 const 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();
88  const G4int nshell = G4AtomicShells::GetNumberOfShells(Z);
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 = (sum > 0.0) ? 1./sum : 0.0;
106  for(G4int i=0; i<nlev; ++i) {
107  sternf[i] *= sum;
108  }
111 }
112 
114 {
115  delete [] sternf;
116  delete [] levE;
117  delete [] sternl;
118  delete [] sternEbar;
119 }
120 
122 {
123  if(fVerbose > 1) {
124  G4cout << "G4DensityEffectCalculator::ComputeDensityCorrection for "
125  << fMaterial->GetName() << ", x= " << x << G4endl;
126  }
128  const G4double exact = FermiDeltaCalculation(x);
129 
130  if(fVerbose > 1) {
131  G4cout << " Delta: computed= " << exact
132  << ", parametrized= " << approx << G4endl;
133  }
134  if(approx > 0. && exact < 0.) {
135  if(fVerbose > 0) {
136  ++fWarnings;
137  if(fWarnings < maxWarnings) {
139  ed << "Sternheimer fit failed for " << fMaterial->GetName()
140  << ", x = " << x << ": Delta exact= "
141  << exact << ", approx= " << approx;
142  G4Exception("G4DensityEffectCalculator::DensityCorrection", "mat008",
143  JustWarning, ed);
144  }
145  }
146  return approx;
147  }
148  // Fall back to approx if exact and approx are very different, under the
149  // assumption that this means the exact calculation has gone haywire
150  // somehow, with the exception of the case where approx is negative. I
151  // have seen this clearly-wrong result occur for substances with extremely
152  // low density (1e-25 g/cc).
153  if(approx >= 0. && std::abs(exact - approx) > 1.) {
154  if(fVerbose > 0) {
155  ++fWarnings;
156  if(fWarnings < maxWarnings) {
158  ed << "Sternheimer exact= " << exact << " and approx= "
159  << approx << " are too different for "
160  << fMaterial->GetName() << ", x = " << x;
161  G4Exception("G4DensityEffectCalculator::DensityCorrection", "mat008",
162  JustWarning, ed);
163  }
164  }
165  return approx;
166  }
167  return exact;
168 }
169 
171 {
172  // Above beta*gamma of 10^10, the exact treatment is within machine
173  // precision of the limiting case, for ordinary solids, at least. The
174  // convergence goes up as the density goes down, but even in a pretty
175  // hard vacuum it converges by 10^20. Also, it's hard to imagine how
176  // this energy is relevant (x = 20 -> 10^19 GeV for muons). So this
177  // is mostly not here for physical reasons, but rather to avoid ugly
178  // discontinuities in the return value.
179  if(x > 20.) { return -1.; }
180 
181  sternx = x;
182  G4double sternrho = Newton(1.5, true);
183 
184  // Negative values, and values much larger than unity are non-physical.
185  // Values between zero and one are also suspect, but not as clearly wrong.
186  if(sternrho <= 0. || sternrho > 100.) {
187  if(fVerbose > 0) {
188  ++fWarnings;
189  if(fWarnings < maxWarnings) {
191  ed << "Sternheimer computation failed for " << fMaterial->GetName()
192  << ", x = " << x << ":\n"
193  << "Could not solve for Sternheimer rho. Probably you have a \n"
194  << "mean ionization energy which is incompatible with your\n"
195  << "distribution of energy levels, or an unusually dense material.\n"
196  << "Number of levels: " << nlev
197  << " Mean ionization energy(eV): " << meanexcite
198  << " Plasma energy(eV): " << plasmaE << "\n";
199  for(G4int i = 0; i < nlev; ++i) {
200  ed << "Level " << i << ": strength " << sternf[i]
201  << ": energy(eV)= " << levE[i] << "\n";
202  }
203  G4Exception("G4DensityEffectCalculator::SetupFermiDeltaCalc", "mat008",
204  JustWarning, ed);
205  }
206  }
207  return -1.;
208  }
209 
210  // Calculate the Sternheimer adjusted energy levels and parameters l_i given
211  // the Sternheimer parameter rho.
212  sternrho /= plasmaE;
213  for(G4int i=0; i<nlev; ++i) {
214  sternEbar[i] = levE[i] * sternrho;
215  sternl[i] = std::sqrt(gpow->powN(sternEbar[i], 2) + 2./3.*sternf[i]);
216  }
217 
218  // Make imphirical initial guess
219  const G4double sternL = Newton(sternrho, false);
220  if(sternL > -1.) {
221  return DeltaOnceSolved(sternL);
222  }
223 
224  return -1.; // Signal the caller to use the Sternheimer approximation,
225  // because we have been unable to solve the exact form.
226 }
227 
228 /* Newton's method for finding roots. Adapted from G4PolynominalSolver, but
229  * without the assumption that the input is a polynomial. Also, here we
230  * always expect the roots to be positive, so return -1 as an error value. */
232 {
233  const G4int maxIter = 100;
234  G4int nbad = 0, ngood = 0;
235 
236  G4double lambda(start), value(0.), dvalue(0.);
237 
238  if(fVerbose > 2) {
239  G4cout << "G4DensityEffectCalculator::Newton: strat= " << start
240  << " type: " << first << G4endl;
241  }
242  while(true) {
243  if(first) {
244  value = FRho(lambda);
245  dvalue = DFRho(lambda);
246  } else {
247  value = Ell(lambda);
248  dvalue = DEll(lambda);
249  }
250  if(dvalue == 0.0) { break; }
251  const G4double del = value/dvalue;
252  lambda -= del;
253 
254  const G4double eps = std::abs(del);
255  if(eps <= 1.e-12) {
256  ++ngood;
257  if(ngood == 2) {
258  if(fVerbose > 2) {
259  G4cout << " Converged with result= " << lambda << G4endl;
260  }
261  return lambda;
262  }
263  } else {
264  ++nbad;
265  }
266  if(nbad > maxIter || eps > 1.) { break; }
267  }
268  if(fVerbose > 2) {
269  G4cout << " Failed to converge last value= " << value
270  << " dvalue= " << dvalue << " lambda= " << lambda << G4endl;
271  }
272  return -1.;
273 }
274 
275 /* Return the derivative of the equation used
276  * to solve for the Sternheimer parameter rho. */
278 {
279  G4double ans = 0.0;
280  for(G4int i = 0; i < nlev; ++i) {
281  if(sternf[i] > 0.) {
282  ans += sternf[i] * gpow->powN(levE[i], 2) * rho /
283  (gpow->powN(levE[i] * rho, 2)
284  + 2./3. * sternf[i] * gpow->powN(plasmaE, 2));
285  }
286  }
287  return ans;
288 }
289 
290 /* Return the functional value for the equation used
291  * to solve for the Sternheimer parameter rho. */
293 {
294  G4double ans = 0.0;
295  for(G4int i = 0; i<nlev; ++i) {
296  if(sternf[i] > 0.) {
297  ans += sternf[i] * G4Log(gpow->powN(levE[i]*rho, 2) +
298  2./3. * sternf[i]*gpow->powN(plasmaE, 2));
299  }
300  }
301  ans *= 0.5; // pulled out of loop for efficiency
302 
303  if(fConductivity > 0.) {
304  ans += fConductivity * G4Log(plasmaE * std::sqrt(fConductivity));
305  }
306  ans -= G4Log(meanexcite);
307  return ans;
308 }
309 
310 /* Return the derivative for the equation used to
311  * solve for the Sternheimer parameter l, called 'L' here. */
313 {
314  G4double ans = 0.;
315  for(G4int i=0; i<nlev; ++i) {
316  if(sternf[i] > 0 && (sternEbar[i] > 0. || L != 0.)) {
317  const G4double y = gpow->powN(sternEbar[i], 2);
318  ans += sternf[i]/gpow->powN(y + L*L, 2);
319  }
320  }
321  ans *= (-2*L); // pulled out of the loop for efficiency
322  return ans;
323 }
324 
325 /* Return the functional value for the equation used to
326  * solve for the Sternheimer parameter l, called 'L' here. */
328 {
329  G4double ans = 0.;
330  for(G4int i=0; i<nlev; ++i) {
331  if(sternf[i] > 0. && (sternEbar[i] > 0. || L != 0.)) {
332  ans += sternf[i]/(gpow->powN(sternEbar[i], 2) + L*L);
333  }
334  }
335  ans -= gpow->powZ(10, -2 * sternx);
336  return ans;
337 }
338 
345 {
346  G4double ans = 0.;
347  for(G4int i=0; i<nlev; ++i) {
348  if(sternf[i] > 0.) {
349  ans += sternf[i] * G4Log((gpow->powN(sternl[i], 2)
350  + gpow->powN(sternL, 2))/gpow->powN(sternl[i], 2));
351  }
352  }
353  ans -= gpow->powN(sternL, 2)/(1 + gpow->powZ(10, 2 * sternx));
354  return ans;
355 }