Methodologies

For the JUNO detector, the PMTs closely packed around the LS sphere form a PMT shell, which is a 2D surface instead of a 3D surface, so the PMT charge hit pattern can be projected onto a 2D plane to perform convenient fitting. Similar to the world map of the Earth, there are several projection methods to obtain a 2D map, such as the longitude–latitude projection method and the spherical projection method, but every method introduces distortion through the projection. For reconstruction, different projection methods yielded similar reconstruction results. Finally, we selected the longitude–latitude projection method. The ϕ angle of the PMT position is for the X-axis on the 2D plane and the θ angle is for the Y-axis.

When a muon passes through the LS, more than 90% of the photons are emitted from the scintillation process, which emits isotropic light. The closest PMTs to the muon incident point or exit point will collect the most PEs, which will create two clusters on the PMT charge-hit pattern (Fig. 3). Figure 5 shows an example of the distribution of the number of photoelectrons (nPE) for a single through-going-muon event and a double through-going-muon event. PMTs with a charge greater than 75 PEs generally correspond to clusters caused by muon tracks. To highlight the clusters that were caused by muon incident and exit points, we needed to eliminate the PMT noise (∼1 PE), radioactivity (several PEs), and PMTs far from clusters (tens of PEs). As a result, PMTs with charges smaller than 75 PEs are discarded. In a simulation study using the JUNO offline software [12, 13], when the maximum discard charge varies from 70 to 80 PEs, the reconstruction result is not changed significantly; more details can be found later (Section 4.2 and Fig. 14).

Fig. 14Full-screen View

(Color online) (a) Reconstruction efficiency for single and double through-going muons. For double through-going muons, the distances to the center are calculated according to their single track. (b) The effect of PE cut (Sect. 3.1) on reconstruction efficiency.

For the PMT charge pattern, a 2D peak searching and fitting algorithm named “TSpectrum2” was investigated and applied. “TSpectrum2” is a ROOT tool [34], which is based on a 2D Gaussian function fitting and it searches for peaks in source spectrum with a deconvolution method. “TSpectrum2” can obtain the maximum local peak by calculating the spectrum of the standard deviations of the smoothed second derivatives in two dimensions [35]. If the fit is successful, the number of identified peaks and their locations are returned. The identified peaks were arranged according to their heights in the spectrum. In our study, the maximum number of clusters was set to seven to list the highest seven peaks only, which is sufficient for searching for single muons (2 peaks), double muons (4 peaks), and triple muons (6 peaks). The tracking precision of triple muons is much lower than that of double muons, so it is not described in this paper, but we still set the maximum to seven to retain the reconstruction potential for triple muons when using this method, which will be studied further in the future. “TSpectrum2” runs very fast (Sect. 4.3) and it can mathematically find the cluster seed candidates.

However, the “TSpectrum2” fitting usually finds more seeds than at the muon incident and exit points. The unmatched seeds were caused by the statistical fluctuation of the PMT charge. In the simulation, the collected charges of adjacent PMTs surrounding the largest charge collection PMT do not continuously drop from the cluster center to the cluster edge. However, some muons have a relatively large energy deposit at a certain position along the track. All of these may result in an incorrect cluster number and fake cluster seeds when matching with the true muon hit points. Figure 6a shows an example of the PMT charge hit pattern of a single through-going-muon event before charge smoothing above the PMT charge threshold. This shows that the fitting finds three cluster seeds near the muon exit point, where we need only one cluster seed. The fitting results for the cluster seeds must be carefully checked and analyzed.

Fig. 6Full-screen View

(Color online) Hit patterns of a single through-going-muon event before and after smoothing. (a) The PMTs with charges smaller than 75 PEs are discarded; as a result, the pixels only with small charges are not shown. The red triangles are the cluster seeds fitted by “TSpectrum2.” Two purple circles indicate the muon true incident and exit points on the pattern. The black stars are the final reconstructed muon track incident and exit points. (b) The smoothed charge cluster. The same labels as in Fig. 3 are used. The cluster seeds (red triangles) were reduced from three to two for the bottom cluster.

To discard fake cluster seeds and obtain true hit points, a charge-weighted algorithm was developed. First, based on each cluster seed, we calculated the charge-weighted position of the fired PMTs in a 4-m radius (corresponding to the size of the normal cluster) using the following formula: $\overline{x}=\frac{{\displaystyle {\sum }_{i=0}^n{q}_i{x}_i}}{{\displaystyle {\sum }_{i=0}^n{q}_i}}$ (1) where xi is the spatial position of PMT_i, n is the number of fired PMTs in a 4-m radius around the cluster seed, and qi is the charge collected by PMT_i. If the total charge in a 4-m radius around the cluster seed is less than 1000 PEs, the cluster seed is discarded.

Second, after the calculation, replace the cluster seed with the charge-weighted position ($\overline{x}$). Then, if two new cluster seeds are too close (<4 m), they are merged into one point according to their charge weights. This step can filter out most of the fake seeds given by “TSpectrum2” fitting and give a correct position. Third, if the gathered cluster charge within a radius of 4 m around the cluster seed is too small (<10% of the total charge of all clusters), the seed is also discarded. This can eliminate the small-charge clusters found by “TSpectrum2” fitting which cannot match with the true incident and exit points of the muons. Finally, the charge-weighted positions of the remaining clusters were considered as muon incident and exit points.

Distortion elimination and performance optimization

To make the muon track reconstruction more precise, we eliminated distortion using the rotation method, removed the statistical fluctuation of the charge cluster using the smoothing method, and transferred the cluster seeds on the PMT charge pattern onto the true muon track.

As mentioned previously, the projection method (normal projection) causes non-negligible distortion, particularly near the polar zone. To avoid projection distortion, we rotated the PMT positions 90° around the X-axis to allow the north and south poles to be placed at the equator (Fig. 7), and then, the new PMT positions are again projected onto θ–ϕ 2D plane (called rotation-X projection). The cluster seeds in the θ range (-50°, 50°) were selected in both the normal projection (before the rotation) and rotation-X projection (after the rotation). Some repeated seeds in the overlap area of the two types of projection are merged during the charge-weighted step.

Fig. 7Full-screen View

(Color online) Rotation method of PMTs to reduce the distortion of projection. The dotted line is the equator and the solid line is the muon track. After rotating 90° around the X-axis, the north and south poles are located at the equator.

Sometimes, some clusters are split from one cluster into two clusters when the incident or exit points are located at the edge of the 2D plane (ϕ > 160° or ϕ < -160°). First, the PMTs were rotated 90° around the Z-axis (rotation-Z projection) to obtain the precise position of the clusters that were located at the edge of the 2D plane. Then, combined with the normal projection and rotation-Z projection, we check whether some fitted points are located at the edge. If a cluster is split into two parts, the fitting result from the rotation-Z projection is used.

In addition, apart from the peaks induced by the incident and exit points of the muon, there were also some relatively small peaks. This is because of the statistical fluctuation in the PE number collected by the PMTs or a relatively large-energy deposit at a certain location along the muon track. Figure 6a shows all the reconstructed peaks. Meanwhile, the “TSpectrum2” tool may supply fitted points at the non-hit position. To handle this effect, each pixel's charge was smoothed (q^s) based on the adjacent pixel charges using the smoothing formula $q^s=\frac{{\displaystyle {\sum }_{i=0}^N\left(\frac{w_0}{w_1+d_i}{q}_i\right)}}{{\displaystyle {\sum }_{i=0}^N\left(\frac{w_0}{w_1+d_i}\right)}}$ (2) where qi is the charge of the ith adjacent pixel, w₀ and w₁ are the weight factors, with the best tuning result w₀/w₁ = 4/9, and di is the distance factor between the smoothing pixel and the ith adjacent pixel. The distance factor of the nearest pixel was one. The further the pixel distance, the larger di is and the lower the weight of qi. Here N is the number of adjacent pixels that are used for charge smoothing and surround the specified pixel, which is set to 24. A comparison of Figs. 6a and b shows that the charge smoothing method can make the cluster peaks protrude further, reduce the number of false fitting peaks, and can also improve the reconstruction efficiency.

In fact, the charge cluster center is not exactly equal to the true incident and exit points of the muon. This is because there is a water buffer with a thickness of 1.8 m between the PMT spherical shell and the LS sphere (Fig. 8). When a muon enters the PMT shell, most photons originate from the LS, and the PMT charge cluster is close to the muon hit point of the LS. As a result, the track reconstructed by the PMT charge cluster center does not correspond to the true muon track. For this reason, the cluster charge centers (yellow clusters in Fig. 8) must be translated into true muon hit points on the LS sphere (red clusters in Fig. 8), according to their spatial positions. Finally, we reconstructed the muon track (green dotted line) using the reconstructed incident and exit points.

Fig. 8Full-screen View

(Color online) Schematic diagram of a single-muon track reconstruction with two clusters. In JUNO detector, there is a water buffer with a thickness of 1.8 m between the PMT spherical shell and LS sphere. The cluster charge centers (yellow clusters) make a fake track and need to be corrected by a geometric effect to obtain the reconstructed track.

Double-muon event reconstruction strategies

As listed in Table 1, a double-muon event typically corresponds to four hit points. Every two points comprise a track candidate. There are three combinations of four hit points, and two of them are in parallel; thus, they can be considered as muon bundle candidates, as shown in Fig. 7(a). The PMT timing information was used to further identify the true muon tracks.

With each reconstructed point in the center, the FHTs of the surrounding PMTs within 4 m were averaged as the cluster time. The two earlier times were treated as two incident points, and the other two points were treated as two exit points. This time strategy can help to determine the correct combination.

If the two tracks are too close, then only two cluster seeds can be found (Fig. 9(b)). In this case, the two cluster seeds are connected and reconstructed into a single track. Although the two tracks are close, the overlapped cluster is larger than the normal cluster. The distance between two muon tracks can be deduced from the overlapped cluster information. The primary requirement of muon track reconstruction is to veto muon-induced isotopes through a cylindrical volume cut along the muon track. Although only one track is reconstructed, a cylindrical cut with a larger radius will be considered and can provide the same background reduction ability.

Fig. 9Full-screen View

(Color online) Plots of different double-muon events hitting the CD and their clusters. (a) Double-muon event which has two parallel tracks hits the detector almost at the same time, resulting in four isolated clusters. (b) Two close tracks in a muon event hit the CD and only resulted in two clusters because of the overlap. (c) Two relatively close (∼4 m) tracks with a small angle (< 2°), the two incident points correspond to one cluster and the two exit points correspond to two clusters. (d) One clipping track in a double-muon event and the resulting three clusters. (e) One stopping-muon event in a double muon resulting in three clusters.

If double muons have only three charge clusters, it is not easy to distinguish between the two muon tracks. In this study, we identified three cases in a three-cluster situation (Figs. 9(c)–(e)). One case occurred when two incident points overlapped (Fig. 9(c)). We calculated the average FHT ($\overline{T}$) of the surrounding PMTs within 4 m of each cluster. Because two tracks in the double muons shoot in the LS at the same time, the $\overline{T}$ of the incident or exit points are close. In three values of $\overline{T}$, if the two $\overline{T}$ differences are less than 5 ns, both clusters are recognized as exit points and the other is the incident point. Then, the two muons can be correctly tracked using the incident point twice. Although the incident point is not very precise, the larger radius cylindrical cut according to the incident cluster size along each track can still provide the same background reduction of the muon-induced isotopes.

The other two cases, clipping (Fig. 9(d)) or stops in the CD (Fig. 9(e)), also produced three clusters. The projection on the 2D charge histogram is shown in Fig. 10, respectively. Although the two clusters' time ($\overline{T}$) difference is larger than 5 ns, we still need to select two seeds to make a track. A direction determination method was developed for the 2D cluster plot. When a muon track in the LS is further from the CD center, such as a clipping muon, the PMT charge cluster is shaped like an ellipsoid. The direction of the elliptical major axis was along the muon track direction. The major axis directions of the ellipse can be expressed by the vectors $\overleftarrow{V_A}$, $\overleftarrow{V_B}$, and $\overleftarrow{V_C}$ for the three clusters (Fig. 10). The angles of the candidate track (LAB, LBC, or LAC) with its corresponding cluster vectors can be calculated using the following equations: $\begin{array}{l}\angle 1=\angle \overleftarrow{V_A}\&L_{AB}+\angle \overleftarrow{V_B}\&L_{AB}\hfill \\ \angle 2=\angle \overleftarrow{V_B}\&L_{BC}+\angle \overleftarrow{V_C}\&L_{BC}\hfill \\ \angle 3=\angle \overleftarrow{V_A}\&L_{AC}+\angle \overleftarrow{V_C}\&L_{AC}\hfill \end{array}$ (3) where $\angle \overleftarrow{V_A}\&L_{AB}$ is the acute angle of $\overleftarrow{V_A}$ and line LAB, as are the other angles. The minimum of three angles $\angle 1$, $\angle 2$, and $\angle 3$ is selected, and the corresponding two clusters are set to one muon track, and the line direction is selected as the muon direction.

Fig. 10Full-screen View

(Color online) An example of the three cluster charge patterns from a double-muon event and single-muon event. (a) An example of the three clusters' (A, B, and C) charge patterns by a double-muon event which have one edge muon track. The black dashed lines are the vectors (${\overrightarrow{V}}_A$, ${\overrightarrow{V}}_B$, and ${\overrightarrow{V}}_C$) for the three clusters. The red dashed lines correspond to the candidate tracks (LAB, LBC, and LAC). (b) An example of the three cluster charge patterns by a single-muon event whose distance to the CD center is about 14 m and near the edge of the CD.

A single-muon event may also have three clusters, which may confuse double-muon reconstruction. When a single-muon track is close to the edge of the LS sphere (the spatial distance between the muon track and LS center is from 13 to 16 m), the photon emissions on the muon track are all close to the PMTs and almost all PMTs are fired (Fig. 10). For this case, we still calculated the corresponding three angle values using Equation (3). The main feature is that all three angles are less than 5°. Therefore, we used this feature to distinguish between single and double muons. The midpoint (Fig. 10) is a false point, and the remaining two points are selected as incident and exit points.

Reconstruction performances of single and double muons

To obtain further exposure for the signal (neutrinos or others), the basic veto strategy is to apply a partial volume veto along the muon track for a period of time. The reconstructed and true track relative positions can be described using two parameters, Δ D and α. Here Δ D is defined as Δ D = D_t - D_r, where D_r and D_t are the distances between the reconstructed and true tracks from the center of the CD. We use α to denote the acute angle between two tracks (Fig. 11).

Fig. 11Full-screen View

(Color online) Schematic of the parameters Δ D and α. Here Δ D = D_t - D_r, where D_r and D_t are the distance of the reconstructed and true track from the CD center, and α is the acute angle between the two tracks.

The distributions of Δ D and α indicate the reconstruction performance. The single and double through-going-muon reconstruction performance at a distance of 16 m from the detector center is shown in Fig. 12. The distance to the center is divided into eight slices with a 2-m interval, which indicates a muon passing through different detector volumes. For a single muon, the average biases of Δ D in all slices (black squares in Fig. 12) were less than 5 cm, and the average biases of α (the black squares in Fig. 12) increase from 0.5° to 1.3° from the detector center to the edge. To reconstruct the double muons, as shown in Figs. 12c and d, the average biases of Δ D increased from 2 to 15 cm, and the average biases of α increased from 0.6° to 1.5°. The reconstruction resolution for each slice is shown in Fig. 13. In the internal volume (R ≤ 10 m), the resolution of Δ D and α is better than 20 cm and 0.5° for a single muon, and better than 30 cm and 0.6° for double muons. In the case of a muon passing through an external volume (10舁m ≤ R < 16 m), the resolution of Δ D and α increases to approximately 22 cm and 0.9° for a single muon and to approximately 35 cm and 1.2° for double muons, respectively. This is because the track length is short near the edge, and the corresponding clusters of muon incident and exit points may overlap; thus, there will be more deviations in the search for cluster seeds based on the charge-weight method. As mentioned in Sect. 2, a clipping muon (16舁m ≤ R < 17.7 m) creates a short track, and the two clusters easily overlap into one cluster; therefore, they are difficult to reconstruct into one track, so they cannot be evaluated by the parameters Δ D and α.

Fig. 12Full-screen View

(Color online) Reconstruction distribution of Δ D and α versus muon distance to the CD center for a single muon in (a) and (b) and for double muons in (c) and (d). Different distances to the center indicate a muon passing through different detector volumes. For double muons, Δ D, α, and the distances to the center are calculated according to their single track. In these figures, the color corresponds to the number of muons in each bin, and the black squares correspond to the average values of Δ D or α in each bin. The error bars of the black squares are very small and cannot be observed.

Fig. 13Full-screen View

(Color online) Resolutions of Δ D and α for single and double through-going muons in different detector volumes.

Reconstruction efficiency and muon bundle veto strategy discussion

As mentioned in Sect. 3, to match two candidate seeds to a muon track, several optimization methods were applied, especially in the case of more than two cluster seeds. As a result, most muon tracks can be reconstructed well. However, there are still some cases in which the muon track cannot be reconstructed (reconstructed track unavailable) or the reconstructed track is far from the true track (Δ D or α is larger than the five standard deviations of their distributions). For a single and double through-going muons, the percentages of these two cases are 1.6% and 1.0% and 1.7% and 3%, respectively. After reconstruction, reconstruction efficiency can be defined as the ratio of successful reconstruction to the total number of events. Figure 14a shows the reconstruction efficiency of single and double through-going muons. A single through-going muon with high reconstruction efficiency has an average efficiency of 98.3% ± 0.1%. The average reconstruction efficiency of the double through-going muons was 95.5% ± 0.2%. In addition, as shown in Fig. 14b, the reconstruction efficiencies of single and double through-going muons are stable when the PE cut varies from 70 PEs to 80 PEs.

A muon veto strategy can be developed based on the reconstructed muon track. A simple approach is to veto a cylindrical volume along the muon trajectory in a given time period [3]. Considering the position dependence of the reconstruction bias and resolution (Sect. 4.1), the radius of the cylinder can be modified by increasing the reconstruction bias and resolution [27]. Furthermore, a distance-dependent veto strategy was introduced for the detection of ⁸B solar neutrinos at JUNO [36], which investigates the time and distance distributions between the isotopes and muon tracks in detail and can significantly improve the signal-to-background-noise ratio. Considering the non-negligible contribution from the muon bundle and the higher estimated muon rate caused by the shallower overburden in JUNO, the veto strategy of JUNO must be optimized to obtain more exposure.

Using the algorithm in this study, we can reconstruct different types of muons and develop their corresponding veto strategies in the future. For a muon event with an available reconstructed track, a distance-dependent cylindrical veto cut along the reconstructed track is applied. The rejected detector volume and veto time can be determined according to the space and time distributions from the muon track to ⁹Li/⁸He [36]. If a double-muon event is reconstructed with two nearby reconstructed tracks (track distance d < 3 m), a larger cylindrical radius is used in the above strategy. For a clipping or stopping muon, their tracks are difficult to reconstruct because of the heavy smearing of clusters; thus, a spherical veto using the cluster seed as the center are used for conservative exclusion. For multiple-muon events, the cluster features of the charge pattern are more complicated, and it is difficult to match them to the correct muon tracks. As a result, the tracking precision of multiple muons was much lower than that of single and double muons, so it is not shown in this paper. In reconstruction, if a muon event is tagged as multiple muons, a whole-volume veto is applied. In addition, for failed reconstruction events whose track or cluster seeds are unavailable, we also veto a whole volume. However, a new method (called neutron veto) is under development that can use the corresponding neutron to precisely determine the position of muon-induced isotopes. A neutron veto uses a spherical veto on neutrons due to isotope and neutron space correlations. This is expected to save more dead volume and allow a longer veto time window for long-lived isotope rejection. On the other hand, the neutron veto is also expected to have a good ability to reject ⁹Li/⁸He from the shower muon. Thus, in addition to the muon veto, the neutron veto can be used to extend the veto strategy in the future.

Reconstruction speed performance

The time distribution of muon reconstruction is shown in Fig. 15, which was approximately 98 ms/event. The CPU (Intel Xeon E5-2660 v4) implementation of the reconstruction method in this study provides a significant processing speed improvement by a factor of 51 compared with the fastest light method (5000 ms/event) [28, 30]. This speed makes it possible to apply the algorithm in online event classification, which requires fast online filtering for possible data rate reduction.

Fig. 15Full-screen View

Reconstruction speed of a muon event.