fix: Regularize GainMatrixSmoother inverse

[andiwand]

issue: #1215

mitigate singular covariance matrix in kalman smoother by adding a regularization term $\lambda I$

[codecov]

Codecov Report

Merging #1298 (14bd7ed) into main (f95a332) will decrease coverage by 0.00%.
The diff coverage is 0.00%.

@@            Coverage Diff             @@
##             main    #1298      +/-   ##
==========================================
- Coverage   47.43%   47.42%   -0.01%     
==========================================
  Files         375      375              
  Lines       19787    19788       +1     
  Branches     9286     9287       +1     
==========================================
  Hits         9385     9385              
  Misses       4021     4021              
- Partials     6381     6382       +1     

Impacted Files	Coverage Δ
Core/src/TrackFitting/GainMatrixSmoother.cpp	10.63% <0.00%> (-0.24%)	⬇️

📣 Codecov can now indicate which changes are the most critical in Pull Requests. Learn more

[asalzburger]

Hi @andiwand - that's a good step, there's also a more robust formulation of the KFilter and Smoother, as outlined in the original paper by R. Fruehwirth ... shall we try to implement that as a version as well?

[andiwand]

@asalzburger I agree we should try other implementations too. the other algorithm I found is this one https://en.wikipedia.org/wiki/Kalman_filter#Modified_Bryson%E2%80%93Frazier_smoother which does not rely on the inverse of the covariance. I tried to lay out different ideas here #1215 (comment)

I also thought about using a toy model/example to compare different implementations at different noise levels to see how stable they are. maybe a nice thing to add to the thesis

[andiwand]

@paulgessinger the hash comparison fails now but I could not spot the error in the log anymore. success?

how should we go about validating the results?

[paulgessinger]

I think this is simple enough to just put in to work around the immediate failures. No objections to more advanced solutions of course.

It would be good to look at the performance plots produced by the physmon job, as well as ideally checking the resolution & pull plots that you ran a while ago.

[paulgessinger]

Not sure why the bridge CI job fails. Might be some software update on the GPU runner broke the nvidia environment for docker. If that's the case, I'll have to look into it.

[andiwand]

also interesting is that linux_physmon seems to be stuck? https://github.com/acts-project/acts/runs/7107138181?check_suite_focus=true

[paulgessinger]

Oh indeed.

[paulgessinger]

The physmon job again seems to hang. This is very strange.

[andiwand]

it does not even tell where it hangs. and I don't see how this could be related to my changes

[paulgessinger]

Maybe we need to let it fail and see if that tells us how it far it gets.

[paulgessinger]

Ok the GPU job should be fixed again.

[andiwand]

I looked at these plots here https://github.com/acts-project/acts/suites/7163935748/artifacts/285350562

I don't understand all of them. some look fine even the tho the statistic tests fail. some look better than the previous version. some look a little funky, especially those with pT. looking into the particle gun configuration it should use [1,10] GeV for p. somehow the reference has pT up to 100 GeV?

we can also play a little with the epsilon scale

@paulgessinger

[andiwand]

weird pT thing

much cleaner track efficiency

less biased pull phi?

more zero centered q/p?

[andiwand]

adopted epsilon to 10^-13 which seems to keep the histcmp green. putting the analysis comparison back in hoping it will not crash again

[asalzburger]

That's very interesting, indeed - the 16fold structure in the efficiency looks like we had a problem with the overlap sections in the innermost layer ... we need to check what's going on there.

[andiwand]

@asalzburger my thought on this was that these dips are the missing tracks where the smoothing failed. since the smoothing is "fixed" now they appear in the efficiency.

maybe the smearing produces too little variance if the particle hits the module with a 90° inclination and that is why the gain matrix blew up?

[asalzburger]

@asalzburger my thought on this was that these dips are the missing tracks where the smoothing failed. since the smoothing is "fixed" now they appear in the efficiency.

maybe the smearing produces too little variance if the particle hits the module with a 90° inclination and that is why the gain matrix blew up?

Yes, indeed - I was wondering, why it would fail at certain values - could be the incident angle or the overlap region of detectors, can you plot the local coordinates of the hits in the two cases (on the first layer)? That should answer the question

[asalzburger]

The pT is indeed weird. I'm not sure why we have values that high in the reference. If the particle gun config didn't change it should be the same range.

As far as I understood - these jobs now have been configured to run 1-10 GeV ... are we sure that the generation in the monitored job is restricted to generate in pT (and not p?)

[andiwand]

I was digging into this if pT or p is sampled from [1,10] GeV

physmon calls into ckf_tracks https://github.com/acts-project/acts/blob/main/CI/physmon/physmon.py#L82-L93 and truth_tracking_kalman https://github.com/acts-project/acts/blob/main/CI/physmon/physmon.py#L55-L61
no momentum param set in ckf_tracks https://github.com/acts-project/acts/blob/main/Examples/Scripts/Python/ckf_tracks.py#L175-L182 or truth_tracking_kalman https://github.com/acts-project/acts/blob/main/Examples/Scripts/Python/truth_tracking_kalman.py#L93-L100
particle gun defaults to None https://github.com/acts-project/acts/blob/main/Examples/Scripts/Python/particle_gun.py#L49
so defaults should come from ParametricParticleGenerator.hpp https://github.com/acts-project/acts/blob/main/Examples/Algorithms/Generators/ActsExamples/Generators/ParametricParticleGenerator.hpp#L50

leaving us with

    /// Low, high (exclusive) for absolute/transverse momentum.
    double pMin = 1 * Acts::UnitConstants::GeV;
    double pMax = 10 * Acts::UnitConstants::GeV;
    /// Indicate if the momentum referse to transverse momentum
    bool pTransverse = false;

[asalzburger]

So, pT should be even lower ?

…

On Fri, Jul 1, 2022 at 3:19 PM Andreas Stefl ***@***.***> wrote: I was digging into this if pT or p is sampled from [1,10] GeV physmon calls into ckf_tracks https://github.com/acts-project/acts/blob/main/CI/physmon/physmon.py#L82-L93 and truth_tracking_kalman https://github.com/acts-project/acts/blob/main/CI/physmon/physmon.py#L55-L61 no momentum param set in ckf_tracks https://github.com/acts-project/acts/blob/main/Examples/Scripts/Python/ckf_tracks.py#L175-L182 or truth_tracking_kalman https://github.com/acts-project/acts/blob/main/Examples/Scripts/Python/truth_tracking_kalman.py#L93-L100 particle gun defaults to None https://github.com/acts-project/acts/blob/main/Examples/Scripts/Python/particle_gun.py#L49 so defaults should come from ParametricParticleGenerator.hpp https://github.com/acts-project/acts/blob/main/Examples/Algorithms/Generators/ActsExamples/Generators/ParametricParticleGenerator.hpp#L50 leaving us with /// Low, high (exclusive) for absolute/transverse momentum. double pMin = 1 * Acts::UnitConstants::GeV; double pMax = 10 * Acts::UnitConstants::GeV; /// Indicate if the momentum referse to transverse momentum bool pTransverse = false; — Reply to this email directly, view it on GitHub <#1298 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AGLD7B4B3KZKST7HLATHRMTVR3V7VANCNFSM52ET4M6A> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[andiwand]

I would think so yes. we should really get a plot of the generated tho because I don't trust it 100%

[andiwand]

norm of p as a histogram for truth_smeared, truth_estimated, seeded

[andiwand]

I extracted these points at https://github.com/acts-project/acts/blob/main/Core/src/TrackFitting/GainMatrixSmoother.cpp#L65 with ts.referenceSurface().localToGlobal(gctx, {ts.filtered()[eBoundLoc0], ts.filtered()[eBoundLoc1]}, {0, 0, 0})

looks like it happens at random?

data.csv

[paulgessinger]

The only thing that still confuses me are the momenta.
Does your plot show total p or pT now?

[andiwand]

The only thing that still confuses me are the momenta. Does your plot show total p or pT now?

the histogram here is p not pT

this confirms the [1-10] GeV gun that I suspected here

[andiwand]

I agree this is all very confusing. the smoothed pT looks more reasonable with my changes but still not perfect as it should be lower than 10 GeV

[paulgessinger]

Ok, so it is in fact sampling total momentum, and the plot at least makes sense. I still don't know why we get these few outliers at very high pT but this could be something with the sampling/conversion, so out-of-scope for this PR.

Let's merge this then.