Effect of passive metallic layers on muon energy estimation by means of deflection angle for muon scattering tomography: A comparative study based on GEANT4 simulations

During the encounter between the primary cosmic rays and the earth atmosphere, a non-negligible number of muons are generated over a wide energy spectrum. The fundamental basis, on which the muon scattering tomography is founded, is to follow the propagation of the cosmic ray muons within the volume-of-interest (VOI) where the entering muons of a certain energy deflect from their initial directions in the wake of the physical processes primarily hinging on the atomic number, the material density, and the material thickness [1, 2, 3]. Among the detection modules existing in the tomographic setups based on the muon scattering are the plastic scintillators that have substantially found their application by accentuating their favorable aspects like fast rise and decay times, high optical transmission, ease of manufacturing, low cost, and large available size [4]. The hodoscope structure for the scattering-based tomography consists of two sections that are installed atop and beneath the VOI under the investigation [5], and each section is composed of two or more distinct detector layers [6, 7] occasionally made out of plastic scintillators usually with a moderate thickness.

Typical muon tomography systems are inherently incapable of directly measuring the kinetic energy of incoming muons, as that would demand the presence of a magnetic field or of Cherenkov detectors or of a very precise time-of-flight measurement, and all those options would make the cost of the apparatus increase by a large factor. However, the capability of assessing, even roughly, the energy of the incoming muons on an event per event basis can significantly enhance the precision of the tomography [8], as the deflection process depends directly on energy. In the course of the muon propagation through the detection system, the hodoscope components slightly contribute to the deviation of the transversing muons up to a certain extent, and this tiny contribution might serve to categorize the detected muons by building a binary relation between the deflection angle and the muon energy [9].

By recalling the intrinsic properties of the plastic scintillators such as low atomic number and low density apart from a low thickness on the basis of the deflection angle, the strategic introduction of the passive metallic layers is expected to numerically ameliorate the angular outcomes for the purpose of the energy estimation thanks to the deflection power of the inserted material, which leads to the visible reduction in the angular uncertainty as well as the drastic augmentation in the tiny angular values. In this study, we categorize the incoming muons into three groups based on their deflection angle within the hodoscope that we use as a proxy for their kinetic energy. We compare the classification performance in this three-group energy structure between our current tomographic prototype that contains three detector layers manufactured from polyvinyl toluene with a thickness of 0.4 cm as well as an accuracy of 1 mrad in both the top section and the bottom section, and an alternative hodoscope scheme including stainless steel layers in addition. In order to quantify the uncertainty propagated through the energy groups, we present a couple of misclassification probabilities characterized by one standard deviation where the first means supposes a linear finite approximation in one dimension, while the latter procedure is governed by the positively defined modified Gaussian distributions in two dimensions. The present study is organized as follows. In section 2, we express the average deflection angle as well as the corresponding standard deviation in terms of the hit positions and we also introduce the corresponding expressions averaged over the top section as well as the bottom section to diminish the width of the deflection angle. We define two forms of misclassification probabilities based on the standard deviations in section 3. While we demonstrate our simulation setup in the GEANT4 code [10] besides the simulation features in section 4, we exhibit our simulation results in section 5. Finally, section 6 narrates our concluding remarks.

By reminding that the crossing muons undergo deviations from the initial trajectory owing to the interactions with matter, the deflection angle of the incoming muons throughout a set of three detector layers is estimated via constructing two distinct vectors where the first vector is defined in accordance with the hit locations on the first two detector layers, whereas the remaining vector is the sequel of the hit positions on the last two detector layers as illustrated in Fig. 1.

Upon these two vectors that suffice the computation of the angular deviation, the deflection angle of a muon crossing the detector layers denoted by $\theta$ is obtained by using the following expression [11, 12, 13]:

$$\theta=\arccos\left(\frac{\vec{v}_{1}\cdot\vec{v}_{2}}{\left|v_{1}\right|\left|v_{2}\right|}\right)$$

(2.1)

Since the detector layers accept a significant number of muons, the average deflection angle at a particular energy is determined by averaging the previously calculated deflection angles over $N$ number of the non-absorbed/non-decayed muons as indicated in

$$\bar{\theta}=\frac{1}{N}\sum_{i=1}^{N}\theta_{i}$$

(2.2)

which refers to a standard deviation as written in

$$\delta\theta=\sqrt{\frac{1}{N}\sum_{i=1}^{N}(\theta_{i}-\bar{\theta})^{2}}$$

(2.3)

In light of the fact that the deflection angle is a consequence of a stochastic process, the standard deviation of the deflection angle is anticipated to be lessened with the intention of a better energy estimation. By harking back to the entire tomographic setup consisting of two separate sections (i.e. the top section and the bottom section), each of which includes three detector layers, the deflection angle determined for each non-absorbed/non-decayed muon in each section is first averaged over the number of sections, then the mean deflection angle of both sections is re-averaged over the number of the non-absorbed/non-decayed muons. Essentially, the following term outlines the average deflection angle of two separate hodoscopes specified by $x$ and $y$ at a certain energy value over $N$ number of the non-absorbed/non-decayed muons:

$$\bar{\theta}_{\frac{x+y}{2}}=\frac{1}{N}\sum_{i=1}^{N}\frac{\theta_{x,i}+\theta_{y,i}}{2}$$

(2.4)

Accordingly, the resulting standard deviation dependent on the contributions from both the top hodoscope as well as the bottom hodoscope is formulated as noted in

$$\delta\theta=\sqrt{\frac{1}{N}\sum_{i=1}^{N}\biggl{(}\frac{\theta_{x,i}+\theta_{y,i}}{2}-\bar{\theta}_{\frac{x+y}{2}}\biggr{)}^{2}}$$

(2.5)

By highlighting the probability that even fairly different muon energies routinely generate a substantial set of similar deflection angles, it is possible to foresee that the standard deviations at discrete energies are supposed to coincide in a considerably common sub-interval. In the first method, by assuming that two adjacent energy groups, i.e. labeled by $A$ and $B$ in a descending order, intersect in terms of the deflection angle by virtue of one standard deviation in one dimension, the misclassification probability is the ratio between the overlapping length $\bar{\theta}_{B}+\delta\theta_{B}-\bar{\theta}_{A}+\delta\theta_{A}$ and the entire length $\bar{\theta}_{A}+\delta\theta_{A}-\bar{\theta}_{B}+\delta\theta_{B}$ as follows

$$P_{\rm Linear}=\begin{cases}\displaystyle\frac{\mbox{Overlapping length}}{\mbox{Entire length}}=\displaystyle\frac{\bar{\theta}_{B}+\delta\theta_{B}-\bar{\theta}_{A}+\delta\theta_{A}}{\bar{\theta}_{A}+\delta\theta_{A}-\bar{\theta}_{B}+\delta\theta_{B}}&\bar{\theta}_{B}+\delta\theta_{B}>\bar{\theta}_{A}-\delta\theta_{A}\\ 0&\bar{\theta}_{B}+\delta\theta_{B}\leq\bar{\theta}_{A}-\delta\theta_{A}\\ \end{cases}$$

(3.1)

In other words, two overlapping error bars defined by one standard deviation at two neighboring energy values are explicitly referred to the finite linear approximation in Eq. (3.1).

In the second procedure, we consider that two positively defined modified Gaussian probability density functions (PDFs) represent the mean deflection angles and its standard deviations at kinetic energies $A$ and $B$, respectively, as follows

$$G^{\prime}(\bar{\theta}_{A},\delta\theta_{A},\theta)=\frac{G(\bar{\theta}_{A},\delta\theta_{A},\theta)}{\displaystyle\int_{0}^{\infty}G(\bar{\theta}_{A},\delta\theta_{A},\theta)d\theta}$$

(3.2)

and

$$G^{\prime}(\bar{\theta}_{B},\delta\theta_{B},\theta)=\frac{G(\bar{\theta}_{B},\delta\theta_{B},\theta)}{\displaystyle\int_{0}^{\infty}G(\bar{\theta}_{B},\delta\theta_{B},\theta)d\theta}$$

(3.3)

Then, the overlapping coefficient (OVL) of these two Gaussian PDFs might be obtained by calculating the intersection area as shown in

$$\mbox{OVL}=\mbox{Intersection area}=\displaystyle\int_{0}^{\infty}\min[G^{\prime}(\bar{\theta}_{A},\delta\theta_{A},\theta),G^{\prime}(\bar{\theta}_{B},\delta\theta_{B},\theta)]d\theta$$

(3.4)

Thus, in the case where the intervals constitute a Gaussian distribution, the misclassification probability is determined by using the following expression:

$$P_{\rm Gaussian}=\frac{\rm OVL}{\displaystyle\int_{0}^{\infty}G^{\prime}(\bar{\theta}_{A},\delta\theta_{A},\theta)d\theta+\displaystyle\int_{0}^{\infty}G^{\prime}(\bar{\theta}_{B},\delta\theta_{B},\theta)d\theta-{\rm OVL}}=\frac{\rm OVL}{\rm 2-OVL}$$

(3.5)

The present study is conducted by means of the GEANT4 simulations, and we initially track the muon hit locations on the detector layers made out of polyvinyl toluene in order to calculate the average deflection angles as well as the corresponding standard deviations for the associated groups of kinetic energies. The simulation geometry for our present prototype without passive metallic layers in addition to the alternative hodoscope scheme containing stainless steel layers is depicted in Fig. 2. As described in Fig. 2(a), the dimensions of each the detector layer fabricated from polyvinyl toluene are $100\times 0.4\times 100$ $\rm cm^{3}$, and the inter-layer spacing at each section is 10 cm. In the proposed hodoscope setup, where we introduced stainless steel layers as illustrated in Fig. 2(b), the dimensions of each metallic layer are $100\times 0.4\times 100$ $\rm cm^{3}$, i.e. identical to those of the detector layers. In both the tomographic schemes, the gap between the top hodoscope and the bottom hodoscope is 100 cm into the bargain.

A central mono-directional beam with the uniform energy distribution that is generated at y=85 cm via G4ParticleGun is employed, and the generated muons are crossing in the vertically downward direction, i.e. from the top edge of the simulation box through the bottom edge as indicated with the black arrows in Fig. 2. Since the current aperture of the entire detection geometry commonly only accepts the narrow angles apart from the very rare entries around the corners, this beam setup is considered significantly reliable by reminding that the distribution of the incident angle ($\alpha$) approximately corresponds to $\cos^{2}\alpha$ for an interval between $-\pi/2$ and $\pi/2$ [14]. We favor a uniform energy distribution lying on an energy interval between 0.25 and 8 GeV so as to achieve the absorption minimization especially in the case of stainless steel layers besides the numerical accuracy optimization [8]. We partition the entire energy interval bounded by 0.25 and 8 GeV into three energy groups such that the first energy group consists of the energy values between 0.25 and 0.75 GeV, whereas the second energy group is composed of the muon energies between 0.75 and 3.75 GeV, and the muon energies between 3.75 and 8 GeV finally constitute the third energy group.

For each of the three energy groups, the number of the simulated muons is 10⁴. All the materials in the current study are defined in accordance with the GEANT4/NIST material database. The reference physics list used in these simulations is FTFP$\_$BERT. The muon tracking is sustained by G4Step, and the registered hit positions are post-processed with the assistance of a Python script where the deflection angle is first determined for every single non-absorbed/non-decayed muon of a given energy group at each hodoscope. Finally, we conclude the post-processing stage by applying Eqs. (2.4) and (2.5) on the datasets of deflection angle attained through the top hodoscope and the bottom hodoscope, thereby obtaining the average deflection angle and the corresponding standard deviation.

Our GEANT4 simulations initially yield the average deflection angles as well as the associated standard deviations for both our current prototype and our proposed new scheme as tabulated in Table LABEL:avgangle. By contrasting the present prototype with the alternative hodoscope in agreement with the average deflection angles in Table LABEL:avgangle, the first positive influence of stainless steel layers in terms of the detection performance is already revealed. To better elucidate, by reminding the existing detector accuracy of 1 mrad, all the average deflection angles generated by the proposed hodoscope containing stainless steel layers exceed the angular value of 1 mrad for the three-group energy structure; on the other hand, the average deflection angles obtained through the present configuration remains below this detector accuracy in two among three energy groups, i.e. $\bar{E}=2.25$ GeV and $\bar{E}>3$ GeV.

On top of the augmented average deflection angles that lead to a better detection efficiency, the presence of stainless steel layers is partially expected to result in the diminished uncertainty. To verify this hypothesis over our GEANT4 simulations, the linear overlap based on the standard deviations is shown in Fig. 3, and we observe that the uncertainty in the deflection angle is apparently decreased in the case of stainless layers compared to the present hodoscope scheme without any passive metallic layers.

In order to further evaluate the effect of stainless steel layers, the areal overlap of the positively defined modified Gaussian distributions as defined in Eqs. (3.2) and (3.3) is displayed in Fig. 4, and we face a similar reduction in the angular uncertainty for the reason of the areal contraction.

In the long run, the misclassification probabilities calculated by means of these methodologies are listed in Table 2, and it is numerically exhibited that the proposed hodocope scheme with stainless steel layers leads to the remarkably lower misclassification probabilities in comparison with the present tomographic setup without stainless steel layers.

In this study, we investigate an alternative tomographic setup by introducing stainless steel layers into the bottom section and the top section with the aim of the uncertainty reduction in the deflection angle. Upon our GEANT4 simulations combined with two distinct forms of the misclassification probabilities, we experience that the stainless steel layers slightly yield more favorable angular values with the diminished uncertainty as opposed to our current prototype apart from the increased average deflection angles that fulfill the present detector accuracy. Within the perspective of an indirect energy classification dependent on the deflection angle, our simulation outcomes qualitatively and quantitatively indicate that the presence of stainless steel layers might be an inexpensive upgrade in the potential applications of the muon scattering tomography.