Detector geometry design

Three components were simulated based on GEANT4.10.6p-03 software, as shown in Fig. 1. Figure 1(a) shows the structural model of a single detector. The material structure of an EJ-200 plastic scintillator was used as a reference. The dimensions of the model were set to 800 mm × 800 mm × 50 mm, and the material was filled with molecules with a chemical formula of C₉H₁₀, density of 1.023 g/cm³, scintillation light yield of 10,000 photons/MeV, and peak emission wavelength of 425 nm. Based on the emission spectrum of EJ-200 [34] and the relationship between the reflectivity of the reflector layer and collection efficiency, a reflector layer with a reflectivity of 0.88 was placed on the surface of the plastic scintillator. Assuming a wavelength of 425 nm for the scintillation photons in the EJ-200 plastic scintillator, the estimated velocity (v) was approximately 0.633c. PMTs were constructed at the four corners of the detector field for response simulation. The effective diameter of the photocathode was set as 46 mm. As shown in Fig. 1(b), two sets of track-detector modules were positioned separately at the top and bottom of the object under investigation. This arrangement enabled muon-scattering imaging based on the varying scattering angles of the tracks.

Fig. 1Full-screen View

(a) Single-layer detector and material parameters. (b) Muon-scattering imaging system. (c) Basic performance simulation

Basic performance simulation

The incident particle source was generated using an independent particle-source library (CRY) [35]. Owing to the inherent incident angles of natural muons, to determine the muon trajectory, an incident muon was recorded, and the data of the response photons collected by PMTs were analyzed only when a signal response was observed from both layers of detectors.

The dimensions and optical parameters of the detector model components considered in the simulation were as follows. EJ-200 had a refractive index of 1.58. Based on the relationship between reflectance and efficiency shown in Fig. 1(a), the reflective foil was modeled as an additional volume located 0.1 mm away from the scintillator surface to enhance the light collection on the PMT photocathode. To simulate the actual conditions of the PMT window, as shown in Fig. 1(c), a thin volume (0.2 mm thick) of the optical-coupling material is defined with silicon grease with a refractive index of 1.47. The area conversion from the square scintillator-end face (area S₁ = 25 cm²) to the PMT sensitive-area surface (area S₂ = 21.16 cm²) and the coupling effects between the scintillator and PMT must be considered. The critical angle for the total internal reflection ${\theta }_{\text{c}}={{\displaystyle \sin }}^{-1}\left(n_{\text{grease}}/n_{\text{scint}}\right)$ was calculated as 68.5 °. The lower half of Fig. 1(c) illustrates the method based on the Geant4 simulation and PMT response-function measurements. This method involves recording the time at which each photon reaches the PMT, obtaining the distribution of photon-arrival times at the PMT photocathode, and using a photon response function to obtain simulated pulse signals [36]. When a single photoelectron arrives at the photocathode, the PMT generates a pulsed signal. The pulse waveform of a single photoelectron was simulated using the following time response function [37]: $v_i\left(t\right)=G{C}_{\text{e}}\frac{t^2{e}^{-\frac{t^2}{{\tau }^2}}}{{\displaystyle \int }{t}^2{e}^{-\frac{t^2}{{\tau }^2}}\text{d}t}$, where G is the gain of the PMT, C_e represents the charge-to-voltage conversion factor, and τ represents the time constant of 2.5 ns, which is the rise time of the PMT detector. Because the arrival time of each photoelectron at the PMT is different, the simulated output signal of the PMT is obtained by summing the pulse waveforms of individual photoelectrons, as follows: $V_{\text{PMT}}\left(t\right)={\displaystyle {\sum }_{i=1}^n{v}_i}\left(t\right)$.

Selection of machine-learning models

To further improve the accuracy of muon-hit-position reconstruction, machine-learning algorithms were introduced for optimization. Four commonly used machine-learning models were investigated: Gaussian process regression (GPR), support vector machines (SVMs), decision trees (DTs), and backpropagation (BP) neural networks. The position coordinates obtained from the localization algorithm described in the previous section were included as features in the feature dataset. The features in this dataset also included the peak and start times of the PMT responses obtained from the simulations. A flowchart of the machine-learning models is shown in Fig. 3(a).

Fig. 3Full-screen View

(Color online) (a) Training process of the machine-learning models. (b) Selection of data points for the training sets. (c) Box plots showing the performance of the four machine-learning models. (d). Average error and residual plots of the validation set for four machine-learning models

When training a model, data points that exhibit significant differences must be selected to improve the generalizability of the model. This prevents the training results from being biased towards frequently occurring data, thereby ensuring a balanced data distribution and preventing the model from becoming ineffective. However, it is crucial to guard against overfitting, which can result in a model that performs well on the training set but poorly on the test set. To mitigate this, a portion of the training data should be randomly reserved as a validation set. Using a validation set to tune the hyperparameters of the machine-learning model, we can avoid overfitting the training set and improve the model’s generalizability. Box plots provide a visual representation of the dataset’s central tendency, dispersion, and outliers, allowing the evaluation of the performance of the machine-learning algorithms. The box length in a box plot represents the dispersion of the data. The median position indicates the central tendency of the data, whereas the upper and lower quartiles (Q1 and Q3) represent the degree of data dispersion. Points outside the box, called outliers, represent data values outside the error range.

During the data-training process, the same algorithm was used to build models for the x- and y-axis coordinate-data values separately [24], with the aim of improving the accuracy of the model predictions. As shown in Fig. 3(b), the simulated muons hit the scintillator plate at the marked positions, and the PMT response data were obtained to create a feature dataset. These hit points were spaced 10 mm apart, and a total of 5776 known positional data points were used to train and validate the model. To mitigate the effects of randomness and augment the training dataset, we simulated each point five times, resulting in a total of 28880 datasets. The test dataset was independent of the training dataset. Following the muon-incidence pattern, the muons were randomly incident on the detector field, and the collected data were stored as the test set.

The performances of the four machine-learning models are shown in Fig. 3(c). The box area for the SVM model is the largest among the four models, indicating the poorest performance. The box area for the BP model falls between those of the GPR and DT models, indicating a moderate performance. Both the GPR and DT models have smaller box areas; however, the GPR model has fewer discrete points. This suggests that among the four models, the GPR model exhibits the best performance. To verify whether the models were overfitting, a validation set was used to assess their performance. Because the dataset was sufficiently large, randomly selecting 350 data points for the validation set did not affect the conclusions. The error results of the four machine-learning models on the validation set are shown in Fig. 3(d). We used the root mean square error as a performance metric for the machine-learning model, which amplified the performance differences of the model under large errors. The mean absolute error was also used to assess the average prediction accuracy of the model. Among the four models, the SVM model exhibited the largest differences in the residual values, indicating the worst generalization ability. However, the GPR model exhibited a smoother trend and the smallest differences in residual values, indicating the best generalization ability. The results of the residual plot were consistent with those of the box plot, indicating that overfitting did not occur and the model performance was reliable. Therefore, subsequent research used the GPR model to improve the positioning accuracy.

Position resolution under the GPR algorithm

This section presents the validation of the effectiveness of the simulation results by comparing the reconstructed muon-hit position images obtained in the experiments with those generated in the simulations. In the experiment, approximately 4×10⁴ valid muon events were obtained through data selection and then divided into training, validation, and test sets in a 2:1:1 ratio. The collected trigger time, time over threshold, and values obtained through the time-difference algorithms were used as the feature library for the GPR algorithm. In the simulation, the corresponding time information obtained through the Geant4 simulation served as the feature library for the GPR algorithm with an equal collection of 4d4 valid events. The reconstructed position results are presented in Fig. 4. The distribution of the reconstructed muon positions for the MWDC is shown in Fig. 4(a), exhibiting a trend of more muons at the center and fewer muons at the edges. This is attributed to the pronounced edge effects of the outermost wires in the MWDC, which resulted in a lower count of valid muon data near the edges after selection. The project also compared the muon reconstruction effects of plastic scintillators with three- and four-corner PMT readouts, as depicted in the third image in Fig. 4(a) and the first image in Fig. 4(b). The results indicate that using a three-corner PMT leads to a significant reduction in the accuracy of position reconstruction at the corner because of missing data from the PMT readout of the corner, causing a noticeable inward shift of the reconstructed points in that region. The position reconstruction results for plastic scintillators with four-corner PMT readouts are shown in Fig. 4(b), where the lower data volume near the edges caused by the edge effects of the MWDC detector led to a training imbalance in the feature values. Consequently, the reconstruction results at the edges of the detector were inferior to those in the central region. The reconstructed position results from the simulations are shown in Fig. 4(c). The results exhibit a more even distribution of data on the detector panel, resulting in better position reconstruction overall. However, the reconstruction results at the edge positions tended to deviate from the actual edges and shifted towards the interior. From Sect. 2, we can infer that we have calculated the critical angle for the total internal reflection on the surface of the PMT to be approximately 68.5 °. Thus, in this situation, if the muon impact point is located at the edge of the detector, the scintillation light travels the shortest path to reach the two adjacent PMTs, exceeding the acceptance angle of the PMT photocathode. However, the scintillation light that reaches the PMT via a longer path owing to refraction remains unaffected. This results in a delay in the arrival time of the collected scintillation photons at these PMTs compared with the actual time. However, the propagation of the scintillation light to distant nonadjacent PMTs was not affected by this phenomenon. Therefore, the reconstructed image tended to bend inward at the edge positions.

Fig. 4Full-screen View

(Color online) Comparison of experimental and simulated results. (a) Muon distribution measured by MWDC with muon-hit-position reconstruction under the three-corner readout for a 400 mm × 400 mm effective area. (b) Muon-hit-position reconstruction and position resolution obtained under the GPR algorithm for a 400 mm × 400 mm effective area. (c) Position reconstruction and position resolution under simulation

The experiment yielded an average error of σx = 20.9 mm on the x-axis, σy = 18.6 mm on the y-axis, and an overall error of σ = 27.9 mm. In the 800 mm × 800 mm position reconstruction simulation, the average error on the x-axis was σx = 9.6 mm, that on the y-axis was σy = 10.9 mm, and the position resolution was σ = 14.5 mm.

Muon-track reconstruction

To investigate the effectiveness of the reconstruction of the muon trajectories in the detector system, as shown in Fig. 5(a), when the distance between the detector panels changes, the corresponding angles vary, where θ_1max<θ_2max. Therefore, we simulated the reconstruction results for 100 muon trajectories to determine the effective angular resolution of the detector in this configuration. The angular distributions of the reconstructed trajectories are presented in Fig. 5(b).

Fig. 5Full-screen View

(Color online) (a) Muon-track reconstruction. (b) Detector receivable angle

The muon flux at sea level followed a cos² θ distribution relative to the zenith angle. Considering the detector-acceptance angles, the detected muon flux relative to the zenith angle exhibited a distribution of $2\pi I_0{\displaystyle \int }I\left(\theta \right)\sin \theta \text{d}\theta$. A comparison of this result with Fig. 5(b) shows that, despite some trajectory deviations due to precision errors, the overall distribution of muons aligned with the expected pattern within the detector structure.