Low-energy pion detection with a Micromegas
In the ASACUSA antihydrogen experiment at CERN, antiprotons and positrons are trapped in an electromagnetic field and some of them recombine into antihydrogen, some of them annihilate. Most of the final state of antiproton annihilation consists of charged pions. The trap is placed inside magnetic coils to confine the particles.
To be able to track the pions a Micromegas detector is being utilised. The image below illustrates the Micromegas structure:
In the example we will simulate the drift and multiplication of electrons produced by a charged particle crossing the Micromegas and the resulting induced signal on the readout strips.
After putting together the code below the result of the simulation will look something like the image below. On the left the 2D x-y projection of the avalanche is shown, while on the right the contour plot of the electric potential is shown.
A working version of the code is available in the directory Examples/Asacusa of the Garfield++ source tree.
Gas mixture
We choose a mixture of 90% Argon and 10% isobutane, at room temperature and atmospheric pressure.
MediumMagboltz gas("ar", 90., "ic4h10", 10.);
Detector geometry
We will model the electric fields in the drift gap and amplification gap using two separate ComponentAnalyticField objects, each with a parallel-plate configuration.
ComponentAnalyticField cmpA;
ComponentAnalyticField cmpD;
cmpA.SetMedium(&gas);
cmpD.SetMedium(&gas);
// Amplification gap.
constexpr double yMesh = 0.01;
constexpr double vMesh = -500.;
cmpA.AddPlaneY(yMesh, vMesh, "mesh");
cmpA.AddPlaneY(0., 0., "bottom");
// Drift gap.
constexpr double yDrift = yMesh + 0.3;
constexpr double vDrift = -3000.;
cmpD.AddPlaneY(yDrift, vDrift, "top");
cmpD.AddPlaneY(yMesh, vMesh, "mesh");
The detector operates in a magnetic field of 2 Tesla perpendicular to the electric field.
constexpr double bfield = 2.;
cmpD.SetMagneticField(0, 0, bfield);
cmpA.SetMagneticField(0, 0, bfield);
We put the two ComponentAnalyticField objects together in a Sensor.
Sensor sensor;
sensor.AddComponent(&cmpD);
sensor.AddComponent(&cmpA);
We define three strips along z (perpendicular to the plane of the screen, as seen on the figure above) on the anode plane.
constexpr double pitch = 0.07;
constexpr double interstrip = 0.01;
constexpr double hw = 0.5 * (pitch - interstrip);
const double xStrip1 = hw;
cmpA.AddStripOnPlaneY('z', 0., xStrip1 - hw, xStrip1 + hw, "strip1");
const double xStrip2 = xStrip1 + pitch;
cmpA.AddStripOnPlaneY('z', 0., xStrip2 - hw, xStrip2 + hw, "strip2");
const double xStrip3 = xStrip2 + pitch;
cmpA.AddStripOnPlaneY('z', 0., xStrip3 - hw, xStrip3 + hw, "strip3");
and request the induced current to be calculated for them.
sensor.AddElectrode(&cmpA, "strip1");
sensor.AddElectrode(&cmpA, "strip2");
sensor.AddElectrode(&cmpA, "strip3");
Simulating a track
We can now start simulating particles. First we construct a TrackHeed object, and configure it to simulate a negatively charged pion with a momentum of 300 MeV/c.
TrackHeed track(&sensor);
track.SetParticle("pi-");
track.SetMomentum(300.e6);
To simulate the drift and multiplication of the electrons generated along the track we use the microscopic tracking method.
AvalancheMicroscopic aval(&sensor);
We now generate a track, retrieve the "clusters" along the track and simulate the drift/avalanche for each of the primary electrons. The drift lines of the electrons released in the drift gap will stop once they hit the mesh plane. This is of course an artifact due to the fact that we have modelled the detector using two parallel-plate chambers. We therefore need to move the electrons "by hand" to the amplification gap to simulate the avalanches.
track.NewTrack(xt, yt, zt, 0, 0, -1, 0);
// Loop over the clusters.
for (const auto& cluster : track.GetClusters()) {
for (const auto& electron : cluster.electrons) {
aval.AvalancheElectron(electron.x, electron.y, electron.z, electron.t,
0.1, 0., 0., 0.);
// Get the end point of the initial electron's trajectory.
const auto& p1 = aval.GetElectrons().front().path.back();
// Skip electrons that didn't hit the mesh plane.
if (fabs(p1.y - yMesh) > 1.e-6) continue;
// Move the electron into the amplification gap.
aval.AvalancheElectron(p1.x, yMesh - 1.e-6, p1.z, p1.t,
p1.energy, 0, -1., 0.);
}
}
Signal
Finally, we retrieve the collected charge on each of the strips and calculate the reconstructed hit position as the weighted average of the strip centre coordinates.
sensor.IntegrateSignals();
const double q1 = sensor.GetSignal("strip1", nTimeBins - 1);
const double q2 = sensor.GetSignal("strip2", nTimeBins - 1);
const double q3 = sensor.GetSignal("strip3", nTimeBins - 1);
const double sum = q1 + q2 + q3;
double mean = (q1 * xStrip1 + q2 * xStrip2 + q3 * xStrip3) / sum;