PARSIFAL: a toolkit for triple-GEM parametrized simulation

PARSIFAL (PARametrized SImulation by Farinelli And Lavezzi) is a fast and reliable software tool that reproduces the complete response of a triple-GEM detector to the passage of a charged particle, taking into account the main physical effects. Starting from the detector configuration and the particle information, PARSIFAL reproduces ionization, spatial and temporal diffusion, effect of magnetic field, if present, and GEM amplification to provide the dependable triple-GEM detector response. In the design and optimization stages of this kind of detectors, simulations play an important role. Accurate and robust software programs, such as GARFIELD++, can simulate the transport of electrons and ions in a gas medium and their interaction with the electric field, but they are CPU-time consuming. The necessity to reduce the processing time while maintaining the precision of a full simulation is the main driver of this work. For a given set of geometrical and electrical settings, GARFIELD++ is run once-and-for-all to provide the input parameters for PARSIFAL. Once PARSIFAL is initialized and run, it produces the detector output, including the signal induction and the output of the electronics. The results of the analysis of the simulated data obtained with PARSIFAL are compared with the results of the experimental data collected during a testbeam: some tuning factors are applied to the simulation to improve the agreement. This paper describes the structure of the code and the methodology used to match the output to the experimental data.


Introduction
In 1997 F. Sauli [1] invented a new technique for signal amplification in gaseous detectors: the Gaseous Electron Multiplier (GEM). It consists of a 50 µm thick Kapton foil, copper coated (5 µm) on both sides, with a high density pattern of holes with a diameter of 50 µm and a pitch of 140 µm ( Fig. 1, left, middle). The application of an electric potential difference of a few hundreds of volts between the two faces creates an electric field of about 100 kV/cm. Such an intense field is enough to start the avalanche multiplication of the electrons, generated by the ionization of the specific gas mixture, drifting inside the holes (Fig. 1, right). Usually, multiple stages of amplification are stacked together in order to obtain the proper gain applying a lower voltage difference on each stage, thus lowering the discharge probability [4]. Figure 2 shows the full design of a triple-GEM, which consists of a cathode, three stages of GEM and an anode, segmented in strips or pads, for the signal readout. For the design and optimization of a detection system, possessing a tool to describe and predict the device behavior and response is of great importance. Moreover, this tool should grow from the experimental measurements and have the capability to go beyond them. The most used and robust existing software to simulate gaseous detectors is GARFIELD++, which The copper thickness t is 5 µm and the Kapton thickness T is 50 µm [2]. A high voltage difference is applied between the copper layers and the generated electric field lines inside the holes are shown in the right plot [3]. "is an object-oriented toolkit for the detailed simulation of particle detectors which use a gas mixture or a semiconductor material as sensitive medium" [5]. It comprises additional packages, such as: • HEED [6], for the primary and secondary ionizations; • MAGBOLTZ [7], to describe the electron diffusion effect and the drift in magnetic and electric fields.
The input map of electric fields and the geometric description can be provided by external tools, e.g. ANSYS [8]. The interactions between electrons and gas molecules are described at microscopic level. A triple-GEM simulation performed by GARFIELD++ is really CPU-time demanding (around one day for one single track on a standard distributed computing farm). For an extended study with high statistics, it is mandatory to reduce the time by several orders of magnitude. Pre-existing studies on this topic approached the triple-GEM simulation taking into account the main physics processes independently [9]. In this paper, a new stand-alone software is described, which, starting from these premises, can fully simulate a triple-GEM in a fast and reliable way. The reconstruction algorithms were applied to the simulation output and the performance of the simulated detector was measured.

Common setup between simulations and tests
While easily generalizable to different geometries and configurations, the characteristics of the simulated triple-GEM, used for this development, resemble the ones used in experimental measurements. As shown in the layout of Fig. 2, four gaps are defined by the positioning of the electrodes: • one drift gap, between the cathode and GEM 1, of 5 mm; • two transfer gaps, between GEM 1 and GEM 2, GEM 2 and GEM 3, of 2 mm; • one induction gap, between GEM 3 and the anode, of 2 mm.
The gas mixture filling the detector is Ar + 10% iC 4 H 10 . The high voltage applied to each GEM foil is 275 V. The electric fields in the different gaps are 1.5/2.75/2.75/5 kV/cm respectively. The used beam of ionizing particles was composed of muons with momentum of 150 GeV/c. Both cases with magnetic field 1T and without were considered, as well as different incident angles between the incident muons and the GEM chamber. The segmented anode has a double view readout, with orthogonal x y strips, having a pitch of 650 µm (x active area width = 580 µm, y active area width = 130 µm). The experimental tests were conducted in the H4 line of SPS at CERN North Area, within the RD51 collaboration [10]. The readout was performed by APV-25 ASIC [11]: it samples the signal 27 times, one every 25 ns and the shaping time is τ = 50 ns. The experimental data were analyzed and the results can be found in [12]. After the validation, PARSIFAL can be extended to further configurations in terms of gas mixtures, geometric and electrical settings, magnetic fields as well as different types of particles and energies.

Simulation Methodology
The variables used in PARSIFAL are extracted from GARFIELD++ simulations, as presented in the following sections. The simulation is composed of four independent parts, which describe ionization, diffusion, amplification and induction. This approximation is proven to be valid in a low rate environment. In order to obtain results comparable to the experimental data, a description of the response of the APV-25 ASIC is included in the code. Moreover, the reconstruction of the particle position is implemented. It is a well known issue in the MPGD (Micro Pattern Gaseous Detector) community that GARFIELD++ evaluation of the (single) GEM gain is around a factor two lower than the real one, so a tuning factor is expected to be necessary, as will be described later.

Ionization
Gas-based detectors exploit ionization to highlight the passage of particles: the interaction of the charged radiation with the gas atoms generates primary electrons which, if they carry enough energy, further ionize the atoms creating a cluster of secondary electrons. The former process follows a Poissonian statistics, thus it can be simulated by extracting the inter-cluster distances from a proper Exponential function. The latter process was studied with GARFIELD++ simulations and the number of secondary electrons in each cluster was tabulated up to 100 e − /cluster; these numbers were cross-checked with published tables in literature [14].
Concerning the GARFIELD++ simulations, ten thousand muons were shot into 5 mm of gas: the mean number of primary ionizations per centimeter was extracted and used as a fundamental parameter in PARSIFAL input along with the number of secondary ionizations per cluster. The distribution of the number of primary ionizations, which follows from the simulation procedure, is shown in Fig. 3: it can be fitted with a Poissonian function, as expected.

Drift of electrons
In this section, the positive direction z coincides with the direction of motion of the electrons, from the cathode to the anode, while the x reference axis is set along the anode strips. The transport of electrons is dominated by the gas characteristics which affect the diffusion, the presence of an electric field and, possibly, of a magnetic field. When the magnetic field is present, it is set orthogonal to the xz plane and its presence affects only the x coordinate, not the y. Studies on these effects were run separately for each gap. Each gap is delimited by two electrodes, as described in Sec. 2. For each gap, ten thousand electrons were drifted through the gap up to an imaginary x y plane, placed 150 µm before the final electrode and parallel to it. The distribution of the arrival positions on this plane with respect to the initial one was considered.
In the drift gap, the starting positions of the electrons were randomly distributed along the whole gap, up to the plane placed 150 µm before GEM 1, since the primary electrons from ionization are generated along the whole drift gap. The diffusion of the electrons in the drift gap depends then on their starting position: for this reason, the drift gap has been divided into several horizontal slices along the distance from the cathode. For each slice, the position distribution on the final plane was considered and fitted to evaluate the mean and sigma values. The behavior of the mean vs distance (as well as the sigma vs distance) has been fitted with a polynomial and the parameters were plugged into PARSIFAL (see Fig. 4).
Conversely, the z coordinate of the entrance point of the electrons in the gaps following each GEM  position is fixed. For the simulation of the two identical transfer gaps, the electrons were generated isotropically on a plane placed 150 µm after the GEM foil and drifted up to the final plane, placed 150 µm before the next GEM. The distribution on the final plane is shown in Fig. 5, left: it has been fitted with a Gaussian function and the mean and sigma values have been used as input parameters of PARSIFAL. In the induction gap, an analogous simulation was run, but the drift from a plane 150 µm after GEM 3 to the anode plane was considered. The resulting distribution is shown in Fig. 5  Similarly, the distributions of the drift time of the electrons were studied and inserted among the PARSIFAL parameters. All the simulations were run both with a magnetic field of 1T and without magnetic field. Figures 4 and 5 are related to the case with magnetic field, as they both present a shift of the mean value of the distribution under the effect of the Lorentz force. In the case without the magnetic field, the transverse diffusion is still present but all the distributions are centered around zero.

GEM properties
A single foil of GEM is fully characterized by its transparency and its gain, as it will be described in the following. In order to have the highest gain, all the electrons should go through the GEM holes where the multiplication happens. Some of the electric field lines which should drive the electrons to the holes, however, fall on the sides of the GEM, hence not all the electrons drifting towards the GEM contribute to the avalanche multiplication. The collection efficiency coll is defined as the number of electrons entering the hole of a GEM, divided by the number of electrons impinging on the GEM. All the electrons which enter the hole undergo multiplication, but not all of the electrons in the avalanche leave the GEM. Also here, some of the electric field lines drive some electrons on the surface of the GEM itself. The extraction efficiency extr is defined as the number of electrons leaving the GEM divided by the number of electrons present in the avalanche. The collection efficiency depends mostly on the field of the gap before the GEM, while the extraction efficiency on the field of the gap after the GEM. The product between the collection and extraction efficiencies is the transparency of the GEM. The intrinsic gain of a GEM is the number of electrons generated in the avalanche per single electron entering the hole. Also an effective gain can be defined, by multiplying the intrinsic gain by the transparency. The effective gain roughly gives the final number of electrons one can expect after the multiplication stage, given the initial number of electrons hitting the GEM. The gain value has an exponential dependence on the high voltage applied to the GEM [2]. The GARFIELD++ simulations to evaluate all these properties were conducted for each GEM separately. Ten thousand electrons have been simulated with GARFIELD++ from a plane placed 150 µm before the GEM to another plane placed 150 µm after the same GEM, with the avalanche microscopic description switched on; all the parameters have been extracted from these results. The fluctuations of the gain of a single GEM are described by the Polya distribution P(G): where G is the mean intrinsic gain, θ is a parameter connected to the variance of the distribution, C 0 is a constant and Γ is the Gamma function. In order to evaluate the full triple-GEM effective gain, one million complete simulations were run. In each simulation, one electron was followed, considering the probability of entering GEM 1 (collection efficiency), sampling the intrinsic gain G from the Polya distribution of GEM 1 as shown in Eq. 3.1, and for each electron in the avalanche applying a probability to be extracted. The same steps have been reproduced on the other two GEMs, simulating a full chain: Eventually, a histogram representing the effective gain of the full triple-GEM is filled and used in PARSIFAL to sample the gain during simulation (Fig. 6).

Induction
The signal in the triple-GEM is readout by the strips on the anode plane, connected to the electronics. The movement of the electrons in the induction gap, from the moment they exit GEM 3 to the one they arrive on the anode, induces a current on the strips which depends on the electron drift velocity and on the distance of the electron from the strip. Once all the electrons have arrived on the anode, the signal is over. The instantaneous current i ind (t) induced on the i th -strip is given by the Shockley-Ramo theorem [15,16] as: where e is the electron charge, v drift is the drift velocity and E w is the weighting field, which is a computational artifice and corresponds to the electric field generated by the electrode under consideration, when kept at 1 V, with all the other electrodes set to 0 V. Its analytical calculation (from Ref. [17]) is reported in Eq. 3.4 where E 1z is the z component of the weighting field, V 1 = 1 V, x and z are the coordinates of the point where the calculation is performed, on and orthogonal to the anode respectively, D is the induction gap thickness and w is the strip size. This calculation is just an approximation since it does not consider the non-active area between the strips (pitch ≡ active area) and it does consider only one view. Correction factors are foreseen to account for both the approximations. Another approximation which was used in PARSIFAL is to consider for simplicity the drift velocity as constant. Moreover, the weighting potential was used instead of the field and the calculation of the induced current from Eq. 3.3 was replaced by the computation of the induced charge Q ind (t): where ∆V w (t) is the gradient of weighting potential and φ is the angle between the weighting field and the drift velocity. The weighting potentials for strip pitches of 130 µm and 580 µm (dimensions One important property of the electronic signal is that it ends once all the electrons are collected by the strip and that the total charge induced on the strip corresponds to the total number of electrons actually falling on the strip. Following this assumption, a fast induction was added to PARSIFAL, in order to obtain a correct resulting signal with a much faster method. In the fast induction, for each electron entering the induction gap, the strip of the arrival of the electron is identified with a simple extrapolation taking into account the starting position at the entrance of the induction gap, the magnetic field and the diffusion effect, and the electron contributes to its charge by one electronic charge, released at the time of arrival of the electron on the anode. The consistency of the inducted/collected charge on the i th -strip is checked by simulating the induction of the same electron with both methods. The optimal matching between the two results shows that the charge can be inferred by the fast induction, with a huge improvement in computing time, being the fast induction around 30 times faster than the full one.

Electronics
As already said, also the electronics plays a role in the measurement and so a simulation of the same electronics used in testbeams, the APV-25 ASIC, was implemented in the code. Each APV-25 channel was simplified as a pre-amplifier and a shaper. The pre-amplifier was sketched as a simple integrator, without any electronic gain implemented. It integrates the signal continuously, so our approximation provides the integration every 1 ns (see Fig. 8(b)). The shaper is simplified as a CR-RC circuit, with shaping time τ = 50 ns. In order to compute the shaped integrated charge Q shaped (t) for each APV-25 time bin (27 time bins, each 25 ns long), the following shaping function is applied: with Q preamp as the integrated charge from the pre-amplifier, t as time bin time, t 0 as time of beginning of the signal. For each time bin t i , the electronic signal Q APV-25 is the the sum of all the shaping functions coming from previous time bins, evaluated at the time t i : For each strip, a random noise is added, sampling a Gaussian function centered at zero and with a sigma evaluated from the data (see figure 8(c)).

Reconstruction of data
As shown in Fig. 8(c) and (d), the simulated signal resembles with an optimum accuracy the real one and can be reconstructed in the same way. The hit charge, i.e. the charge measured by each strip, comes directly from the maximum value of the Q APV-25 histogram (see Fig. 8(c)), while the hit time comes from a fit of the rising edge of the signal with a Fermi-Dirac function [12]. A strip is firing if the measured charge is above a certain threshold Q thr . When more adjacent strips are firing, they can be grouped to form a cluster. The total cluster charge as well as the multiplicity of strips in one cluster, defined as cluster size, are two of the variables used to evaluate the goodness of the simulation when compared to real data (see section 5). The other two important quantities are the position reconstructed via the charge centroid (CC) and via the micro-Time Projection Chamber (µ-TPC) methods. The CC method computes the cluster position x CC as the average of the positions x i of the firing strips, weighted by the charge q i measured by each of them, as following: The µ-TPC method uses the drift gap as a time projection chamber of few millimeters: by multiplying the drift velocity, known from GARFIELD++ simulations, and the measured signal time on each strip, from the Fermi-Dirac fit, the bi-dimensional position of each primary ionization in the drift gap can be computed. By fitting all these points, the position x µ-TPC , where the charged incident particle passed in the middle of the drift gap, can be computed from the line fit parameters, as following: where gap is the drift gap thickness, a is the slope and b the intercept of the fitting line (see Fig. 9).

Tuning to experimental data
The first requirement from a detector simulation is to replicate the experimental results with a very good accuracy: this makes the validation of the code a necessity to guarantee a reliable prediction.
Since some approximations are used in the simulation, their impact is discussed in Sec. 5.1, before the actual description of the tuning procedure (Sec. 5.2).

General comments on simulation
PARSIFAL, in its current layout, can simulate one track outcome in less than two seconds, which is a great improvement with respect to GARFIELD++ full simulation, which takes about one day to simulate the same event. Further improvements can be added, for example to exploit the parallel computing about the drift and avalanche formation of different electrons, hence enhancing the CPU-time gain even more. It must be stressed that this good time performance is the outcome of some approximations: • the origin of secondary electrons is the same, in position and time, of the primaries; • the signal is created only by the electrons from ionization in the drift gap: in reality, about 2% of the signal is due to electrons created in the first transfer gap, which undergo only two multiplication stages [18]. For now, this contribution has been neglected; • in full induction, the drift velocity is considered as constant (35 µm/ns, in drift gap); • in full induction, the weighting field/potential is computed analytically with only one view while a double view anode is considered in this paper. To cope with this discrepancy a multiplication factor to account for the charge sharing was introduced; • in the fast induction, the time of arrival of each electron on the strip is considered as the time associated with each charge deposition, while the induction of the signal happens during the entire drift to the induction gap.
It is worth to mention here an issue which is known in literature: the mismatch between the simulated and the experimental gain in a GEM. Simulations performed with GARFIELD++ show a gain which is systematically lower than the real one by a factor ranging from 1.5 to 2.5. This leads to the necessity to add a tuning factor to the simulation to obtain a perfect matching between the simulated and measured charge. This issue will be treated in section 5.2.

Tuning procedure
In order to check and tune the results obtained by the simulation, the experimental data collected with the setup described in section 2 were used and four sentinel variables were chosen as indicators of the goodness of the accordance: the cluster charge, the cluster size, the spatial resolution on the position reconstructed with CC and with µ-TPC methods. A first version of the tuning is focused on two main elements: the gain value and the spatial diffusion. The gain, as already said, needs a correction since it is based on GARFIELD++ simulations which are affected by the factor two underestimation issue. Moreover, since the cluster size of the simulated data is lower w.r.t. the experimental one and a since there is a mismatch in the µ-TPC reconstructed position resolutions, also a tuning factor to the spatial diffusion was considered.
The conversion factor, to convert the ASIC ADC steps to fC, is fixed to the testbeam value, i.e. 30 ADC = 1 fC [19]; in the same way, the charge threshold was not tuned (45 ADC). The noise level, estimated from random trigger runs to be around 15 ADC, was not tuned as well. The optimization of all these values could enter a more detailed version of the tuning procedure in the future.
The tuning procedure consists of spanning over a set of parameter value combinations among which the one which better matches the experimental data is selected. The simulations are run for each combination, on a set of seven incident angles [0, 5,10,15,20,30,40] degrees and the global χ 2 is computed to finally find the parameters which provide the minimum value. The global χ 2 is the sum of the χ 2 p on the single parameters in order to perform the minimization on the four sentinel variables simultaneously: where, for each parameter p, the corresponding χ 2 p defined as in Eq. 5.2 In Eq. 5.2, only the error on experimental data is taken into account since the statistical error on the simulation is negligible. The samples are in fact of ten thousand tracks each and this keeps the statistical error around 1%, while the error on experimental data was evaluated to be around 5% for the spatial resolution and 10% for cluster size and charge. The simulations were run with the fast induction and without magnetic field. The range [1,8] was scanned for the tuning factor of the gain, with steps of 0.1, while the range from [1,3] was scanned for the tuning factor of the diffusion, with steps of 0.1. The best combination of parameters scored a χ 2 r ed ∼ 3, for a gain tuning factor of 6.8 and a diffusion tuning factor of 1.5. The comparisons between the simulated and experimental data are shown in Fig. 10.

Double-view readout
The charge collection on both views was developed adding the diffusion of the electrons on the y view and duplicating the induction of the charge on the second readout dimension. The magnetic field is orthogonal to the xz plane, thus the diffusion in the y direction is not affected by the presence of the magnetic field. The total amount of electrons is split between the views: one electron is assumed to be collected only by one view and for large numbers the charge sharing is granted. Several charge measurements, acquired at different gain values, have been used to evaluate the charge sharing between x and y views from the experimental data. In Fig. 11, the charge on the x view as a function of the total charge is reported and a linear fit extracts the ratio between the charge collected by this view and the total one. The value of 0.44 was included in the gain tuning factor since in the section 5 the charge was assumed to be collected entirely on one view. From the experimental data the charge sharing results Q y Q x = 1.13.

Validation of the tuning with double view readout in magnetic field
A validation of the simulation was obtained with a different dataset in order not to bias the results. In the validation procedure also the magnetic field effect was taken into account, to further prove the potential of the simulation with PARSIFAL. The accordance between the experimental and simulated values for this new dataset is shown in Fig. 12. Figure 12 highlights that the tuning factor evaluated with one view and without magnetic field grants an optimal matching between experimental and simulated data also on both views and with magnetic field. Table 1 shows the level of agreement defined as the ratio Figure 10. Comparison between the simulated (blue squares) and experimental (red circles) data for the four sentinel variables vs the incident angle: on top the cluster size (left) and charge (right), on bottom the spatial resolution for CC (left) and µ-TPC (right). The µ−TPC at 5 • is in between the most performing region and the less performing one, thus the non optimal agreement in correspondence of this incident angle is mostly due to this. The experimental data analysis can be found in [12,13]. No magnetic field is used in these measurement.

Conclusion
A parametrized simulation of a Micro Pattern Gas Detector has been successfully implemented: the quality of PARSIFAL results, exemplified by its validation against testbeam data, allows to simulate the full avalanche development and the signal formation in a triple-GEM, in a small amount of time (less than two seconds). The level of agreement with the results obtained on real data is about 1%. The parametrization starts from the detector description in ANSYS and GARFIELD++. The simulated signal is extracted from the ionization, amplification, diffusion and induction; it reproduces the real signal shape thanks to tuning factors taken into account in order to overcome the limits of the introduced approximations. Among the tuning factors, particularly interesting is the one on the gain: the value of 2.41 has to be applied to the GEM gain simulation from GARFIELD++ in order to reproduce the real data. This can extend the knowledge on the performance of a triple-GEM to a wide range of ionizing particles at different kinetic energies. The approach on which PARSIFAL is based can be extended to different configurations of geometry, high voltage settings and gas mixtures. Moreover, since some of the physics involved is shared among different technologies, the code can be extended to simulate also other gaseous detectors. In fact, with some modifications, it can be adapted directly to other MPGDs, e.g. MicroMegas and µ-RWELL. The code will be released to the MPGD community and on the GIT repository.