# Parallelization of Garfield++ and neBEM to simulate space-charge effects in
RPCs

Tanay Dey<EMAIL_ADDRESS>Purba Bhattacharya  Supratik Mukhopadhyay
Nayana Majumdar  Abhishek Seal  Subhasis Chattopadhyay

###### Abstract

Numerical simulation of avalanches, saturated avalanches, and streamers can
help us understand the detector physics of Resistive Plate Chambers (RPC). 3D
Monte Carlo simulation of an avalanche inside an RPC, the transition from
avalanche to saturated avalanche to streamer may help the search for the
optimum voltage and alternate gas mixtures. This task is dauntingly resource
hungry, especially when space-charge effects become important, which often
coincides with important regimes of operation of these devices. By modifying
the electric field inside the RPC dynamically, the space charge plays a
crucial role in determining the response of the detector. In this work, a
numerical model has been proposed to calculate the dynamic space-charge field
inside an RPC and the same has been implemented in the Garfield++ framework.
By modeling space charge as large number of line charges and using the
multithreading technique OpenMP to calculate electric field, drift line,
electron gain, and space-charge field, it has been possible to maintain time
consumption within reasonable limits. For this purpose, a new class,
pAvalancheMC has been introduced in Garfield++. The calculations have been
successfully verified with those from existing solvers and an example is
provided to show the performance of pAvalancheMC. Moreover, the details of the
transition of an avalanche into a saturated avalanche have been discussed. The
induced-charge distribution is calculated for a timing RPC and results are
verified with the experiment.

###### keywords:

Avalanche-induced secondary effects , Charge transport , Detector modelling
and simulations , Electric fields , Multiplication and induction , OpenMP ,
Parallel computation , Resistive-plate chambers , Space charge , Physics of
Gases, Plasmas and Electric Discharges.

††journal: Journal of Computer Physics Communications

[a]organization=Homi Bhabha National Institute,city=Mumbai, country=India
[b]organization=Variable Energy Cyclotron centre, city=Kolkata, country=India

[c]organisation=Department of Physics, School of Basic and Applied Sciences,
Adamas University, city=Kolkata, country=India

[d]organization=Saha Institute of Nuclear Physics, city=Kolkata, country=India
[e]Regent Education and Research Foundation, city=Kolkata, country=India

##  1 Introduction

Development of gaseous ionization detectors like Gas Electron Multipliers
(GEM)[1], Resistive Plate Chambers (RPC)[2, 3], Time Projection Chambers
(TPC)[4], etc., has been crucial for progress in various fields like high
energy physics, astronomy, and medical physics. Even though the geometry of an
RPC is simple in comparison to other gaseous detectors, the physics processes
associated with it (e.g., primary ionization, charge transport, electron
multiplication, and signal formation) are as complex as any other. They are
also conisdered reliable enough to be regularly used as timing, triggering,
and tracking devices in numerous experiments [5, 6, 7, 8]. A precise
simulation model is needed to model RPCs and explain the experimental results
adequately. A one-dimensional hydrodynamic model is beneficial for quick
results at small computational time [9, 10, 11]. This model works well when
the number of space charges is significant, as the collective behavior of
space charges meets the criteria of considering the system as a continuum
[12]. A detailed Monte Carlo simulation model, where particles are traced
individually, is the other option [13, 14, 15]. However, the large number of
charges in the avalanches slows down the particle model significantly; hence
it becomes computationally expensive.

Garfield++ [16, 17, 18] is a C++-based simulation tool to simulate ionization
detectors. It has its own geometry tools, which can be used to make simple
virtual detectors. However, for complex geometries, one can borrow field maps
and geometry files from external field solvers like COMSOL [19], neBEM [20,
21] etc. Several methods to generate avalanches inside gaseous detectors in
Garfield++, for example, microscopic tracking, Monte Carlo tracking, etc., are
based on particle models. These tracking methods are very detailed but slow in
a computational sense because all the algorithms run serially. One more
drawback of tracking methods in Garfield++ is the absence of the dynamic
space-charge effect, which is crucial when the number of space charges becomes
significantly large ($\geq 10^{6}$) so as to modify the applied field.

To address the above issues of Garfield++, we have written an algorithm and
introduced an extra class, pAvalancheMC, to generate avalanches with
multithreading techniques using OpenMP [22]. In this class, we have
incorporated the space charge effect using the line charge model discussed in
[23, 24]. Earlier a 2-D model was developed in [13], in which cylindrical co-
ordinates $r$, $z$ and $\phi$ were used under the assumption of rotational
symmetry of the avalanche around the $z$-axis. Consequently, the present model
can be considered as a 3D Cartesian extension of the model developed in [13,
14, 15]. In our model we have not considered rotational symmetry.

Another approach to parallelize Garfield++ using MPI can be found in [25, 26],
named pGARFIELD. However, pGARFIELD does not include the space-charge effect.

In Section 2, we describe the method of generation of uniform and uncorrelated
parallel random numbers using the TRandom3 class of ROOT [27], which is the
heart of the Monte Carlo simulations. Similarly, use of OpenMP and FastVol to
accelerate electric field computation are described in Section 3. Section 4
describes the steps or algorithm of generation of a Monte Carlo avalanche,
including the dynamic space-charge effect. In Section 5, we provide an example
of avalanche saturation with a gas mixture of 80% Ar and 20% $\ce{CO_{2}}$ in
an 2 mm gas gap RPC, generated through pAvalancheMC. In Section 7 we provide
another example of avalanche saturation and streamer formation with a gas
mixture of 97% $\ce{C_{2}H_{2}F_{4}}$, 2.5% $\ce{i-C_{4}H_{10}}$, and 0.5%
$\ce{SF_{6}}$. in an 2 mm gas-gap RPC, generated through pAvalancheMC. It may
be mentioned that the results obtained in this work are quite similar to those
in [13, 14, 15]. In all these computations similar physics processes such as
radial blow-up and contraction of the avalanches have been observed. Despite
the differences between the two models, it is encouraging to note that the
order of magnitude of space-charge effect is comparable. In Section 8, the
performance of the multithreaded code has been discussed. Section 9 compares
induced-charge distributions for three different voltages 1720 V, 1730 V and
1735 V with space-charge effect, calculated using neBEM and the proposed
version of Garfield++.

##  2 Uncorrelated and thread safe parallel Random number generation using
Trandom3 and OpenMP

We know that the outcome of any Monte Carlo simulation depends on the quality
of the generated random numbers. It is important that the generated random
numbers are uncorrelated, and the corresponding random engine has a
considerable period. TRandom3 is a class of ROOT [28], that generates uniform
random numbers. It is based on the Mersenne Twister (MT) algorithm and
provides a period of $2^{19937}-1$ and a 623-diminsionally equidistributed
uniform pseudorandom number generator up to 32-bit accuracy [29].

We have chosen TRandom3 as the random engine in our simulation code. The
random numbers generated serially using TRandom3 are uncorrelated. However, we
have seen that the former engine leaves some correlation between random
numbers generated in parallel in different threads (see Fig. 1) when all
threads use the same random engine and the same seed.

The random numbers corresponding to different threads, namely 0, 1, and 2,
generated from the same object of the random engine are strongly correlated
and mostly lying on three planes (see Fig. 1).

Figure 1: Correlation between random numbers generated in parallel using the
same seed by threads 0,1,2.

To solve this problem, objects of the TRandom3 class have been assigned for
three threads with three different seed values. The 3D correlation plot using
random numbers corresponding to threads 0, 1 and 2 (see Fig. 2) is uniformly
distributed over the volume.

Figure 2: Correlation between random numbers generated in parallel using
different seeds by threads 0,1,2.

The method outlined above is implemented in class pAvalancheMC. Using the
function SetNumberOfThreads(n), one can specify the number of threads used for
random number generation. The different objects of the TRandom3 class for
different threads can be assigned using the same function. Using a unique
thread ID created by OpenMP, all threads can generate random numbers from
their corresponding objects. An example of C++ code for the usage of the
function SetNumberOfThreads(n) is given below:

⬇

void SetNumberOfThreads(int nThread)

{

rnd.resize(nThread);

for(int ThreadId=0; ThreadId<nThread; ThreadId++)

rnd[ThreadId]= new TRandom3(0);

}

##  3 Accelerated electric field solution using OpenMP and FastVol

The Open Multi-Processing (OpenMP) is an Application Programming Interface
(API) that supports multi-platform shared memory multiprocessor programming in
C, C++ and FORTRAN on most processor architectures and operating systems [22].
It consists of a set of compiler directives, library routines and environment
variables that influence runtime behavior and uses a simple, scalable API for
developing parallel applications on platforms ranging from the standard
desktop to supercomputers. We have successfully implemented OpenMP for the
neBEM field solver [21]. The parallelization has been implemented in several
sub-functions of the toolkit, such as computation of the influence coefficient
matrix, matrix inversion and evaluation of the field and potential at desired
locations. These routines are computation intensive since there can be
thousands of elements where the charge densities need to be evaluated and the
influence of all these elements need to be taken into account. The task is
further complicated due to the use of repetitive structures, in order to
conform to the real geometry of a detector. This has proved to be very
important in improving the computational efficiency of the solver. We have
tested these implementations on up to 24 cores. The observed reduction in the
computational time is significant while the precision of the solution is
preserved.

Even after adopting OpenMP for solving for the charge distribution on all the
material interfaces of a given device, the time needed to calculate potentials
and fields for a complex device can become prohibitive. This is especially so
if the device (not usually true for an RPC) is composed of hundreds of
primitives, thousands of elements and several tens of repetitions. Reduction
of time taken to compute electrostatic properties becomes increasingly
important when complex processes such as avalanche, Monte-Carlo tracking and
micro-tracking are being modelled using Garfield++.

In order to model these phenomena within a reasonable time, we have
implemented the concept of using pre-computed values of potential and field at
large number of nodal points in a set of suitable volumes. These rectangular
volumes are chosen such that they can be repeated to represent any region of a
given device and simple trilinear interpolation is used to find the properties
at non-nodal points. The associated volume is named as the Fast Volume
(FastVol) which may be staggered, if necessary. In order to preserve the
accuracy despite the use of trilinear interpolation, it is natural that the
nodes should be chosen such that they are sparse in regions where the
potential and fields are changing slowly and closely packed where these
properties are changing fast.

##  4 Steps of an avalanche simulation with space-charge effect in Garfield++

This section will discuss techniques for implementing the space-charge effect
with multithread in the pAvalancheMC class of Garfield++. In further
discussions, we will use the short name $E_{a}$ as the applied field, $E_{s}$
as the space-charge field, and $E_{t}$ as the sum of applied and space-charge
fields.

The necessary steps to simulate an avalanche are given below:

  1. (i)

Selection of primaries to start avalanche according to their generation time.

  2. (ii)

Parallel calculations of drift line of avalanche charges.

  3. (iii)

Calculation of gain in each step.

  4. (iv)

Generation of cylindrical grid of space charge region.

  5. (v)

Calculation of dynamic space-charge field.

###  4.1 Selection of primaries to start avalanche according to their
generation time

When a charged particle passes through a gas detector, it can interact with
gas molecules and leave some primary ionizing particles along its path. A
track of a charged particle through the gas gap and corresponding primary
ionizations can be simulated using the class TrackHeed [30] of Garfield++. The
locations and time of generation of primary clusters are then stored in an
array. The primaries (here electrons) can be generated at different instants
of time along the track. Hence, the starting time of the avalanche generation
from different primaries can be different. Therefore, if we arrange primaries
in ascending order based on the generation time, then depending on the
simulation step $\delta t$, they will start the avalanche process. For
example, if $t$ is the generation of time of any primary electron, it will
start an avalanche when $n\delta t\geq t$, where $n$ is the number of steps
completed in the simulation.

###  4.2 Parallel calculations of drift line of avalanche charges

The AvalancheMC class of Garfield++ can simulate the drift path of charges
serially. Hence, it becomes slow when the number of charges is large. To solve
this issue, we assign the secondary particles to the available threads using
OpenMP and propagate them in the gaseous medium simultaneously at every step
of time to calculate further drift points. The diffusion part of pAvalancheMC
is kept the same as in-class AvalancheMC, where the thermal diffusion is
considered Gaussian and anisotropic when an electric field is present [13].
Therefore, the diffusion will be in two directions, longitudinal and
transverse. For total electric field $E_{t}$, the drift velocity
($v_{D}(E_{t})$), longitudinal ($D_{L}(E_{t})$) and transverse diffusion
($D_{T}(E_{t})$) constants have been calculated by using MAGBOLTZ [31, 32].
The variance of longitudinal Gaussian distribution is
$\sigma_{L}=D_{L}\sqrt{v_{D}\delta t}$ and that of transverse diffusion is
$\sigma_{T}=D_{T}\sqrt{v_{D}\delta t}$.

###  4.3 Calculation of Gain in each step

The calculation of gain in the pAvalancheMC class is the same as in the
AvalancheMC class. It is based on the Yule-Furry model and includes the effect
of electron attachment [33, 34]. At each step, the probability of ionisation
and attachment is calculated using the Townsend ionisation ($\alpha$) and
attachment ($\eta$) coefficients, where $\alpha$ and $\eta$ are computed using
MAGBOLTZ [31, 32].

###  4.4 Calculation of dynamic space-charge field

The space-charge effect is turned on when the number of electrons crosses a
threshold. This threshold is determined by the user and, for the computations
presented in this paper, the value has been set to $10^{4}$. The generation of
grid elements and finding locations of the space charges in the grids have
been divided into four steps. The steps are as follows:

Step 1: Let us consider at a particular time step of a growing avalanche, the
space charges are distributed as shown in Fig. 3a. To calculate the field for
this charge distribution, the $z$-directional space is divided into $S$
elements of size $\delta z$, as shown in Fig. 3a.

Step 2: For a particular $z$ position, the avalanche region is divided into R
co-centric rings each having a thickness $\delta r$ as shown in Fig. 3b. For
all $z$, the thickness $\delta z$ is considered as small as 0.001 cm.

Step 3: In the third step, each ring can be divided into $L$ curved segments
(see Fig. 3c). Hence, each curved segment will lie between angles of $\phi$
and $\phi+\delta\phi$, where $\phi$ is the azimuthal angle subtended to the
center of the circle or ring (see Fig. 3d). Therefore, the $\phi$ directional
space will also be divided into $L$ segments of size $\delta\phi$, so that
$L\delta\phi=360$° (for more details see [23]).

Step 4: In this calculation, $\delta\phi$ is considered very small ($\approx
1$°) so that each curved segment can be considered as a straight line of
length $r\delta\phi$ (see Fig. 3d).

(a)

(b)

(c)

(d)

Figure 3: Steps to generate cylindrical grid, (a) dividing z space into
elements of size $\delta z$, (b) radial space divided into several co-centric
rings of increasing radius r and thickness $\delta r$, (c) divide rings into
several small curved segments, (d) decrease size of curved elements and
approximation by straight line.

Therefore, the entire avalanche volume is segmented into small volumes
(voxel). The charges within a specific voxel is represented by a line charge
with constant charge density situated within that volume. The total number of
voxels $T=S\times R\times L$.

The space-charge field at any voxel is the sum of the electric fields due to
other voxels containing non-zero charges. The space-charge field is considered
to be uniform within a voxel. The total field at any location ($r,\phi,z$) is
the sum of the applied field and space-charge field. Since voxels are
considered as line segments, the line charge formula has been used to
calculate the space-charge field, as discussed in [23, 24].

###  4.5 Implementation of positive and negative ions

Positive ions are considered to be generated at the ionization locations. The
attachment positions of the electrons within the gas gap are considered as the
positions of the singly charged negative ions. It may be noted that, since the
electron drift velocity is much higher than that of the ions, the ions are
approximated to be static in this paper. However, implementation of ion drift
is possible, if necessary, by providing additional information about ion
mobility.

##  5 Instance of a simulated avalanche inside an RPC with space-charge effect
in Ar and $\ce{CO_{2}}$ gas-mixture

###  5.1 Avalanche in absence of the negative ion effect

This section will discuss an example of an avalanche simulated in an RPC using
pAvalancheMC with and without the space-charge effect. A gas mixture of 80% Ar
and 20% $\ce{CO_{2}}$ is used. The electrodes of the RPC are made of a 2 mm
thick Bakelite with an area of $30\times 30~{}\text{cm}^{2}$. The gas gap is
also taken as 2 mm. A constant electric field along $+z$ ($E_{a|z}$) of 23.5
kV/cm has been applied across the gas gap with the help of ComponentConstant
class of Garfield++. The time-step of the simulation is taken as 20 ps. Muon
tracks and primary ionisation inside the gas gap are generated using HEED
[30]. It is noted that the contribution of the negative ions are not
considered while performing this avalanche.

In Fig. 4a, the electron gain at each step of time with and without the space-
charge effect is compared. The blue and pink curves specify the electron gain
without and with the space-charge effect, respectively. For avalanche times
below 10 ns, the two curves almost overlap. After 10 ns, the avalanche charge
without space-charge effect continuously grows, while the avalanche charge
with space-charge effect (pink curve) grows slowly, reaches a peak at around
15.46 ns, and then starts showing a saturation effect after 17.96 ns, at an
electron count of approximately $10^{8}$.

(a)

(b)

Figure 4: (a) Variation of number of electrons and signal with time, (b)
variation of total electric field $E_{t}^{max}$ and $E_{t}^{min}$ with time in
absence of negative ions.

(a)

(b)

(c)

(d)

Figure 5: (a) Variation of number of electrons with time in presence of
negative ions, (b) Variation of total electric field $E_{t}^{max}$ and
$E_{t}^{min}$ with time in presence of negative ions, (c) variation of ratio
between number of positive ions ($n^{+}$) and negative ions ($n^{-}$) with
time for a gas mixture 80% of Ar and 20% of CO2, (d) variation of ionisation
($\alpha$) and attachment ($\eta$) coefficients with electric field for a gas
mixture 80% of Ar and 20% of CO2.

###  5.2 Avalanche in presence of the negative ion effect

The configurations of the simulation are the same as those in subsection 5.1.
The variation in the number of electrons with time in the presence of the
negative ions is shown in Fig. 5a, and the variation in the maximum of the
total electric field $E_{t}^{max}$ is shown in Fig. 5b. A result from the
comparison of the time variation in the number of electrons with the space
charge effect of Figs. 4a (pink curve) and 5a indicates that the gain has been
reduced by approximately a factor of four.

As in Fig. 4a, saturation in the variation of the number of electrons with
time can be observed in Fig. 5a. However, the knee point or starting point of
the saturation region for the pink curve in Figs. 4a and 5a is at
approximately 17.6 ns and 21 ns, respectively. The number of electrons at the
knee point of Fig. 4a is approximately 8.7 times the number of electrons at
the knee point of Fig. 5a. Moreover, it is found that the saturation curve of
Fig. 5a is steeper than the pink curve of Fig. 4a. This is because, in the
presence of negative ions, the maximum electric field in the saturation region
lies between 17-20 kV/cm and minimum electric field is close to zero, as shown
in Fig. 5b, whereas in the absence of negative ions, it was beyond 30 kV/cm as
shown in Fig. 4b. Hence, the gain factor is reduced.

The variation in the ratio of the number of positive ions ($n^{+}$) to
negative ions ($n^{-}$) is illustrated in Fig. 5c. From the previous figure,
it becomes evident that prior to 10 ns, the ratio fluctuates between 20 and
52. After 10 ns, the ratio gradually decreases, and it stabilizes at
approximately 11 after 18 ns. This substantial imbalance between negative and
positive ions can be attributed to the difference of magnitude between
attachment ($\eta$) and ionization ($\alpha$) coefficients, as depicted in
Fig. 5d. It is apparent from Fig. 5d that at around 12 kV/cm, $\alpha$ and
$\eta$ attain equal values. As the electric field increases, the difference
between $\alpha$ and $\eta$ increases rapidly. Therefore, below 12 kV/cm,
attachment dominates over ionization. As previously mentioned, in the
saturation region, the maximum electric field in certain parts of the
avalanche region reaches around 20 kV/cm. At this point, the difference
between $\alpha$ and $\eta$ is approximately 55 1/cm, indicating that the
ionisation process prevails. Consequently, in the saturation region, the
number of $n^{+}$ and $n^{-}$ becomes balanced, resulting in a constant ratio.

##  6 Visualisation of transport of electrons under space charge effect

In order to understand the transport of charged particles through the gas
volume, the distribution of $z$-positions vs. radial positions of the
electrons at different instants of avalanche growth have been analyzed ( Fig.
6). The following discussion contains the result of the analysis.

It is evident from subsection 5.2 that the Ar-CO2 gas mixture does not exhibit
a strong attachment coefficient, resulting in a higher concentration of
positive ions compared to negative ions. Therefore, for the sake of
expeditious computation and simplicity, we neglect the effect of negative ions
and focus solely on the transport of electrons.

We can divide the pink curve of Fig. 4a into three regions: a) before 10 ns,
where space-charge effects are negligible, b) from 10 ns to 15.46 ns, when
gain increases rapidly, and c) after 15.46 ns, where the gain drops and
maintains a saturated value.

####  6.0.1  From 0 ns to 10 ns, overlap region of gain of Fig. 4a

  1. i.

At time 0 ns (see Fig. 6a) there are in total 58 electrons distributed in 7
clusters. A small number of primary ions are generated along with the primary
electrons. They are neglected since they do not participate strongly in the
rest of the evolution that is being followed in the present study. It is noted
that $z$=0.1 is the anode plane and $z$=-0.1 is the cathode plane. The color
bar represents the number of electrons in a ring of radius r and thickness
$\delta r$.

  2. ii.

In Fig. 6b, it is shown that at time 4.5 ns electrons have already started
participating in the ionization process, and a total of 15,760 electrons and
16,218 ions are generated. In this figure, two clusters near the cathode plane
at $z$=-0.1 merge and form a big cluster. Therefore, the total number of
electron clusters at this stage is 4. The maximum charge density is in the
middle of the distribution (at radius $r=0$). The electron distribution is
roughly symmetric.

  3. iii.

In Fig. 6c at time 7.46 ns, there are only two electron clusters. The total
number of electrons and ions is 612,558 and 666,083, respectively. Also, the
charged center of the merged distribution has reached the anode, and the
electron density contour is smooth and symmetric.

(a)

(b)

(c)

(d)

Figure 6: Location of electrons in $z-r$ plane at time (a) 0. ns, (b) 4.5 ns,
(c) 7.46 ns and (d) 10.46 ns. The inlet shows a zoom into the highlighted
electron cluster.

(a)

(b)

(c)

(d)

Figure 7: (a) $z$-component of space-charge field ($E_{s|z}$) at 10.46 ns, (b)
Radial component of space-charge field ($E_{s|r}$) at 10.46 ns, (c)
$\phi$-component of space-charge field ($E_{s|\phi}$) at 10.46 ns, (d) Shape
of the electron cloud and total electric field ( $E_{t}$) at different voxels
at time 10.46 ns.

####  6.0.2  From 10 ns to 15.46 ns, till peak region of Fig. 4a

  1. i.

At 10.46 ns, the blue and pink curves of Fig. 4a start to separate, which
signifies the relative importance of the space-charge effect. There are
5,897,046 electrons and 7,768,542 ions (see Fig. 6d). At this stage, the last
two clusters also merge and form a single big cluster. A very close look at
Fig. 6d shows that the charge center of the last cluster from the anode has
not remained at $r=0$; rather, it shifted to the right. This is because
$z$-directional space-charge field ($E_{s|z}$) at $z=-0.04$ and $r=0$ is more
negative (opposite to the applied field) in comparison to other regions (see
Fig. 7a). Therefore, the total field ($E_{t}$) and the ionization probability
at the center are smaller compared to the other regions. The radial component
of the space-charge field ( $E_{s|r}$) plays an important role in spreading
the avalanche along the transverse direction. The radial field at 10.46 ns at
different $z$ and radius $r$ has been shown in Fig. 7b. The radial field
$E_{s|r}$ of ions is positive, and of electrons is negative. A positive sign
indicates an outward radial field from the center, and a negative sign
signifies an inward radial field to the center. It is noted that near the
region $z$=-0.04 and $r$=0.005 (see Fig. 7b) the radial field is negative,
meaning the number of electrons dominates the number of ions at those places.
The phi ($E_{s|\phi}$) component of the field measures the axial symmetry of
the avalanche charge distribution. The charge distribution will be axially
symmetric when the value of $E_{s|\phi}$ is close to zero. A nonzero value of
$E_{s|\phi}$ signifies axial asymmetry. In Fig. 7c the variation of
$E_{s|\phi}$ with $z$ position and radius $r$ has been shown. At 10.46 ns,
$E_{s|\phi}$ is not as strong as $E_{s|r}$ and $E_{s|z}$ (maximum $E_{s|\phi}$
is only 26% of $E_{s|z}$). However, $E_{s|\phi}$ can still act as a
perturbation to the applied field.

The shape of the electron cloud corresponding to time 10.46 ns is shown in
Fig. 7d, and the magnitude of the total field ($E_{t}$) in V/cm at different
voxels are shown in the color bar of the same figure. The maximum positive and
negative deviation of the magnitude of the total field from the initially
applied field are 5.62% and 8.0%, respectively.

(a)

(b)

(c)

Figure 8: (a) Location of electrons in $z-r$ plane at time 12.46 ns, (b)
Location of electrons in $z-r$ plane at time 12.96 ns, (c) Location of
electrons in $z-r$ plane at time 13.46 ns.

(a)

(b)

(c)

Figure 9: (a) $z$-component of space-charge field at 12.46 ns, (b) radial
component of space-charge field at 12.46 ns, (c) $\phi$-component of space-
charge field at 12.46 ns.

  2. ii.

From the collective view of Figs. 8a, 8b and 8c it can be said that from 12.46
ns to 13.46 ns the charge center again gradually move back to the center
$r=0$. The $z$-directional space-charge field ($E_{s|z}$) at 12.46 ns near
$z=-0.08$ and $r=0$ becomes positive and higher than in other regions (see
Figure Fig. 9a). Therefore, the ionization rate becomes high at that location.
This is the reason the charge center again moves towards $r=0$. In Fig. 9b,
the radial field at a different location at time 12.46 ns is shown. At time
12.46 ns, the radial field increases significantly (absolute maximum along
positive and negative directions around 5 kV/cm and 3 kV/cm, respectively) at
some places, which are reflected in the transverse or radial spread
($r_{max}$=0.044 cm) of the avalanche (see Fig. 8a). The variation of
$E_{s|\phi}$ at different $z$ location and radius has been shown in Fig. 9c.
The maximum magnitude of $E_{s|\phi}$ at time 12.46 ns is 0.8 kV/cm, which is
double the value of $E_{s|\phi}$ at 10.46 ns.

It is clear from Fig. 10a (color bar represents the number of electrons at a
given voxel) that the range of the electron cloud along the $z$ direction is
-0.1 cm to -0.04 cm. As shown in Fig. 10b (color bar represents the number of
ions at a given voxel), the range of the ion cloud is -0.1 cm to 0.05 cm.
Therefore, along the $z$-axis, ions are on both sides of the electron cloud.
Hence, at the tip and tail of the electron cloud, ions attract the nearest
electrons toward them. As a result, the charge distribution becomes asymmetric
along the $z$-axis.

(a)

(b)

Figure 10: (a) Shape of the electron cloud and number of electrons at
different voxels (color bar) at time 12.46 ns, (b) Shape of the ion cloud and
number of ions at different voxels (color bar) at time 12.46 ns.

In Figs. 11a, 11b and 11c, the shape of the electron cloud and magnitude of
total field at different places of the cloud and at different times 12.46
ns,12.96 ns and 13.46 ns have been shown. It is found from the former figures
that though the electron density at different places inside the cloud is
changing with time, the overall shape of the distribution is not changed much.
In this period of transition, the maximum positive and negative deviation of
the total electric field from the applied field is 36.02% and 27.5%,
respectively.

(a)

(b)

(c)

Figure 11: (a) Shape of the electron cloud and electric field magnitude
($E_{t}$ in V/cm) at different voxels (color bar) at time 12.46 ns, (b) Shape
of the electron cloud and electric field magnitude ($E_{t}$ in V/cm) at
different voxels (color bar) at time 12.96 ns, (d) Shape of the electron cloud
and electric field magnitude ($E_{t}$ in V/cm) at different voxels (color bar)
at time 13.46 ns.

####  6.0.3  After 15.46 ns, peak to saturated region of Figure Fig. 4a.

  1. i.

At 15.46 ns, the electron gain reaches its maximum value (see Fig. 4a). All
the electron clusters are merged together (see Fig. 12a). Since already many
electrons have left the gas gap, the total number of ions (253,862,910) is far
greater than the total number of electrons (160,516,288). From Fig. 13a, it is
confirmed that due to negative $E_{s|z}$, the net z-directional field has
reduced by a large amount at the tip and the center ($r=0$) of the avalanche.
Also, due to the significant difference in the electron and ion numbers, the
radial field of the ions dominates; hence, only a positive radial field exists
everywhere. As a result, the remaining electrons feel a strong radial force
towards the center, as shown in Fig. 13b. Therefore, the transverse spread is
stopped, and all the electrons tend to converge toward the center. The
component $E_{s|\phi}$ attains a maximum value of 2 kV/cm approximately (see
Fig. 13c) which indicates significant axial asymmetry of the space-charge
distribution. In Fig. 12b, the shape of the electron cloud and electric field
at different voxels has been shown. The maximum positive and negative
deviation of a total electric field from the applied field at 15.46 ns is
84.4% and 50.5%, respectively.

(a)

(b)

(c)

(d)

Figure 12: (a) Location of electrons in $z-r$ plane at time 15.46 ns, (b)
Shape of the electron cloud and electric field magnitude ($E_{t}$ in V/cm) at
different voxels (color bar) at time 15.46 ns, (c) Location of electrons in
$z-r$ plane at time 16.96 ns, (d) Shape of the electron cloud and electric
field magnitude ($E_{t}$ in V/cm) at different voxels (color bar) at time
16.96 ns.

(a)

(b)

(c)

Figure 13: (a) $z$-component of space-charge field at 15.46 ns, (b) radial
component of space-charge field at 15.46 ns, (c) $\phi$-component of space-
charge field at 15.46 ns.

(a)

(b)

(c)

(d)

Figure 14: (a) Location of electrons in $z-r$ plane at time 18.46 ns, (b)
Shape of the electron cloud and electric field magnitude ($E_{t}$ in V/cm)
(color bar) at different voxels at time 18.46 ns, (c) Location of electrons in
$z-r$ plane at time 23.46 ns, (d) Shape of the electron cloud and electric
field magnitude ($E_{t}$ in V/cm) (color bar) at different voxels at time
23.46 ns.

  2. ii.

At time 16.96 ns, there are 126,594,692 electrons and 282,514,583 ions left
inside the gas gap (see Fig. 12c). This is the time when the space-charge
field attains its maximum value. The maximum positive and negative deviation
of the magnitude of the total field (space-charge + applied field) from the
initially applied field are 132.8% and 86.2% respectively. The shape of the
electron cloud becomes conical (see Fig. 12d), where at the tip, the field is
reduced drastically, and at the tail, it becomes more than double the initial
value of the field.

  3. iii.

At 18.46 ns, the gain curve with space-charge effect of Fig. 4a contains a
second knee point after which the saturation region starts. At this stage,
there are 100,950,517 electrons and 288,968,288 ions. The density contours are
not smooth near the center (see Fig. 14a). There are three possibilities for
the reduction of the number of outermost electrons (i) attachment, (ii)
leaving the gas gap, and (iii) at some places, the inward radial field is much
stronger so that it attracts the electrons towards the center. Therefore, the
conical shape gradually modifies into cylindrical (see Fig. 14b). The maximum
positive and negative deviation of the magnitude of the total field ($E_{t}$)
from the initially applied field are 85.8% and 99.9%.

  4. iv.

At 23.46 ns, the gain of the avalanche is saturated. In most of the voxels,
the total field is less than 1 kV/cm (see Fig. 14d). Therefore, at those
places, the drift velocity of the electrons is less than 0.008 cm/ns (see Fig.
15a). Also, in these low-field regions, the attachment coefficient ($\eta$)
dominates over the first Townsend coefficient (alpha) ($\alpha$). Hence the
value of $\alpha-\eta$ becomes negative (see Fig. 15b). At higher field
regions, there are some electrons capable of producing ionisations. However,
the number of electrons created in ionisations is compensated by the number of
electrons attached; that’s why the electron gain saturates, as shown in Fig.
4a. The maximum positive and negative deviation of the magnitude of the total
field ($E_{t}$) from the initially applied field are 41.5% and 99.8%.

(a)

(b)

Figure 15: (a) Drift velocity as a function of electric field, (b) Effective
Townsend coefficient as a function of electric field.

##  7  Avalanche simulations inside the RPC with $\ce{C_{2}H_{2}F_{4}}$,
$\ce{i-C_{4}H_{10}}$, and $\ce{SF_{6}}$ gas mixture

###  7.1 Avalanche in absence of negative ion

(a)

(b)

(c)

Figure 16: Variation of electron gain with time steps of the avalanche for
electric fields of 49.85 kV/cm, and 50 kV/cm, (a) without space-charge effect,
and (b) with space-charge effect, (c) variation of maximum and minimum values
of total electric field with time.

In this particular section, we present another example of the avalanches
occurring inside an RPC with a 2 mm gas-gap. For our investigation, we used a
gas-mixture containing 97% $\ce{C_{2}H_{2}F_{4}}$, 2.5% $\ce{i-C_{4}H_{10}}$,
and 0.5% $\ce{SF_{6}}$. A muon track and primary ionisations inside the gas
gap are simulated using HEED. We performed avalanche simulations, first
without taking into account the space-charge effect, and then with the space-
charge effect. Two different electric fields 49.85 kV/cm, and 50 kV/cm along
+z direction are considerd for these simulations. Our findings show the
variation in the number of electrons with time, as shown in Figs. 16a (without
space-charge effect) and 16b (with space-charge effect). From Fig. 16a it is
found that in absence of the space-charge effect the number of electrons grow
continuously with time. Similar to the case discussed in Section 5, the number
of electrons in the presence of the space-charge effect initially grows with
time, reaches a peak, and then decreases, as shown in Fig. 16b. In the same
figure, the green curve exhibits a brief saturation region after 16.5 ns,
followed by a steep increase at 17.5 ns. In the case of an applied field of 50
kV/cm (red curve in Fig. 16b), a saturation region is not observed. However,
in this case, the number of electrons reaches a peak, starts decreasing, and
suddenly increases without having reached a saturation period. The variation
of the maximum and minimum of the total field $E_{t}^{max/min}$ with time for
both applied fields is shown in Fig. 16c.

###  7.2 Avalanche in presence of negative ion

(a)

(b)

(c)

(d)

Figure 17: (a) Variation of number of electrons with time in presence of
negative ions. Variation of total electric field, (b) $E_{t}^{max}$ and (c)
$E_{t}^{min}$ with time. (d) Variation of ratio between number of positive
ions ($n^{+}$) and negative ions ($n^{-}$) with time for a gas mixture 97% of
$\ce{C_{2}H_{2}F_{4}}$, 2.5% of $\ce{i-C_{4}H_{10}}$, and 0.5% of
$\ce{SF_{6}}$.

Avalanche processes take place within the RPC with the space charge effect,
under three distinct applied electric fields: 49.8 kV/cm, 51 kV/cm, and 51.8
kV/cm. The geometric configuration and gas mixture remain consistent with
those described in subsection 7.1. The variation in the number of electrons
over time for the three different applied electric fields is depicted in Fig.
17a. As the applied electric field increases, the peaks of the curves in Fig.
17a shift towards smaller time values (from 13.5 ns at 49.8 kV/cm to 12.5 ns
at 51 kV/cm). At 51.8 kV/cm, the number of electrons exhibits a significant
increase beyond 10 ns. However, unlike that observed in subsection 7.1, no
saturation region or abrupt steep rise in electron count is observed here. The
fluctuations in $E^{t}_{max}$ and $E^{t}_{min}$ corresponding to the three
applied voltages are illustrated in Figs. 17b and 17c, respectively. With the
growth of electrons, positive and negative ions, $E^{t}_{max}$ surpasses the
applied field, while $E^{t}_{min}$ falls below it. Following the peak (refer
to the red and black lines in Fig. 17a), the electron count begins to decline.
Consequently, for the applied fields of 49.85 kV/cm and 51 kV/cm (see Figs.
17b and 17c), $E^{t}_{max}$ decreases and eventually becomes nearly equivalent
to $E^{t}_{min}$. Nonetheless, for the scenario at 51.8 kV/cm, the electron
count never decreases, maintaining an upward trend in $E^{t}_{max}$ and
$E^{t}_{min}$, exhibiting an opposite trend as shown in Figs. 17b and 17c.

The time variation of the ratio between the number of positive ions ($n^{+}$)
and negative ions ($n^{-}$) is depicted in Fig. 17d. The presence of SF6
enhances the electronegativity of the gas mixture used here, resulting in a
high attachment coefficient. Consequently, the probability of negative ion
generation is high. From Fig. 17d, it can be concluded that for an applied
field of 49.8 kV/cm at the initial stage of the avalanche, there are more
positive ions than negative ions. This leads to a maximum ratio between
$n^{+}$ and $n^{-}$ exceeding 22. As time progresses, this ratio gradually
decreases and approaches 1 at the end of the avalanche. The same trend is
observed for the case of an applied field of 51 kV/cm. However, when the
applied field is 51.8 kV/cm, the avalanche diverges, and at 12 ns, the final
ratio is approximately 6.

##  8 Timing performance of pAvalancheMC using OpenMp

We have used the OpenMP multithreading technique to parallelize the class
pAvalancheMC of Garfield++. In the following, a timing RPC of area 30 cm
$\times$ 30 cm, 0.3 mm gas gap, and 2 mm thick bakelite electrode has been
studied in order to assess the timing performance of the parallelized code. A
voltage of $\pm$1720 V has been applied on the graphite surface. As a result,
the average electric field inside the gas gap is 43 kV/cm. A gas mixture of
$\ce{C_{2}H_{2}F_{4}}$ (85%), $\ce{i-C_{4}H_{1}0}$ (5%), $\ce{SF_{6}}$ (10%)
is used.

A single primary electron is placed near the negative electrode to generate
avalanches. A set of $10^{4}$ avalanches has been generated repeatedly by
varying the number of thread values ($N$) in OpenMP. Let $T_{N}$ represent the
time taken to generate $10^{4}$ events with a total $N$ number of threads. The
variation in total time ($T_{N}$) with and without space-charge effect as a
function of $N$ threads has been shown in Figs. 18a and 18c, respectively. It
is clear from the preliminary figures that in both cases, with and without the
space-charge effect, $T_{N}$ decreases with the increase in the total number
of threads $N$. The speed-up performance of the code can be calculated by
taking the ratios between sequential time ($T_{1}$) with $T_{N}$. Therefore,
speed-up factor $v_{p}$ can be written as:

$v_{p}=\frac{T_{1}}{T_{N}}.$ (1)

The speed-up performance $v_{p}$ as a function of number of threads ($N$) has
been shown in Figs. 18b and 18d with and without space-charge effect,
respectively. The $v_{p}$ increases with the number of threads ($N$). The data
points of both Figs. 18b and 18d is fitted with a non linear function:

$f(N)=p_{0}-p_{1}\,e^{-(p_{2}\,N^{p_{3}})},\,(N\geq 1)$ (2)

where the fit parameters are shown in the same figures. From equation 2, it
can be said that for large $N$, the contribution of the second term is
significantly less; hence the parameter $p_{0}$ signifies the maximum speed-up
which can be achieved. Therefore, based on parameter $p_{0}$, the maximum
speed up of 5.46 and 7.2 have been observed without and with space-charge
effect, respectively.

The saturation in speed-up can be explained by overhead issues. The overhead
arises during the communication and synchronization between threads. More
details about the overhead in OpenMP can be found in [35, 36]. It is also
noted that the time plotted in Figs. 18a and 18c is the total time to complete
all avalanches, which is the combination of time of parallel and non-parallel
computations. Therefore, $T_{N}$ is not the time purely taken by $N$ threads;
instead, it is the sum of time by $N$ threads for parallel computations and a
single thread for non-parallel or sequential computations.

In Figs. 19a and 19b, the distributions of electron gain of $10^{4}$
avalanches have been compared with and without space-charge effect,
respectively. As a result, it is found that the mean gain is reduced by order
of 10 when the space-charge effect is considered, and the shape of the
distribution is also modified significantly. Again, it is expected that the
average electron gain should not vary much with the number of threads ($N$).
In the same Figs. 19a and 19b, the gain calculations have been carried out
with single thread and a total of 18 threads, where it can be seen that they
are matched within some statistical uncertainties.

(a)

(b)

(c)

(d)

Figure 18:  Variation of execution time to complete $10^{4}$ avalanches with
number of threads (a) without space-charge effect and (c) with space-charge
effect; speed up performance with number of threads (b) without space-charge
effect, (d) with space-charge effect.

(a)

(b)

Figure 19: Distribution of electron gain of $10^{4}$ avalanches (a) without
space-charge effect (b) with space-charge effect.

##  9 Induced-charge distribution of avalanches

Charge can be induced on the readout electrodes due to the movement of the
electrons and ions. The induced charge $Q^{ind}$ due to $N_{av}$ number of $q$
point charge moving inside the RPC can be calculated by using Ramo’s equation
as follows [37]:

$Q^{ind}=\int_{0}^{t}dt\sum_{n=0}^{N_{av}}\,\,q\,({\phi}_{w}^{n}({r_{f}}(t))-{\phi}_{w}^{n}({r_{i}}(t))),$
(3)

where ${\phi}_{w}^{n}({r_{f}}(t))$ and $\,{\phi}_{w}^{n}({r_{i}}(t)$ are the
weighting potential at initial ($r_{i}$) and final position ($r_{f}$) of the
step calculated by using neBEM.

The induced-charge distributions of a timing RPC for three different applied
voltages, 1720 V, 1730 V, and 1735 V, are shown in Fig. 20. The RPC’s geometry
and other configurations remain the same as in Section 8. The three
distributions of Fig. 20 have been fitted with the Polya function given as
follows [38, 39] :

$f(Q^{ind})=a(\frac{Q^{ind}\,b}{c})^{b-1}\,e^{-\frac{b}{c}Q^{ind}},$ (4)

where parameter $a$ is the scaling factor, $b$ is a free parameter determining
the shape of the distribution, and $c$ is the mean charge. In all fitting
processes with Polya function, the left inefficiency peak of the charge
distributions of the Fig. 20 has not been considered. The fitted values of
$a$, $b$, and $c$ for different voltages given in the Table 1. From the values
of fit parameters of Fig. 20, it can be said that the value of the parameter
$c$ or mean charge increases with the increase in applied voltage. Moreover,
the value of $b$ also shifted towards a higher value with the voltage
increment, determining the broadness of the charge spectrum. From experimental
results such as [40, 13], the shifting of the mean and broadening of the shape
of induced-charge distribution with increasing applied voltage are also
observed.

Voltage (V) | a | b | c (fc)
---|---|---|---
1720 | $\begin{array}[]{c}139.7\\\ \pm\\\ 68.7\end{array}$ | $\begin{array}[]{c}3.8\\\ \pm\\\ 0.5\end{array}$ | $\begin{array}[]{c}0.34\\\ \pm\\\ 0.01\end{array}$
1730 | $\begin{array}[]{c}19.8\\\ \pm\\\ 7.9\end{array}$ | $\begin{array}[]{c}4.9\\\ \pm\\\ 0.3\end{array}$ | $\begin{array}[]{c}0.8\\\ \pm\\\ 0.01\end{array}$
1735 | $\begin{array}[]{c}0.06\\\ \pm\\\ 0.07\end{array}$ | $\begin{array}[]{c}8.2\\\ \pm\\\ 0.7\end{array}$ | $\begin{array}[]{c}1.3\\\ \pm\\\ 0.02\end{array}$
Table 1: Fit parameters of induced-charge distribution of Fig. 20. Figure 20:
Comparison between induced-charge distributions corresponding to three applied
voltages 1720 V, 1730 V and 1735 V.

##  10 Summary

In this paper, we have discussed the method of parallel generation of
uncorrelated and uniform random numbers. The problem of correlation between
the random numbers generated by multiple threads has been successfully solved
using different seed values and objects of TRandom3.

The steps of implementing the dynamic space-charge field inside Garfield++
have been discussed. In this context, the detailed modeling of a saturated
avalanche has been examined to understand the space-charge effect. The
changing of the dynamic space-charge field and the evolution of the 3D shape
of the avalanche with time have also been shown. As a result, we have found
that at the time of saturation, the rate of attachment is dominant in the low-
field regions, while in the high-field regions the rate of ionization is
dominant. Hence, overall, the number of charges remains the same, which brings
this saturation.

The influence of negative ions has also been carefully examined. In both gas
mixtures containing (a) Ar and CO2 and (b) C2H2F4, i-C4H${}_{1}0$, and SF6,
electron gain is significantly reduced due to the presence of negative ions.
The simulation shows instability in the absence of negative ions for the gas
mixture containing SF6, which is resolved when they are included in the
simulation.

Parallelization has been acheived by implementing multithreading techniques in
Garfield++. With OpenMP implemented in neBEM and the avalanche generation
process has made it possible to achieve approximately 5.4 times and 7.2 times
speed up without and with the space-charge effect, respectively.

Finally, the induced-charge distribution for a timing RPC has been calculated
for three different voltages. The shape of the spectra of the induced-charge
distribution predicted by the present simulations are similar to those
observed in other experimental results [40].

In future, we plan to implement a photon transport model in our code to study
more effectively the avalanche-to-streamer transition. Moreover, the study of
other gaseous detectors will also be possible with the model described in the
present paper.

## Acknowledgement

All authors are grateful to the INO collaboration, the HEP experiment division
of VECC, Applied Nuclear Physics Division of SINP and Physics Department,
School of Basic and Applied Sciences, Adamas University for providing the
resources and help. The authors thank the reviewer for constructive and
meticulous review that helped improvement in both content and form of the
manuscript.

## References

  * [1] F. Sauli, The gas electron multiplier (GEM): Operating principles and applications, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 805 (2016) 2–24, special Issue in memory of Glenn F. Knoll.  doi:https://doi.org/10.1016/j.nima.2015.07.060.
  * [2] R. Santonico, R. Cardarelli, Development of resistive plate counters, Nuclear Instruments and Methods in Physics Research 187 (2) (1981) 377–380.  doi:https://doi.org/10.1016/0029-554X(81)90363-3.
  * [3] R. Cardarelli and R. Santonico and A.Di Biagio and A. Lucci, Progress in resistive plate counters, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 263 (1) (1988) 20 – 25.  doi:https://doi.org/10.1016/0168-9002(88)91011-X.
  * [4] H. J. Hilke, Time projection chambers, Rept. Prog. Phys. 73 (2010) 116201.  doi:10.1088/0034-4885/73/11/116201.
  * [5] S. Goswami, The status of INO, Journal of Physics: Conference Series 888 (2017) 012025.  doi:10.1088/1742-6596/888/1/012025.
  * [6] P. Kumari et al., Improved-RPC for the CMS muon system upgrade for the HL-LHC, Journal of Instrumentation 15 (11) (2020) C11012–C11012.  doi:10.1088/1748-0221/15/11/c11012.
  * [7] F. Bossù et al., Performance of the RPC-based ALICE muon trigger system at the LHC, Journal of Instrumentation 7 (12) (2012) T12002–T12002.  doi:{10.1088/1748-0221/7/12/t12002}.
  * [8] M. Mondal, T. Dey, S. Chattopadhyay, J. Saini, Z. Ahammed, Performance of a prototype bakelite RPC at GIF++ using self-triggered electronics for the CBM experiment at FAIR, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment (2021) 166042doi:{https://doi.org/10.1016/j.nima.2021.166042}.
  * [9] A. Moshaii, L. K. Khorashad, M. Eskandari, S. Hosseini, RPC simulation in avalanche and streamer modes using transport equations for electrons and ions, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 661 (2012) S168–S171, x. Workshop on Resistive Plate Chambers and Related Detectors (RPC 2010).  doi:https://doi.org/10.1016/j.nima.2010.09.133.
  * [10] J. Datta, S. Tripathy, N. Majumdar, S. Mukhopadhyay, Study of streamer development in resistive plate chamber, Journal of Instrumentation 15 (12) (2020) C12006.  doi:10.1088/1748-0221/15/12/C12006.
URL https://dx.doi.org/10.1088/1748-0221/15/12/C12006

  * [11] P. Fonte, Survey of physical modelling in resistive plate chambers, Journal of Instrumentation 8 (11) (2013) P11001.  doi:10.1088/1748-0221/8/11/P11001.
URL https://dx.doi.org/10.1088/1748-0221/8/11/P11001

  * [12] P. Rout, J. Datta, P. Roy, P. Bhattacharya, S. Mukhopadhyay, N. Majumdar, S. Sarkar, Fast simulation of avalanche and streamer in GEM detector using hydrodynamic approach, Journal of Instrumentation 16 (02) (2021) P02018.  doi:10.1088/1748-0221/16/02/P02018.
URL https://dx.doi.org/10.1088/1748-0221/16/02/P02018

  * [13] C. Lippmann, Detector physics of resistive plate chambers (cern-thesis-2003-035) (2003).
URL https://cds.cern.ch/record/1303626/files/CERN-THESIS-2003-035.pdf

  * [14] C. Lippmann, W. Riegler, Space charge effects in resistive plate chambers, Nucl. Instrum. Meth. A 517 (2004) 54–76.  doi:10.1016/j.nima.2003.08.174.
  * [15] C. Lippmann, W. Riegler, B. Schnizer, Space charge effects and induced signals in resistive plate chambers, Nucl. Instrum. Meth. A 508 (2003) 19–22.  doi:10.1016/S0168-9002(03)01270-1.
  * [16] H. Schindler, Garfield++ User′s Guide, https://garfieldpp.web.cern.ch/garfieldpp (April 2020).
  * [17] R. Veenhof, Garfield, a drift chamber simulation program, Conf. Proc. C 9306149 (1993) 66–71.
  * [18] R. Veenhof, Garfield, recent developments, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 419 (2) (1998) 726–730.  doi:https://doi.org/10.1016/S0168-9002(98)00851-1.
  * [19] https://www.comsol.co.in/.
  * [20] N. Majumdar, S. Mukhopadhyay, S. Bhattacharya, Computation of 3D electrostatic weighting field in Resistive Plate Chambers, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 595 (2) (2008) 346–352.  doi:https://doi.org/10.1016/j.nima.2008.07.033.
  * [21] N. Majumdar, S. Mukhopadhyay, S. Bhattacharya, Three-dimensional electrostatic field simulation of a resistive plate chamber, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 602 (3) (2009) 719–722, proceedings of the 9th International Workshop on Resistive Plate Chambers and Related Detectors.  doi:https://doi.org/10.1016/j.nima.2008.12.098.
  * [22] https://www.openmp.org/.
  * [23] T. Dey, S. Mukhopadhyay, S. Chattopadhyay, J. Sadukhan, Numerical study of space charge electric field inside resistive plate chamber, Journal of Instrumentation 15 (11) (2020) C11005–C11005.  doi:10.1088/1748-0221/15/11/c11005.
  * [24] T. Dey, S. Mukhopadhyay, S. Chattopadhyay, Numerical study of effects of electrode parameters and image charge on the electric field configuration of RPCs, Journal of Instrumentation 17 (04) (2022) P04015.  doi:10.1088/1748-0221/17/04/p04015.
  * [25] Sheharyar, A., Bouhali, O., Castaneda, A., Speeding up and parallelizing the garfield++, EPJ Web Conf. 174 (2018) 06004.  doi:10.1051/epjconf/201817406004.
  * [26] O. Bouhali, A. Sheharyar, T. Mohamed, Accelerating avalanche simulation in gas based charged particle detectors, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 901 (2018) 92–98.  doi:https://doi.org/10.1016/j.nima.2018.05.061.
  * [27] ROOT Data analysis framework User′s Guide, https://root.cern.ch/root/htmldoc/guides/users-guide/ROOTUsersGuideA4.pdf (May 2018).
  * [28] R. Brun, F. Rademakers, ROOT: An object oriented data analysis framework, Nucl. Instrum. Meth. A 389 (1997) 81–86.  doi:10.1016/S0168-9002(97)00048-X.
  * [29] M. M. T. Nishimura, Mersenne twister: A 623-dimensionally equidistributed uniform pseudo-random number generator (Jan 1998).  doi:https://doi.org/10.1145/272991.272995.
  * [30] I. Smirnov, Modeling of ionization produced by fast charged particles in gases, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 554 (1) (2005) 474–493.  doi:https://doi.org/10.1016/j.nima.2005.08.064.
  * [31] A multiterm boltzmann analysis of drift velocity, diffusion, gain and magnetic-field effects in argon-methane-water-vapour mixtures, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 283 (3) (1989) 716–722.  doi:https://doi.org/10.1016/0168-9002(89)91446-0.
  * [32] S. Biagi, Monte carlo simulation of electron drift and diffusion in counting gases under the influence of electric and magnetic fields, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 421 (1) (1999) 234–240.  doi:https://doi.org/10.1016/S0168-9002(98)01233-9.
  * [33] W. H. Furry, On fluctuation phenomena in the passage of high energy electrons through lead, Phys. Rev. 52 (1937) 569–581.  doi:10.1103/PhysRev.52.569.
URL https://link.aps.org/doi/10.1103/PhysRev.52.569

  * [34] H. Schindler, Microscopic Simulation of Particle Detectors, presented 13 Dec 2012 (2012).
URL https://cds.cern.ch/record/1500583

  * [35] J. M. Bull, Measuring synchronisation and scheduling overheads in openmp, in: In Proceedings of First European Workshop on OpenMP, 1999, pp. 99–105.
  * [36] M. K. Bane and G. D. Riley, Extended Overhead Analysis for OpenMP, Springer Berlin Heidelberg, Berlin, Heidelberg, 2002, pp. 162–166.  doi:https://doi.org/10.1007/3-540-45706-2\\_20.
  * [37] Z. He, Review of the shockley ramo theorem and its application in semiconductor gamma-ray detectors, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 463 (1) (2001) 250–267.  doi:https://doi.org/10.1016/S0168-9002(01)00223-6.
  * [38] M Abbrescia et al., The simulation of resistive plate chambers in avalanche mode: charge spectra and efficiency, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 431 (3) (1999) 413–427.  doi:https://doi.org/10.1016/S0168-9002(99)00374-5.
  * [39] M. Kobayashi, An estimation of the effective number of electrons contributing to the coordinate measurement with a TPC, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 562 (1) (2006) 136–140.  doi:https://doi.org/10.1016/j.nima.2006.03.001.
  * [40] P. Fonte, V. Peskov, High-resolution TOF with RPCs (2001) 6 pdoi:10.1016/S0168-9002(01)01914-3.