Hello, It's been 6 months since I submitted this fix. The NOvA experiment depends on this code, so it would be helpful if it worked in mainline Geant4. Please let me know if there are any problems accepting this patch. Thanks, Matthew Created attachment 713 [details]
Plot showing delta as a function of x=log(p/m) with and without this fix
This plot shows an example of the effect of this patch. If SetDensityEffectCalculatorFlag() is not used, then with or without this patch, the approximation from
R. M. Sternheimer and R. F. Peierls. General expression for the density effect for the ionization loss of charged particles. Phys. Rev., B3:36813692, 1971
is used.
In the current version of Geant4, if SetDensityEffectCalculatorFlag() is used, then the calculation fails for the reasons described in the bug report. At high x, the Sternheimer 1971 approximation is used as a fall-back. At low x, the unphysical value -1 is returned.
With the patch in this bug report, the calculation from
R. M. Sternheimer, M. J. Berger, and S. M. Seltzer. Density effect for the ionization loss of charged particles in various substances. Atom. Data Nucl. Data Tabl., 30:261, 1984
is successfully performed, resulting in a smoother evolution of delta around x = 0 to x = 1.
Created attachment 714 [details]
Plot showing muon range with and without this fix between 500 MeV and 1.5 GeV
This plot shows, using the same color scheme as the plot of delta, the effect of using applying the patch in this bug report.
Without SetDensityEffectCalculatorFlag, the Sternheimer 1971 approximation is used, and it does not matter whether this patch is applied or not.
With SetDensityEffectCalculatorFlag, the current version of Geant4 returns an unphysical negative value of delta for low x, which results in muon range being slightly reduced.
With this patch applied and SetDensityEffectCalculatorFlag, the calculation is successful, resulting in larger values of delta around x = 0 to x = 1, and ranges that are somewhat over 1% longer.
Hello Matthew, sorry for the late answer. I was extremely busy. In 2020 it was not possible to take you code directly because there were crashes and arithmetic exceptions in various tests. Solution, which I proposed you have approved. I was not sure if that solution is valid and agree to introduce fixes you propose, but cannot take your propose blindly. It is needed to have a test for various materials. VI Hi. Sounds good. Could you propose a list of test materials? Hello Matthew, I would think that the test should be able to run for any material. By default, it may be run for G4_H, G4_C, G4_Al, G4_Fe, G4_Pb, G4_Galactic, G4_WATER, G4_PbWO4. Vladimir Created attachment 715 [details]
Delta values, approximate and calculated, for 8 test materials
For the suggested materials (G4_H, G4_C, G4_Al, G4_Fe, G4_Pb, G4_Galactic, G4_WATER, G4_PbWO4), here is a plot of delta with and without SetDensityEffectCalculatorFlag().
The dashed black line is without SetDensityEffectCalculatorFlag(), so it uses the Sternheimer parameterization, either as tabulated in G4DensityEffectData, if available, or else from the "general expression" of Sternheimer and Peierls (PRB 3, 11 (1971)).
The solid green line gives the result with SetDensityEffectCalculatorFlag(), so it uses the full calculation from Sternheimer 1984.
Most calculated values are within 0.1 or 0.2 of the approximate ones. However, note the case of PbWO4, which differs by up to 0.474. This difference is reasonable, however. Sternheimer and Peierls note in their 1971 "general expression" paper that the largest deviation in their tests, 0.381, was for the H2-Ne mixture. In both cases, there is a combination of elements with very different I-values. Since the "general expression" only uses the effective combined I-value for the mixture, it does not do a good job of evaluating this sort of substance. The full calculation instead uses the full list of atomic shells.
Created attachment 716 [details]
Difference of delta values, approximate and calculated, for 8 test materials
For the same 8 substances, the difference of delta values for the two treatments. The discontinuities you see here are expected, because the Sternheimer parameterization is defined piecewise.
Another comment (which applies to these plots and the previous ones). As it stands, the user is responsible for specifying whether a material is a conductor. By default, no materials are conductors. I have set G4_Al, G4_Pb and G4_Fe as conductors for these tests.
It is straightforward to run tests for a larger set of materials.
Created attachment 717 [details]
Incremental patch fixing treatment of a few conductors
1) Testing revealed a few cases where conductors would incorrectly trigger the check for there being no solutions to find. This was most evident for lithium and beryllium, but had small effects for some other materials as well. Now only skip trying to find solutions in the case that the material is not a conductor, because conductors *must* have solutions.
2) Adjust the fallback procedure so that it will trigger if the approximate result for delta is >= 0 instead of > 0. This avoids cases where delta = -1 gets returned.
This patch should be applied after the first patch above. If desired, I can (of course) produce a single patch which does both.
Created attachment 718 [details]
Delta values, approximate and calculated, for 308 test materials
After the second patch above, this is a test for all of the materials in G4NistMaterialBuilder. For purposes of this test, I set the following as conductors:
G4_Ac, G4_Ag, G4_Al, G4_Am, G4_Au, G4_Ba, G4_Be, G4_Bi, G4_Bk, G4_BRASS, G4_BRONZE, G4_Ca, G4_Cd, G4_Ce, G4_Cf, G4_Cm, G4_Co, G4_Cr, G4_Cs, G4_Cu, G4_Dy, G4_Er, G4_Eu, G4_Fe, G4_Fr, G4_Ga, G4_Gd, G4_Hf, G4_Hg, G4_Ho, G4_In, G4_Ir, G4_K, G4_La, G4_Li, G4_Lu, G4_Mg, G4_Mn, G4_Mo, G4_Na, G4_Nb, G4_Nd, G4_Ni, G4_Np, G4_Os, G4_Pa, G4_Pb, G4_Pd, G4_Pm, G4_Pr, G4_Pt, G4_Pu, G4_Ra, G4_Rb, G4_Re, G4_Rh, G4_Ru, G4_Sc, G4_Sm, G4_Sn, G4_Sr, G4_STAINLESS-STEEL, G4_Ta, G4_Tb, G4_Tc, G4_Te, G4_Th, G4_Ti, G4_Tl, G4_Tm, G4_U, G4_V, G4_W, G4_Y, G4_Yb, G4_Zn, G4_Zr
Note discontinuities in the *approximate* delta values for G4_Mg and G4_Tm. I believe this is caused by typos in G4DensityEffectData.cc. I will address that in a separate bug report.
The biggest difference found between the approximate and exact values of delta is for G4_Fr, presumably related to the fact that it's impossible to actually have a macroscopic amount of francium. The next attachment will have plots of the differences.
Created attachment 719 [details]
Difference of delta values, approximate and calculated, for 308 test materials
As above, but for the difference between the approximate treatment and the full calculation.
Created attachment 724 [details]
Difference of delta values, approximate and calculated, for tabulated materials only
This is the same as the above set of plots showing the difference in calculated and parameterized delta values, but this set includes only the materials that have tabulated parameters, and not those that use the "general expression" of Sternheimer and Peierls 1971.
In this case, the differences are rather small as expected.
Created attachment 725 [details]
Difference of delta values, approximate and calculated, for untabulated materials only
These are the remainder of the materials, i.e. the ones that use the "general expression". As expected, the differences between the full calculation and the parameterized values are larger for this set of materials.
Created attachment 726 [details]
Comparison of Delta_max values
Delta_max is the greatest difference reported in Sternheimer, Berger and Seltzer between their calculated values of delta and the approximate values obtained from their recommended parameter values for the Sternheimer approximation.
For all tabulated materials, this attachment gives:
Column 1: Material name
Column 2: Delta_max as reported in Sternheimer, Berger and Seltzer 1984
Column 3: Greatest difference between deltas with and without SetDensityEffectCalculator flag for a scan over x in increments of 0.01. If my calculation and SBS 1984's calculation are the same, this would be identical to column 2.
Column 4: Difference of columns 2 and 3.
For 93% of materials, the difference between SBS's Delta_max and mine are very small, less than 0.05. We would not expect the difference to be exactly zero, because SBS 1984 does not specify in what granularity they scanned over x. There's also some ambiguity to how SBS 1984 handled conductors, as noted in my code.
Some of the larger differences can be readily explained. G4_C, graphite, and porous graphite are different because the I-value has been updated since 1984. SBS used 78eV, but now we use 81eV. I suspect the other differences are similarly related to updated I-values and/or shell binding energies since 1984.
Created attachment 730 [details]
New variant of the class
Hello Matthew,
can you, please, verify the class, which is now in the Geant4 master branch. Is it what you expect?
Vladimir
Yes, look good. Thanks. Hello Matthew, the bug is fixed and the fix will be available with the next public release 11.0 (December, 2021). It will be backported to the next patch to 10.7. VI |
Created attachment 672 [details] Patch correcting the calculation In #2121 I submitted a patch that added an implementation of the exact Sternheimer calculation of the Fermi density effect. As adapted into the released version of Geant4, however, it always returns -1 for delta. This is because of several changes made after I submitted the code: * The criteria for convergence in Newton's method was changed, and the new criteria often it to report no convergence. The attached patch reverts to the criteria in my original code. * My original patch tried several starting points for Newton's method. This was replaced with a single starting point, which usually fails. The attached patch reverts to my original method. * A check that determines whether the density effect is zero because all values of the variable L are negative was removed. This caused the routine to fail when the correct result is to return delta = 0. I have reinstated this check. Additionally, the calculation was wrong in some cases even with these problems fixed because: * In my submitted version, conduction electrons were handled by adding one more element to the arrays representing energy levels. This was changed to having a separate variable hold the number of conduction electrons, but in several places in the code, they were then ignored since they were no longer part of the array. I have added the necessary code to handle them. Finally: * I have corrected a comment I wrote before which said "l^2", but should say "l". * I have changed new code that used "sum" to mean "1/sum" and "sternrho" to mean "sternrho/plasmaE" to instead always mean what they are called. I have tested the behavior after applying my new patch, and the results are as expected. Thanks, Matthew