Track Error Propagation

The error propagation package serves to propagate one particle together with its error from a given trajectory state until a user-defined target is reached (a surface, a volume, a given track length,…).

Physics

The error propagator package computes the average trajectory that a particle would follow. This means that the physics list must have the following characteristics:

  • No multiple scattering

  • No random fluctuations for energy loss

  • No creation of secondary tracks

  • No hadronic processes

It has also to be taken into account that when the propagation is done backwards (in the direction opposed to the one the original track traveled) the energy loss has to be changed into an energy gain.

All this is done in the G4ErrorPhysicsList class, that is automatically set by G4ErrorPropagatorManager as the Geant4 physics list. It sets G4ErrorEnergyLoss as unique electromagnetic process. This process uses the Geant4 class G4EnergyLossForExtrapolator to compute the average energy loss for forwards or backwards propagation. To avoid getting too different energy loss calculation when the propagation is done forwards (when the energy at the beginning of the step is used) or backwards (when the energy at the end of the step is used, always smaller than at the beginning) G4ErrorEnergyLoss computes once the energy loss and then replaces the original energy loss by subtracting/adding half of this value (what is approximately the same as computing the energy loss with the energy at the middle of the step). In this way, a better calculation of the energy loss is obtained with a minimal impact on the total CPU time.

The user may use his/her own physics list instead of G4ErrorPhysicsList. As it is not needed to define a physics list when running this package, the user may have not realized that somewhere else in his/her application it has been defined; therefore a warning will be sent to advert the user that he is using a physics list different to G4ErrorPhysicsList. If a new physics list is used, it should also initialize the G4ErrorMessenger with the classes that serve to limit the step:

G4ErrorEnergyLoss* eLossProcess = new G4ErrorEnergyLoss;
G4ErrorStepLengthLimitProcess* stepLengthLimitProcess = new G4ErrorStepLengthLimitProcess;
G4ErrorMagFieldLimitProcess* magFieldLimitProcess = new G4ErrorMagFieldLimitProcess;
new G4ErrorMessenger( stepLengthLimitProcess, magFieldLimitProcess, eLossProcess );

To ease the use of this package in the reconstruction code, the physics list, whether G4ErrorPhysicsList or the user’s one, will be automatically initialized before starting the track propagation if it has not been done by the user.

Trajectory state

The user has to provide the particle trajectory state at the initial point. To do this it has to create an object of one of the children classes of G4ErrorTrajState, providing:

  • Particle type

  • Position

  • Momentum

  • Trajectory error matrix

G4ErrorTrajState( const G4String& partType,
                  const G4Point3D& pos,
                  const G4Vector3D& mom,
                  const G4ErrorTrajErr& errmat = G4ErrorTrajErr(5,0) );

A particle trajectory is characterized by five independent variables as a function of one parameter (e.g. the path length). Among the five variables, one is related to the curvature (to the absolute value of the momentum), two are related to the direction of the particle and the other two are related to the spatial location.

There are two possible representations of these five parameters in the error propagator package: as a free trajectory state, class G4ErrorTrajStateFree, or as a trajectory state on a surface, class G4ErrorTrajStateonSurface.

Free trajectory state

In the free trajectory state representation the five trajectory parameters are

  • G4double fInvP

  • G4double fLambda

  • G4double fPhi

  • G4double fYPerp

  • G4double fZPerp

where fInvP is the inverse of the momentum. fLambda and fPhi are the dip and azimuthal angles related to the momentum components in the following way:

p\_x = p cos(lambda) cos(phi) p\_y = p cos(lambda) sin(phi) p\_z = p sin(lambda), that is, lambda = 90 - theta, where theta is the usual angle with respect to the Z axis.

fYperp and fZperp are the coordinates of the trajectory in a local orthonormal reference frame with the X axis along the particle direction, the Y axis being parallel to the X-Y plane (obtained by the vectorial product of the global Z axis and the momentum).

Trajectory state on a surface

In the trajectory state on a surface representation the five trajectory parameters are

  • G4double fInvP

  • G4double fPV

  • G4double fPW

  • G4double fV

  • G4double fW

where fInvP is the inverse of the momentum; fPV and fPW are the momentum components in an orthonormal coordinate system with axis U, V and W; fV and fW are the position components on this coordinate system.

For this representation the user has to provide the plane where the parameters are calculated. This can be done by providing two vectors, V and W, contained in the plane:

G4ErrorSurfaceTrajState( const G4String& partType,
                         const G4Point3D& pos,
                         const G4Vector3D& mom,
                         const G4Vector3D& vecV,
                         const G4Vector3D& vecW,
                         const G4ErrorTrajErr& errmat = G4ErrorTrajErr(5,0) );

or by providing a plane

G4ErrorSurfaceTrajState( const G4String& partType,
                         const G4Point3D& pos,
                         const G4Vector3D& mom,
                         const G4Plane3D& plane,
                         const G4ErrorTrajErr& errmat = G4ErrorTrajErr(5,0) );

In this second case the vector V is calculated as the vector in the plane perpendicular to the global vector X (if the plane normal is equal to X, Z is used instead) and W is calculated as the vector in the plane perpendicular to V.

Trajectory state error

The 5X5 error matrix should also be provided at the creation of the trajectory state as a G4ErrorTrajErr object. If it is not provided a default object will be created filled with null values.

Currently the G4ErrorTrajErr is a G4ErrorSymMatrix, a simplified version of CLHEP HepSymMatrix.

The error matrix is given in units of GeV and cm. Therefore you should do the conversion if your code is using other units.

Targets

The user has to define up to where the propagation must be done: the target. The target can be a surface G4ErrorSurfaceTarget, which is not part of the Geant4 geometry. It can also be the surface of a Geant4 volume G4ErrorGeomVolumeTarget, so that the particle will be stopped when it enters this volume. Or it can be that the particle is stopped when a certain track length is reached, by implementing a G4ErrorTrackLengthTarget.

Surface target

When the user chooses a G4ErrorSurfaceTarget as target, the track is propagated until the surface is reached. This surface is not part of Geant4 geometry, but usually traverses many Geant4 volumes. The class G4ErrorNavigator takes care of the double navigation: for each step the step length is calculated as the minimum of the step length in the full geometry (up to a Geant4 volume surface) and the distance to the user-defined surface. To do it, G4ErrorNavigator inherits from G4Navigator and overwrites the methods ComputeStep() and ComputeSafety(). Two types of surface are currently supported (more types could be easily implemented at user request): plane and cylindrical.

Plane surface target

G4ErrorPlaneSurfaceTarget implements an infinite plane surface. The surface can be given as the four coefficients of the plane equation ax+by+cz+d = 0:

G4ErrorPlaneSurfaceTarget(G4double a=0,
                          G4double b=0,
                          G4double c=0,
                          G4double d=0);

or as the normal to the plane and a point contained in it:

G4ErrorPlaneSurfaceTarget(const G4Normal3D &n,
                          const G4Point3D &p);

or as three points contained in it:

G4ErrorPlaneSurfaceTarget(const G4Point3D &p1,
                          const G4Point3D &p2,
                          const G4Point3D &p3);

Cylindrical surface target

G4ErrorCylSurfaceTarget implements an infinite-length cylindrical surface (a cylinder without end-caps). The surface can be given as the radius, the translation and the rotation

G4ErrorCylSurfaceTarget( const G4double& radius,
                         const G4ThreeVector& trans=G4ThreeVector(),
                         const G4RotationMatrix& rotm=G4RotationMatrix() );

or as the radius and the affine transformation

G4ErrorCylSurfaceTarget( const G4double& radius,
                         const G4AffineTransform& trans );

Geometry volume target

When the user chooses a G4ErrorGeomVolumeTarget as target, the track is propagated until the surface of a Geant4 volume is reached. User can choose if the track will be stopped only when the track enters the volume, only when the track exits the volume or in both cases.

The object has to be instantiated giving the name of a logical volume existing in the geometry:

G4ErrorGeomVolumeTarget( const G4String& name );

Track Length target

When the user chooses a G4ErrorTrackLengthTarget as target, the track is propagated until the given track length is reached.

The object has to be instantiated giving the value of the track length:

G4ErrorTrackLengthTarget(const G4double maxTrkLength );

It is implemented as a G4VDiscreteProcess and it limits the step in PostStepGetPhysicalInteractionLength. To ease its use, the process is registered to all particles in the constructor.

Managing the track propagation

The user needs to propagate just one track, so there is no need of run and events. neither of G4VPrimaryGeneratorAction. G4ErrorPropagator creates a track from the information given in the G4ErrorTrajState and manages the step propagation. The propagation is done by the standard Geant4 methods, invoking G4SteppingManager::Stepping() to propagate each step.

After one step is propagated, G4ErrorPropagator takes cares of propagating the track errors for this step, what is done by G4ErrorTrajStateFree::PropagateError(). The equations of error propagation are only implemented in the representation of G4ErrorTrajStateFree. Therefore if the user has provided instead a G4ErrorTrajStateOnSurface object, it will be transformed into a G4ErrorTrajStateFree at the beginning of tracking, and at the end it is converted back into G4ErrorTrajStateOnSurface on the target surface (on the normal plane to the surface at the final point).

The user G4VUserTrackingAction::PreUserTrackingAction( const G4Track* ) and G4VUserTrackingAction::PreUserTrackingAction( const G4Track* ) are also invoked at the beginning and at the end of the track propagation.

G4ErrorPropagator stops the tracking when one of the three conditions is true:

  • Energy is exhausted

  • World boundary is reached

  • User-defined target is reached

In case the defined target is not reached, G4ErrorPropagator::Propagate() returns a negative value.

The propagation of a trajectory state until a user defined target can be done by invoking the method of G4ErrorPropagatorManager

G4int Propagate( G4ErrorTrajState* currentTS, const G4ErrorTarget* target,
                 G4ErrorMode mode = G4ErrorMode_PropForwards );

You can get the pointer to the only instance of G4ErrorPropagatorManager with

G4ErrorPropagatorManager* g4emgr = G4ErrorPropagatorManager::GetErrorPropagatorManager();

Another possibility is to invoke the propagation step by step, returning control to the user after each step. This can be done with the method

G4int PropagateOneStep( G4ErrorTrajState* currentTS,
                        G4ErrorMode mode = G4ErrorMode_PropForwards );

In this case you should register the target first with the command

G4ErrorPropagatorData::GetG4ErrorPropagatorData()->SetTarget( theG4eTarget );

Error propagation

As in the GEANT3-based GEANE package, the error propagation is based on the equations of the European Muon Collaboration, that take into account:

  • Error from curved trajectory in magnetic field

  • Error from multiple scattering

  • Error from ionization

The formulas assume propagation along an helix. This means that it is necessary to make steps small enough to assure magnetic field constantness and not too big energy loss.

Limiting the step

There are three ways to limit the step. The first one is by using a fixed length value. This can be set by invoking the user command:

G4UImanager::GetUIpointer()->ApplyCommand("/geant4e/limits/stepLength MY_VALUE MY_UNIT");

The second one is by setting the maximum percentage of energy loss in the step (or energy gain is propagation is backwards). This can be set by invoking the user command:

G4UImanager::GetUIpointer()->ApplyCommand("/geant4e/limits/energyLoss MY_VALUE");

The last one is by setting the maximum difference between the value of the magnetic field at the beginning and at the end of the step. Indeed what is limited is the curvature, or exactly the value of the magnetic field divided by the value of the momentum transverse to the field. This can be set by invoking the user command:

G4UImanager::GetUIpointer()->ApplyCommand("/geant4e/limits/magField MY_VALUE");

The classes that limit the step are implemented as Geant4 processes. Therefore, the invocation of the above-mentioned commands should only be done after the initialization (for example after G4ErrorPropagatorManager::InitGeant4e().