Event vertex and time reconstruction in large-volume liquid scintillator detectors

Zi-Yuan Li· Yu-Mei Zhang· Guo-Fu Cao · Zi-Yan Deng · Gui-Hong Huang ·Wei-Dong Li · Tao Lin· Liang-Jian Wen · Miao Yu · Jia-Heng Zou ·Wu-Ming Luo · Zheng-Yun You

Abstract Large-volume liquid scintillator detectors with ultra-low background levels have been widely used to study neutrino physics and search for dark matter. Event vertex and event time are not only useful for event selection but also essential for the reconstruction of event energy. In this study, four event vertex and event time reconstruction algorithms using charge and time information collected by photomultiplier tubes were analyzed comprehensively. The effects of photomultiplier tube properties were also investigated. The results indicate that the transit time spread is the main effect degrading the vertex reconstruction, while the effect of dark noise is limited.In addition,when the event is close to the detector boundary, the charge information provides better performance for vertex reconstruction than the time information.

Keywords JUNO · Liquid scintillator detector · Neutrino experiment · Vertex reconstruction · Time reconstruction

1 Introduction

Liquid scintillators (LSs) have widely been used as detection medium for neutrinos in experiments such as Kamioka Liquid Scintillator Antineutrino Detector (Kam-LAND) [1], Borexino [2], Double Chooz [3], Daya Bay[4], and Reactor Experiment for Neutrino Oscillation(RENO) [5]. KamLAND revealed a large mixing angle(LMA) solution for solar neutrino oscillations. Borexino confirmed the Mikheyev–Smirnov–Wolfenstein (MSW)LMA [6] model in the sub-MeV region for solar neutrino oscillations.Double Chooz,Daya Bay,and RENO reported nonzero measurements for the mixing angle θ13. The size of such detectors varies from hundreds to thousands of cubic meters.Large-volume liquid scintillator detectors are widely used in the next generation of neutrino experiments,aiming to solve problems such as neutrinoless double-beta decay (SNO+ at the Sudbury Neutrino Observatory [7])and neutrino mass ordering (Jiangmen Underground Neutrino Observatory (JUNO) [8]).

The sensitivity of these experiments is limited by the energy resolution, detector volume, and detector background. These detectors typically contain a fiducial volume, where the signal-to-noise ratio is maximal. To distinguish between events occurring in the fiducial and non-fiducial regions,the event vertex is reconstructed using the charge and time distribution of photons collected by the photomultiplier tubes (PMTs). Most importantly, to achieve a high energy resolution, an accurate vertex is essential to correct for energy non-uniformity. In addition,the event vertex and event time information can also be used for particle identification, direction reconstruction,event classification, and other purposes.

A previous study [9] investigated vertex reconstruction with time information in JUNO, without discussing the event time reconstruction, dark noise effect, and the improvement based on the charge information when the event is close to the detector boundary. The aim of this study was to analyze vertex and time reconstruction for point-like events in JUNO under more realistic conditions.The main contributions of this study are as follows:

* The event time was reconstructed to provide the start time information of the event,which was important for event alignment, event correlation, etc.

* The dark noise from PMTs was considered and its effect on the vertex reconstruction was properly controlled.

* Two types of large PMTs were considered and handled separately mainly because of the difference in transit time spread (TTS).

* An algorithm to provide a more accurate initial vertex value was developed to improve the performance,especially at the detector boundary region.

* An algorithm employing charge information to reconstruct the event vertex at the detector boundary was developed.

The remainder of this paper is organized as follows. In Sect. 2, a brief introduction of the JUNO detector and the configurations of the PMTs used in this study is provided.In Sect. 3,the optical processes are described and a simple optical model is introduced. In Sect. 4, two simple algorithms that can quickly provide initial values are compared.In Sects. 5 and 6, two complex algorithms that provide relatively good performance are introduced in detail.Finally, a performance summary, the discussion, and the conclusions are provided in Sects. 7, 8,and 9,respectively.

2 The JUNO detector

Fig. 1 (Color online) Schematic of the JUNO detector

A schematic of the JUNO detector is shown in Fig. 1.The central detector (CD) is the main part of the JUNO detector assembly, with an acrylic sphere with a diameter of 35.4 m as inner layer.The acrylic sphere is supported by a stainless steel latticed shell with a diameter of 40.1 m and filled with approximately 20,000 tons of LS as target for neutrino detection. The composition of the LS includes linear alkylbenzene(LAB)as solvent,2,5-diphenyloxazole(PPO) as fluor, and p-bis-(o-methylstyryl)-benzene (bis-MSB) as wavelength shifter.To collect photons,the CD is surrounded by approximately 17,600 20-inch PMTs,including 5000 Hamamatsu R12860 PMTs, 12,600 North Night Vision Technology (NNVT) GDG-6201 PMTs, and approximately 25,600 3-inch XP72B22 PMTs. Around the CD, a 2.5-m-thick water pool is used to shield external radioactivity from the surrounding rocks and is combined with approximately 2,000 20-inch GDG-6201 PMTs to serve as a water Cherenkov detector to veto cosmic ray muons. A top tracker detector, consisting of plastic scintillators, located above the water pool is used for the identification and veto of muon tracks. A more detailed description of the JUNO detector can be found in Refs.[8, 11].

The main factors affecting the reconstruction of the event vertex and time include the TTS and dark noise of the PMTs. In this study,only the 20-inch PMTs of the CD were used for reconstruction. The number and parameters of the Hamamatsu and NNVT PMTs are summarized in Table 1 [12]. In principle, the 3-inch PMTs with σ ?1.6 ns TTS could also be included to improve the reconstruction performance; however, these PMTs [13]were not considered in this study because of their small photon detection coverage (～3%) [14, 15].

Table 1 Number and parameters of the PMTs used in the reconstruction. The TTS and dark noise rate are the mean values of the distribution measured during the mass testing.However,these are not the final values for JUNO

3 Optical processes

When a charged particle deposits energy in the scintillator, the solvent enters an excited state and transfers energy to the fluor in a non-radiative manner. Scintillation photons are then emitted along the particle track through the radiative de-excitation of the excited fluor within a limited time.The emitted scintillation photons can undergo several different processes while propagating through a large LS detector. At short wavelengths (＜410 nm),photons are mostly absorbed and then re-emitted at longer wavelengths, which maximizes the detection efficiency of the PMTs. At long wavelengths (＞410 nm), photons mainly undergo Rayleigh scattering.A more detailed study of the wavelength-dependent absorption and re-emission can be found in Ref. [16]. Additionally, the refractive indices at 420 nm are 1.50 and 1.34 for the LS and water,respectively.The difference in the refractive indices results in refraction and total reflection at the boundary of the two media, which affects the time-of-flight of the photons.When using the time information in the reconstruction,the time-of-flight is crucial, and it is calculated using the equation

where tof, dm, and vmare the time-of-flight, optical path length, and effective light speed, respectively, and m represents different media,in this case the LS and water,in the JUNO experiment.The acrylic sphere(thickness of 12 cm)and acrylic cover(thickness of 1 cm)in front of each PMT were ignored in this study because their refractive indices were similar to that of the LS and their thickness was small compared to that of the LS and water.

3.1 Optical path length

The optical path length can be characterized by the start and end positions in the detector, which are the vertex of the event and the position of the PMTs, respectively. The bold cyan curve in Fig. 2 shows a typical example of the optical path of photons detected by the PMT in the JUNO simulation using the event display [17, 18], and further examples are shown by the thin green curves. There are multiple physically possible paths between these two positions,each of which has a different optical path length,as follows:

* owing to absorption and re-emission, the re-emitted photon is not in the same absorption position and the propagation direction also changes;

Fig.2 (Color online)Event display of the optical path from the event vertex to the PMT in the JUNO simulation. The red circle ring is the event vertex and the gray bulbs with blue caps represent the PMTs

* owing to scattering, the photon changes the original direction of the propagation; and

* owing to refraction and total reflection,the photon does not travel in a straight line.

As shown in Fig. 2, owing to the various aforementioned optical processes, it is difficult to predict the actual optical path length for each photon.In this paper,a simple optical model is proposed, which uses a straight line connecting the vertex and the PMTs to calculate the optical path length (Fig. 3), and combines with the effective light speed to correct for the time-of-flight. Using this simple optical model reasonable results can be obtained, as discussed in Sect. 5.

In Fig. 3, {r→0，t0} represents the event vertex and start time, {r→i，ti} is the position of the ith PMT and the time of the earliest arriving photon detected by it. The angle between the normal direction of the ith PMT and the vector of the position of the ith PMT pointing to the event vertex is θiand αi=arccos（^r→0·^r→i）.The optical path length of the photon arriving at the ith PMT is dpathlength，i=|r→i-r→0|=dLS，i+dwater，iand the corresponding time-of-flight is tofi.The optical path length in the LS and water can be calculated by simply solving the trigonometric equation.

Fig.3 (Color online)Optical path length from the event vertex to the ith PMT. O denotes the center of the detector

3.2 Effective light speed

According to Ref. [16], the emission spectrum of scintillation photons is in the range of approximately 300–600 nm. Typically, the group velocity of the wave packet is used to describe the photon propagation in the medium,which is given by the equation

where vg is the group velocity, c is the speed of light in vacuum, n is the refractive index,and λ is the wavelength.

By fitting the Sellmeier equation [19], which describes the dispersion of the measurement in Refs. [20] and [21],the refractive index of the LS and water at different wavelengths is shown in the upper panel of Fig. 4. The group velocity of the LS and water can be calculated using Eq. 2 at different wavelengths,as shown in the lower panel of Fig. 4.

The propagation speed of photons in water(vwater)was determined as the average speed weighted by the probability density function of the photon wavelength, which was obtained from a Monte Carlo(MC)simulation.As the absorption and re-emission change the initial wavelength,determining the propagation speed of photons in the LS(vLS) is more complicated. To consider all wavelengthdependent effects that affect the propagation speed of photons, the effective light speed veff is introduced. In addition, veff also mitigates the effects by the simplified optical model,which,for example,ignores the refraction at the interface between the LS and water, as well as the change in the optical path length due to Rayleigh scattering. The exact value for veff can be determined using a data-driven method based on the calibration data as follows: place γ sources along the Z-axis, use vLS at 420 nm as the initial value of veff in the reconstruction algorithm and then, calibrate veff such that the source positions can be appropriately reconstructed. As no calibration data was available for JUNO, in this study, simulated calibration data were used, and the optimized values for the effective refractive index(c／veff)were 1.546 in the LS and 1.373 in water.In the future,the same method can be applied to the experimental calibration data.

4 Initial value for vertex and time

The TMinuit package [22] was used for the minimization procedure in the time likelihood and in the charge likelihood algorithm introduced in Sects. 5 and 6. When there are multiple local minima in the parameter space, an inaccurate initial value results in local instead of global minima, resulting in a lower reconstruction efficiency. For detectors such as JUNO, the initial value needs to be treated carefully because of the total reflection, as discussed in the following subsections.

4.1 Charge-based algorithm

The charge-based algorithm is essentially based on the charge-weighted average of the positions of the PMTs in an event, and the event vertex can be determined using the equation

where qiis the charge of the pulses detected by the ith PMT and r→0and r→iare defined in Fig. 3. A scale factor a is introduced because the charge-based algorithm is inherently biased and an ideal point-like event in a spherical detector is covered by a uniform photocathode. Even if all propagation-related effects, such as absorption and scattering are ignored, the result of a simple integral of the intersections of all photons with the sphere surface shows that the reconstructed position of the event is 2/3 of the true position. The value of a can be tuned based on the calibration data along the Z-axis. In this study, a=1.3 was used,which was sufficient to provide an initial estimate for the event vertex.

As can be seen in Fig. 5, even with the scale factor,owing to total reflection, the reconstructed vertex deviates up to 3 m near the detector boundary. According to Ref.[19], total reflection occurs only when the event vertex is located at an R larger than RLS·nwater／nLS ≈15.9 m,where RLS is the radius of the acrylic sphere, nLS and nwater are the refractive indices in the LS and water,respectively. The total reflection region is defined as R ＞15.9 m while R＜15.9 m is the central region. If the result from the charge-based algorithm is used as the initial value for the time likelihood algorithm, approximately 18% of events is reconstructed at a local minimum position. In addition, it should be noted that the charge-based algorithm is not able to provide an initial value for the event generation time t0. Therefore, a fast time-based algorithm needs to be introduced, which can provide more accurate initial values.

Fig. 5 (Color online) Heatmap of Rrec (upper panel) and Rrec-Rtrue (lower panel) as a function of Rtrue for 4-MeV e+ uniformly distributed in space reconstructed by the charge-based algorithm

4.2 Time-based algorithm

The time-based algorithm uses the distribution of the time-of-flight correction time Δt (defined in Eq. 4) of an event to reconstruct its vertex and t0. In practice, the algorithm finds the reconstructed vertex and t0using the following iterations:

1. Apply the charge-based algorithm to obtain the initial vertex.

2. Calculate time-of-flight correction time Δt for the ith PMT as

where j is the iteration step and ti, tofiare defined in Fig. 3. Plot the Δt distribution for all triggered PMTs,and label the peak position as Δtpeak.

where r→0,and r→iare defined in Fig. 3.To minimize the effect of scattering, reflection, and dark noise on the bias of the reconstructed vertex, only the pulses appearing in the （-10ns，+5ns） window around Δtpeak are included. The time cut also suppresses the effect of the late scintillation photons. The number of triggered PMTs in the window is Npeak.

The distribution of Δt at different iteration steps is shown in Fig. 6. At the beginning of the iteration, the Δt distribution is wide because the initial vertex is far from the true vertex. As the number of iterations increases, the Δt distribution becomes more concentrated. Finally, when the requirement in step 4 is satisfied,the iteration stops.In the final step, r→0is the reconstructed vertex and Δtpeak is the reconstructed time t0.

After the time-of-flight correction, the Δt distribution is independent of the event vertex. However, because the earliest arrival time is used, according to the first-order statistic, as discussed in Refs. [23–25], tiis related to the number of photoelectrons Nipe detected by ith PMT. To reduce the bias of the vertex reconstruction, the following form of the time–Npe correction is applied,and in Eq. 4 tiis replaced by t′i:

The parameters(p0, p1, p2)with the corresponding values of (9.42, 0.74, -4.60) for Hamamatsu PMTs and (41.31,-12.04,-20.02)for NNVT PMTs were found to minimize the bias and energy dependence of the reconstruction in this study. The difference in the parameters is mainly due to the difference in the TTSs of the PMTs. Following the correction, the times of different PMTs with different values of Npe are aligned.

Fig. 6 (Color online) Δt distribution at different iteration steps j

As shown in Fig. 7,the time-based algorithm provided a more accurate reconstructed vertex than the charge-based algorithm (Fig. 5). In addition, after the time–Npe correction, the reconstruction shows no obvious bias within the entire detector, even in the total reflection region. The reconstructed result was used as the initial value for the time likelihood algorithm.

5 Time likelihood algorithm

5.1 Principle of the algorithm

The time likelihood algorithm uses the scintillator response function to reconstruct the event vertex. The variable residual time tres（r→0，t0） for the ith PMT can be described as

where tires is the residual time of the ith PMT and r→0,t0,ti,and tofiare defined in Fig. 3.

The scintillator response function mainly consists of the emission time profile of the scintillation photons and the TTS and the dark noise of PMTs. In principle, the additional delays introduced by the absorption, re-emission,scattering,and total reflection of the photon arriving to the PMT depend on the distance between the emission position and the individual PMTs. However, the differences are only noticeable for the late arrival hits, which are largely suppressed by the requirement for the earliest arriving photons in the time likelihood algorithm. Therefore, in the first-order approximation,the scintillator response function can be considered to be the same for all positions inside the scintillator. The scintillator response function can be described as follows.

Fig. 7 (Color online) Heatmap of Rrec (upper panel) and Rrec-Rtrue (lower panel) as a function of Rtrue for 4-MeV e+ uniformly distributed in space reconstructed by the time-based algorithm

As described in Sect. 3, when a charged particle interacts with a scintillator molecule, the molecule is excited,then de-excites, and emits photons. Typically, the scintillator has more than one component;thus,the emission time profile of the scintillation photons,f（tres）,can be described as

where each k component is characterized by its decay time τkand intensity ρk. The different components result from the different excited states of the scintillator molecules.To consider the spread in the arrival time of photons at the PMTs, a convolution with a Gaussian function is applied, given by

where σ is the TTS of PMTs and ν is the average transit time.

The dark noise, which occurs without incident photons in the PMTs,is not correlated with any physical event.The fraction of the dark noise in the total number of photoelectrons εdn can be calculated based on the data acquisition(DAQ)windows,dark noise rate,and light yield of the LS. The probability of dark noise ε（tres） is constant over time,where∫DAQ ε（tres）dtres =εdn.By adding ε（tres）to g（tres） and renormalizing its integral to 1, the probability density function(PDF)of the scintillator response function can be written as

The distribution of the residual time tres of an event for a hypothetical vertex can be compared with p（tres）. The best-fitting vertex and t0are chosen by minimizing the negative log-likelihood function

The parameters in Eq. 10 can be measured experimentally[26–29].In this work,the PDF from the MC simulation for the methodology study was employed.

5.2 Probability density function

The PDF of the scintillator response function for PMTs detecting a single photoelectron was obtained from the MC simulation, using a 4.4-MeV γ source located at the center of the detector, such that the distance to all PMTs is the same. For PMTs detecting multiple photoelectrons, the time of the earliest arriving photon is biased toward an earlier time. Therefore, the PDF need to be modified according to the first-order statistic of p（tres） or the socalled first photoelectron timing technique [23–25] as

where pNpe（tres） is the PDF of the scintillator response function when the PMTs detect Npe hits.

The PDF of two kinds of PMTs is shown in Fig. 8: the upper panel is for Hamamatsu while the lower panel is for NNVT PMTs. As the PDF is affected by the time resolution of the PMTs, the PDF of the NNVT is wider because of its inferior TTS. The inset in the lower panel shows the PDF on a logarithmic scale, and the time constant contribution of the dark noise ε（tres） is clearly visible.

5.3 Reconstruction performance

Fig. 8 (Color online) PDF of the scintillator response function for PMTs detecting different numbers of photoelectrons.The upper panel shows the response function for Hamamatsu, the lower panel for the NNVT PMTs

The reconstructed vertex was compared with the true vertex in spherical coordinates (R，θ，φ) for the MC e+samples and fitted with a Gaussian function to analyze the bias and resolution.The bias of the reconstruction is shown in Fig. 9, where different colors represent events with different energies. As can be seen in the left panel of Fig. 9,the reconstructed R is consistent with the true value in the central region, while an energy-dependent bias behavior is noticeable near the detector boundary.Given its regular bias behavior, the bias can be corrected with an energy-dependent correction. Moreover, although the reconstructed R is biased, there is no bias in θ and φ, as shown in the middle and right panels of Fig. 9,respectively.

Owing to the low time resolution of the NNVT PMTs,in Fig. 10 reconstruction using only Hamamatsu PMTs is shown (green circles). In this study, we found that the vertex resolution with Hamamatsu PMTs was similar to that of using all PMTs. The reconstruction speed was 3.5 times faster, because the fraction of the Hamamatsu PMTs was approximately 28% of all PMTs in the CD.

The reconstructed event time t0is shown in Fig. 11.The effect of t0is essentially a global shift of an event to match the scintillator response function PDF; in reality, t0is also affected by the trigger time and the time delay from the cable. The absolute value of t0can be neglected; only the relative difference of different events is important for the alignment of events. The small bump near -1.6 ns is correlated with the R bias,and the long tail on the right side results from positronium formation. The variation in the reconstructed t0is within a few nanoseconds.

Fig.10 (Color online)Resolution of the reconstructed R(left panel),θ(middle panel),and φ(right panel)as a function of energy reconstructed by the time likelihood algorithm. Different colors represent different PMT configurations

Table 2 Vertex resolution for different PMT configurations at 1.022 MeV (detection of～1328 Npe in total,corresponding to ～370 Npe detected by Hamamatsu PMTs)

Table 3 Vertex resolution for different PMT configurations at 10.022 MeV (detection of～13280 Npe in total,corresponding to ～3700 Npe detected by Hamamatsu PMTs)

6 Total reflection region reconstructed by the charge likelihood algorithm

The time likelihood method described in Sec. 5 introduces a bias in the R direction when the reconstructing events are close to the acrylic sphere.As mentioned in Ref.[30], using a charge signal with the maximum likelihood method can provide better spatial resolution than the time likelihood algorithm when an event occurs near the detector boundary. In this section, we discuss the charge likelihood algorithm to reconstruct the event vertex in the total reflection region only, while the reconstruction result in the central region is omitted.

The charge likelihood algorithm is based on the distribution of the number of photoelectrons in each PMT.With the mean expected number of photoelectrons μ（r0→，E）detected by each PMT at a given vertex and energy, the probability of observing Npe on a PMT follows a Poisson distribution. Furthermore,

Fig. 11 (Color online) Reconstructed event time t0 at different energies

where Y is the energy scale factor, Ω（r0→，ri） is the solid angle of the ith PMT,εiis the detection efficiency of the ith PMT, f（θi） is the angular response of the ith PMT, θiis defined in Fig. 3, ζmis the attenuation length [31] in materials,and δiis the expected number of dark noise.This equation is based on the assumption that the scintillation light yield is linearly proportional to the energy.

However, Eq. 15 cannot describe properly the contribution of the indirect light, the effect of light shadows because of the geometric structure, and the effect of the total reflection. Another solution is to use the model-independent method described in Ref. [32]: the mean expected number of photoelectrons can be obtained by placing gamma sources at 29 specific positions along the Zaxis,which can be performed using a calibration procedure[33].In this study,different from Ref. [32],we focused on the performance of the vertex reconstruction. The mean expected number of the photoelectron distributions as a function of radius R and angle α is shown in Fig. 12, and the definition of angle α is shown in Fig. 3.

The mean expected number of photoelectrons μ obtained from Fig. 12 was used to calculate the hit probability.Instead of reconstructing(R，θ，φ)at the same time,θ and φ were fixed at the reconstructed values provided by the time likelihood algorithm, and only the event radius R was reconstructed using the charge likelihood algorithm.Therefore, the probability in Eq. 13 can be rewritten as

The reconstruction performance, focusing on the total reflection region, is shown in Figs. 13 and 14. In the total reflection region,the mean value of the reconstructed R was consistentwiththetrueR,andtheresolutionsintheRdirection were 81 mm at 1.022 MeV and 30 mm at 10.022 MeV.

As the charge distribution provides good radial discrimination ability, this algorithm can provide better resolution and a significantly smaller bias compared with those of the time likelihood algorithm in the total reflection region.

7 Performance summary

Fig. 12 (Color online) Mean expected number of photoelectron distribution as a function of radius R and angle α. This map is obtained by placing gamma sources at 29 specific positions along the Z-axis, which can be performed using a calibration procedure [32]

Fig. 13 (Color online) Bias of the reconstructed R in the total reflection region at different energies reconstructed by the charge likelihood algorithm

Fig.14 (Color online)Resolution of reconstructed R as a function of energy calculated by the time likelihood and charge likelihood algorithms (R3 ＞4000m3)

The execution time of the reconstruction for each event was tested on a computing cluster with Intel Xeon Gold 6238R CPUs (2.2 GHz), as shown in Fig. 15. The execution time of the charge-based algorithm was in the order of O（10-4） s per event,which cannot be presented in the figure. The execution time of the time-based and the time likelihood algorithm was proportional to the event energy and could be reduced by using only the Hamamatsu PMTs for the reconstruction. The execution time of the charge likelihood algorithm was independent of the event energy.

The resolutions of the four algorithms in the R direction are shown in Fig. 16.Owing to the large bias of the chargebased algorithm, a correction to remove the position-dependent bias was applied before the analysis of the resolution.

The charged-based algorithm is suitable for online reconstruction tasks that require high speed but do not require high resolution.The time-based algorithm does not rely on MC; it can be used as a data-driven reconstruction method. The time likelihood and charge likelihood algorithms are relatively accurate and each has its own advantages in a specific detector region.

Fig. 15 (Color online) Execution time for the reconstruction for different algorithms

Fig.16 (Color online)Resolution of the reconstructed R as a function of energy for different algorithms

8 Discussion

The vertex resolutions of KamLAND and Borexino are approximately 12 cm and 10 cm at 1 MeV, respectively,and for JUNO it is approximately 10.5 cm.The diameter of JUNO (35.4 m) is several times larger than that of Kam-LAND (13 m) and Borexino (8.5 m). Despite its larger size, JUNO is still able to achieve a similar vertex resolution based only on the PMT time information. In this study, various effects on the vertex reconstruction for JUNO were comprehensively analyzed. As expected, the TTS of the PMT is the dominant factor. The vertex reconstruction capability of JUNO is mainly based on the Hamamatsu PMTs. Although the number of NNVT PMTs is more than twice, their time information is not useful in the vertex reconstruction because of their significantly inferior TTS. After considering the light yield of the LS and the PMT coverage, the number of photons detected by the Hamamatsu PMTs is in the same range as those of the three detectors mentioned above. This provides an explanation of their similar vertex resolutions based only on the PMT time information. To fully exploit the large PMT coverage and large number of PMTs in JUNO, the charge information of PMTs also need to be utilized in addition to the time information to constrain the event vertex, especially near the detector boundary. At the same time, the effect of the dark noise of the PMT can be mitigated with appropriate treatment. A more accurate initial value also improves the performance of the vertex reconstruction. In addition to the event vertex, the event time can also be reconstructed simultaneously,which is a useful variable for downstream analyses.

In LS detectors, in addition to the scintillation photons there are also Cherenkov photons,whose effects need to be studied in the future. The R bias near the detector edge in the current results also indicates that a more accurate PDF of the scintillator response function is needed to include its dependence on the position as well as the particle type.All vertex reconstruction methods based on PMT time information use the time of the first arrival photons only; in principle,later photons might be useful as well.Therefore,novel methods also need to be explored. Our preliminary studies on algorithms based on machine learning [34]showed comparable vertex reconstruction performance,and they need to be further investigated and optimized.

9 Conclusion

In this study, four algorithms for the reconstruction of the event vertex and event time were investigated in detail and verified using MC samples generated by the offline software of JUNO. Considering the TTS and dark noise effects from the PMTs, a vertex resolution of 10 cm or higher can be achieved in the energy range of reactor neutrinos. The TTS has a major effect on the vertex resolution,whereas the effect of dark noise is limited.Near the detector boundary, charge information can constrain the event vertex better compared with time information. The algorithms discussed in this paper are applicable to current and future experiments using similar detection techniques.

Author ContributionsAll authors contributed to the study conception and design. Material preparation, data collection, and analysis were performed by Zi-Yuan Li,Yu-Mei Zhang,Guo-Fu Cao,Zi-Yan Deng, Gui-Hong Huang, Wei-Dong Li, Tao Lin, Liang-Jian Wen,Miao Yu, Jia-Heng Zou, Wu-Ming Luo, and Zheng-Yun You. The first draft of the manuscript was written by Zi-Yuan Li, Yu-Mei Zhang, Wu-Ming Luo, and Zheng-Yun You, and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.