Event vertex and time reconstruction in large volume liquid scintillator detectors

NUCLEAR PHYSICS AND INTERDISCIPLINARY RESEARCH

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

Nuclear Science and Techniques

Vol.32, No.5

Article number 49

Published in print 01 May 2021

Available online 12 May 2021

DOI：10.1007/s41365-021-00885-z

46200

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.

JUNOLiquid scintillator detectorNeutrino experimentVertex reconstructionTime reconstruction

1 Introduction

Liquid scintillators (LSs) have widely been used as detection medium for neutrinos in experiments such as Kamioka Liquid Scintillator Antineutrino Detector (KamLAND) [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-–Wolfenstei (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 θ₁₃. 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

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 Ref. [8, 11].

Fig. 1

(Color online) Schematic of the JUNO detector.

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].

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.

Company	Number	TTS (σ)	Dark noise rate
Hamamatsu	5,000	1.15 ns	15 kHz
NNVT	12,600	7.65 ns	32 kHz

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

t o f = \sum_{m} \frac{d_{m}}{v_{m}},

(1)

where tof, dm, and vm are 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:

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 absorption and re-emission, the re-emitted photon is not in the same absorption position and the propagation direction also changes;

* 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.

Fig. 3

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

In Fig. 3, { ${\vec{r}}_{0}, t_{0}$ } represents the event vertex and start time, { ${\vec{r}}_{i}, t_{i}$ } 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 θi and $α_{i} = \arccos ({\hat{\vec{r}}}_{0} \cdot {\hat{\vec{r}}}_{i})$ . The optical path length of the photon arriving at the ith PMT is $d_{pathlength, i} = | {\vec{r}}_{i} - {\vec{r}}_{0} | = d_{LS, i} + d_{water, i}$ and 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.

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

v_{g} (λ) = \frac{c}{n (λ) - λ \frac{\partial n (λ)}{\partial λ}},

(2)

where v_g 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.

Fig. 4

(Color online) Dependence of the refractive index (upper panel) and group velocity (lower panel) on the wavelength in the LS and water.

The propagation speed of photons in water (v_water) 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 (v_LS) is more complicated. To consider all wavelength-dependent effects that affect the propagation speed of photons, the effective light speed v_eff is introduced. In addition, v_eff 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 v_eff can be determined using a data-driven method based on the calibration data as follows: place γ sources along the Z-axis, use v_LS at 420 nm as the initial value of v_eff in the reconstruction algorithm and then, calibrate v_eff 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/v_eff) 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

{\vec{r}}_{0} = a \cdot \frac{\sum_{i} q_{i} \cdot {\vec{r}}_{i}}{\sum_{i} q_{i}},

(3)

where qi is the charge of the pulses detected by the ith PMT and ${\vec{r}}_{0}$ and ${\vec{r}}_{i}$ are 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 R_LS⋅n_water/n_LS≈15.9 m, where R_LS is the radius of the acrylic sphere, n_LS and n_water 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 t₀. Therefore, a fast time-based algorithm needs to be introduced, which can provides more accurate initial values.

Fig. 5

(Color online) Heatmap of R_rec (upper panel) and R_rec-R_true (lower panel) as a function of R_true 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 t₀. In practice, the algorithm finds the reconstructed vertex and t₀ using 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

Δ t_{i} (j) = t_{i} - t o f_{i} (j),

(4)

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

3. Calculate the correction vector $[\vec{r} (j)]$ as

\vec{δ} [\vec{r} (j)] = \frac{\sum_{i} (\frac{Δ t_{i} (j) - Δ t^{peak} (j)}{t o f_{i} (j)}) \cdot ({\vec{r}}_{0} (j) - {\vec{r}}_{i})}{N^{peak} (j)},

(5)

where ${\vec{r}}_{0}$ , and ${\vec{r}}_{i}$ are 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 (-10 ns, +5 ns) window around Δt^peak are included. The time cut also suppresses the effect of the late scintillation photons. The number of triggered PMTs in the window is N^peak.

4. If $\vec{δ} [\vec{r} (j)] < 1 mm$ or j=100, stop the iteration; otherwise, update the vertex with ${\vec{r}}_{0} (j + 1) = {\vec{r}}_{0} (j) + \vec{δ} [\vec{r} (j)]$ and go to step 2 to start a new round of iteration.

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, ${\vec{r}}_{0}$ is the reconstructed vertex and Δt^peak is the reconstructed time t₀.

Fig. 6

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

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 Ref. [23-25], ti is related to the number of photoelectrons $N_{pe}^{i}$ detected by ith PMT. To reduce the bias of the vertex reconstruction, the following form of the time–N_pe correction is applied, and in Eq. 4 ti is replaced by ${t^{'}}_{i}$ :

{t^{'}}_{i} = t_{i} - p 0 / \sqrt{N_{pe}^{i}} - p 1 - p 2 / N_{pe}^{i} .

(6)

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 N_pe are aligned.

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–N_pe 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.

Fig. 7

(Color online) Heatmap of R_rec (upper panel) and R_rec-R_true (lower panel) as a function of R_true for 4-MeV e⁺ uniformly distributed in space reconstructed by the time-based 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 $t_{res} ({\vec{r}}_{0}, t_{0})$ for the ith PMT can be described as

t_{res}^{i} ({\vec{r}}_{0}, t_{0}) = t_{i} - t o f_{i} - t_{0},

(7)

where $t_{res}^{i}$ is the residual time of the ith PMT and ${\vec{r}}_{0}$ , t₀, ti, and tofi are 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.

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(t_res), can be described as

f (t_{res}) = \sum_{k} \frac{ρ_{k}}{τ_{k}} e^{\frac{- t_{res}}{τ_{k}}}, \sum_{k} ρ_{k} = 1,

(8)

where each k component is characterized by its decay time τk and 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

g (t_{res}) = \frac{1}{\sqrt{2 π} σ} e^{- \frac{{(t_{res} - ν)}^{2}}{2 σ^{2}}} \cdot f (t_{res}) .

(9)

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 ε(t_res) is constant over time, where ∫_DAQε(t_res)dt_res=ε_dn. By adding ε(t_res) to g(t_res) and renormalizing its integral to 1, the probability density function (PDF) of the scintillator response function can be written as

p (t_{res}) = (1 - ε_{dn}) \cdot g (t_{res}) + ε (t_{res}) .

(10)

The distribution of the residual time t_res of an event for a hypothetical vertex can be compared with p(t_res). The best-fitting vertex and t₀ are chosen by minimizing the negative log-likelihood function

L ({\vec{r}}_{0}, t_{0}) = - \ln (\prod_{i} p (t_{res}^{i})) .

(11)

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(t_res) or the so-called first photoelectron timing technique [23-25] as

p_{N_{pe}} (t_{res}) = N_{pe} p (t_{res}) {(\int_{t_{res}}^{\infty} p (x) d x)}^{N_{pe} - 1},

(12)

where pNpe(t_res) is the PDF of the scintillator response function when the PMTs detect N_pe 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 ε(t_res) is clearly visible.

Fig. 8

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

5.3 Reconstruction performance

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.

Fig. 9

(Color online) Bias of the reconstructed R (left panel), θ (middle panel), and ϕ (right panel) for different energies calculated by the time likelihood algorithm.

The spatial resolution of the vertex reconstruction as a function of energy is shown in Fig. 10. The R bias was corrected before the analysis of the resolution. To study the individual effect of the TTS and dark noise on the vertex reconstruction, different MC samples were produced with and without these effects. The vertex reconstruction results are shown in Fig. 10. The magenta circles represent the default PMT configuration, as described in Sect. 2. The red triangles represent an ideal configuration, which assumes perfect PMTs without the effects of the TTS and dark noise. The black squares represent the configuration of PMTs including only the dark noise effect, while the blue inverted triangles represent the PMT configuration including only the TTS effect. The exact values of the vertex resolution at 1.022 MeV and 10.022 MeV are summarized in Tables 2 and 3, respectively. The energy E_true includes the energy of the annihilation gamma rays. The light yield was approximately 1300 detected Npe per 1 MeV of deposited energy in JUNO, and the energy nonlinearities on the light yield were ignored in the approximation. As can be seen in Tables 2 and 3, the dark noise has no effect at high energy and its effect at low energy is also highly limited. The largest effect results from the TTS in the time likelihood algorithm. The energy-dependent vertex resolution is approximately proportional to $1 / \sqrt{N_{pe}}$ [24].

Vertex resolution for different PMTs configurations at 1.022 MeV (detection of ∼1328 N_pe in total, corresponding to ∼370 N_pe detected by Hamamatsu PMTs)

PMT configuration	R (mm)	θ (degrees)	ϕ (degrees)
Ideal	60	0.25	0.31
With dark noise only	62	0.27	0.34
With TTS only	89	0.37	0.44
With TTS and dark noise	103	0.40	0.47
With TTS and dark noise (Hamamatsu PMTs only)	105	0.42	0.49

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

PMT configuration	R (mm)	θ (degrees)	ϕ (degrees)
Ideal	19	0.08	0.11
With dark noise only	19	0.08	0.11
With TTS only	31	0.13	0.16
With TTS and dark noise	31	0.13	0.16
With TTS and dark noise (Hamamatsu PMTs only)	32	0.14	0.17

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.

Owing to the low time resolution of the NNVT PMTs, in Fig. 10 reconstuction 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 t₀ is shown in Fig. 11. The effect of t₀ is essentially a global shift of an event to match the scintillator response function PDF; in reality, t₀ is also affected by the trigger time and the time delay from the cable. The absolute value of t₀ can 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 t₀ is within a few nanoseconds.

Fig. 11

(Color online) Reconstructed event time t₀ at different energies.

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 $μ ({\vec{r}}_{0}, E)$ detected by each PMT at a given vertex and energy, the probability of observing N_pe on a PMT follows a Poisson distribution. Furthermore,

* Probability for the jth PMT with no hits: $P_{nohit}^{j} ({\vec{r}}_{0}, E) = e^{- μ_{j}}$ ,

* Probability for the ith PMT with $N_{pe}^{i}$ hits: $P_{hit}^{i} ({\vec{r}}_{0}, E) = \frac{μ_{i}^{N_{pe}^{i}} e^{- μ_{i}}}{N_{pe}^{i}!}$ .

Therefore, the probability of observing a hit pattern for an event can be written as

p ({\vec{r}}_{0}, E) = \prod_{j} P_{nohit}^{j} ({\vec{r}}_{0}, E) \cdot \prod_{i} P_{hit}^{i} ({\vec{r}}_{0}, E) .

(13)

The best-fit values of ${\vec{r}}_{0}$ and E can be obtained by minimizing the negative log-likelihood function

L ({\vec{r}}_{0}, E) = - \ln (p ({\vec{r}}_{0}, E)) .

(14)

In principle, $μ ({\vec{r}}_{0}, E)$ can be expressed by the equation

μ_{i} ({\vec{r}}_{0}, E) = Y \cdot \frac{Ω ({\vec{r}}_{0}, r_{i})}{4 π} \cdot ε_{i} \cdot f (θ_{i}) \cdot e^{- \sum_{m} \frac{d_{m}}{ζ_{m}}} \cdot E + δ_{i},

(15)

where Y is the energy scale factor, $Ω ({\vec{r}}_{0}, r_{i})$ is the solid angle of the ith PMT, εi is the detection efficiency of the ith PMT, f(θi) is the angular response of the ith PMT, θi is defined in Fig. 3, ζm is the attenuation length [31] in materials, and δi is 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 Z-axis, 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.

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].

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

p (R, E) = \prod_{j} P_{nohit}^{j} (R, E) \cdot \prod_{i} P_{hit}^{i} (R, E) .

(16)

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 consistent with the true R, and the resolutions in the R direction were 81 mm at 1.022 MeV and 30 mm at 10.022 MeV.

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 (R³>4000m³)

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

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.

Fig. 15

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

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

Fig. 16

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

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.

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 KamLAND (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 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 Conclusions

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.

References

[1]

KamLAND Collaboration,

First Results from KamLAND: Evidence for Reactor Antineutrino Disappearance

. Phys. Rev. Lett. 90, 021802 (2003). doi: 10.1103/PhysRevLett.90.021802