Modelling a GEM
In this example, we will simulate the drift of electrons and ions in a standard GEM, which consists of a 50 μm thick kapton foil coated on both sides with a 5μm layer of copper, with a hexagonal pattern of etched holes (outer hole diameter: 70 μm, inner hole diameter: 50 μm, pitch: 140 μm).
The source code of this example application can be found in the folder Examples/Gem of the Garfield++ source tree.
Creating a field map
Using ANSYS at CERN
To have access to ANSYS on lxplus, you have to add yourself to the parc-service-users e-group. The executable is located in /afs/cern.ch/project/parc/bin/. You may wish to include this directory in your PATH environment variable.
As a first step, we need to calculate the electric field in the GEM. Currently, with Garfield++, this calculation is usually be performed using an external (finite-element) field solver. In this example, we use Ansys, but the steps for importing and using field maps from other field solvers are very similar.
The technique for calculating electric fields using the finite element program Ansys is explained in the following set of slides:
The Ansys run should produce four text files:
- ELIST.lis: the list of elements with pointers to the material property table and to the node list;
- NLIST.lis: the list of nodes, with their position in space;
- MPLIST.lis: material property table
- PRNSOL.lis: estimated potentials at each of the nodes.
These files can be moved to any convenient location.
Importing the field map
In the following we assume that the field has been calculated using Ansys and that the resulting *.lis files are located in the current working directory.
Importing three-dimensional Ansys field maps and evaluating the electric field and potential by interpolation in these field maps is dealt with by the class ComponentAnsys123. The initialisation of ComponentAnys123 consists of
- loading the mesh (ELIST.lis, NLIST.lis), the list of nodal solutions (PRNSOL.lis), and the material properties (MPLIST.lis);
- specifying the unit in which distances are supposed to be given in these files;
- setting the appropriate periodicities/symmetries.
ComponentAnsys123 fm;
fm.Initialise("ELIST.lis", "NLIST.lis", "MPLIST.lis", "PRNSOL.lis", "mm");
// Set the periodicities.
fm.EnableMirrorPeriodicityX();
fm.EnableMirrorPeriodicityY();
// Print some information about the cell dimensions.
fm.PrintRange();
fm = ROOT.Garfield.ComponentAnsys123()
fm.Initialise("ELIST.lis", "NLIST.lis", "MPLIST.lis", "PRNSOL.lis", "mm")
# Set the periodicities.
fm.EnableMirrorPeriodicityX()
fm.EnableMirrorPeriodicityY()
# Print some information about the cell dimensions.
fm.PrintRange()
We need to apply mirror periodicity in x and y in ComponentAnsys123 since we had exploited the symmetry of the geometry and modelled only half a hole in Ansys.
Using the class ViewField, we do a visual inspection of the field map to make sure it makes sense. We first make a plot of the potential along the hole axis (z axis).
ViewField fieldView(&fm);
// Plot the potential along the hole axis.
fieldView.PlotProfile(0., 0., 0.02, 0., 0., -0.02);
fieldView = ROOT.Garfield.ViewField(fm)
# Plot the potential along the hole axis.
fieldView.PlotProfile(0., 0., 0.02, 0., 0., -0.02)
We also make a contour plot of the potential in the x − z plane.
constexpr double pitch = 0.014;
// Set the normal vector (0, -1, 0) of the viewing plane.
fieldView.SetPlane(0, -1, 0, 0, 0, 0);
// Set the plot limits in the viewing plane.
fieldView.SetArea(-0.5 * pitch, -0.02, 0.5 * pitch, 0.02);
fieldView.SetVoltageRange(-160., 160.);
fieldView.PlotContour();
pitch = 0.014
# Set the normal vector (0, -1, 0) of the viewing plane.
fieldView.SetPlane(0, -1, 0, 0, 0, 0)
# Set the plot limits in the viewing plane.
fieldView.SetArea(-0.5 * pitch, -0.02, 0.5 * pitch, 0.02)
fieldView.SetVoltageRange(-160., 160.)
fieldView.PlotContour()
Gas
We use a gas mixture of 80% argon and 20% CO2 at room temperature and atmospheric pressure (760 Torr).
MediumMagboltz gas("ar", 80., "co2", 20.);
// Set temperature [K] and pressure [Torr].
gas.SetTemperature(293.15);
gas.SetPressure(760.);
gas = ROOT.Garfield.MediumMagboltz("ar", 80., "co2", 20.)
gas.SetTemperature(293.15)
gas.SetPressure(760.)
Electron collision rates
In this example, we will track the electrons in the gas using a "microscopic" Monte Carlo method, based on the electron-atom/molecule cross sections in the Magboltz program. As discussed in more detail in the user guide, the algorithm takes as input the collision rates (as function of the electron's kinetic energy) for each scattering mechanism that can take place in the gas. The preparation of the tables of collision rates and the interpolation in these tables is done by the class MediumMagboltz, which - as the name suggests - provides an interface to Magboltz.
In MediumMagboltz the collision rates are stored on an evenly spaced energy grid. The max. energy can be set by the user. For avalanche calculations, 50-200 eV is usually a reasonable choice.
gas.SetMaxElectronEnergy(200.);
gas.Initialise(true);
gas.SetMaxElectronEnergy(200.)
gas.Initialise(True)
Penning transfer
Argon includes a number of excitation levels with an excitation energy exceeding the ionisation potential of CO2 (13.78 eV). These excited states can contribute to the gain, since (part of) the excitation energy can be transferred to a CO2 molecule through collisions or by photo-ionisation.
In the simulation, this so-called Penning effect can be described in terms of a probability r that an excitation is converted to an ionising collision.
// Penning transfer probability.
constexpr double rPenning = 0.57;
// Mean distance from the point of excitation.
constexpr double lambdaPenning = 0.;
gas.EnablePenningTransfer(rPenning, lambdaPenning, "ar");
# Penning transfer probability.
rPenning = 0.57
# Mean distance from the point of excitation.
lambdaPenning = 0.
gas.EnablePenningTransfer(rPenning, lambdaPenning, "ar")
Ion transport properties
Unlike electrons, ions cannnot be tracked microscopically in Garfield++. Moreover, there is no program like Magboltz that can compute the macroscopic transport properties (such as the drift velocity), so we have to provide these data ourselves. In this example, we use the mobility of Ar+ ions in Ar as an approximation because there is no literature data for drift in the mixture. Example files with mobility data for various pure gases are located in the Data directory of the source tree and are copied to the install directory during installation. If you don't provide a full path or if the mobility file is not in your current working directory, MediumMagboltz will look for it in the subdirectory share/Garfield/Data of the install folder.
gas.LoadIonMobility("IonMobility_Ar+_Ar.txt");
gas.LoadIonMobility('IonMobility_Ar+_Ar.txt')
Associating the gas to a field map region
In order to track a particle through the detector, we have to tell ComponentAnsys123 which field map material corresponds to which Medium object. The gas can be distinguished from the other materials (here: kapton and copper) by its dielectric constant, in our case ε= 1. The function
fm.SetGas(&gas);
fm.SetGas(gas)
sets the drift medium for all materials with ε = 1.
Electron and ion transport
Finally, we have to create a Sensor class, which is basically an assembly of components. In our case, the sensor has only one component.
Sensor sensor(&fm);
sensor = ROOT.Garfield.Sensor(fm)
The Sensor object acts as an interface to the transport classes.
Normally, particles are transported until they exit the mesh. To speed up the calculation we restrict the drift region to −100 μm < z < +250 μm.
sensor.SetArea(-5 * pitch, -5 * pitch, -0.01, 5 * pitch, 5 * pitch, 0.025)
sensor.SetArea(-5 * pitch, -5 * pitch, -0.01, 5 * pitch, 5 * pitch, 0.025)
Microscopic tracking of electrons is dealt with by the class AvalancheMicroscopic.
AvalancheMicroscopic aval(&sensor);
aval = ROOT.Garfield.AvalancheMicroscopic(sensor)
For tracking the ions, we use the class AvalancheMC, which takes as input the macroscopic drift velocity as function of the electric field and simulates drift lines using a Monte Carlo technique.
AvalancheMC drift(&sensor);
drift.SetDistanceSteps(2.e-4)`
drift = ROOT.Garfield.AvalancheMC(sensor)
drift.SetDistanceSteps(2.e-4)
We are now ready to simulate an electron avalanche in the GEM. We place the initial electron 200 μm above the centre of the GEM hole and set its initial energy to a typical energy in the electric field of the drift gap.
// Set the initial position [cm] and starting time [ns].
double x0 = 0., y0 = 0., z0 = 0.02, t0 = 0.;
// Set the initial energy [eV].
double e0 = 0.1;
// Set the initial direction (x, y, z).
// In case of a null vector, the direction is randomized.
double dx0 = 0., dy0 = 0., dz0 = 0.;
// Calculate an electron avalanche.
aval.AvalancheElectron(x0, y0, z0, t0, e0, dx0, dy0, dz0);
# Set the initial position [cm] and starting time [ns].
x0 = 0.
y0 = 0.
z0 = 0.02
t0 = 0.
# Set the initial energy [eV].
e0 = 0.1
# Set the initial direction (x, y, z).
# In case of a null vector, the direction is randomized.
dx0 = 0.
dy0 = 0.
dz0 = 0.
# Calculate an electron avalanche.
aval.AvalancheElectron(x0, y0, z0, t0, e0, dx0, dy0, dz0)
After the calculation, we can extract information such as the number of electrons/ions produced in the avalanche. We loop over all the electron trajectories, and calculate an ion drift line starting from the same initial point as the electron.
int ne, ni;
// Get the number of electrons and ions in the avalanche.
aval.GetAvalancheSize(ne, ni);
// Loop over the electrons in the avalanche.
for (const auto& electron : aval.GetElectrons()) {
// Simulate an ion drift line from the starting point of the electron.
const auto& p0 = electron.path[0];
drift.DriftIon(p0.x, p0.y, p0.z, p0.t);
# Get the number of electrons and ions in the avalanche.
ne, ni = aval.GetAvalancheSize()
# Loop over the electrons in the avalanche.
for electron in aval.GetElectrons():
p0 = electron.path[0]
drift.DriftIon(p0.x, p0.y, p0.z, p0.t)
We can subsequently retrieve the endpoint of the ion drift line using
// ...
drift.DriftIon(p0.x, p0.y, p0.z, p0.t);
const auto& p2 = drift.GetIons().front().path.back();
std::cout << "Ion stopped at (" << p2.x << ", " << p2.y << ", " << p2.z << ").\n";
// ...
drift.DriftIon(p0.x, p0.y, p0.z, p0.t)
p2 = drift.GetIons().front().path.back()
Visualizing the drift lines
To plot the electron and ion drift lines together with the geometry, we use the classes ViewDrift and ViewFEMesh. When setting up the AvalancheMicroscopic and AvalancheMC objects, we need to switch on the storage of the drift line points and attach a pointer to a ViewDrift object.
ViewDrift driftView;
aval.EnablePlotting(&driftView);
drift.EnablePlotting(&driftView);
driftView = ROOT.Garfield.ViewDrift()
aval.EnablePlotting(driftView)
drift.EnablePlotting(driftView)
After having calculated the electron avalanche and ion drift lines, we plot the trajectories together with the mesh.
ViewFEMesh* meshView = new ViewFEMesh(&fm);
meshView->SetArea(-2 * pitch, -2 * pitch, -0.02, 2 * pitch, 2 * pitch, 0.02);
// x-z projection
meshView->SetPlane(0, -1, 0, 0, 0, 0);
meshView->SetFillMesh(true);
// Set the color of the kapton.
meshView->SetColor(2, kYellow + 3);
meshView->EnableAxes();
meshView->SetViewDrift(&driftView);
meshView->Plot()`
meshView = ROOT.Garfield.ViewFEMesh(fm)
meshView.SetArea(-2 * pitch, -0.02, 2 * pitch, 0.02)
# x-z projection.
meshView.SetPlane(0, -1, 0, 0, 0, 0)
meshView.SetFillMesh(True)
# Set the color of the kapton.
meshView.SetColor(2, ROOT.kYellow + 3)
meshView.EnableAxes()
meshView.SetViewDrift(driftView)
meshView.Plot()
Animation
To visualize the evolution of an avalanche, it can be instructive to perform the simulation in time slices. The animation below starts with a single electron at 25 μm above the upper metal. The strong electric field in the GEM hole causes the electron to ionise gas molecules and start an avalanche. The electrons created in the avalanche (represented by orange dots) drift towards the readout pad (or the next GEM layer), while the ions (red dots) drift back through the GEM hole.
The source code for creating the above animation can be found in the directory Examples/GemMovie of the Garfield++ source tree. The program can either create a gif or a sequence of png files which can subsequently be converted to a video format using e.g. ImageMagick or MEncoder.
Note that in the Ansys script used for this animation, the metal volumes area meshed to show the full geometry.