ALICE TPC Readout Chamber (ROC)
I - The ALICE Time Projection Chamber (TPC)
In this example, we simulate the behavior of the Multi-Wire Proportional Chambers used during Run 1 and Run 2 of the LHC in the ALICE Time Projection Chamber (TPC) to detect fast ionizing particles passing through the TPC gas. The particle ionizes the gas along its track, thus creating electrons that are drifted along the axis of the TPC due to an electric field applied by a high voltage electrode at the middle of the TPC. The drifting electrons eventually reach after a few microseconds the Multi-Wire Proportional Chambers (MWPC) where they are amplified. These MWPC consist of a plane of anode wires, a plane of cathode wires and a plane of gating wires, to which a gating voltage can be applied (alternating voltages) to prevent the ions from drifting back to the drift volume. Under these layers of wires, the pad plane is where the actual signal is read out.
In what follows we first build the geometry of the MWPC and add the fields (sections II and III), then we present two scenarios: one for the plotting of electron and ion drift lines for a charged pion generating electrons which are drifted to the MWPC and create avalanches at the wires (sections IV-VI); and one for simulating the signal resulting from a single avalanche at one of the wires (section VII).
Working versions of the codes discussed below along with the Magboltz gas file can be found in the directory Examples/AliceTPC of the Garfield++ source tree.
II - Loading Magboltz gas files
First, we create a new gas medium at the desired temperature and pressure. The gas mixture used in the ALICE TPC is Ne-CO2-N2 (ratio 90/10/5).
MediumMagboltz gas("ne", 85.72, "co2", 9.52, "n2", 4.76);
// Set the temperature [K] and pressure [Torr].
gas.SetTemperature(293.15);
gas.SetPressure(750);
// Read the electron transport coefficients from a .gas file.
gas.LoadGasFile("my_gas_file.gas");
// Read the ion mobility table.
gas.LoadIonMobility("IonMobility_Ne+_Ne.txt");
III - Defining the MWPC layout
We adopt a coordinate system where the x-axis is the radial direction of the TPC (perpendicular to the wires), the y-axis is the axis of the TPC cylinder and the z-axis runs along the MWPC wires. We start with the component describing the electric field used for the drifting of electrons, which pass through the open gating. Later on, we will use another component for ions, which drift under closed gating. The layout has a repetitive structure along the x-axis with a period of 0.25 cm (corresponding to the spacing of the anode and cathode wires). The anode (sense) wires are 2 cm above the pad plane and the cathode wires are 2 cm above the anode wires.
ComponentAnalyticField cmpe;
cmpe.SetMedium(&gas);
// Distance between rows of wires [cm].
constexpr double gap = iroc ? 0.2 : 0.3;
// Periodicity (wire spacing) [cm].
constexpr double period = 0.25;
cmpe.SetPeriodicityX(period);
Next we add the wires with their voltages.
// Add the sense wires.
constexpr double xs = 0.;
constexpr double ys = gap;
constexpr double ds = 0.0020; // wire diameter [cm]
constexpr double vs = 1570.; // wire potential [V]
cmpe.AddWire(xs, ys, ds, vs, "s", 100., 50., 19.3, 1);
// Add the cathode wires.
constexpr double xc = 0.5 * period;
constexpr double yc = 2 * gap;
cmpe.AddWire(xc, yc, dc, 0., "c", 100., 50., 19.3, 1);
// Add the gate wires.
constexpr double xg1 = 0.25 * period;
constexpr double xg2 = 0.75 * period;
constexpr double yg = 2. * gap + 0.3;
cmp.AddWire(xg1, yg, dg, vg, "g", 100., 50., 19.3, 1);
cmp.AddWire(xg2, yg, dg, vg, "g", 100., 50., 19.3, 1);
Note that the trap radius (last parameter in the AddWire function) is a virtual radius within which the electron cannot escape. It is used to make sure a drift line does not "miss" a wire whose radius is too small. However, in certain cases like for the gating wires, it is preferable to keep the trap radius to a single multiple of the wire radius (i.e. to the actual wire radius), to avoid stopping the drift line prematurely at the gating grid.
Next, the pad plane and high-voltage electrode are added with respective voltages 0 and vHV.
cmpe.AddPlaneY(0., 0., "pad_plane");
cmpe.AddPlaneY(yHV, vHV);
The magnetic field (0.5 T) is parallel to the 400 V/cm drift field created by the HV electrode.
cmpe.SetMagneticField(0, 0.5, 0);
We repeat the above steps for the ions (creating another ComponentAnalyticField object named cmpi), using closed gating this time.
Finally we create a Sensor object and add the two components to it.
Sensor sensor;
sensor.AddComponent(&cmpe);
sensor.AddComponent(&cmpi)`
cmpe when drifting electrons, cmpi when drifting ions) and deactivate the other.
To simulate the drift lines we will use a DriftLineRKF object.
DriftLineRKF drift(&sensor);
In what follows, two scenarios are presented. In the first one, we simulate an ionizing particle and drift the electrons to the MWPCs. Then we introduce the ions formed by the avalanche process and plot their drift lines on the same graph.
The second scenario is a single avalanche occurring on one of the sensing wires. We are interested in plotting the signal on the pad plane. To go to this scenario, jump to section VII.
IV - Plotting electric field and drift lines
We plot the contours of the potential using the class ViewField.
// Plot isopotential contours.
ViewField fieldView(&cmpe);
fieldView.SetArea(xmin, 0., xmax, 5 * gap);
fieldView.SetVoltageRange(-400., 1000.);
fieldView.PlotContour();
We also create a ViewDrift object and pass it to DriftLineRKF.
// Plot the drift line.
ViewDrift driftView;
drift.EnablePlotting(&driftView);
Once the drift lines have been computed, we visualize them together with the geometry of the chamber.
ViewCell cellView(&cmpe);
cellView.SetArea(xmin, 0., xmax, 5 * gap);
cellView.Plot2d();
driftView.SetArea(xmin, 0., xmax, 5 * gap);
driftView.SetCanvas(cellView.GetCanvas());
driftView.Plot(true, false);
V - Shooting an ionizing particle
After all the setup we can finally start simulating particles. Below we show how to shoot an ionizing particle (here a 1 GeV / c pion), using the interface with HEED.
TrackHeed track(&sensor);
track.SetParticle("pi");
track.SetMomentum(1.e9);
track.NewTrack(xt, yt, 0., 0., 1., 0., 0.);
VI - Drifting Electrons and Ions
The code for electron drifting and avalanche creation along with ion drifting is shown below. The outer for loop gets the clusters formed by the ionization of the gas by the fast particle. Each cluster in turn contains several electrons (a primary and potentially some secondaries). The second for loop iterates over all these electrons within the same cluster. Next each electron is drifted to the MWPC. A "trick" is used to take into account the transverse diffusion, and at the same time reduce the running time of the program. If a particle is shot far away from the ROCs, it is brought to a y-coordinate of 1.2 cm (just above the gating grid) smeared in a Gaussian manner with σ related to the diffusion in space (dT represents the diffusion coefficient at the drift field value of the TPC [0.0198 cm1/2]). Once the electrons attain the sensing wires, ions are introduced with a Gaussian smearing in φ (angle around the wire). For practical purposes and for the sake of plotting only, a small gain is used.
for (const auto& cluster : track.GetClusters()) {
// Retrieve the electrons of the cluster.
for (const auto& electron : cluster.electrons) {
const double y0 = electron.y;
// Smear the coordinates to account for the drift up to the ROC.
constexpr double dT = 0.0198;
const double sigma = dT * sqrt(std::max(y0, 1.2) - 1.2);
const double x0 = electron.x + RndmGaussian() * sigma;
const double z0 = electron.z + RndmGaussian() * sigma;
sensor.EnableComponent(0, true);
sensor.EnableComponent(1, false);
drift.DriftElectron(x0, 1.2, z0, 0.);
double x1 = 0., y1 = gap, z1 = 0., t1 = 0.;
int status = 0;
drift.GetEndPoint(x1, y1, z1, t1, status);
sensor.EnableComponent(1, true);
sensor.EnableComponent(0, false);
if (y1 > 0.5 * gap && y1 < 1.5 * gap) {
// Electron drifted to a sensing wire.
// Sample the gain from a Polya distribution.
const int nIons = static_cast<int>(std::round(RndmPolya(theta) * gain));
for (int j = 0; j < nIons; ++j) {
constexpr double r = 0.01;
const double angle = RndmGaussian(0, 1.4);
drift.DriftIon(x1 + r * sin(angle), gap + r * cos(angle), 0, 0);
}
}
}
}
VII - Signal from a single avalanche
Here we show how to get both the raw and convoluted signal from the drifting ions (the electron signal is neglected here because it is very fast in comparison). First, we need to define the pad plane as an electrode and associate it with a readout to be able to collect the signal. To do this, we need to make sure to give the plane a non-empty string as a label when defining the cell layout (recall that we had labelled it "pad_plane"). Then after creating the Sensor and adding cmpi to it we need to add a line:
sensor.AddElectrode(&cmpi, "pad_plane");
Next, drift ions as in the previous section as follows:
for (int i = 0; i < nIons; i++) {
// Sample the starting point of the ion around the sense wire.
constexpr double r = 0.003;
const double angle = RndmGaussian(0, 1.4);
driftline.DriftIon(r * sin(angle), gap + r * cos(angle), 0, 0);
}
The signal is viewed using the ViewSignal class. Below, we show the code for the raw current signal generated by the drifting ions.
ViewSignal signalView(&sensor);
signalView.PlotSignal("pad_plane");
Now we convolute the signal with the delta response of the amplifier to get a voltage signal. To do this, we first define the impulse response function outside the "main" function as follows:
double transfer(double t) {
const double tau = 160; // peaking time in ns
const double fC_to_mV = 12.7; // mV per fC
// Impulse response of the amplifier.
return fC_to_mV * exp(4) * pow((t / tau), 4) * exp(-4 * t / tau);
}
Then we do the convolution and display the result in a plot:
sensor.SetTransferFunction(transfer);
sensor.ConvoluteSignal();
signalView.PlotSignal("pad_plane");