Problem 1021

Summary: Major error tracking in E field
Product: Geant4 Reporter: Tom Roberts <tjrob>
Component: processes/transportationAssignee: John Apostolakis <John.Apostolakis>
Status: RESOLVED FIXED    
Severity: major    
Priority: P5    
Version: 9.1   
Hardware: All   
OS: All   
URL: http://muonsinc.com/Geant4/G4TrackingBugInEfield.tar.gz
Attachments: tarball containing code that displays the bug

Description Tom Roberts 2008-08-19 21:51:37 CEST 
When tracking a particle in an E field such that the particle turns around,
wildly incorrect values of its momentum and energy can be obtained. This
is most noticeable when the particle stops (i.e. its trajectory remains
parallel to the E field); small amounts of transverse motion reduce the error
considerably. The amount of error depends sensitively on the relationship
between the point the particle stops and the location of the previous step,
and also on the step size.

A conspicuous example: a 1.0169 MeV anti_proton enters a uniform 1 MV/m
E field, and reflects back out with a kinetic energy of 8488.2 MeV (!!!).
Energy conservation of course implies its outcoming energy should equal
its ingoing energy. The entire problem occurs in a single step, and appears
to be an artifact of a fixed step length and a wildly varying velocity
within the step.

The URL and the attachment contain the source of a program that displays this bug.
Comment 1 Tom Roberts 2008-08-19 22:04:46 CEST 
Created attachment 28 [details]
tarball containing code that displays the bug

Original attach failed (I tried to set type to tar.gz manually; this is just application/octet-stream).
Comment 2 John Apostolakis 2009-05-13 02:06:07 CEST 
Tom, 
Thank you for the problem report, and clear example illustrating the issue.  The issue it exposes is the lack of provision in the field propagation code for going through a point at which the velocity is zero. 

The basic issues as I currently understand them are two:
  - we are integrating the equation of motion dP/ds = (1/v) * F  
  (where P=3-momentum, s=distance along trajectory, v=velocity, F=force), so when v approaches zero this is a challenge;
  - variable size steps is used in integration, adjusted so that both the sagita (called d_chord = chord distance) and integration errors are adequately small.  However to avoid endless integration there is a minimum step size, which by default is 0.01 mm.

It is the second issue (minimum step size = 0.01 mm) which dictates the size of the error currently in this case.  When the trial step size for integration of the equation of motion falls below this size, the integration is carried out in one step and the error is accepted whatever its size.

Simply by reducing the minimum step size it should be possible to get adequate accuracy.  I have confirmed this using the value of 1.0e-6 mm, in which case the value of the energy is within reason after the difficult step:
 112     10.000   1015.354      0.000      0.000   159.9426   0.002323

The minimum step size is a parameter of the ChordFinder
        G4ChordFinder( G4MagneticField* itsMagField,
                     G4double         stepMinimum = 1.0e-2 * mm, 
                     G4MagIntegratorStepper* pItsStepper = 0 );  

and can also be set in a variety of other ways.  This should serve as at least a workaround for nor for this case.

A full critical review of a few key classes is likely needed to address the general case of zero velocity, in particular to avoid potential divisions by zero.
Comment 3 John Apostolakis 2015-03-09 16:55:29 CET 
Hi Tom,
Did you try to use my suggestions?
This item has not had any comments since my recipe / suggestions.
John
Comment 4 John Apostolakis 2015-06-04 17:04:40 CEST 
Closing this item.
The current best solution is a small minimum step ( as little as 1.0 e-7 mm ), which should reduce the errors seen.
Comment 5 John Apostolakis 2015-06-04 17:05:12 CEST 
( This can and should be made the default for stepping in an electric field. )