In the old implementation it was tested whether the particle is inside at R_i + \Delta t \cdot \beta_i \cdot c whereas in the new implementation we test at R_i. This could explain the difference in the results.

I have to correct the above statement. R_i + \Delta t \cdot \beta_i \cdot c is used only to determine whether a particle is inside the element:

const double z = R(2) + P(2) * recpgamma;

if ((z > 0.0) && (z <= getElementLength())) {

I don't know what the motivation was to implement it this way. It's not comprehensible to me.

If I again use z = R(2) + P(2) * recpgamma then the results match again. This suggests that the difference is only caused by the difference of the above check.

Ok I am lost (sorry way too many things going on): what are we doing now? Will we there be an update in 1.6.x and 2.0.0 or everything finde now?

I just explained where the differences in the results stem from and that it was justified to update the reference.

We could discuss whether it makes sense to use z = R(2) + P(2) + recpgamma to test z > 0 && z <= getElementLength() or, as we do now, z = R(2). I don't see any reason to use the former.

Ok I understand and agree with the assessment.