Skip to content

Geant4 interface

Introduction

This page describes a way to create an interface between Geant4 and Garfield++. The interface uses the Geant4 physics parameterisation feature. Good starting points to understand this concept are the Geant4 physics process documentation and the parametrisation example code.

Before starting to describe the interface, a few principal points should be discussed. The normal use case of the interface is the simulation of primary ionization, charge transport and/or avalanches in a gas with Garfield++, but within a Geant4 simulation. Depending on the individual simulation, the point at which Garfield++ "takes over" from Geant4 varies. The major options concerning gas detectors are shown in the graph here and explained below.

Task sharing between Geant4 and Garfield

  • Relativistic charged particle enters thin absorber

The Heed interface of Garfield++ can be used to create the primary ionization clusters. Heed uses the Photo-Absorption-Ionisation Model to create the ionisation clusters. The Photo-Absorption-Ionisation model also exists in Geant4 as G4PAIPhotModel for both light and heavy charged particles. So the alternative approach is to let Geant4 create the primary ionization clusters, and then take Garfield++ for the charge transport and the avalanches. How to use the G4PAIPhotModel is shown in example TestEm8.

  • Non-relativistic electron and/or very thick absorbers

Here it get's a little more tricky. Heed tracks the electron without Coulomb scattering. The stopping power of the electron is continuous and based on its initial energy. This works fine for Heed's intended purpose, the tracking of relativistic charged particles in thin absorbers. For slower electrons in thicker absorber (a use case not foreseen), it leads to unrealistically straight electron tracks and an incorrect energy loss. The Geant4 G4PAIPhotModel model, although not intended for this use case either, works better here and shows in combination with the Geant4 Coulomb scattering realistic tracks.

  • Non-relativistic heavy charged particle (α, ion) in gas: Either the Garfield++ SRIM interface or in Geant4 the G4PAIModel can be tried.

  • Very low energy (< 1 keV) electron in gas: The microscopic tracking of Garfield++ can be tried.

  • Gamma particle in gas: Heed can be used to transport the photon.

Please refer to the Garfield++ documentation for details.

To summarize, depending on which type of particle and which energy the user wants to simulate, different forms of interaction between Geant4 and Garfield++ are possible. Either the particle that should be detected, is killed in Geant4 as soon as it enters the gas or is produced in the gas. The same particle is then recreated in Garfield++ using the NewTrack function of TrackHeed. This approach works well for relativistic particles.

The second option is to use a low energy physics model like the G4PAIPhotModel in the gas region in Geant4. The production of secondary electrons in the PAI and PAI photon model heavily depends on the production threshold. The figure below shows for a 1 MeV electron in Ar/CO2 70/30 for the PAI model the influence on the lower production cut on deposited energy and produced secondaries. For values greater than the minimum ionization potential, the deposited energy in Geant4 stays stable. To preserve the correct balance of energy, the production threshold should be set to a value between the minimum ionization potential and the W value.

Production cut

The comparison plot of the deposited energies in Geant4, Heed (energy lost in collisions and number of electrons × W) and the Geant4/Garfield++ interface (number of electrons × W) shows that agreement between Geant4 and the different simulation programs is achieved with a production cut between 21 eV and 24 eV. The exact threshold has to be determined for each gas mixture, particle type and kinetic energy by the user.

Threshold

The production threshold for electrons is set to around 100 eV. The particle that the user wants to detect thus creates electrons in the gas. These electrons are subsequently killed in Geant4, and recreated as delta electrons in Garfield++ using the function TransportDeltaElectron of TrackHeed.

Example

To illustrate the interface between Geant4 and Garfield++, we consider a simple example, based on the Geant4 basic example B4. The setup consists of a drift tube located behind an absorber, with a Garfield++-based user implementation of the detector response to be used inside the gas volume of the drift tube.

Detector description

We start with the Geant4 part of the interface. To describe our setup, we have to implement a class (called GarfieldDetectorConstruction in our example) derived from the abstract base class G4VUserDetectorConstruction. The method Construct() of this class is called during initialization and returns a pointer to the world physical volume.

Inside this function, we define a G4Region, corresponding to the logical volume of the gas-filled tube. In this G4Region, we want to use a Garfield++-based "parameterised model", while for the rest of the geometry (like the target), the regular Geant4 physics remains valid.

class GarfieldDetectorConstruction: public G4VUserDetectorConstruction {
  G4VPhysicalVolume* Construct() {
    // ...
    G4Region* regionGarfield = new G4Region("RegionGarfield");
    regionGarfield->AddRootLogicalVolume(fLogicalVolumeGas);
    // ...
    fGarfieldG4FastSimulationModel = new GarfieldG4FastSimulationModel(
        "GarfieldG4FastSimulationModel", regionGarfield);
    // ...
  }
};

The parameterised model is implemented in the class GarfieldG4FastSimulationModel, which is derived from G4VFastSimulationModel. The name G4VFastSimulationModel reflects the fact that parameterised models are usually simpler and thus faster than the full Geant4 tracking. A Garfield++ model will hardly be "fast" though, so this name should not be taken too literally in this context. The G4VFastSimulationModel has three pure virtual methods, which must be overriden:

class GarfieldG4FastSimulationModel: public G4VFastSimulationModel {
  // Decide (based on the particle name) if the Garfield model is applicable.
  G4bool IsApplicable(const G4ParticleDefinition& particleType);
  // Trigger the Garfield model or not depending on the kinetic energy. 
  G4bool ModelTrigger(const G4FastTrack& fastTrack);
  // The details of the Garfield model are implemented here.
  void DoIt(const G4FastTrack& fastTrack, G4FastStep& fastStep);
};

We will come back to the implementation of the DoIt method later.

As part of the initialization of the GarfieldG4FastSimulationModel class, we also set up the MediumMagboltz, ComponentAnalyticField and Sensor objects used in Garfield++ to describe the gas properties, detector layout, and electric field. This description needs of course to be consistent with the geometry and materials used by Geant4 (GarfieldDetectorConstruction).

  // ...
  // Define the gas mixture.
  fMediumMagboltz = new Garfield::MediumMagboltz("ar", 70., "co2", 30.);
  fMediumMagboltz->SetTemperature(293.15);
  fMediumMagboltz->SetPressure(760.);
  fMediumMagboltz->Initialise(true);
  // Set the Penning transfer efficiency.
  const double rPenning = 0.57;
  fMediumMagboltz->EnablePenningTransfer(rPenning, 0., "ar");
  fMediumMagboltz->LoadGasFile("ar_70_co2_30_1000mbar.gas");

  fComponentAnalyticField = new Garfield::ComponentAnalyticField();
  fComponentAnalyticField->SetMedium(fMediumMagboltz);
  // Wire radius [cm]
  const double rWire = 25.e-4;
  // Tube radius [cm]
  const double rTube = 1.451;
  // Voltages
  const double vWire = 1000.;
  const double vTube = 0.;
  // Add the wire in the center.
  fComponentAnalyticField->AddWire(0., 0., 2 * rWire, vWire, "w");
  // Add the tube.
  fComponentAnalyticField->AddTube(rTube, vTube, 0);
  fSensor = new Garfield::Sensor(fComponentAnalyticField);
  // ...

Physics list

As a next step we have to set up the physics list to be used. In our example, we create a class GarfieldPhysicsList derived from the class G4VModularPhysicsList. In the method ConstructProcess() of this class, we add a G4FastSimulationManagerProcess to the process list of the particles for which the parameterised model shall apply. If, for instance, we want the parameterised model to be applicable to muons, the implementation of the ConstructProcess() method could look like this:

class GarfieldPhysicsList: public G4VModularPhysicsList {

  void ConstructProcess() {
    // Call the base class method.
    G4VModularPhysicsList::ConstructProcess();
    // Add parameterisation.
    // Create an instance of the G4FastSimulationManagerProcess
    auto fastSimProcess_garfield = new G4FastSimulationManagerProcess("G4FSMP_garfield");

    auto theParticleIterator = GetParticleIterator(); 
    theParticleIterator->reset();
    while ((*theParticleIterator)()) { 
      G4ParticleDefinition* particle = theParticleIterator->value();
      G4ProcessManager* pmanager = particle->GetProcessManager();
      // Select particles for which the parameterisation applies.
      if (particle->GetParticleName() == "mu-" ||
          particle->GetParticleName() == "mu+") {
        // Add the G4FastSimulationManagerProcess as discrete process
        // (only PostStep actions).
        pmanager->AddDiscreteProcess(fastSimProcess_garfield);
      }
    }
  }
};

The G4FastSimulationManagerProcess is the physics process of the parametrisation. It serves as interface between the Geant4 tracking and the user parametrisation. At tracking time, the G4FastSimulationManagerProcess collaborates with the G4VFastSimulationManager of the current volume, and allows the user model to be triggered. As you might recall, the G4VFastSimulationManager was attached to the volume by defining a G4Region that contains the volume.

Garfield part of the interface

The Garfield++ part of the interface is contained in the method DoIt of the class GarfieldG4FastSimulationModel. In our simple example, we simulate simulate the track using Heed and calculate the deposited energy.

void GarfieldG4FastSimulationModel::DoIt(const G4FastTrack& fastTrack, 
                                         G4FastStep& fastStep) {
  // Get the local direction of the track (components of the momentum vector).
  G4ThreeVector localdir = fastTrack.GetPrimaryTrackLocalDirection();
  double dx = localdir.x();
  double dy = localdir.y();
  double dz = localdir.z();
  // Get the position of the track.
  G4ThreeVector localPosition = fastTrack.GetPrimaryTrackLocalPosition();
  // Convert the coordinates to cm since Garfield uses cm as standard unit.
  double x0 = localPosition.x() / CLHEP::cm;
  double y0 = localPosition.y() / CLHEP::cm;
  double z0 = localPosition.z() / CLHEP::cm;
  // Get the global time of the track.
  G4double t0 = fastTrack.GetPrimaryTrack()->GetGlobalTime();
  // Get the kinetic energy of the track in eV. 
  // Garfield uses eV as standard unit!
  double ekin_eV = fastTrack.GetPrimaryTrack()->GetKineticEnergy() / eV;
  // Get the particle name.
  auto particleName =
    fastTrack.GetPrimaryTrack()->GetParticleDefinition()->GetParticleName();    
  // Kill the track in Geant4 - from this point on Garfield takes over.
  fastStep.KillPrimaryTrack();

  // Total energy loss along the track
  double esum = 0.;

  // Initialize the Heed track with the information from the Geant4 fast track.
  Garfield::TrackHeed heed; 
  heed.SetParticle(particleName);
  heed.SetKineticEnergy(ekin_eV);
  heed.NewTrack(x0, y0, z0, t0, dx, dy, dz);
  for (const auto& cluster : heed.GetClusters()) {
    esum += cluster.energy;
  }
  G4cout << "Deposited energy: E=" << esum << "eV \n" << G4endl;

}

Running and configuring the example

An extended version of the above example is available in the directory Examples/Geant4GarfieldInterface of the Garfield++ source tree. In this version, the Garfield++ part of the interface is encapsulated in a class GarfieldPhysics. It also allows the user to configure the ionization model (PAIPhot or PAI model in Geant4 or Heed in Garfield++) to be used for the gas volume and the energy ranges for which the models are valid. This configuration is done in the macro physics.mac, which is called from the main macro run.mac. Different particles can be fired with the particle gun in run.mac, either from inside the gas tube or from outside. The material of the absorber can also be set in the macro file. The fired particle or one of its secondaries then ionizes the gas. The electrons are drifted towards the wire, where an avalanche is created when the electrons are about 50 μm away from the wire. The deposited energy, the sum of all avalanches and the gain are scored.

Advanced topic: How to get the electrons created in Garfield++ back into Geant4

In my opinion, there is actually no real reason do bring the electrons that were created in Garfield++ back into Geant4. Normally the scoring and filling of your ROOT histograms can directly be done in the Garfield++ part. But if for whatever reason you would like to do that in Geant4, the way to do it is by creating secondaries in the DoIt method of your class that was derived from the G4VFastSimulationModel. The secondary electrons are created using the CreateSecondaryTrack method of G4FastStep.

Keep in mind that the number of electrons created in avalanches will be very large, so that the creation of secondaries will further slow down the simulation. In the example below we create secondary electrons from the electrons of the primary ionization track.

void GarfieldG4FastSimulationModel::DoIt(const G4FastTrack& fastTrack, 
                                         G4FastStep& fastStep) {
  // ...
  heed.DisableDeltaElectronTransport();
  for (const auto& cluster : heed.GetClusters()) {
    // If delta electron transport is disabled, the cluster contains the 
    // "primary" ionisation electrons, i. e. the photo-electrons and 
    // Auger electrons. 
    for (const auto& electron : cluster.electrons) {
      // Position of the Garfield++ electron
      // Careful with units, Garfield++ uses cm, Geant4 mm
      G4ThreeVector position(10 * electron.x, 10 * electron.y, 10 * electron.z);
      // Momentum direction of the Garfield++ electron    
      G4ThreeVector momentumDirection(electron.dx, electron.dy, electron.dz);
      // Careful with units, Garfield++ uses eV, Geant4 MeV
      G4double eKin_MeV = electron.e * 1.e-6;
      // Time unit is ns in Garfield++ and Geant4
      G4double time = electron.t;
      // Check whether time returned by Garfield++ is larger than starting time of track
      if (te < globaltime) time += globaltime;
      // Create new dynamic particle
      G4DynamicParticle electron(G4Electron::ElectronDefinition(), 
                                 momentumDirection, eKin_MeV);
      // Create secondary electron
      fastStep.CreateSecondaryTrack(electron, position, time, true);
    }
  }
}

After the creation of the secondary electrons, these will become accessible in the SteppingAction or the SensitiveDetector of the Geant4 program. But careful, do not specify in the ModelTriggered method of the G4VFastSimulationModel that the model applies to very low energy electrons. Otherwise after the creation of the secondary electron, the G4VFastSimulationModel is triggered again, possibly resulting in an endless loop.

Contact

Dorothea Pfeiffer