Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Implement exact density effect calculation. #5

Closed
wants to merge 4 commits into from

Conversation

straitm
Copy link

@straitm straitm commented Jan 18, 2019

[I don't know whether this is an official place to submit patches, so I am also
opening a Bugzilla issue.]

Added calculation of the Fermi density effect using atomic data as
described in R.M. Sternheimer et al. "Density Effect For The Ionization Loss
of Charged Particles in Various Substances" Atom. Data Nucl. Data Tabl. 30
(1984) 261-271.

This calculation gives values for the density effect 'delta' which are much
more accurate than that provided by the Sternheimer 3-part approximation,
introduced in his 1952 paper, Phys.Rev. 88 851-859. This is particularly
true when using the general method for producing parameter
values given in his 1971 paper, Phys.Rev. B3 3681-3692, but doing the
full calculation is still a significant improvement over using the
tabulated values given in the 1984 paper for elements and a selected
set of mixtures and compounds.

The purpose of the 3-part approximation was to avoid heavy computation from
the perspective of sliderules (1952) or minicomputers (1971). It
is no longer necessary. I ran tests using the NOvA simulation framework and
could measure no difference in performance between using the approximation
and the full treatment.

A major advantage of this method is that the stopping power of a substance
does not abruptly change when the user modifies it so that it is no
longer a tabulated material. In our case, we noticed this effect when
we changed iron to steel (98% iron, 2% other). Previously, iron got
the 1984 tabulated values while steel got the 1971 general form, which
gives up to a 1.3% difference in mean muon range, with the largest
difference at about 200MeV.

I have tested the code with muons in iron and verified that the mean range
changes as expected between the approximate and exact forms.

Implementation notes:

  • The calculation involves finding roots and can fail. In this case, the
    code falls back to the approximation and prints a message saying it has
    done so. I have not found any cases that prompt a failure, but have only
    tested with materials found in and around the NOvA detector.

  • Attempted to make clear in places that refer to the parameterization,
    such as G4IonisParamMat::SetDensityEffectParameters(), that these calls only
    have any effect if the full calculation fails.

  • I found a separate implementation of the Sternheimer parameterization in
    source/processes/electromagnetic/highenergy/src/G4mplIonisationModel.cc
    I don't know why this parallel implementation exists, so I simply changed
    a comment to make it more clear that it was the approximate form.

  • G4IonisParamMat::DensityCorrection() is no longer inlined.

  • Removed debugging print statements in G4IonisParamMat::ComputeDensityEffect().

Added calculation of the Fermi density effect using atomic data as
described in R.M. Sternheimer et al. "Density Effect For The Ionization Loss
of Charged Particles in Various Substances" Atom. Data Nucl. Data Tabl. 30
(1984) 261-271.

This calculation gives values for the density effect 'delta' which are much
more accurate than that provided by the Sternheimer 3-part approximation,
introduced in his 1952 paper, Phys.Rev. 88 851-859.  This is particularly
true when using the general method for producing parameter
values given in his 1971 paper, Phys.Rev. B3 3681-3692, but doing the
full calculation is still a significant improvement over using the
tabulated values given in the 1984 paper for elements and a selected
set of mixtures and compounds.

The purpose of the 3-part approximation was to avoid heavy computation from
the perspective of sliderules (1952) or minicomputers (1971).  It
is no longer necessary.  I ran tests using the NOvA simulation framework and
could measure no difference in performance between using the approximation
and the full treatment.

A major advantage of this method is that the stopping power of a substance
does not abruptly change when the user modifies it so that it is no
longer a tabulated material.  In our case, we noticed this effect when
we changed iron to steel (98% iron, 2% other).  Previously, iron got
the 1984 tabulated values while steel got the 1971 general form, which
gives up to a 1.3% difference in mean muon range, with the largest
difference at about 200MeV.

I have tested the code with muons in iron and verified that the mean range
changes as expected between the approximate and exact forms.

Implementation notes:

* The calculation involves finding roots and can fail.  In this case, the
code falls back to the approximation and prints a message saying it has
done so.  I have not found any cases that prompt a failure, but have only
tested with materials found in and around the NOvA detector.

* Attempted to make clear in places that refer to the parameterization,
such as G4IonisParamMat::SetDensityEffectParameters(), that these calls only
have any effect if the full calculation fails.

* I found a separate implementation of the Sternheimer parameterization in
source/processes/electromagnetic/highenergy/src/G4mplIonisationModel.cc
I don't know why this parallel implementation exists, so I simply changed
a comment to make it more clear that it was the approximate form.

* G4IonisParamMat::DensityCorrection() is no longer inlined.

* Removed debugging print statements in G4IonisParamMat::ComputeDensityEffect().
Matthew Strait added 2 commits January 29, 2019 17:19
* Zero out the sternf array before using it.  Otherwise we're adding
to uninitialized values.

* Don't divide by Z when counting the number of electrons.  This is already
taken into account so was leading to wrong values (sometimes subtle, sometimes
obvious -- my earlier test cases were in the subtle category).

* Give a better message when calculation fails.  One failure mode that can be
reproduced by hand is if the density of a substance is arbitrarily increased.
This raises the plasma energy until Sternheimer 1984 Eq (8) has no solution.
Pure hydrogen at several g/cc is such a case.
Wrap the warnings that the new exact calculation of the Fermi
density effect can print in checks for the value of
G4NistManager::Instance()->GetVerbose() so that it is
possible for the user to suppress them.
* Removed stdio.h, stdlib.h, math.h headers, use instead "globals.hh".

* Switched to G4double, G4int, G4bool.

* Switched from abs() and fabs() to std::abs().

* Use G4Pow::powN() and G4Pow::powZ() instead of pow().
@civanch
Copy link

civanch commented Sep 18, 2019

This merge request cannot be taken directly, because Geant4 master was modified.

The proposal was discussed at Geant4 EM group meeting https://indico.cern.ch/event/825436/ , see Matthew Strait talk. At the meeting it was agreed, that this development is integrated to Geant4.

During integration process it was realised that the optimal implementation would be to introduce G4DensityEffectCalculator where this PR is implemented. This calculator is enabled by user request. We plan for 10.6 release (December, 2019) add UI command which enables the calculator by material name. The main argument for this choice: the density effect is most important for tracking gaseous detectors where precise dEdx simulation is required. At the same time, it is safe do not enable this computation for all users just now. We need accumulate experience and get user feedback.

@civanch
Copy link

civanch commented Sep 18, 2019

So, this PR cannot be merged directly but may be closed.

@civanch civanch closed this Sep 18, 2019
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
2 participants