fix: Energy loss function

[beomki-yeo]

This PR fixes the Bethe energy loss function in Interactions.cpp

Acts::computeEnergyLossBethe and Acts::computeEnergyLossLandau is implemented from equation 33.5 and 33.11 of Review in Particle Physics

And their first term is computed with computeEpsilon function, following the notation of equation 33.11:

/// Compute epsilon energy pre-factor for RPP2018 eq. 33.11.
///
/// Defined as
///
///     (K/2) * (Z/A)*rho * x * (q²/beta²)
///
/// where (Z/A)*rho is the electron density in the material and x is the
/// traversed length (thickness) of the material.
inline float computeEpsilon(float molarElectronDensity, float thickness,
                            const RelativisticQuantities& rq) {
  return 0.5f * K * molarElectronDensity * thickness * rq.q2OverBeta2;
}

Problem is that the epsilon term of equation 33.5 for Bethe function doesnt have the 0.5f term

Eq. 33.5

Eq. 33.11

If we are going to use the same epsilon term, we need to multiply 2 in the next term of Bethe energy loss function

After the change, the output of Acts::computeEnergyLossBethe gets consistent with Figure 33.2 of the reference.

Update July 18th

Looks like there is a literature error in Landau energy loss function (eq. 33.11). It used the mass term of the incident particle (m), but I guess it should be the mass of electron (m_e) as equation 33.5. I need a cross-check from someone who can access Straggling in thin silicon detectors, which I cannot (Update: It is confirmed that it should be the rest mass of electron)

Finally, to get a consistent result of energy loss function, we need to fix the computeDeltaHalf term as well. But this is already mentioned in the comment which says: "Should we use RPP2018 eq. 33.7 instead w/ tabulated constants?"

/// Compute the density correction factor delta/2.
///
/// Uses RPP2018 eq. 33.6 which is only valid for high energies.
///
/// @todo Should we use RPP2018 eq. 33.7 instead w/ tabulated constants?
inline float computeDeltaHalf(float meanExitationPotential,
                              float molarElectronDensity,
                              const RelativisticQuantities& rq) {

Yeah it does make big difference between 33.6 (valid for very high energy particles) and 33.7. I am not going to include delatHalf term correction in this PR, but it is included in acts-project/detray#282

After fixing codes (including the deltaHalf term), the Bethe energy loss function matches to the value in (PDG)

In the following figure, Bethe energy loss (dE/dx) in Silicon is compared for Acts main, this PR, acts-project/detray#282, and PDG value. The value from this PR is not accurate yet due to the incorrect deltaHalf term

For the most probable energy loss from Landau equation under the same condition of Figure 33.7 of Review in Particle Physics , detray#282 got 0.526 MeV, which is very close to the value in the figure. This PR got 0.739 MeV due to the incorrect delatHalf term

[codecov[bot]]

Codecov Report

Merging #1323 (7bc345b) into main (1d6aeb8) will increase coverage by 0.00%. The diff coverage is 87.50%.

@@           Coverage Diff           @@
##             main    #1323   +/-   ##
=======================================
  Coverage   49.81%   49.81%           
=======================================
  Files         420      421    +1     
  Lines       23874    23876    +2     
  Branches    10835    10836    +1     
=======================================
+ Hits        11894    11895    +1     
  Misses       4366     4366           
- Partials     7614     7615    +1     

Impacted Files	Coverage Δ
Core/src/Material/Interactions.cpp	74.19% <83.33%> (-1.33%)	:arrow_down:
Core/include/Acts/Material/Interactions.hpp	100.00% <100.00%> (ø)

:mega: We’re building smart automated test selection to slash your CI/CD build times. Learn more

[paulgessinger]

Looking at eq. 33.5 and 33.11 this seems indeed like an inconsistency to me. @asalzburger can you have a look at this?

Of course this changes the output of everything we run that runs the fast sim.

[asalzburger]

Hi @beomki-yeo - great catch, I assigned that to me & will cross-check the details.

[beomki-yeo]

OK. BTW I still see the inconsistency in energy loss function (from both Bethe and Landau) compared to the value in the figures of the same reference. I thought the discrepancies are from the approximated mean excitation energy, but it seems not the main reason apparently (the approximation is not that bad). I will provide the details later when I have time but do you guys have any idea what the reason could be

[paulgessinger]

Re your update: this is indeed very interesting. If we're now sort of mitigating literature errors, I guess we should thoroughly document the actual equations being used. Do you think it makes sense to have this as a page on the documentation, which references the source publications, but also lists explicitly the combined equations being used now? (in a future PR)

[AJPfleger]

m

I am not sure if we need a dedicated place in the documentation, since mostly everything is explained in this PR. I am afraid, that we might just mirror GH to the docs and lose some information (e.g. the discussions). I think, adding @paulgessinger 's suggestions as a comment in this PR could be sufficient.

[andiwand]

If we want to track the documentation in this PR we should at least link it in the code since git history is a bit loose after a following change IMO

but it would not hurt to have an MD documenting the formulas used in the code and referencing where they come from. the code can then point to this doc

[paulgessinger]

I agree with @andiwand. What I was suggesting was to have a document which explicitly writes out the equations as they are implemented after this PR, rather than having something like

// like eq 33.5 from Ref 1 but with epsilon from eq 35.22 and adjusted m to account for ...

because I think that will make it much harder to get up to speed on what's happening.

[andiwand]

this seems to have some impact on https://github.com/acts-project/acts/issues/1385#issuecomment-1217968896

do you think we can get this in in the near future? @beomki-yeo @paulgessinger

it looks like some reference tests are failing because the output changed. but otherwise the fix is complete?

[beomki-yeo]

The fix is not complete yet. I still need to modify the density correction effect.

[stale[bot]]

This issue/PR has been automatically marked as stale because it has not had recent activity. The stale label will be removed if any interaction occurs.

[beomki-yeo]

@asalzburger

Here are the key changes:

• Multiply x2 to computeEnergyLossBethe as explained above
• Add densityEffectData member variable to Material: Unless it is specified in the constructor (e.g. fromMassDensity or fromMolarDensity, it will be the default value which doesn't make any changes. I actually put the real data in makeSilicon() and makeBeryllium() defined in PredefinedMaterials.hpp

After the updates, Integration test of PropagationDenseConstant (which uses makeBeryllium() function!) failed in covariance check due to the low tolerance value (0.05) so I had to increase it to 0.07. I am not also sure if the updates in the PR are meant to increase the covariance during the propagation :p

Let me also ping @niermann999 as she might be interested in synchronizing detray and Acts material.

[github-actions[bot]]

:bar_chart: Physics performance monitoring for 7bc345b2cf21803ba4c03a0678aca1d19bb37e8d

Full report Seeding: seeded, truth estimated, orthogonal CKF: seeded, truth smeared, truth estimated, orthogonal IVF: seeded, truth smeared, truth estimated, orthogonal Ambiguity resolution: seeded, orthogonal Truth tracking Truth tracking (GSF)

Vertexing

Vertexing vs. mu

IVF seeded

IVF truth_smeared

IVF truth_estimated

IVF orthogonal

Seeding

Seeding seeded

Seeding truth_estimated

Seeding orthogonal

CKF

CKF seeded

CKF truth_smeared

CKF truth_estimated

CKF orthogonal

Ambiguity resolution

seeded

Truth tracking (Kalman Filter)

Truth tracking

Truth tracking (GSF)

Truth tracking

[beomki-yeo]

Hmm truth tracking test fails :/... @paulgessinger Is there any possibility that physmon uses data simulated from the reference commit?

[asalzburger]

I saw - let's discuss this tomorrow, I have an idea how to validate this.

[andiwand]

Hmm truth tracking test fails :/... @paulgessinger Is there any possibility that physmon uses data simulated from the reference commit?

simulation is done on the fly but we compare to reference. so I guess a failure is expected if we fix something in the energy loss?

[beomki-yeo]

simulation is done on the fly but we compare to reference.

"compare to reference" means "compare to the data simulated from reference commit"? If that is the case, the shift in qoverp pull value distribution can be explained. If not, I made a stupid mistake in the PR

[andiwand]

yes - the reference is just an output generated by a previous run of physmon with a commit from the main branch. so if you change the physics this will reflected in your physmon run and will fail the comparison

if we are happy with the new results you can copy the new output files from the CI and put the into the reference folder before merging. that way we make sure nobody accidentally breaks things again

[stale[bot]]

This issue/PR has been automatically marked as stale because it has not had recent activity. The stale label will be removed if any interaction occurs.

[paulgessinger]

Should we discuss what to do with this PR at some point?

[beomki-yeo]

Yes for sure, would you like to meet in person later today? - the qop pull distribution is distored with the change and I have no clue on it right now

[beomki-yeo]

Hopefully CI will be OK after this run :pray:

[paulgessinger]

Discussed in the meeting that this might have overlap with what was seen in: https://github.com/acts-project/acts/pull/1631

[beomki-yeo]

Actually I was surprised by #1631. If I understood correctly, the simulation is not using the interactions but reconstruction is using them in physmon?

[benjaminhuth]

Actually I was surprised by #1631. If I understood correctly, the simulation is not using the interactions but reconstruction is using them in physmon?

I think exactely thats the case at the moment. I was surprised as well when I saw it.

[beomki-yeo]

There have been a couple of updates...

1. The material interactions are turned on in simulation.py. This makes the pull value distributions better indeed. For example, the current pull value distribution of qop is quite narrow where the sigma is much smaller than 1.

2. After turning on the interactions, the propagation failed sometimes because the momentum of a particle drops to zero. Technically qop is not calculated properly in update function EigenStepper.ipp so I set the qop to the infinity in case of the zero momentum.

Meanwhile I also see many of physmon plots failing the test with the changes You can see the physmon plots here (I don't know how to retrieve the plots generated from the CI directly) : https://drive.google.com/drive/folders/11MNgL9gHVOcJtBkfvwrm-eALcE0VE8D6?usp=share_link

[stale[bot]]

This issue/PR has been automatically marked as stale because it has not had recent activity. The stale label will be removed if any interaction occurs.

[beomki-yeo]

@andiwand @paulgessinger @asalzburger @niermann999 This PR is ready to be reviewed.

Let me summarize the changes:

1. Fix the Bethe energy loss function by multiplying 2
2. Fix the Landau energy loss by replacing the incident particle mass with the electron mass
3. Add the density effect data into the material
4. Add unit tests which compare the PDG and ACTS energy loss

~~There is still a suspicious result in physmon - This PR is not meant to improve the pull distribution because simulation and reconstruction shares the same algorithm. But physmon shows that the pullmeans gets stabilized compared to the reference, which I don't clearly understand:~~

https://herald.harbinger.gessinger.dev/view/acts-project/acts/634576160/truth_tracking.html

Update: I guess having the right most probable energy from the Landau function (only used for the reconstruction) could contribute to the stabilization of pullmeans

[andiwand]

good to go from my side. any last comments before I merge @asalzburger @paulgessinger

[beomki-yeo]

Please don't tell me Athena has been using incident particle mass in Landau energy loss: https://gitlab.cern.ch/atlas/athena/-/blob/master/Tracking/TrkExtrapolation/TrkExTools/src/EnergyLossUpdator.cxx#L348

Ha ha