Problem 177 - Integration of Equation of motion in E.M. Field, relativistic case
Summary: Integration of Equation of motion in E.M. Field, relativistic case
Status: RESOLVED FIXED
Alias: None
Product: Geant4
Classification: Unclassified
Component: geometry/magneticfield (show other problems)
Version: 2.0
Hardware: PC All
: P2 critical
Assignee: John Apostolakis
URL:
Depends on:
Blocks:
 
Reported: 2000-12-08 09:31 CET by lebrun
Modified: 2001-06-20 12:26 CEST (History)
0 users

See Also:


Attachments

Note You need to log in before you can comment on or make changes to this problem.
Description lebrun 2000-12-08 09:31:24 CET
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.
Comment 1 John Apostolakis 2001-06-20 12:25:59 CEST
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.