`EigenStepper` (and extensions) cannot be easily extended for dense environment covariance transport

[andiwand]

Our implementation of the EigenStepper has an extension mechanism which is supposed to automatically switch between "normal" and dense environments if necessary. This was most likely done to have CPU optimal code and not waste resources with the dense extension if we are not in a dense environment.

While implementing the covariance transport for the dense environment I realized that we never calculate a free covariance matrix but only accumulate the jacobian which will be used on demand to transport the bound covariance from the bound covariance at the last surface. Using the free representation is important for the multiple scattering contribution as it will affect the spatial covariance and the direction covariance. The bound spatial coordinates are not necessarily cartesian and in mm so we can only rely on the coordinates and units of the free representation.

Since the stepper extension can be switched any time without notice the free covariance transport needs to be done in all extensions and with some shared state. So shared state can only live in the stepper for now. Both is not optimal.

My proposed solution would be to replace the extension list with a single extension which can relay on its state to be consistent and a mechanism to get a free covariance from the extension by the stepper. This would allow the user to still modify the RK4 integration and to efficiently and reliably handle dense environments.

[andiwand]

started implementing the proposed solution here https://github.com/acts-project/acts/pull/2865

[github-actions[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]

Just a question for my understanding. As I also need to think about it for detray implementation. Isn't it enough to calculate the multiple scattering jacobian only for angles (3X3 matrix for free parameter) and multiply it to the corresponding block of transport jacobian in the stepper?

[andiwand]

The complication I encountered in ACTS is that we do not have a covariance matrix while stepping but only the accumulated jaccobian. We can either accumulate a free param covariance during stepping and then bring it into the end result by transforming it accordingly. Or we can generate the additional covariance term (free or bound) at the end and add it to the transported covariance.

In ACTS this would have to be done by the extension but there are no callbacks wired for this yet and the multiple extension scenario makes this very complicated.

[beomki-yeo]

How does the Multiple scattering work in the dense environment? I have never looked into the code yet

🤔 Naively, I was meaning multiplying the multiple scattering jacobian to jacTransport in Stepper as long as the extension can access this variable: state.jacTransport = FreeMatrix::Identity();

// Per step
jacTransport *= D // from RKN step
block<3,3>(jacTransport,4,4) *= M // from Multiple scattering 

I think such Jacobian is provided in Eq. 4.87 - 4.102 in https://link.springer.com/book/10.1007/978-3-030-65771-0 Even though those are in 2X2 matrix for (phi,theta) but it easy almost joke to convert them into 3x3 one

[beomki-yeo]

And I have mathematically derive the ACTS pointwise multiple scattering free covaricne update from those Jacobians so they are already consistent with our code

[andiwand]

I have an implementation here for the MSC covariance https://github.com/acts-project/acts/blob/862d739e5929eda2f96fc7f38a6f8378269455e3/Core/include/Acts/Propagator/detail/GenericDenseEnvironmentExtension.hpp#L548-L557

I used the free representation here but I think bound should be also possible.

In general I don't think this should effect the transport jacobian but only the final bound covariance. The trajectory should stay the same with and without MSC only the uncertainty in position and direction should go up.

One problem is how to do this in a continuous way. I think this is properly solved in ATHENA by transporting individual pieces of material depending on the amount. I have not implemented this yet in ACTS. Also we need more wiring into the extension to properly handle the update.

[beomki-yeo]

The equations I showed do not touch the track parameter.

I will let you know once the idea is implemented and validated in the detray side - I believe it will much easier than the current implementation :D

[beomki-yeo]

Hmm I tried to find the MSC jacobian ($J_M$) that satisfies the following equation

$$ C_f = C_i + \sigma^2 \begin{pmatrix} \frac{1}{t_x^2 + tz^2} & 0 \newline 0 & 1 \end{pmatrix} = J{M} Ci J^T{M} $$

and it is found that the solution does not exist :/ It's very interesting...

[andiwand]

I fear it is not possible to replace the addition with a jacobian transform - also I am not sure if this would play nice with transporting this additional term correctly.

I remembered not why I used the free representation. The problem is that we do not have a common coordinate system for bound coordinates so the additional term for position would have to be specialized for each surface type. Using the free representation surrounds this problem as the transform to bound will handle the shift in coordinates.

[beomki-yeo]

I understand. But I still don't get how you get the free covariance matrix during the RK4 stepper. Are you planning to obtain it for every step? (Just to avoid the confusion, the free covariance matrix refers the free covariance matrix of track parameters not the additional term)

[andiwand]

Right now I have to do it at every step but that seems very wasteful in terms of CPU. I plan to accumulate the material during the stepping and when the bound covariance is requested I will calculate the addition term in free representation and sneak it into the final bound covariance. I think this is also done in Athena with the slight modification of doing the material in chunks as this seems to be more physical accurate. But I have to check that.

[beomki-yeo]

I believe the thin scatterer method is not supposed to be performed on the accumulated material - Please take a look into the same chapter of book that I referred above. There is a thick scatterer method as well with a reasonable amount of calculation.

Whatever the method is, it might be OK as long as the pull distribution is OK

[andiwand]

I don't mean to necessarily combine all the material into one thin slab but to calculate the additional term in one go and just collecting the material in a vector while stepping.

I will check the reference thanks! I already wondered if there is some formalism for thick material. Athena seems to approximate with thin layers.

[beomki-yeo]

calculate the additional term in one go and just collecting the material in a vector while stepping.

Hmm not sure if the RKN error propagation and MSC cov update are commutative :thinking:

[andiwand]

I was also worried about that and I didn't prove it but I convinced myself by the fact that the MSC is symmetric along the momentum and it does not alter the track so it should really just be an additional term to the final covariance.

[beomki-yeo]

We might need to be careful because the following equations are not the same. The first is what I think normal for two-RKN-step covariance update and second is what you are proposing:

$$ C_f = J_2 (J_1 C_i J^T1 + C{M1}) J^T2 + C{M2} $$

$$ C_f = J_2 J_1 C_i J^T_1 J^T2 + C{M1} + C_{M2} $$

where $Jn$ is the RKN transport jacobian per step and $C{Mn}$ is the MSC addition covariance per step

[andiwand]

Right that would lead to different results but if we put the transport effect of the jacobian into the total MSC covariance (combining multiple steps into one) it should give the same, I would think.

[github-actions[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.