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

fix: Energy loss function #1323

Merged
merged 3 commits into from
Apr 12, 2023
Merged

Conversation

beomki-yeo
Copy link
Contributor

@beomki-yeo beomki-yeo commented Jul 12, 2022

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
image

Eq. 33.11
image

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

@beomki-yeo beomki-yeo changed the title Fix Bethe Energy loss function fix: Bethe Energy loss function Jul 12, 2022
@codecov
Copy link

codecov bot commented Jul 12, 2022

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%) ⬇️
Core/include/Acts/Material/Interactions.hpp 100.00% <100.00%> (ø)

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

@paulgessinger
Copy link
Member

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 asalzburger self-requested a review July 13, 2022 08:54
@asalzburger
Copy link
Contributor

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

@beomki-yeo
Copy link
Contributor Author

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

@beomki-yeo beomki-yeo changed the title fix: Bethe Energy loss function fix: Energy loss function Jul 18, 2022
@paulgessinger
Copy link
Member

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
Copy link
Contributor

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
Copy link
Contributor

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
Copy link
Member

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.

@beomki-yeo beomki-yeo marked this pull request as draft July 19, 2022 15:29
@paulgessinger paulgessinger added the 🚧 WIP Work-in-progress label Jul 21, 2022
@andiwand
Copy link
Contributor

this seems to have some impact on #1385 (comment)

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
Copy link
Contributor Author

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

@stale
Copy link

stale bot commented Sep 20, 2022

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.

@stale stale bot added the Stale label Sep 20, 2022
@stale stale bot removed the Stale label Dec 13, 2022
@beomki-yeo beomki-yeo marked this pull request as ready for review December 13, 2022 12:05
@beomki-yeo beomki-yeo removed the 🚧 WIP Work-in-progress label Dec 13, 2022
@beomki-yeo
Copy link
Contributor Author

beomki-yeo commented Dec 13, 2022

@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
Copy link

github-actions bot commented Dec 13, 2022

📊 Physics performance monitoring for 7bc345b

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
Copy link
Contributor Author

beomki-yeo commented Dec 13, 2022

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

@asalzburger
Copy link
Contributor

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

@andiwand
Copy link
Contributor

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
Copy link
Contributor Author

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
Copy link
Contributor

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

@andiwand
Copy link
Contributor

andiwand commented Apr 6, 2023

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

lets hope not 🤔 it would be pions in that case right?

@beomki-yeo
Copy link
Contributor Author

@andiwand
Copy link
Contributor

andiwand commented Apr 6, 2023

the iPot * iPot does also not align with the equation you posted no?

      // calculate the most probable value of the Landau distribution
      double mpv = m_mpvScale * xi / (beta * beta) *
                   (std::log(m * eta2 * kaz / (iPot * iPot)) + 0.2 -
                    beta * beta - delta); //
                                          // 12.325);

@beomki-yeo
Copy link
Contributor Author

beomki-yeo commented Apr 6, 2023

TBH, I didn't look into other terms seriously; I will crosscheck - but why do you think it is not aligned?

@andiwand
Copy link
Contributor

andiwand commented Apr 6, 2023

I the equation you posted there is only / I not / I^2. maybe some terms cancel but I could not see that

@beomki-yeo
Copy link
Contributor Author

Nah the equation I posted is in the form of log(A/I) + log(B/I) = log(AB/I^2) so I think it is OK

@beomki-yeo
Copy link
Contributor Author

Just a remark on the pure impact of Landau function mass input

monitored(incident particle: muon I guess)
reference(electron)

@andiwand
Copy link
Contributor

andiwand commented Apr 7, 2023

Nah the equation I posted is in the form of log(A/I) + log(B/I) = log(AB/I^2) so I think it is OK

haha I fell for the same trick twice - thanks 😄

lets bring this PR one more time into the round table on tuesday and then merge it

@beomki-yeo
Copy link
Contributor Author

Of course

@beomki-yeo beomki-yeo marked this pull request as draft April 7, 2023 20:59
@beomki-yeo beomki-yeo force-pushed the Fix-Bethe-Function branch 2 times, most recently from 4b7f52a to 89a2962 Compare April 11, 2023 11:14
@beomki-yeo beomki-yeo marked this pull request as ready for review April 11, 2023 13:08
Copy link
Contributor

@asalzburger asalzburger left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@andiwand has ironed that out anyway - again, @beomki-yeo great work!

@asalzburger
Copy link
Contributor

One pending discussion to be resolved, then we can merge that in.

@beomki-yeo
Copy link
Contributor Author

I also need to update the CI reference and root hash tag - I will do that promptly

@beomki-yeo
Copy link
Contributor Author

It is good to go now

@kodiakhq kodiakhq bot merged commit fcdc820 into acts-project:main Apr 12, 2023
37 checks passed
@paulgessinger paulgessinger modified the milestones: next, v25.0.0 Apr 21, 2023
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Component - Core Affects the Core module Impact - Major Significant bug and/or affects a lot of modules
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

6 participants