Simulating a Resistive Plate Chamber

With its foundations laid in the 1980s by R. Santonico, R. Cardarelli et al., Resistive Plate Chambers (RPCs)¹² have been integrated into many of today's HEP experiments, e.g. in the ATLAS trigger system³. In this example, a step-by-step methodology will be suggested to obtain the induced current on a grounded electrode in a single-gap RPC in order to study the relevant physics behind it. A schematic representation of this geometry can be found in the figure below. For this, an extended form of the Ramo-Shockley theorem will be used, which will take the resistive layer's finite resistivity into account. Besides this, a grid-based avalanche statistics method efficiently simulates the Townsend avalanches in the uniform electric field present in the gas gap.

The dynamic weighting potential

Most particle detectors can, by approximation, be said to consist of perfectly conducting electrodes and insulating materials. In these geometries, the reciprocity theorem allows one to calculate the current induced by a moving charge carrier on any electrode by placing the electrode at a unit potential starting from a time, say, t=0. In addition, the charge carrier needs to be removed from the system, and all other electrodes are grounded. The resulting solution to the field equations is a static potential, called the "prompt" weighting potential, which, in turn, can be used to calculate the induced signal on the electrode under study. This is called the Ramo-Shockley theorem⁴⁵. This has been included in the class ComponentParallelPlate which can give the static weighting potential and applied electric field for a general N-layer parallel plate geometry as schematically represented by the figure below. Here the time-dependent dependent "delayed" weighting potential as described in [6,7] is not taken into account.

To set up the ComponentParallelPlate class, the relevant parameter values need to be given to its Setup function. These parameters are the total number of layers, the thickness of each layer, its relative permittivity, and the voltage applied on the top electrode. This is done by passing a series of vectors to the Setup function as shown, for example for an MRPC with six gas gaps, in the code below.

// The geometry of the RPC.
const int N = 15;  // Total number of layers inside the geometry

// Relative permitivity of the layers
const double epMylar = 3.1;     // [1]
const double epGlaverbel = 8.;  // [1]
const double epWindow = 6.;     // [1]
const double epGas = 1.;        // [1]

std::vector eps = {epMylar, epWindow, epGas, epGlaverbel, epGas, epGlaverbel, epGas, epGlaverbel, epGas, epGlaverbel, epGas, epGlaverbel, epGas,   epWindow, epMylar};

// Thickness of the layers
const double dMylar = 0.035;      // [cm]
const double dGlaverbell = 0.07;  // [cm]
const double dWindow = 0.12;      // [cm]
const double dGas = 0.025;        // [cm]

std::vector thickness = {dMylar, dWindow, dGas,  dGlaverbell, dGas, dGlaverbell, dGas, dGlaverbell, dGas, dGlaverbell, dGas, dGlaverbell, dGas, dWindow, dMylar};

// Applied potential
const double voltage = -15e3;  // [V]

ComponentParallelPlate *RPC = new ComponentParallelPlate();
RPC->Setup(N, eps, thickness, voltage);

With their relevant parameters, different readout structure can be added below the resistive layer:

• a plane using AddPlane,
• a strip using AddStrip,
• a pixel using AddPixel.

For this example we will focus on the induced signal of a plane readout labeled "Readout" and add it to the RPC object.

// Adding a readout structure.
const std::string label = "Readout";
RPC->AddPlane(label);

Sensor

Within a Sensor class, an assembly of components can be made, which can be used to calculate the induction of current on an electrode using the dynamic weighting potential from the components. In the context of this example, we need to point it to our RPC and its readout electrode.

// Create the sensor.
Sensor sensor;
sensor.AddComponent(RPC);
sensor.AddElectrode(RPC, label);

The resulting signal will be plotted as a histogram with the time-axis binned.

// Set the time bins.
const unsigned int nTimeBins = 200;
const double tmin = 0.;
const double tmax = 2;
const double tstep = (tmax - tmin) / nTimeBins;
sensor.SetTimeWindow(tmin, tstep, nTimeBins);

Setting the medium.

A gas mixture of 93.5% C₂F₄H₂, 5% iC₄H₁₀ and 1.5% SF₆ is used at room temperature and a pressure of 1 atm.

// Setup the gas, but one can also use a gasfile.
MediumMagboltz gas;
gas.LoadGasFile("c2h2f4_93-5_iso_5_sf6_1-5_bis.gas");

With the medium created one defines the geometry of the gas gap using SolidBox and pair the two with GeometrySimple. Here, the medium fills the gap's full height and is centered at the origin of the xy-plane with a width and depth of 10 cm.

// Setting the drift medium.
double totalThickness = 0.81; // [cm]
SolidBox box(0., totalThickness / 2, 0., 5., totalThickness / 2, 5.);
GeometrySimple geo;
geo.AddSolid(&box, &gas);

Subsequently this information can be passed on to the RPC, giving it the ability to let charge cariers drift in the gap.

RPC->SetGeometry(&geo);
RPC->SetMedium(&gas);

Grid based avalanche statistics calculation.

The electric field's uniformity inside the RPC's gas gap allows for a grid-based Monte Carlo simulation for the avalanche dynamics as the primary electrons propagate through the nodes of a lattice that represents the gas gap. Based on the constant Townsend and attachment coefficient, the AvalancheGrid class calculates the number of electrons that are added by ionization or subtracted from the avalanche by attachment as it moves from node to node. As the avalanches grow, the importance of space-charge effects grows with it. Modeling the resulting suppression of the growth of the avalanche, the total number of electrons in each gas gap will be bounded from above by 1.6 × 10⁷ by default. The saturation value can be modified using the SetMaxAvalancheSize function. This methodology is explained in great detail in the work of W. Riegler and C. Lippmann, and R. Veenhof⁸.

The AvalancheGrid class needs some initial information before it can function: the dimensions of the 3D lattice, the positions of the primary electrons and a reference to the Sensor class (to which we will come back later). The former can be achieved by using the SetGrid function giving the extrema and the number of nodes on the lattice in each dimension.

// Create the AvalancheGrid for grid-based avalanche calculations
int steps = totalThickness * 1e4;
AvalancheGrid avalgrid;
avalgrid.SetGrid(-0.05, 0.05, 5, 0.0, totalThickness, steps, -0.05, 0.05, 5);

Secondly, the initial positions of the primary electrons can be set by the function CreateAvalanche which will snap the particle to the nearest lattice node. Alternatively, one can use the AvalancheMicroscopic class to propagate the electrons throught the gas gap and letting the resulting Townsend avalanches evolve for a predefined amount of time and subsequently extract the relevant data from the class using the GetElectronsFromAvalancheMicroscopic function.

For this example we will opt for the AvalancheMicroscopic class aproach for the electron propagations and avalanche dynamics for the first 0.5 ns, set by the SetTimeWindow function, after which they will be imported into the AvalancheGrid lattice.

// Create the AvalancheMicroscopic.
AvalancheMicroscopic aval;
aval.UseWeightingPotential();

// Set time window where the calculations will be done microscopically.
const double tMaxWindow = 0.15;
aval.SetTimeWindow(0., tMaxWindow);

To interface with the AvalancheMicroscopic class, AvalancheGrid needs the pointer to it. This can be done by using the ImportElectronsFromAvalancheMicroscopic function.

avalgrid.ImportElectronsFromAvalancheMicroscopic(&aval);

Thirdly, a pointer ought to be given to the Sensor object in order to compute the induced signal and retrieve the Townsend coefficient, attachment coefficient and electron drift veleocity of the medium.

avalgrid.SetSensor(&sensor);
aval.SetSensor(&sensor);

Creating events.

The induced signal will be sourced by a high-energy pion tracking through the detector, ionizing the gas and letting the resulting primary electrons from a Townsend avalanche. This will be done using the recipe described in the Examples/Silicon. For this we will use the TrackHeed class to simulate a pion with a momentum of 10¹¹ eV/c.

// Set up pion track.
TrackHeed track(&sensor);
track.SetParticle("pion");
track.SetMomentum(1.e11);

As for its position equation, the pion will traverse the z-axis in the negative direction starting from the top of the gas gap. One such event is shown below; the green line indicates the pion trajectory and the yellow dots the ionization of the gas. The information of the resulting primary charges need to be passed to the AvalancheMicroscopic object, which will evaluate the avalanche dynamics within the aforementioned time window. After this the AvalancheGrid object takes the positions of the electrons in the gap at time t = tMaxWindow.

// Simulate pion track.
track.NewTrack(0., 0., gap, 0., 0., 0., -1);
// Retrieve the clusters along the track.
for (const auto& cluster : track.GetClusters()) {
  // Loop over the electrons in the cluster. 
  for (const auto& electron : cluster) {}
    // Propagate avalanche though gap within the time window set.
    aval.AvalancheElectron(electron.x, electron.y, electron.z, electron.t, 0.1);
    // Stops calculation after tMaxWindow ns.
    avalgrid.GetElectronsFromAvalancheMicroscopic();
  }
}

In order to keep track of in which gas gap each electron is the AsignLayerIndex function is used. This allows the saturation of the total amount of charge in each gas gap due to space-charge effects to be applied to each gas gap separately.

 // Start grid based avalanche calculations starting from where the microsocpic
 // calculations stoped.
avalgrid.AsignLayerIndex(RPC);

Lastly, the AvalancheGrid object can start its calculation on the evolution of the avalanches as it traverses the remaining distance to the resistive layer.

 avalgrid.StartGridAvalanche();
 ```

## The induced signal.

After the event has been simulated and the induced signal sourced in the `Sensor` object under its label it can be plotted using the Plot function of the signalView.

```cpp
 // Preparing the plotting of the induced signal of the electrode readout.
 TCanvas* cSignal = new TCanvas("cSignal", "", 600, 600);
 ViewSignal* signalView = new ViewSignal(&sensor);
 signalView->SetCanvas(cSignal); 
 signalView->Plot(label, true);
 cSignal->Update();
 gSystem->ProcessEvents();

The resulting signal, shown below, is plotted in blue along with its two constituite components: the prompt component in red and the delayed component in blue. This result can subsequently be exported as a .cvs file for further analysis using the function ExportSignal.

sensor.ExportSignal(label, "Signal");

Integrating the signal reveals the total induced charge on the electrode as a function of time and can be plotted and exported in a similar manner after using the ConvoluteSignals function.

sensor.ConvoluteSignals();

A working example can be found in Examples/RPC.

Contact

Djunes Janssens

---

1. R. Santonico , R. Cardarelli, "Development of Resistive Plate Counters", Nucl. Instrum. Meth. 187 (1981). ↩

2. R. Santonico , R. Cardarelli, A. Di Biagio and A. Lucci,"Progress in Resistive Plate Counters", Nucl. Instrum. Meth. A 263 (1988). ↩

3. ATLAS, "ATLAS muon spectrometer: Technical design report", CERN-LHCC-97-22, ATLAS-TDR-10, (1997). ↩

4. S. Ramo, "Currents induced in electron motion", PROC. IRE 27, 584 (1939). ↩

5. W. Shockley, "Currents to Conductors Induced by a Moving Point Charge"., Journal of Applied Physics. 9 (10): 635 (1938). ↩

6. W. Riegler, "Extended theorems for signal induction in particle detectors VCI 2004", Nucl. Instrum. Meth.535 (2004). ↩

7. W. Riegler, "Electric fields, weighting fields, signals and charge diffusion in detectors including resistive materials", JINST 11 (2016). ↩

8. W. Riegler and C. Lippmann and R. Veenhof, "Detector physics and simulation of resistive plate chambers", Nucl. Instrum. Meth. A (2003). ↩