Transportation in Magnetic Field - Further Details

The challenge of integrating all tracks

What leads us to discard tracks looping in a magnetic field

The integration of charged particle tracks in magnetic field is an important part of the computational cost (CPU time). Part of this cost is due to integration of low-energy particles in a volume with low density and strong magnetic field.

In HEP applications the most important type of tracks causing such problems are electrons in the vacuum of beam pipes. Charged particles in volumes near decay volumes and muons in large volumes of air are other examples.

To limit this CPU cost, a type of tracking cut for charged particles was introduced in Geant4 release 7.0 in G4Transportation and G4CoupledTransportation. Tracks which require more than a threshold number of integration steps [maxLoopCount] (currently 1,000) during a physics/tracking step are marked as ‘looping’ and are considered candidates for being killed - i.e. they can potentially be abandoned after the current step, and have their energy deposited locally.

Enhancements introduced in release 10.6 provide more comprehensive information about the tracks killed, in the form of G4Exception warning messages.

This section describes this policy, the parameters which the user is able to set to tune it, and recent refinements implemented in Geant4 10.5.

Cost of integration

Occasionally tracks ‘looping’ in a strong magnetic field, making little progress even over thousands of integration steps. This is due to a combination of a strong magnetic field and a thin material (gas or vacuum) in which the size of a physics step is substantially larger than the radius of curvature of the track.

The preferred integration method for tracks in an EM field is the Runge-Kutta method. This and other similar methods are well suited to variations in magnetic fields and step sizes up to a few times the radius of curvature of the charged tracks.

However when the step sizes are hundreds or thousands of times larger than the curvature of the track, these methods are expensive as they do not progress the integration of a track adequately.

The amount of CPU time which can be consumed by one or few such tracks can very large, sometimes contributing per cent increases to the simulation of some primary particles. Some tracks with a very small drift velocity (projection of the velocity along the vector of the magnetic field) can stop the progress of a simulation if they are not limited or integrated using alternative means.

So it is important to limit the number of integration steps spent on these tracks. The module for propagation in field in Geant4 flags tracks which take more than a certain number (default 1,000) integration steps without reaching the requested end of the step size, which was determined by the physics and geometry.

Parameters for eliminating or controling which particles are killed

The Geant4 G4Transportation and G4CoupledTransportation processes are tasked to select which of the tracks flagged as looping are killed and which survive. To balance the potential significant cost of integrating looping particles, three thresholds exist: - the ‘Warning’ Energy: a track with energy below this value that is found to loop is killed silently (no warning.) Above the ‘Warning Energy’, if a track is selected for killing then a warning is generated. - the ‘Important’ Energy: the threshold energy above which a track will survive for multiple steps if found looping. - the number of extra ‘tracking’ steps for important particles. These tracks will be only be killed only if they still loop more than this number of trial steps. ( So in effect the number of integration steps will be this number times the maximum number of steps allowed in G4PropagatorInField.

In versions of Geant4 from 7.0 up to release 10.4, Transportation did not examine the types of a charged particle - all types of particles were killed if they fulfilled the same criteria. A short message was written in the G4cout output that gave the energy and location of the killed track. This printout was under the G4VERBOSE flag, so it was suppressed if the G4_NO_VERBOSE configuration option was chosen at installation.

In Geant4 10.5 several changes have been implemented:

  • only stable particles are killed. ( Re-enabling the killing of unstable particles as an option is envisioned. )

  • each particle with energy above the warning energy which is killed generates a detailed warning (using G4Exception) with the full information about the particle location, the current volume and its material, and the particle momentum and energy.

  • for the first 5 tracks killed a detailed description is printed that describes the criteria and parameters which are used to decide what tracks are killed, and provides a first guidance regarding how to ‘save’ tracks by chaning the values of thresholds or by adopting different integration methods.

Below we discuss the different way in which a user can change the thresholds for killing ‘looping’ tracks, which criteria can be used to ensure that a track continues to propagate and for how many steps an ‘important’ track that is ‘looping’ can survive.

Two techniques are demonstrated below. An example of using them is available in the extended example field01, in the directory examples/extended/field/field01.

Using preset thresholds for killing loopers

This method is new in Geant4 release 10.5, and uses the G4PhysicsListHelper which has methods to choose a pre-selected set of parameter values. The choices are between a set each of low and high thresholds. Either one can be enabled by calling correspondingly method.

It is possible to select a set of pre-selected values of the parameters for looping particles using

New functionality in G4PhysicsListHelper, introduced in Geant4 release 10.5, enables the Transportation process chosen to be provided with this set of parameters. This reuses the AddTransportation method, which is called in each thread.

To configure with low values of the thresholds, appropriate for typical applications using low-energy physics, choose

#include "G4PhysicsListHelper.hh"

int main( int, char**)
{
  auto plHelper = G4PhysicsListHelper::GetPhysicsListHelper();

  plHelper->UseLowLooperThresholds();
  //
  // Use low values for the thresholds
  //    Warning  1 keV, Important 1 MeV, trials 10

  auto physList = new FTFP_BERT();
  // ...
}

The original high values of the parameters can be selected with a similar call

plHelper->UseHighLooperThresholds();
//
// Configures with the original (high) values of parameters. Currently:
//    Warning 100 MeV, Important 250 MeV, trials 10

and are chosen as starting points for energy-frontier HEP experiments.

The above sequence is demonstrated in the main() of field01.cc, part of the extended field examples ( examples/extended/field/field01 .)

This configuration method works only if modular physics lists are used, or if the AddTransportation() method is used to construct the transportation in a user physics list.

These calls must be done before a physics list is instantiated, in particular before G4PhysicsListHelper::AddTransportation() is called during the construction of a physics list. Else the configuration of the parameters does not occur.

These methods must be called before the physics is constructed - i.e. typically before G4RunManager ‘s Initialise method is called.

In order for this method to work the physics list must be constructed in one of two ways:

  • a preconstructed physics lists, from the list of recommended physics lists, or

  • the list must be constructed using the G4ModularPhysicsList and its AddTransportation method.

Note that in each thread the AddTransportation instantiates a single common transportation process which is then used by all particles types.

Users who build a physics list without making use of G4ModularPhysicsList and its AddTransportation method, are responsible to register a Transportation process to each particle type, and to set its parameters appropriately. This would allow the most finer grained control, and would also allow different thresholds to be chosen for different particle types.

Finer-grain control of the parameters for killing looping particles

A new feature to set any value to these parameters is introduced in Geant4 release 11.1 using the class G4TransportationParameters. If an instance of G4TransportationParameters exists, the constructor of G4Transportation will utilise the values it stored to inititalise its own parameters.

auto transportParams= G4TransportationParameters::Instance();

transportParams->SetWarningEnergy(  warningE );
transportParams->SetImportantEnergy( importantE );
transportParams->SetNumberOfTrials( numTrials );
G4cout << "Using G4TransportationParameters to set looper parameters."  << G4endl;

A couple of caveats exist. First is that its values will be used by default for all instances of G4Transporation and its derived classes: whether it is

  • the single instance typically registered for all particles (as is done in the modular physics lists), or

  • an instance created separately and registered by a user to one or more (charged) particles as a replacement.

Secondly, if it exists, the values of all the parameters that G4TransportationParameters stores are currently used to overwrite the existing default values in G4Transportation.

So, if you create G4TransportationParameters it is your responsibility to set the values of all of its parameters.

Full control of the parameters for killing looping particles

Whether you use one of the previous methods or not, it is possible to exercise full fine-grained control over each values for each type of particle separately.

The user can choose arbitrary values for the different parametes related to killing loopers and also refine the integration of charged particle propagation in particular volumes in order to eliminate or reduce the incidence of looping tracks.

This is also the only method which will work in all Geant4 versions since 7.0.

To obtain reliable configuration of the G4Transportation (or G4CoupledTransportation ) process in a potentially multi-threaded application, we configure it using a G4VUserRunAction. In particular such configuration can be undertaken in the BeginOfRunAction methods.

For example, to ensure that only looping particles with energy 10 keV are killed silently we change the value of the ‘Warning’ Energy

transport->SetThresholdWarningEnergy( 1.0 * CLHEP::keV  );

After this each time a (stable) looping track with energy over 1.0 keV is killed by this transportation process, it will generate a warning.

The second configurable energy threshold is labelled the ‘important’ energy and it enables tracks above its value to survive a chosen number of ‘tracking’ steps. They will be only be killed only if they are still looping after the given number of tracking steps.

These are demonstrated also in the F01RunAction’s ChangeLooperParameters method, which is called by the BeginOfRunAction.

To obtain the appropriate Transportation object for a particular particle type G4ParticleDefinition *particleDef; either obtain it manually obtain directly if we know its type

G4VProcess* partclTransport =
 particleDef->GetProcessManager()->GetProcess("G4Transportation");

auto transport= dynamic_cast<G4Transportation*>(partclTransport);

or write code which can adapt to different types to different types which inherit from G4Transportation

Listing 104 A method to find the transportation object for a particle type. It can also find the ordinary G4Transportation or the derived classes G4CoupledTransportation and the newest G4TransportationWithMsc.
G4Transportation* FindTransport( G4ParticleDefinition* particleDef )
{
  G4Transportation *transport;
  G4VProcess* cplTransport= nullptr, *transportMsc= nullptr;

  auto pm= particleDef->GetProcessManager();

  G4VProcess* ordTransport = pm->GetProcess("G4Transportation");
  transport= dynamic_cast<G4Transportation*>(ordTransport);

  if( !transport ) {
    // Maybe it is G4CoupledTransportation ...
    cplTransport= pm->GetProcess("G4CoupledTransportation");
    if ( vpCoupledTransport) {
      transport= dynamic_cast<G4Transportation*>(cplTransport);
    }
  }
  return transport;
}

or else use (or copy) the helper method F01RunAction::FindTransportation

auto transport= FindTransportation(G4Electron::Definition(), true);

This example method returns returns a `G4Transportation *`. ( whereas in release 11.0 it returned a pair `std::pair<G4Transportation*, G4CoupledTransportation*>`. )

Since Geant4 11.1, G4CoupledTransportation and the new G4TransportationWithMsc classes inherit from G4Transportation you can use common code to configure them (as shown in Listing 104).

In case a different Transportation type is used which does not inherit from G4Transportation, such as G4ITTransportation, and G4DNABrownianTransportation (both relevant for Geant4 DNA) similar code is required for each such class.

auto pm= particleDef->GetProcessManager();
G4VProcess* vpItTransport= pm->GetProcess("G4ITTransportation");
auto itTransport= dynamic_cast<G4ITTransportation*>(vpItTransport);

NOTE: Up to (and including) release 11.0 G4CoupledTransportation was an independent class, not inheriting from G4Transportation. If your application needs to be backward compatible with previous releases of Geant4 (including release 10.1 through 10.7 and 11.0) you must ignore this new inheritance relationship.

When using earlier version of Geant4 it was necessary to treat instances of G4CoupledTransportation separately:

Listing 105 Demonstration of applying cuts to an instance of G4CoupledTransportation.
G4VProcess* vProcCoupled= pm->GetProcess("G4CoupledTransportation");
G4CoupledTransportation* coupledTransport=
     dynamic_cast<G4CoupledTransportation*>(vProcCoupled);

coupledTransport->SetThresholdWarningEnergy( 1.0 * CLHEP::keV  );
coupledTransport->SetThresholdImportantEnergy( 1.0 * CLHEP::MeV );
coupledTransport->SetNumberOfTrials( 20 );

It is still possible to use this approach (listing Listing 105) if you are maintaining an application which must work both older versions and the current release 11.1.

However moving forward the code can now be simplified, as demonstrated in the next two code excerpts. First by obtaining a G4Transportation`, e.g. as in listing Listing 104.

Then using common code, as in listing Listing 106 to overwrite the thresholds in a G4Transportation (or derived) class as found in F01RunAction ‘s ChangeLooperParameters method

Listing 106 Adapted extract of the method ChangeLooperParameters from F01RunAction
void ChangeLooperParameters(const G4ParticleDefinition* particleDef )
{
  auto transport= FindTransportation( particleDef );

  // Since Geant4 11.1 the following code works for several transportation
  //  classes:
  //   - ordinary G4Transporation,
  //   - G4CoupledTransportation used for parallel worlds, and
  //   - G4TransportationWithMsc used to speed up charged particles.

  if( transport != nullptr )
  {
    // Change the values of the looping particle parameters of Transportation
    transport->SetThresholdWarningEnergy(  warningEnergy );
    transport->SetThresholdImportantEnergy(  importantEnergy );
    transport->SetThresholdTrials( numberOfTrials );
  }
}

Note that for all pre-configured and modular physics lists share a single Transportation process for all types of particles. So the parameters for killing loopers will be shared by all particle types in this case.

If this is not the desired behaviour, it is necessary to register a separate instance of the Transporation process for a particular type of particle. See the subsection ReplacingTransportation about how this can be done.

F01RunAction plays the role of a helper object, which holds the proposed (new) values of parameters, and which can allow them to be set, e.g., in the main() function

runAction->SetWarningEnergy(   10.0 * CLHEP::keV );

F01RunAction then forwards them to the Transportation object of each thread at the start of each run.

Using a helper object to forward parameter changes

Since the type of the transportation it can be useful to use a helper object to hold the desired values for the parameters (thresholds, number of iterations), and to forward them to the Transportation class.

This is demonstrated in the class F01RunAction and its ChangeLooperParameters method of the field01 extended example.

It copes with either transportation class, G4Transportation or a G4CoupledTransportation, and passes new values of parameters as needed.

In field01 the methods of F01RunAction

runAction->SetImportantEnergy( 0.1  * CLHEP::MeV );
runAction->SetNumberOfTrials( 30 );

which the run action passes to the G4Transportation or G4CoupledTransportation object registered for the electron in F01RunAction ‘s method ChangeLooperParameters.

How to replace the Transportation Process of a particle type

The most advanced use case of controling requires a separate instance of a G4Transporation (or G4CoupledTransporation or other process).

If you have configured your application to use G4TransporationWithMsc only for electrons

Currently in order to undertake this it is necessary to use the property of the transportation of being first in the process list and interact directly with the process manager:

Listing 107 Replacing the G4Transportation process for one particle type - electrons
G4ParticleDefinition* particleDef= G4Electron::Definition();
G4int verboseLevel= 0;

G4ProcessManager* procManager = particleDef->GetProcessManager();
auto plist = procManager->GetProcessList();

procManager->RemoveProcess(0);  // Remove the current Transport

auto transport = new G4Transportation(verboseLevel);
// Here we can adjust the parameters for this instance/particle, e.g.
transport->SetThresholdWarningEnergy(  30.0 * CLHEP::keV  );
transport->SetThresholdImportantEnergy( 3.0 * CLHEP::MeV );
transport->SetNumberOfTrials( 50 );

procManager->AddProcess(transport, -1, 0, 0);
// Add the new type of Transport(ation)

You can repeat this for positrons - but we recommend much lower thresholds for positrons as their annihilation will produce two 0.5 MeV gammas.

Avoiding loopers or reducing the incidence of looping particles

There are different ways to reduce the occurence of looping particles. This section will provide an overview, and refer the user to the detailed information on particle propagation in a magnetic field for details.

Volumes which have a strong field and contain vacuum, large air cavities or large volumes of gases are prime candidates for causing integration difficulties for low energy charged particles, which result in looping particles.

  • A very simple way to reduce the incidence of looping particles is to reduce the maximum step size which particles that interact very infrequently can travel. Geant4 attempts to estimate an effective maximum using the diameter of the world volume, and frequently the maximum step size is large if the experimental hall is used as the world volume and has large dimentions. A smaller value can be impose by using the method. G4PropagatorInField::SetMaximumStepSize()

  • Another ways is to change the maximum number of integration substeps. The default value is 1000, but it can be obtained from G4PropagatorInField

    auto *transportMgr = G4TransportationManager::GetTransportationManager();
G4PropagatorInField* propFld= transportMgr->GetPropagatorInField();
G4cout << " The maximum number of substeps for integration is "
       << propFld->GetMaxLoopCount() << G4endl;

    It can be also be changed simply by calling the corresponding set method, e.g.

    propFld->SetMaxLoopCount(2500);

  • Assign a separate G4FieldManager class to each such volume. It can use adapted methods for the integration of the ODEs of motion.

  • One solution is to use a helical stepper, such as G4HelixImplicitEuler or G4HelixHeum which are inherently for steps over multiple ‘turns’ of a helix-like track. Their reason d’ etre is the ability to use the helix solution as baseline ‘first-order’ like solution and treat deviations from this as something like pertubations.

  • A new integration driver G4BFieldIntegrationDriver, introduced in Geant4 10.6 samples the value of the magnetic field at the start of each step and using the estimated track curvature determines whether the current step will traverse an angler smaller or larger than \(2 * \pi\). For larger steps the hybrid-helical stepper G4HelixHeum is used, and for smaller steps the G4DormandPrince745 stepper is used. This driver is the default driver in Geant4 10.6, created by G4ChordFinder for magnetic fields, and when G4FieldManager ‘s CreateChordFinder method is called. Note that it is applicable only for charged particles in pure magnetic field.

  • An older approach, usable before Geant4 10.6, is to use the G4HelixMixedStepper. This also combines a helix stepper for large steps with a Runge-Kutta stepper for small and intermediate step sizes. It does this by checking the value of the field at the start of every integration. As a result it is less efficient than the new method G4BFieldIntegrationDriver.

Ensuring progress for particles

To allow better progress for looping particles, the default behaviour was changed in Geant4 10.6, so that a track’s parameters are changed within a Geant4 step after it has undertaken 100 integration substeps.

For the remainder of the current step the track ‘tightness’ parameter delta chord is relaxed, first by a factor of 2 at the hundredth step, and again by another factor of 2 after every 100 subsequent steps. The original value of delta chord is restored at the end of the current Geant4 step.

This will allow tracks which are in a tight spiral with a radius of curvature less than the user-defined delta chord to make substantial progress.

Note that this applies only to particles below the ‘Important’ energy threshold, which would be killed if their integration is not completed within a single Geant4 step.