ECCE @ EIC Software
Reference for
ECCE @ EIC
simulation and reconstruction software on GitHub
Home page
Related Pages
Modules
Namespaces
Classes
Files
External Links
File List
File Members
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:851859, 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
47
#include "
G4DensityEffectCalculator.hh
"
48
#include "
G4AtomicShells.hh
"
49
#include "
G4NistManager.hh
"
50
#include "
G4Pow.hh
"
51
52
static
G4Pow
*
gpow
=
G4Pow::GetInstance
();
53
54
const
G4int
maxWarnings
= 20;
55
56
G4DensityEffectCalculator::G4DensityEffectCalculator
(
const
G4Material
*
mat
,
G4int
n
)
57
: fMaterial(mat), fVerbose(0), fWarnings(0), nlev(n)
58
{
59
fVerbose
=
std::max
(
fVerbose
,
G4NistManager::Instance
()->GetVerbose());
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.;
77
const
G4double
tot =
fMaterial
->
GetTotNbOfAtomsPerVolume
();
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.
92
const
G4double
xx
= frac*
G4AtomicShells::GetNumberOfElectrons
(Z, i);
93
if
(i < nshell-1 || !conductor) {
94
sternf
[sh] +=
xx
;
95
}
else
{
96
fConductivity
+=
xx
;
97
}
98
levE
[sh] =
G4AtomicShells::GetBindingEnergy
(Z, i)/
CLHEP::eV
;
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
}
109
plasmaE
=
fMaterial
->
GetIonisation
()->
GetPlasmaEnergy
()/
CLHEP::eV
;
110
meanexcite
=
fMaterial
->
GetIonisation
()->
GetMeanExcitationEnergy
()/
CLHEP::eV
;
111
}
112
113
G4DensityEffectCalculator::~G4DensityEffectCalculator
()
114
{
115
delete
[]
sternf
;
116
delete
[]
levE
;
117
delete
[]
sternl
;
118
delete
[]
sternEbar
;
119
}
120
121
G4double
G4DensityEffectCalculator::ComputeDensityCorrection
(
G4double
x
)
122
{
123
if
(
fVerbose
> 1) {
124
G4cout
<<
"G4DensityEffectCalculator::ComputeDensityCorrection for "
125
<<
fMaterial
->
GetName
() <<
", x= "
<< x <<
G4endl
;
126
}
127
const
G4double
approx
=
fMaterial
->
GetIonisation
()->
GetDensityCorrection
(x);
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
) {
138
G4ExceptionDescription
ed;
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
) {
157
G4ExceptionDescription
ed;
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
170
G4double
G4DensityEffectCalculator::FermiDeltaCalculation
(
G4double
x
)
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
) {
190
G4ExceptionDescription
ed;
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. */
231
G4double
G4DensityEffectCalculator::Newton
(
G4double
start
,
G4bool
first)
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. */
277
G4double
G4DensityEffectCalculator::DFRho
(
G4double
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. */
292
G4double
G4DensityEffectCalculator::FRho
(
G4double
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. */
312
G4double
G4DensityEffectCalculator::DEll
(
G4double
L
)
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. */
327
G4double
G4DensityEffectCalculator::Ell
(
G4double
L
)
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
344
G4double
G4DensityEffectCalculator::DeltaOnceSolved
(
G4double
sternL)
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
}
geant4
tree
geant4-10.6-release
source
materials
src
G4DensityEffectCalculator.cc
Built by
Jin Huang
. updated:
Wed Jun 29 2022 17:25:23
using
1.8.2 with
ECCE GitHub integration