Step length is wrongly computed in field propagation

[hahnjo]

As described in #145:

PR https://github.com/apt-sim/AdePT/pull/95 introduced and used momentumXYMag as follows: https://github.com/apt-sim/AdePT/blob/fba5bd2373a67ffcf2478aec1d93a0d9e07171ba/magneticfield/inc/fieldPropagatorConstBz.h#L78-L91 If there is no z-direction, ie direction[2] = 0, then momentumXYMag = momentumMag. Otherwise, and especially as direction[2] approaches 1 or -1, momentumXYMag will get smaller and smaller. This leads to an increasing value for curv and a decreasing value of safeLength. So in essence, the more a track is flying in z-direction and the less it is influenced by the constant field in z-direction, the shorter safeLength will get and the more chord steps magnetic field propagation will do.

Contrary to what I wrote back then, this is not only a problem for performance, but also correctness since we kill particles that cannot be propagated in max_iterations chord steps.

[agheata]

I'll take care of adding a patch for this

[agheata]

The correct approach for estimating the safeLength value is shown below:

Using this approach with example13 using cms2018 and 3.8T field, I see a reduction of the number of killed tracks and a small speedup:

master:
24: Run time: 14.8801
24: 
24: 
24: Mean energy deposit          9.97193 GeV
24: Mean number of gamma         3638.32
24: Mean number of e-            11932.6
24: Mean number of e+            316.15
24: Mean number of charged steps 25869.4
24: Mean number of neutral steps 33551.8
24: Mean number of hits          24431.2
24: Killed in propagation:       926

fixed version:

24: Run time: 14.2466
24: 
24: 
24: Mean energy deposit          9.97475 GeV
24: Mean number of gamma         3639.57
24: Mean number of e-            11952.5
24: Mean number of e+            316.49
24: Mean number of charged steps 26036.8
24: Mean number of neutral steps 33579.2
24: Mean number of hits          24530.2
24: Killed in propagation:       787

[hahnjo]

Hm, my understanding of the error the propagation is bounding has been a different one, shouldn't we look at the "miss distance" as described in https://geant4-userdoc.web.cern.ch/UsersGuides/ForApplicationDeveloper/html/Detector/electroMagneticField.html? That one would be in the inside of the chord, not in the straight line propagation of the original direction. In particular, I don't think point P is ever used in the propagation, we're doing a linear propagation (and boundary check) from A to B.

[agheata]

Yes we propagate on chords and I agree that the "miss" distance calculated as the maximum difference between the trajectory and the chord better estimates the error. Getting the dependency of the "miss error" on the step on the helix may be more difficult however than the estimate obtained as above, which for "small displacements" over-estimates the "miss error" by about a factor of 2. So this approach may give an underestimate of the safe_length, but easier to compute and well-behaved in (hopefully) all cases. This does not exclude of course a more accurate possible approach.

[jonapost]

Indeed in Geant4 we consider the sagitta - the difference between the radius the height of the triangle define by the endpoints and the centre of the circle. The distance is estimated numerically using an estimate of the position of the midpoint of the curved trajectory. But I prefer (for a pure magnetic field that is not too anisotropic) to use the 'helix' construction. I am not sure that you are taking into consideration that the values are being projected onto to plane. This is the reason that the (1-dirz) * (1+dirz) factor is used.

The crude operation of abandoning the integration of tracks after a fixed maximum number of integration sub-steps (which I acknowledge introducing) is what is causing the bias, not calculating the 'miss distance' correctly. I think that this crude operation should be kept only as a short-term measure, and a better alternative sought as soon as possible. Once we are able to keep propagating the particles by passing the integration work (of unfinished elements) onto subsequent 'physical' steps we can dispense with this altogether - there can be a limit to the total length a particle can fly instead of the arbitrary number of integration steps.

By the way, in Geant4 after a fixed number of integration sub-steps that don't exhaust the "distance to physical interaction", the value of the 'miss-distance' is increased by a given factor - this may be anathema with respect to the need to carry a small amount of extra state per track ... Is there not a way to carry such information only for charged tracks which are still in flight ? ( Or some similar filter for other attributes ? )

[agheata]

I am not sure that you are taking into consideration that the values are being projected onto to plane. This is the reason that the (1-dirz) * (1+dirz) factor is used.

This is what the fix #172 is about.

[hahnjo]

The crude operation of abandoning the integration of tracks after a fixed maximum number of integration sub-steps (which I acknowledge introducing) is what is causing the bias, not calculating the 'miss distance' correctly.

This issue is not about the bias that I've been reporting for months, this was #145. And there the reason wasn't anything about abandoning tracks, but that particles got hard stuck on boundaries which I fixed in #170