Proposed improvement of the integration of motion in Geant4. About six months ago, a need to fix the change of the energy of a particle in presence of an electric field arised. While the equation of motion (G4EqMagElectricField.cc) was analytically correct, the variation of energy in a given step was never propagated up to the track itself. A set of changes to the core code has been proposed. This fix was rigged up to make the least amount of modifications, a mere patch, not a clean fix. It is clear that this patch has a serious limitation: the Runge Integration becomes highly inaccurate if the relativist gamma factor (E/m) get above 5 to 10, or so.. (This depends on the step size, the geometrical configuration or other factor limiting the G4 step, such that the Runge-Kutta steps are also limited.). This unacceptable loss of accuracy comes from a relatively straightforward numerical property of the integration of motion themselves: The variable that are integrated, (see array y[] and dydx[] in G4MagIntegratorDriver) are the position in physical space and the velocity: the first three elements of this y array are the X,Y,Z position of the particle and the next three elements are the derivatives of this position with respect to the lab. time, e.g., the velocity in the lab frame. dydx are the derivative of y with respect to the step length, measured along the path of the particle. This means that dydx[3] -> [5] are the change of velocity in the step, proportional to the acceleration. As we all know, the acceleration that the particle receives when the gamma factor is high is negligeable, and is not directly propotional to the step length, even if the the electric field is constant and tangent to the particle's path. However, in this case, the increment in energy is linear with the step size, we'll come back to that. Although the equations is G4EqMagElectricField are analytically correct, because of this non-linearity in the velocity increment: (i) Over the course of one integration step, the particle velocity can exceed the speed of light. A patch to prevent this has been installed, but what do we do to give a correct estimation of the energy? (ii) The error in the energy estimated by the Runge-Kutta algorithm in the step becomes very inaccurate: The energy must be compute from the change of velocity, which is very tiny, since the velocity has to change only by a negligeable amount. The quantity yError[3] -> [5] computed in G4MagErrorStepper:Stepper is not a good indicator of the error in the energy change of the particle. The real fix is clear: we should not compute the change in velocity along one step, but the change in the vector beta.gamma, which is directly proportional to the change in momentum, not velocity. Unfortunately, this implies changing the meaning of the y and dydx array in the stepper. Each Stepper should have in fact it's own definition for the left-hand side portion of the equation of motion should be: Imagine a Stern-Gerlach experiment, where the particle velocity does not changes, but it's spin orientation does. The error in one step should be based on the spin, evidently. The codes allows for this, and it will work well as long as the relative change of orientation of the spin is linear with the step size. Thus, this "linear relative change" dydx is a property of the combination of the equation of motion and the stepper rather than the G4FieldTrack itself, since it depends upon which equation of motion we are trying to integrate. I therefore propose to make the y and dydx array true members of the each stepper, and have different physical meaning for each stepper. Currently the member functions DumpToArray and LoadFromArray belong to the G4FieldTrack class. This is wrong, since their meaning is than unique across the entire system. I propose to make them virtual member functions of the inteface class G4MagIntegratorStepper, and upgrade all the implemented Stepper (3 of them) so that they are defined according to what we are trying to integrate. Also, add a specific Stepper for the E.M. fields, and make corresponding change to the equation of motion in G4EqMagElectricField. The change in energy in G4FieldTrack will than be guaranteed to remain linear with step length, no matter what the relativistic boost is... This is a somewhat of a big change. Hopefully, in only a few directories. I would appreciate feed-back. I have not tried them, I simply write this note as a preliminary brief to the Geant4 Supreme Court. I do apologize for sending it, whitout due process, to both the "state" Fermilab Court and to the "federal" one, at CERN. I will though try to implement these changes in the next few weeks, unless these changes have already been made at Cern, or elsewhere... Paul.
Thank you for your report. Release 3.1 of Geant4 now expresses the equations of motion in terms of the momentum (in place of the velocity), in order to correct the issue of numerical stability in the relativistic limit.