# Reconstruction of Point Events in Liquid-Scintillator Detectors Subjected to
Total Reflection

Wei Dou(窦威)  Department of Engineering Physics, Tsinghua University, Beijing,
China  Center for High Energy Physics, Tsinghua University, Beijing, China
Key Laboratory of Particle & Radiation Imaging (Tsinghua University), Ministry
of Education, Beijing, China  Benda Xu(续本达)111Corresponding author.
<EMAIL_ADDRESS>Department of Engineering Physics, Tsinghua University,
Beijing, China  Center for High Energy Physics, Tsinghua University, Beijing,
China  Key Laboratory of Particle & Radiation Imaging (Tsinghua University),
Ministry of Education, Beijing, China  Kavli Institute for the Physics and
Mathematics of the Universe, UTIAS, the University of Tokyo, Tokyo, Japan
Jianfeng Zhou(周建锋)  Xingfan Information Technology Co., Ningbo, China  Zhe
Wang(王 吉 吉 )  Department of Engineering Physics, Tsinghua University, Beijing,
China  Center for High Energy Physics, Tsinghua University, Beijing, China
Key Laboratory of Particle & Radiation Imaging (Tsinghua University), Ministry
of Education, Beijing, China  Shaomin Chen(陈少敏)  Department of Engineering
Physics, Tsinghua University, Beijing, China  Center for High Energy Physics,
Tsinghua University, Beijing, China  Key Laboratory of Particle & Radiation
Imaging (Tsinghua University), Ministry of Education, Beijing, China

###### Abstract

The outer water buffer is an economic option to shield the external radiative
backgrounds for liquid-scintillator neutrino detectors. However, the
consequential total reflection of scintillation light at the media boundary
introduces extra complexity to the detector optics. This paper develops a
precise detector-response model by investigating how total reflection
complicates photon propagation and degrades reconstruction. We first
parameterize the detector response by regression, providing an unbiased energy
and vertex reconstruction in the total reflection region while keeping the
number of parameters under control. From the experience of event degeneracy at
the Jinping prototype, we then identify the root cause as the multimodality in
the reconstruction likelihood function, determined by the refractive index of
the buffer, detector scale and PMT coverage. To avoid multimodality, we
propose a straightforward criterion based on the expected photo-electron-count
ratios between neighboring PMTs. The criterion will be used to ensure success
in future liquid-scintillator detectors by guaranteeing the effectiveness of
event reconstruction.

Keywords: event reconstruction, liquid scintillator, spherical harmonics,
total reflection

##  1 Introduction

Liquid-scintillator (LS) detectors, such as Borexino [1], KamLAND [2], SNO+
[3] and JUNO [4], obtain lower energy thresholds and better energy resolutions
than water Cherenkov detectors by enhanced light yields. It contains rich
contemporary physics topics, including neutrino mass ordering [5, 6],
neutrinoless double beta decay [7, 8, 9] as well as terrestrial, solar and
supernova neutrinos [10].

The crucial step of energy reconstruction in LS detectors is to predict the
_photon electron_ (PE) count and timing on each _photomultiplier tube_ (PMT)
for a given vertex. Vertex reconstruction thus influences energy resolution.
As illustrated in figure 1, if an LS detector is surrounded by a water buffer,
each PMT can be blind to photons in a specific region due to _total
reflection_ (TR), causing a great challenge to vertex reconstruction.

Figure 1: The particle deposits energy at position $\mathbf{r}$ on time
$t_{0}$. The red lines are two typical photon tracks. Photons can reach the
blue PMT directly but not the grey one due to TR. $r_{\mathrm{LS}}$ is the
radius of the LS. The $i$-th PMT is located at position
$\mathbf{r}_{\mathrm{PMT},i}$ with
$|\mathbf{r}_{\mathrm{PMT},i}|=r_{\mathrm{PMT}}$. A function of
$(\mathbf{r},\mathbf{r}_{\mathrm{PMT},i})$ can be simplified to that of
$(r,\theta)$ by spherical symmetry, where $r=|\mathbf{r}|$ and $\theta$ is the
central angle defined by the vertex and the PMT.

The prevalent event reconstruction methods fall into 3 groups. _Barycenter_
(BC) averages all PMT positions weighted by PE. It is biased but fast [11, 12,
13], often employed as the initial values to more advanced algorithms.
_Maximum likelihood estimation_ (MLE) uses predicted PE [14] or timing [11,
15, 16] or both [17, 18]. Its theoretical uncertainty is discussed by C.
Galbiati et al [17]. The vertex reconstruction results are biased at the TR
region if only use timing information [13]. Li [13] and Huang [19] predict PE
using interpolation with a map to achieve unbiased results but introduce a
high degree of freedom. Machine learning methods [20] such as neural networks
and decision trees have good modeling power to describe TR optics but need
high-quality labeled training datasets only available from detector
simulations.

In this paper, we develop an improved MLE-based event reconstruction resilient
to TR. Section 2 formulates a detector-response model and tests it with Monte
Carlo (MC) simulation data. Section 3 analyzes the event degeneracy caused by
the multimodality in the likelihood function. Section 4 shows the
reconstruction performance. Finally, section 5 discusses the extensibility and
limitation of our method. The symbol conventions for this paper are in table
1.

Table 1: Definitions of symbols variable | meaning (r.v. for random variable) | first appearance in section
---|---|---
$E,t_{0}$ | visible energy and event start time | 2 Detector response by regression
$\mathbf{r},\mathbf{r}_{\mathrm{PMT}}$ | vertex and PMT’s position
$i,j$ | PMT index and hit index
$r_{\mathrm{LS}},r_{\mathrm{PMT}}$ | radii of the LS and PMT’s position
$r=|\mathbf{r}|$ | vertex radius | 2.1 PE and timing prediction
$\theta$ | central angle defined by PMT and vertex
$\hat{E},\hat{\mathbf{r}}$ | reconstructed energy and position |  4.1 Results from simulations
$\theta_{\mathrm{acr}}$ | central angle defined by vertex and |  5.1 Other reflections by acrylic shell
| point of incidence to the acrylic shell |
$n_{i}$ | observed PE on $i$th PMT | 2 Detector response by regression
$t_{ij}$ | observed $j$th hit’s timing on $i$th PMT
$q_{ij}$ | reconstructed $j$th hit’s charge on $i$th PMT | 4.2 Results from raw data
$N_{\mathrm{PMT}}$ | number of PMTs |  3.1 Cosine distance
$N_{\mathrm{hit}}$ | total hit number |  4.1 Results from simulations
$\lambda_{i},\lambda_{i}(r,\theta)$ | predicted PE on $i$th PMT | 2 Detector response by regression
$T_{i},T_{i}(r,\theta)$ | predicted timing on $i$th PMT
$c^{\lambda}_{l},c^{\lambda}_{l}(r)$ |  $l$th PE coefficient | 2.1 PE and timing prediction
$c^{T}_{l},c^{T}_{l}(r)$ |  $l$th timing coefficient
$c^{\lambda}_{l,m},c^{T}_{l,m}$ | PE and timing coefficient
$\mu_{i},U_{i}$ | PE and timing PMT-specific offset
$P_{l}$ | Legendre polynomial
$\psi(t)$ | scintillation time profile |  2 Detector response by regression
$R^{0}(t)$ | timing PDF by $\psi(t)$ and TTS, etc. |
$\mathscr{L}$ | likelihood function | 2.1 PE and timing prediction
$\mathscr{R}_{\tau},\mathscr{R}_{\tau}(t)$ | loss function using $\tau$-quantile
$R(t)$ | quantile regression approximation to $R^{0}(t)$
$\tau,t_{\mathrm{s}}$ | quantile value and time scale
$\tau_{r},\tau_{d}$ | rise and decay time constant |  2.2 Training and validation
$D_{\mathrm{KL}}$ | Kullback-Leibler divergence |  4.2 Results from raw data
$\Lambda$ | PE pattern | 3.1 Cosine distance
$d_{\cos}$ | cosine distance
$C_{1r}(\mathbf{r})$ | contour function | 3.2 Cosine distance for 3-PMT case
$\Omega$ | solid angle
$\beta$ | angle of incidence on a PMT

##  2 Detector response by regression

This section develops a model to predict the PE and timing in an LS detector
that is suitable for TR. The radius of the LS is $r_{\mathrm{LS}}$. Figure 1
shows how a typical LS detector works. An ionizing particle begins to deposit
energy $E$ at the position $\mathbf{r}$ on time $t_{0}$. It produces
scintillation photons, each obeying a scintillation _time profile_
$\mathrm{\psi}(t)$. PMTs are located at position $\mathbf{r}_{\mathrm{PMT}}$
with $|\mathbf{r}_{\mathrm{PMT}}|=r_{\mathrm{PMT}}$. A photon travels to a PMT
with _time of flight_ (TOF) and induces a PE with the probability of _quantum
efficiency_ (QE). The PE gets amplified by a series of dynodes inside the PMT
with _transit time_ (TT), whose _spread_ is TTS. The whole process is called
an _event_.

The event is characterized by PE counts $n_{i}$ and timings $t_{ij}$ extracted
from the waveforms on the $i$th PMT, where $j\in\mathbb{N}$ and $j<n_{i}$.
Throughout this paper, we use PE to be the short for “PE counts”. _Detector
response model_ predicts the PE $\lambda_{i}$ and timing $T_{i}$.
$\lambda_{i}$ is the average of the Poissonian $n_{i}$. $T_{i}$ is a shift to
the $t_{ij}$’s _probability distribution function_ (PDF) $R^{0}(t-T_{i})$,
which is determined by $t_{0}$, TOF, TT and $\psi(t)$.

###  2.1 PE and timing prediction

$\lambda_{i}$ and $T_{i}$ are the functions of
$(\mathbf{r},\mathbf{r}_{\mathrm{PMT}})$. If the detector is spherically
symmetric, the function can be rewritten in $(r,\theta)$. Here $r$ is the
vertex radius and $\theta$ is the central angle defined by PMT and vertex. By
_varying coefficient model_ [21], we fit the conditional distribution of
$\theta$ given $r$, then fit the coefficients with $r$ by the method of least
squares.

In the first step, $n_{i}$ is modeled by Poisson regression [22] as count
data. $\lambda_{i}$ is connected to the observation $n_{i}$ by the Poisson
log-likelihood

$\displaystyle\log\mathscr{L}=-\sum_{i}\lambda_{i}+\sum_{i}n_{i}\log\lambda_{i}+\mathrm{Const.}$
(1)

A logarithm, called a _link function_ in generalized linear model [22]
terminology, connects $\lambda_{i}$ with the predictor variable as a linear
combination of Legendre polynomials $P_{l}(\cos\theta)$. $l$ is the order and
$c^{\lambda}_{l}(r)$’s are the regression coefficients.

$\lambda_{i}(r,\theta_{i})=\mathbb{E}\left[n_{i}\mid
r,\theta_{i}\right]=\exp\left\\{\sum_{l}c^{\lambda}_{l}(r)P_{l}(\cos\theta_{i})+\log
E+\mu_{i}\right\\}$ (2)

Here $E$ is the energy, $\mu_{i}$ is the correction term to match the PMT-
specific differences dominated by QE. We restrict $\sum_{i}\mu_{i}=0$, $E$ and
$\mu_{i}$ enter the regression as offsets. Eqs. (1) and (2) models the
dependence on $\theta_{i}$ and outputs a set of $r$-dependent coefficients
$c^{\lambda}_{l}(r)$.

Similarly, quantile regression [23] is used to model timing, reducing the
influence of the heavy tail of $R^{0}(t-T_{i})$. The predicted $T_{i}$ for the
$\tau$-quantile is

$T_{i}(r,\theta_{i})=\arg\underset{c}{\min}\left[\sum_{j}\mathscr{R}_{\tau}(t_{ij}-c)\right]=\sum_{l}c^{T}_{l}(r)P_{l}(\cos\theta_{i})+U_{i}+t_{0},$
(3)

where $\mathscr{R}_{\tau}(t)$ is the loss function defined by

$\mathscr{R}_{\tau}(t)=\left\\{\begin{array}[]{ll}-t(1-\tau)&\mbox{, if
$t<0$;}\\\ t\tau&\mbox{, if $t\geq 0$.}\end{array}\right.$ (4)

$t_{0}$ is the time offset and $U_{i}$ is the PMT-specific offset dominated by
TT with $\sum_{i}U_{i}=0$. Like $c^{\lambda}_{l}(r)$, $c^{T}_{l}(r)$ encodes
the $\theta$ dependence of timing.

Minimizing $\mathscr{R}_{\tau}(t)$ in eqs. (3) and (4) is equivalent to
maximizing a likelihood function
$\mathscr{L}\propto\exp(-\mathscr{R}_{\tau}/t_{\mathrm{s}})$, where
$t_{\mathrm{s}}$ is an arbitrary positive real number encoding a _time scale_.
The normalizing constant is

$\displaystyle\int_{-\infty}^{\infty}\exp(-\frac{\mathscr{R}_{\tau}(t)}{t_{\mathrm{s}}})\mathrm{d}t=\frac{t_{\mathrm{s}}}{1-\tau}+\frac{t_{\mathrm{s}}}{\tau}=\frac{t_{\mathrm{s}}}{\tau(1-\tau)}$
(5)

Therefore, quantile regression approximates the timing PDF $R^{0}(t)$ by

$R(t)=\frac{\tau(1-\tau)}{t_{\mathrm{s}}}\exp(-\frac{\mathscr{R}_{\tau}(t)}{t_{\mathrm{s}}}).$
(6)

The examples of the loss function $\mathscr{R}_{\tau}(t)$ and timing PDF
$R(t)$ are shown in figure 2.

(a) The shape of $\mathscr{R}_{\tau}$(top) and $R$(bottom)

(b) The CDF of $R(t-T_{i})$

Figure 2: (2) shows the $\mathscr{R}_{\tau}(t-T_{i})$ (top) and $R(t-T_{i})$
(bottom). $T_{i}$ is $27.5\text{\,}\mathrm{n}\mathrm{s}$ and $\tau$ is 0.1.
“MC” represents $R^{0}(t)$. It is the histogram of simulated event timings at
the detector center, see section 2.2. The best $t_{\mathrm{s}}$ is
$3.0\text{\,}\mathrm{n}\mathrm{s}$. (2) is the _cululative distribution
function_ (CDF) of $R^{0}(t)$ as “MC” and $R(t-T_{i})$ as “Fit”. The steepest
slope is at the peak of $R^{0}(t)$, leading to $\tau\approx 0.1$.

In the second step, we use another set of Legendre polynomials to fit
coefficients $c^{\lambda}_{l}(r)$ and $c^{T}_{l}(r)$. $r$ is scaled to
$[-1,1]$ by dividing the LS radius $r_{\mathrm{LS}}$. Due to the symmetry, we
only use the even orders to guarantee the $c^{\lambda}_{l}(r)$ and
$c^{T}_{l}(r)$’s derivatives are 0 at the detector center,

$c^{\lambda}_{l}(r)=\sum_{m=0}^{\infty}{c^{\lambda}_{l,2m}}P_{2m}(r/r_{\mathrm{LS}}),\quad\quad
c^{T}_{l}(r)=\sum_{m=0}^{\infty}{c^{T}_{l,2m}}P_{2m}(r/r_{\mathrm{LS}}).$ (7)

The above varying coefficient model requires simulation or calibration at
fixed radii. Alternatively, if the simulated events are uniformly distributed
in the detector, the above 2-step requirement can be relaxed with a one-step
regression,

$\lambda_{i}(r,\theta_{i})=\exp\left[\sum_{m}c^{\lambda}_{m}F_{m}(r,\cos\theta_{i})\right],\quad\quad
T_{i}(r,\theta_{i})=\sum_{m}c^{T}_{m}F_{m}(r,\cos\theta_{i}).$ (8)

$F_{m}(r,\cos\theta)$ is a binary basis function of $r$ and $\theta$ at the
$m$th order such as Zernike polynomials [24]. Another way to construct a
binary basis function is to product the two Legendre polynomials from the
varying coefficient model, called _double Legendre_.

$F_{l,m}(r,\cos\theta)=P_{m}(r/r_{\mathrm{LS}})\times P_{l}(\cos\theta)$ (9)

which is indexed by two subscripts $l$ and $m$. Due to the memory constraints
of our computing system, a regression cannot handle more than 800 parameters
in one pass. Consequently binary basis models are more restricted than the
varying coefficient one, although the former models more symmetrically handle
$r$ and $\theta$. Their best order selections are discussed in section 2.2.

###  2.2 Training and validation

The training dataset is generated by a custom piece of software, _Jinping
Simulation and Analysis Package_ (JSAP), based on GEANT4 [25]. Two geometries
are defined as in figure 3: the first utilizes the _Jinping prototype_ [26],
featuring 30 PMTs, and the second is an _ideal detector_ upgraded to 120 PMTs,
following the Fibonacci arrangement [27]. The number of PMTs in the ideal
detector follows the criterion in section 3.3. The central LS [28] is in an
acrylic shell with a water buffer. The Jinping prototype has an outlet and a
support base, which are removed in the ideal detector. Table 2 lists some
essential parameters. In reality, the QE and TT for each PMT are different,
and the TTS affects $R^{0}(t)$. We ignore such difference since $\psi(t)$
dominates the $R^{0}(t)$. $\mu_{i}$ in eq. (1) and $U_{i}$ in eq. (3) can be
obtained by calibration. Without loss of generality, we assume PMTs are
identical by setting $\mu_{i}=U_{i}=0$.

Table 2: Important parameters in simulation parameter | value (Jinping prototype) | value (ideal detector) | parameter | value
---|---|---|---|---
$r_{\mathrm{LS}}$ | $650\text{\,}\mathrm{m}\mathrm{m}$ | $650\text{\,}\mathrm{m}\mathrm{m}$ | QE | 0.2
$r_{\mathrm{PMT}}$ | $832\text{\,}\mathrm{m}\mathrm{m}$ | $900\text{\,}\mathrm{m}\mathrm{m}$ | TTS | $2.2\text{\,}\mathrm{n}\mathrm{s}$
number of PMTs | 30 | 120 | $t_{0}$ | $0\text{\,}\mathrm{n}\mathrm{s}$
Figure 3: Geometries of the Jinping prototype (left) and the ideal detector
(right).

$\psi(t)$ is parameterized in eq. (10) by the _decay time constant_
$\tau_{d}=$ $26.0\text{\,}\mathrm{n}\mathrm{s}$ and the _rise time constant_
$\tau_{r}=$ $1.6\text{\,}\mathrm{n}\mathrm{s}$,

$\mathrm{\psi}(t)\sim\exp\left(-\frac{t}{\tau_{d}}\right)\left[1-\exp\left(-\frac{t}{\tau_{r}}\right)\right].$
(10)

The average refraction index of LS is 1.48, causing TR to occur at
$r>0.88r_{\mathrm{LS}}$ which is $570\text{\,}\mathrm{m}\mathrm{m}$ for the
Jinping prototype.

(a) First six $c^{\lambda}_{l}(r)$ of Jinping prototype under different
energies.

(b) First six $c^{T}_{l}(r)$ of Jinping prototype under different energies.

(c) Predicted PE

(d) Ratio of simulation over predicted PEs

(e) Predicted timing

(f) Inhomogeneous Poisson process

Figure 4: (4) and (4) are regression results. TR occurs at the
$0.88r_{\mathrm{LS}}$ (dashed line). The error bars are $\pm 1$ standard
deviations. (4) – (4) uses $2\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$ $e^{-}$
in Jinping prototype. (4) and (4) are the PE and timing responses,
respectively. In (4), we shall look into the focus structures in the slices of
the dashed line $r=0.4r_{\mathrm{LS}}$ and dotted line $\theta=5\pi/6$ in
section 5.1. (4) is the ratio between the predicted PE and its MC truth. (4)
shows the average function of the inhomogeneous Poisson process compared with
the MC at $r=0.99r_{\mathrm{LS}}$ and $\theta=0$.

We focus on $e^{-}$ in this work. The training dataset includes multiple
batches. Vertices in each batch have the same $r$, ranging from
$0\text{\,}\mathrm{m}\mathrm{m}$ to $550\text{\,}\mathrm{m}\mathrm{m}$ with
step $10\text{\,}\mathrm{m}\mathrm{m}$, and from
$550\text{\,}\mathrm{m}\mathrm{m}$ to $640\text{\,}\mathrm{m}\mathrm{m}$ with
step $2\text{\,}\mathrm{m}\mathrm{m}$. Figure 4 subtracts $\log E$ in the 0th
order and the first six orders of $c^{\lambda}_{l}$ under 1, 2 and 3
$\mathrm{M}\mathrm{e}\mathrm{V}$. Figure 4 shows the first six orders of
$c^{T}_{l}$. The coefficient changes dramatically at $0.88r_{\mathrm{LS}}$,
where the TR just happens. It is evident from figure 4 and 4 that
$c^{\lambda}_{l}$ and $c^{T}_{l}$ depend on energy only when $l=0$. Therefore
one energy sample alone captures the detector response in any other energies.
The Jinping prototype utilizes a 25-PMT threshold [29] and the trigger
efficiency drops for events near the boundary. We take
$2\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$ events that are energetic enough to
guarantee $\sim 100\text{\,}\mathrm{\char 37\relax}$ trigger efficiency.

We can express $\lambda_{i}(r,\theta_{i})$ and $T_{i}(r,\theta_{i})$ in a
single heat map of a disk. Figure 4 and 4 shows the predicted PE and timing,
respectively. The ratio of the predicted PE and its truth is almost 1 in
figure 4, showing a good fit. The models of PE and timing at a specific
$r,\theta$ can be combined into an inhomogeneous Poisson process with average
function $\lambda_{i}(r,\theta_{i})R(t-T_{i}(r,\theta_{i}))$, shown in figure
4.

To determine the optimal number of parameters, we simulate another 15000
events as the validation dataset. A better goodness of fit for PE is
manifested by a higher score in log-likelihood given by eqs. (1) and (7).
Figure 5 shows the highest score requires the number of parameters to be
$35\times 35$ of $(\theta,r)$. The binary polynomials are shown for comparison
in figure 5. The number of parameters is up to 441 for Zernike and 750 for
double Legendre polynomials. Their scores are limited by the number of
parameters.

(a) Relative scores for different sub-models

(b) Relative scores vs the degree of freedom

Figure 5: Relative scores to the best model. (5) indicates the number of
parameters $35\times 35$ of $(\theta,r)$ to be the best. (5) adds the binary
polynomials. The numbers are the number of parameters on $\theta$ in the
varying coefficient models. “ _Zernike_ ” and “ _dLeg_ ” are Zernike and
double Legendre polynomials, respectively.

It is difficult to make a similar selection for $c^{T}_{l,m}$ because their
scores rely on $\tau$ and $t_{\mathrm{s}}$. Since PE is dominant in the small
detectors [17], we choose $10\times 35$ on $(\theta,r)$ for $c^{T}_{l,m}$ to
balance speed and accuracy.

###  2.3 Likelihood for reconstruction

We use MLE for reconstruction with the likelihood function satisfying

$\displaystyle\log\mathscr{L}(\mathbf{r},E,t_{0})$
$\displaystyle\sim\overbrace{\sum_{i}\sum_{j}\frac{\mathscr{R}_{\tau}(t_{ij}-T_{i}-t_{0})}{t_{\mathrm{s}}}}^{\mathrm{timing\
part}}+\overbrace{\sum_{i}\left(-\frac{E\lambda_{i}}{E_{0}}+n_{i}\log\frac{E\lambda_{i}}{E_{0}}\right)}^{\mathrm{PE\
part}}.$ (11)

The former is the timing part and the latter is the PE part. $E$ only
contributes to the PE part. $E_{0}$ is energy
($2\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$) in the training dataset. The
unbiased energy estimation of eq. (11) is

$\displaystyle\frac{\partial\log\mathscr{L}}{\partial
E}=0\Rightarrow\hat{E}=\frac{\sum_{i}n_{i}}{\sum_{i}{\lambda}_{i}(\mathbf{r})}E_{0}$
(12)

We use _Sequential Least Squares Programming_ [30] to maximize eq. (11). The
energy is calculated in each iteration by eq. (12) to reduce the time costs.

Most gradient-based methods are local optimizers. The origin of the local
maxima will be discussed in section 3. Our strategy to obtain the global
maximum is as follows. We generate two grids to calculate the predicted PE.
Each grid is $30\times 50\times 50$ equally spaced on $r$,
$\cos\theta_{\mathrm{v}}$ and $\phi_{\mathrm{v}}$. The inner grid covers the
radius of $0\text{\,}\mathrm{m}\mathrm{m}570\text{\,}\mathrm{m}\mathrm{m}$ and
the outer covers the rest. We choose the best points of the inner and outer
grids as initial values for the gradient optimizer. The larger of the two
$\log\mathscr{L}$’s is recorded.

##  3 Multimodality of the likelihood function

The convexity of the minus-log-likelihood in eq. (11) is crucial for
reconstruction. The reconstruction results are sensitive to the initial values
of the gradient optimizer if the minus-log-likelihood is not convex, or
equivalently, the likelihood function is _multimodal_. Using a poor vertex
undermines the energy resolution and the background rejection by fiducial
volume cuts.

In this section, we concentrate on PE (but not timing) since it is dominant in
the small detectors [17]. We replace the likelihood function in eq. (11) with
_cosine distance_ to eliminate the influence of fluctuations in the observed
PEs. We further narrow our focus to the 3 closest PMTs, giving a criterion of
multimodality in the likelihood function.

###  3.1 Cosine distance

We define _pattern_ $\Lambda$ as a vector containing the predicted PE on the
PMT space.

$\Lambda=\\{\lambda_{1},\lambda_{2},\cdots,\lambda_{N_{\mathrm{PMT}}}\\}$ (13)

If $\mathbf{r}_{1}$ and $\mathbf{r}_{2}$ are both the solution for an event
and $\mathbf{r}_{1}\neq\mathbf{r}_{2}$, we have
$\Lambda(\mathbf{r}_{1})=\alpha\Lambda(\mathbf{r}_{2})$, and $\alpha$ is an
arbitrary positive number. The _cosine distance_

$d_{\mathrm{cos}}\left[\Lambda(\mathbf{r}_{1})\mid\Lambda(\mathbf{r}_{2})\right]=1-\frac{\Lambda(\mathbf{r}_{1})\cdot\Lambda(\mathbf{r}_{2})}{\norm{\Lambda(\mathbf{r}_{1})}\norm{\Lambda(\mathbf{r}_{2})}}$
(14)

is zero.

$d_{\mathrm{cos}}$ matches well with the $-\log\mathscr{L}$. Figure 6 shows an
event at $(0,0,100)$ $\mathrm{m}\mathrm{m}$ of the Jinping prototype. Two
local solutions with different initial values are shown in black dashed lines.
The green line scans the $-\log\mathscr{L}$ between the two. The blue line
shows the $d_{\cos}$ of these points by the true expected PEs. Figure 6 and 6
are at $(540,0,0)$ $\mathrm{m}\mathrm{m}$ in the ideal detector. Figure 6
demonstrates a good unimodal event, in which both $d_{\cos}$ and
$-\log\mathscr{L}$ are convex. Figure 6 is an example of multimodal event
where $-\log\mathscr{L}$ is not convex, but the $d_{\cos}$ is the same as
figure 6. Observation fluctuations lead to the $-\log\mathscr{L}$ difference
between the two. Of all events at $(540,0,0)$ $\mathrm{m}\mathrm{m}$, the
fraction of bad events like figure 6 is less than $1\text{\,}\%$.

(a) Vertex at $(0,0,100)$ $\mathrm{m}\mathrm{m}$

(b) vertex at $(540,0,0)$ $\mathrm{m}\mathrm{m}$

(c) Vertex at $(540,0,0)$ $\mathrm{m}\mathrm{m}$

Figure 6: Comparison of $d_{\cos}$ with $-\log\mathscr{L}$. All events at
$(0,0,100)$ $\mathrm{m}\mathrm{m}$ degenerate with $r=$
$600\text{\,}\mathrm{m}\mathrm{m}$ on the $z$ axis. (6) is one typical event
and the two local solutions are in black dashed lines. The true vertex is in
the red line. The green line scans the $-\log\mathscr{L}$ and the blue line
scans the $d_{\cos}$, showing that they both have two local minima. (6) and
(6) use true events at $(540,0,0)$ $\mathrm{m}\mathrm{m}$, comparing one
unimodal with one multimodal event. $d_{\cos}$ itself is convex. The
difference in $-\log\mathscr{L}$ between (6) and (6) is from fluctuations.

(a) cosine distance on $x$ axis

(b) Cosine distance on $z$ axis

(c) Similarity for (400,0,0) $\mathrm{m}\mathrm{m}$ to $xOy$ plane

(d) zoom for (7)

Figure 7: Cosine distance of different cases. The color shows the value of
$d_{\cos}[\Lambda(\mathbf{r}_{1})\mid\Lambda(\mathbf{r}_{2})]$. (7)
$\mathbf{r}_{1}$ and $\mathbf{r}_{2}$ are both on the $x$ axis. (7)
$\mathbf{r}_{1}$ and $\mathbf{r}_{2}$ are both on the $z$ axis. The red dashed
line shows the center and radius at $600\text{\,}\mathrm{m}\mathrm{m}$
degenerate. (7) $\mathbf{r}_{1}$ is $(400,0,0)\ \mathrm{mm}$ and
$\mathbf{r}_{2}$ scans the $xOy$ plane. (7) is the zoom for (7). The red star
and the red circles degenerate.

Figure 7 maps out multimodality of $d_{\cos}$. $\mathbf{r}_{1}$ and
$\mathbf{r}_{2}$ are all on the $x$ axis in figure 7, and on the $z$ axis in
figure 7. The slice in figure 7 shows $z=$ $0\text{\,}\mathrm{m}\mathrm{m}$
and $600\text{\,}\mathrm{m}\mathrm{m}$ degenerate. In figure 7,
$\mathbf{r}_{1}$ is $400\text{\,}\mathrm{m}\mathrm{m}$ on the $x$ axis and
$\mathbf{r}_{2}$ scans the $xOy$ plane. Figure 7 zooms around the true
position. Local minima are found near the detector boundary. The shape of
$d_{\cos}$ map is remarkably consistent with the reconstruction results in
figure 10, which will be discussed in section 4.1.

###  3.2 Cosine distance for 3-PMT case

Observing that closer PMTs contribute more to the likelihood function, we only
consider the three closest PMTs to the vertex for simplicity. Considering
symmetry, in figure 8 the vertices are constrained in the shadow region
without loss of generality. We define the _contour function_ by

$C_{1i}(\mathbf{r})=\frac{\lambda_{1}(\mathbf{r})}{\lambda_{i}(\mathbf{r})},$
(15)

where $\lambda_{1}(\mathbf{r})$ is the expected number of PEs at the closest
PMT to $\mathbf{r}$, and $\lambda_{i}(\mathbf{r})$ is that of the 2nd or 3rd
closest PMT with $i=2$ or $3$ respectively. In a two-dimensional slice of the
detector, a certain value of $C_{1i}(\mathbf{r})$ defines an _equivalent
ratio_ (ER) line. If one ER line of $C_{12}$ and another in $C_{13}$ has two
intersections $\mathbf{r}_{1}$ and $\mathbf{r}_{2}$,
$d_{\mathrm{cos}}\left[\Lambda(\mathbf{r}_{1})\mid\Lambda(\mathbf{r}_{2})\right]\approx
0$. Because the 3 dominanting PMTs satisfies
$\lambda_{1}(\mathbf{r}_{1}):\lambda_{2}(\mathbf{r}_{1}):\lambda_{3}(\mathbf{r}_{1})=\lambda_{1}(\mathbf{r}_{2}):\lambda_{2}(\mathbf{r}_{2}):\lambda_{3}(\mathbf{r}_{2})$.
Therefore it is the key to count the intersections between any pair of
$C_{12}$ and $C_{13}$ ER lines.

We study the multimodality affected by angle of incidence $\beta$ for
_homogeneous materials_ and TR for _inhomogeneous materials_ respectively. In
the former case, the LS and the buffer have the same refractive index.
Assuming the PMTs are small enough and their surface is flat, the predicted PE
$\lambda_{i}$ is proportional to the solid angle $\Omega_{i}$. $\Omega_{i}$ is
subtended by the $i$th PMT from the vertex,

$\lambda_{i}(\mathbf{r})\propto\Omega_{i}(\mathbf{r})\propto\frac{1}{{\lVert\mathbf{r}-\mathbf{r}_{\mathrm{PMT},i}\rVert}_{2}^{2}}\cos\beta_{i},$
(16)

where $\mathbf{r}$, $\mathbf{r}_{\mathrm{PMT},i}$ are defined in figure 1.
$\beta$ is the angle of incidence on a PMT. Eq. (16) well describes TAO [31]
to be equipped with SiPM and XMASS [32] with flat-photocathode R10789 PMTs
[33]. To comparatively study the effect by $\beta_{i}$ by leaving it out,

$\lambda_{i}(\mathbf{r})\propto\Omega_{i}(\mathbf{r})\propto\frac{1}{{\lVert\mathbf{r}-\mathbf{r}_{\mathrm{PMT},i}\rVert}_{2}^{2}}.$
(17)

We set $r_{\mathrm{PMT}}=r_{\mathrm{LS}}$, namely no buffer is installed
between LS and PMTs. The PMTs located at $0$, $2\pi/3$ and $4\pi/3$ directions
are indexed by 1, 2 and 3. $C_{12}$ and $C_{13}$ are symmetric. Figure 8
demonstrates the ER lines with $\beta$ and figure 8 demonstrates those
without. The two intersections circled red in figure 8 show that $\beta$ leads
to multimodality.

(a) ER lines of eq. (17)

(b) ER lines of eq. (16)

(c) ER lines of eq. (16)

(d) Jinping outlet arrangement

(e) Jinping general arrangement

(f) Ideal detector arrangement

Figure 8: ER lines of different detector response models. PMTs are in black
triangles. (8) – (8) are the homogeneous materials with
$r_{\mathrm{PMT}}=r_{\mathrm{LS}}$. For (8) and (8), PMTs are arranged at 0,
$2\pi/3$, $4\pi/3$. For (8), PMTs are arranged at 0, $\pm\pi/100$. (8) – (8)
are the inhomogeneous materials, using the coefficients shown in figure 4. (8)
is the outlet direction, PMTs are arranged at $2\pi/5$ and $\pm\pi/5$, lacking
at 0. The intersections of $C_{12}=1$ and $C_{13}=1$ are at the central and
top regions. (8) is the general direction. PMTs are arranged at $0$ and
$\pm\pi/5$. More than one intersections of $C_{12}=1$ and $C_{13}=4.5$ are at
$350\text{\,}\mathrm{m}\mathrm{m}$ and $550\text{\,}\mathrm{m}\mathrm{m}$. (8)
is the ideal detector without multimodality. PMTs are arranged at 0,
$\pm\pi/18$. For all subfigures, the vertices are constrained in the shadow
region. We circle the multiple intersections of ER lines in (8), (8), (8) and
(8).

To account for the suspicion that multimodality is induced by sparse PMT
arrangement in figure 8, consider two intersections with a more compact
arrangement under eq. (16) in figure 8. Regardless of the PMT arrangement, the
shadow boundary is not only an ER line but also the bisector of the two
neighboring PMTs. The point on the bisector at $r=r_{\mathrm{LS}}$ is always
one of the intersections under eq. (16). PMT is almost blind to this point due
to a large $\beta$. To avoid this, a non-fluorescent buffer is necessary to
isolate scintillation events away from the PMTs.

We use _inhomogeneous materials_ to study the TR effect. The buffer’s
refractive index is different from the central detector. For the Jinping
prototype, we see from the fitted model that TR makes ER lines segmented. A
two-dimensional slice of the Jinping prototype includes 10 PMTs. We pick the 3
PMTs at $(0,\pm{\pi/5})$ to represent the _general direction_ in figure 8,
matching figure 7. It is evident from figure 8 that events at
$350\text{\,}\mathrm{m}\mathrm{m}$ and $550\text{\,}\mathrm{m}\mathrm{m}$
degenerate. We also pick $(\pm{\pi/5},2\pi/5)$ as the _outlet direction_ in
figure 8, matching figure 7. Vertices at $z=$
$600\text{\,}\mathrm{m}\mathrm{m}$ and the center degenerate. Such 3-PMT plots
are reasonable simplifications to reproduce the trends of $d_{\cos}$ in all-
PMT cases.

The ER lines near the TR region are distorted compared to the homogeneous
materials. The rich functional structures in this region are a hotbed to
multimodality. A solution is to reduce the gaps between PMTs, making them
close enough so as not to fall into each other’s TR regions. For example, the
ideal detector has no multiple intersections in ER lines in figure 8.

###  3.3 A criterion against multimodality

In figure 8, the predicted PE is approximately proportional to the inverse
square of the distance from the vertex to PMT. The incident angle and the TR
break that trend. These effects are so strong that the closest PMTs could
receive fewer PEs, making the ER lines distorted and intersect multiple times.
They are the seeds to multimodality. Therefore a perfect inverse-squared
detector in figure 8 is always free from multimodality.

Taking $\lambda_{1}$ and $\lambda_{2}$ to be the same meaning as eq. (15), in
figures 8 and 8, adding buffer is equivalent to reducing $\lambda_{1}$.
Conditions of figures 8 and 8 improve by reducing the gap between neighboring
PMTs, making $\lambda_{2}$ larger. The two phenomena can be unified by
$\lambda_{1}/\lambda_{2}$: for any event, $\lambda_{1}/\lambda_{2}<10$ is
required to avoid multimodality.

It is sufficient to examine if an extreme vertex with the biggest
$\lambda_{1}/\lambda_{2}$ is less than 10. Such an extreme vertex is
approximately realized in the center-PMT direction at $r=r_{\mathrm{LS}}$,
illustrated in figure 9. Note that $\lambda_{1}/\lambda_{2}<10$ is a necessary
condition, we might construct some case where the TR region is fully contained
in the gap between 2 PMTs to embed a lot of degeneracy but still satisfies
$\lambda_{1}/\lambda_{2}<10$. $\lambda_{1}$ and $\lambda_{2}$ is related to
the $r_{\mathrm{LS}}$, $r_{\mathrm{PMT}}$, $N_{\mathrm{PMT}}$ and the buffer’s
refractive index.

(a) The criterion $\lambda_{1}/\lambda_{2}>10$.

(b) Least needed number of the PMTs based on the refractive index of the outer
buffer and $r_{\mathrm{PMT}}$

Figure 9: (9) shows the criterion to avoid multimodality.
$\lambda_{1}/\lambda_{2}$ should be less than 10 and the TR region should be
excluded. (9) shows the least number of needed PMTs. $r_{\mathrm{LS}}$ is
fixed to $645\text{\,}\mathrm{m}\mathrm{m}$. For the Jinping prototype, the
average refractive index of the water is 1.33 and $r_{\mathrm{PMT}}$ is
$832\text{\,}\mathrm{m}\mathrm{m}$, showing that at least 120 PMTs are
necessary. The red star is the ideal detector with $r_{\mathrm{PMT}}=$
$900\text{\,}\mathrm{m}\mathrm{m}$, requiring at least 90 PMTs. The larger the
$r_{\mathrm{PMT}}$, the refractive index and the number of PMTs help to avoid
multimodality.

Keeping an identical distance between neighboring PMTs to extend our criterion
to three dimensions, $N_{\mathrm{PMT,3D}}$ scales approximately as
$N_{\mathrm{PMT,2D}}^{2}/3$. Fixing
$r_{\mathrm{LS}}=$650\text{\,}\mathrm{m}\mathrm{m}$$, figure 9 shows the least
needed number of PMTs under different $r_{\mathrm{PMT}}$ and the buffer
refractive index. For a certain $N_{\mathrm{PMT}}$, larger $r_{\mathrm{PMT}}$
or the refractive index gives better resilience against degeneracy. For
$r_{\mathrm{PMT}}=$ $832\text{\,}\mathrm{m}\mathrm{m}$ and the refractive
index 1.33, 120 PMTs are needed. Using 8-inch PMTs, the PMT coverage,
calculated by
$N_{\mathrm{PMT}}\times(10\mathrm{mm})^{2}/4r^{2}_{\mathrm{PMT}}$, should
exceed $40\text{\,}\%$. Be aware that the larger $r_{\mathrm{PMT}}$ leads to
poorer PMT coverage, thus $r_{\mathrm{PMT}}$ should be as close to the lower
limit in figure 9 as possible. $r_{\mathrm{PMT}}$’s other influences will be
discussed in section 5.1.

The ideal detector with $r_{\mathrm{PMT}}=$ $832\text{\,}\mathrm{m}\mathrm{m}$
requires at least 90 PMTs. To reduce the influences of fluctuations, we use
120 PMTs to guarantee perfect reconstruction performance.

##  4 Reconstruction results

We verify our vertex reconstruction on MC and raw data. For MC, we use the
Jinping prototype to verify the conditions of multimodality in section 3.3,
and the ideal detector to study the reconstruction bias and resolution. The
simulations use $2\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$ $e^{-}$ on $x$ and
$z$ axes: $x$ represents the general direction, and $z$ the outlet direction,
aligned with the definitions in section 3.2. For raw data, we analyze the
^214Bi-^214Po cascade signals by the model fitted from the Jinping-prototype
simulation.

###  4.1 Results from simulations

We compare the results using MLE and BC. For BC, the reconstructed vertex
$\hat{\mathbf{r}}$ is

$\hat{\mathbf{r}}=1.5\times\frac{\sum_{i}q_{i}\mathbf{r}_{\mathrm{PMT},i}}{\sum_{i}q_{i}}$
(18)

and the energy is scaled from the total number of hits,

$\frac{\hat{E}}{[\mathrm{MeV}]}=\frac{N_{\mathrm{hit}}}{65}$ (19)

where $1.5$ is a correction factor and 65 is the average total PE at
$1\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$ at the Jinping prototype. The
simulation data ranges from
$0\text{\,}\mathrm{m}\mathrm{m}650\text{\,}\mathrm{m}\mathrm{m}$ with steps of
$10\text{\,}\mathrm{m}\mathrm{m}$. Each step contains 5000 events.

(a) MLE

(b) BC

(c) MLE

(d) BC

(e) Energy reconstruction on $x$ axis

(f) Energy reconstruction on $z$ axis

(g) MLE

(h) BC

(i) MLE

(j) BC

(k) Energy reconstruction on $x$ axis

(l) Energy reconstruction on $z$ axis

Figure 10: Vertex and energy reconstruction performance of
$2\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$ $e^{-}$. (10) – (10) are from the
Jinping prototype. The red dashed lines in (10) and (10) show the
degeneracies. (10) – (10) are from the ideal detector. In (10) – (10) and (10)
– (10), the light grey, grey and black represent bins with 1–10, 10–100 and
$>100$ events, respectively. (10), (10) and (10), (10) are the energy
reconstruction results of $2\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$ $e^{-}$
on the Jinping prototype and the ideal detector. “Fit(ReconV)” and
“Fit(TrueV)” are energies with reconstructed and true vertices, respectively.
“Scaled total PE” is by eq. (19).

Figure 10 shows a reasonable vertex reconstruction in the general direction.
That is in contrast to the outlet direction in figure 10, where events at the
detector center and around $(0,0,600)$ $\mathrm{m}\mathrm{m}$ are
indistinguishable. Similarly, events at $350\text{\,}\mathrm{m}\mathrm{m}$ and
$550\text{\,}\mathrm{m}\mathrm{m}$ degenerate. Such degeneracy is a
consequence of multimodality in the likelihood function as discussed in
section 3.3. Notice that the $d_{\cos}$ relations in figures 7 and 7 predict
the major patterns of MLE reconstructions in figures 10 and 10.

In figures 10 and 10, similar worsen trends appear for BC from the general to
the outlet direction. They both have severe biases.

For the ideal detector, no outlet is considered. $x$ and $z$ axes are
equivalent. MLE performs well without degeneracy in figures 10 and 10, proving
the effectiveness of the $\lambda_{1}/\lambda_{2}<10$ criteria in section 3.3.
BC, on the contrary, still reconstructs badly in the TR region. Our detector
model in section 2 describes the TR region well and consequently MLE provides
a big improvement upon the BC method.

The energy reconstruction for the Jinping detector is shown in figure 10 and
10. The scaled total PE decreases rapidly in the TR region, which introduces a
big bias. The reconstructed energy by MLE is even worse, due to the wrong
vertices. The true vertices give an unbiased energy estimation, indicating
that the detector response model is accurate. For the ideal detector of figure
10 and 10, BC has similar trends. The MLE is almost unbiased and improved by
eliminating the vertex degeneracy.

Excluding the regions $|\hat{\mathbf{r}}|>$ $600\text{\,}\mathrm{m}\mathrm{m}$
and $\sqrt{\hat{x}^{2}+\hat{y}^{2}}<$ $50\text{\,}\mathrm{m}\mathrm{m}$, the
standard deviation of Jinping prototype is $60\text{\,}\mathrm{m}\mathrm{m}$
at the center and $20\text{\,}\mathrm{m}\mathrm{m}$ in the TR region. The
energy reconstruction heavily depends on the vertex. When using the true
vertex, the results are almost unbiased even in the TR region. The energy
resolution of the Jinping prototype is approximately
$11\text{\,}\%$$/\sqrt{E\mathrm{[MeV]}}$.

###  4.2 Results from raw data

The raw data of the Jinping prototype is the PMT waveform $w_{i}(t)$, which is
a convolution of hits $\sum_{j}\delta(t-t_{ij})$ and the _single PE response_
$V_{\mathrm{PE}}(t)$,

$w_{i}(t)=\sum_{j}\delta(t-t_{ij})\otimes
V_{\mathrm{PE}}(t)+\epsilon(t)=\sum_{j}V_{\mathrm{PE}}(t-t_{ij})+\epsilon(t),$
(20)

$\epsilon(t)$ is Gaussian white noise. We use _Richardson-Lucy direct
demodulation_ (LucyDDM) for deconvolution [34]. The input is the timing series
$t_{ij}$, returning the gain modified charge $q_{ij}$. Traditional LucyDDM is
biased in photon density due to an artificial threshold. Xu et al. [35] use a
rescaling factor to significantly reduce the bias. In terms of the photon
density resolution, LucyDDM matches the fit method [35] and consumes less
time.

Xu et al. [35] also utilize the non-normalized Kullback-Leibler (KL)
divergence [36, 37] for reconstruction, which is a special case of _density
power divergence_ [37]. It joins the waveform analysis results
$\sum_{j}q_{ij}\delta(t-t_{ij})$ and the inhomogeneous Poisson process. We
take the non-normalized KL divergence as the pseudo-minus-log-likelihood in
this work.

$\displaystyle-\log\mathscr{L}(E,\mathbf{r},t_{0})$
$\displaystyle=\sum_{i}D_{\mathrm{KL}}\left[\sum_{j}q_{ij}\delta(t-t_{ij})\;\middle|\;E\lambda_{i}(r,\theta_{i})R(t-t_{0}-T_{i}(r,\theta_{i}))\right]$
(21) $\displaystyle=-\sum_{i}\left[\sum_{j}q_{ij}\log
R(t_{ij}-t_{0}-T_{i})+\left(\sum_{j}q_{ij}\right)\log(E\lambda_{i})-E\lambda_{i}\right]+\mathrm{Const.}$
$\displaystyle\sim\overbrace{\sum_{i}\sum_{j}q_{ij}\frac{-\mathscr{R}_{\tau}(t_{ij}-T_{i}-t_{0})}{t_{\mathrm{s}}}}^{\mathrm{timing\
part}}+\overbrace{\sum_{i}\left[-q_{i}\log(E\lambda_{i})+E\lambda_{i}\right]}^{\mathrm{PE\
part}},$

where $q_{i}=\sum_{j}q_{ij}$.

The detector response model is based on simulation since there is no dedicated
calibration runs [29]. We check the ^214Bi-^214Po cascade signal of $e^{-}$
prompt and $\alpha$ delayed. The cuts are listed below:

  1. 1.

$|\hat{\mathbf{r}}_{\mathrm{prompt}}|<$600\text{\,}\mathrm{m}\mathrm{m}$$ and
$|\hat{\mathbf{r}}_{\mathrm{delay}}|<$600\text{\,}\mathrm{m}\mathrm{m}$$,

  2. 2.

$\sqrt{\hat{x}_{\mathrm{prompt}}^{2}+\hat{y}_{\mathrm{prompt}}^{2}}>$50\text{\,}\mathrm{m}\mathrm{m}$$
and
$\sqrt{\hat{x}_{\mathrm{delay}}^{2}+\hat{y}_{\mathrm{delay}}^{2}}>$50\text{\,}\mathrm{m}\mathrm{m}$$,

  3. 3.

Visible energy of $e^{-}$ is less than
$3.5\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$,

  4. 4.

Visible energy of $\alpha$ is in
$0.7\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}1\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$,

  5. 5.

Delayed time between the prompt and delayed signal in
$10\text{\,}\mathrm{u}\mathrm{s}1000\text{\,}\mathrm{u}\mathrm{s}$,

  6. 6.

$|\hat{\mathbf{r}}_{\mathrm{prompt}}-\hat{\mathbf{r}}_{\mathrm{delay}}|<$300\text{\,}\mathrm{m}\mathrm{m}$$.

Figures 11, 11, 11 and 11 show data distribution of the six cuts except the
3rd, 4th, 5th and 6th, respectively. The cut is colored red in each subfigure.

(a) Prompt energy cut

(b) Delayed energy cut

(c) Delayed time cut

(d) Distance cut

Figure 11: Reconstructions of ^214Bi-^214Po cascades in raw data.
$\mathscr{N}(\mu,\sigma^{2})$ is a Gaussian distribution with the mean $\mu$
and variance $\sigma^{2}$. (11) and (11) is fitted by Gaussian. (11) is fitted
by an exponential distribution.

Most backgrounds are gammas from the PMTs. For simplicity, we treat it as a
constant. The 25-PMT threshold distorts the beta decay energy spectrum. We use
Gaussian to fit the prompt and delayed signals. The average energy of the
$e^{-}$ is approximately $\sim$
$1.879\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$, and the $\alpha$ peaks at
$0.765\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$ due to the ionization
quenching. The fitted half-life is close to the true value
($164.3\text{\,}\mathrm{u}\mathrm{s}$). The results show that most selected
events are cascade signals and our reconstruction works well.

##  5 Discussion

###  5.1 Other reflections by acrylic shell

The closer PMT often has a larger gradient in the likelihood function. We
select the three closest PMTs in section 3.2 since they are more sensitive to
vertex positions. In inhomogeneous materials, normal reflections other than TR
at the media boundary focus light, making predicted PEs peak at some regions
that are otherwise as dim as their neighbors. Such focal points and lines are
sources to multimodality in the likelihood function.

When the vertices $r/r_{\mathrm{LS}}\in[0.4,0.5]$, figure 12 (top) shows the
photons reflected by the acrylic shell are focused at $\cos\theta\simeq-1$. We
record $\theta_{\mathrm{acr}}$, which is the central angle defined by the
vertex and the point of incidence to the acrylic shell, and compare the
$\theta_{\mathrm{acr}}$’s distribution of blue PMT at $\theta\simeq 0$ and
black PMT at $\theta\simeq\pi$. In figure 12 (bottom), most photons hit the
blue PMT directly, but quite a lot of photons bounce to the black PMT after
reflection by the acrylic shell. $r=0.4r_{\mathrm{LS}}$ is therefore denoted
as the _focus region_. In figure 12, the slice of the predicted PE at
$r=0.4r_{\mathrm{LS}}$ peaks at $\cos\theta\simeq-1$, introducing degeneracy
of PE predictions. It is also observed at the Borexino CTF [38]. The dashed
horizontal line drawn at the abnormal $\cos\theta\simeq-1$ peak intersects
with the PE prediction line, which defines a _degenerate region_.

(a) Photon tracks at $r=$ $260\text{\,}\mathrm{m}\mathrm{m}$

(b) Predicted PE at $r=$ $260\text{\,}\mathrm{m}\mathrm{m}$

(c) Predicted PE at $\theta=\frac{5}{6}\pi$

Figure 12: Focus effect at $r=0.4r_{\mathrm{LS}}=$
$260\text{\,}\mathrm{m}\mathrm{m}$. The tracks of detected photons are
sketched in (12). The top of (12) defines $\theta_{\mathrm{acr}}$ of the
tracks. The vertex is relative to PMT1 at $\cos\theta_{1}=-1$ and PMT2 at
$\cos\theta_{2}=1$. The bottom of (12) is the log-scale distribution of
$\cos\theta_{\mathrm{acr}}$. (12) is a slice of $r=0.4r_{\mathrm{LS}}$ shown
in figure 4. “MC” is the histogram of the simulated events, and “Fit”
represents the fitted response in section 2. The inlet zooms in the peak at
$\cos\theta=-1$. (12) is the predicted PE at $\theta=5\pi/6$ conditioned to
different $r_{\mathrm{PMT}}$. Each dashed line in (12) and (12) covers a
degenerate region.

The opposite PMT plays an important role in the focus region. Vertex
reconstruction at $r$\sim 260\text{\,}\mathrm{m}\mathrm{m}$$ of simulated
events from the Jinping prototype is shown in figure 13. The hotspots in the
zenith $\theta_{\mathrm{v}}$ and azimuth $\phi_{\mathrm{v}}$ map coincide with
the PMT-center directions. Our speculation is verified by masking out one PMT
in eq. (11), resulting in figure 13. The hotspot disappears around the masked
PMT.

(a) Recon with all PMTs

(b) Recon with a closed PMT

Figure 13: Vertex reconstruction results at $r=$
$260\text{\,}\mathrm{m}\mathrm{m}$ subjected to the focus effect. The color is
the count of the reconstructed vertices. (13) is the
$\theta_{\mathrm{v}}$-$\phi_{\mathrm{v}}$ map of the vertex reconstruction
results in spherical coordinates. The circles indicate the PMT-center
directions. (13) is the reconstruction masking out the triangle PMT, noted as
“masked”.

For $r>0.8r_{\mathrm{LS}}$, photons can reflect multiple times. Figure 12
takes the slice of $\theta=5\pi/6$ as an example. The superposition of
multiple reflections makes the predicted PE peak at around
$r>0.8r_{\mathrm{LS}}$, causing degeneracy. The two local solutions of figure
6 coincides with the degenerate region of the
$r_{\mathrm{PMT}}=$900\text{\,}\mathrm{m}\mathrm{m}$$ in figure 12. Evident in
figure 4, such a multiple-reflection region forms a belt, unfortunately it is
no longer possible to select out 3 PMTs to derive a simple criterion.

Generally speaking, degeneracy is inevitable for inhomogeneous materials. Each
PMT contributes a specific degenerate region to the likelihood function. For
example, vertices around $r=0.5r_{\mathrm{LS}}$ could degenerate with both the
focus region and the TR region simultaneously. The situation is further
complicated by the fluctuations on every PMT. Luckily, if the degeneracy scale
is less than the vertex resolution, the reconstruction result is effectively
not ambiguous anymore. Higher PMT coverage gives a larger slope in the
predicted PEs, narrowing the degenerate region, as illustrated by the
$r_{\mathrm{PMT}}=$850\text{\,}\mathrm{m}\mathrm{m}$$ line in figure 12.
Therefore we recommend a small $r_{\mathrm{PMT}}$ at the lower limit in figure
9.

###  5.2 Breaking of spherical symmetry

In figure 10, energy reconstruction using true vertices “Fit(TrueV)” is still
biased when $r=$ $620\text{\,}\mathrm{m}\mathrm{m}$ at $z$ axis. Due to the
outlet in figure 3, only a small fraction of the photons can be detected in
the outlet direction compared to the general direction.

(a) A schematic of the outlet shielding

(b) Photon tracks of general and outlet direction

Figure 14: (14) Vertices are at $r=$ $620\text{\,}\mathrm{m}\mathrm{m}$. The
_shielded tracks_ are the detected tracks in the general direction but not in
the outlet direction. The numbers in the zoom axis show the
$\theta_{\mathrm{acr}}$’s value in radians. (14) The outlet region is on the
left side of the grey dashed line. The distribution of the
$\theta_{\mathrm{acr}}$ of the outlet and general directions are noted as
“Outlet” and “General”.

In figure 14, the outlet breaks spherical symmetry. Figure 14 compares the
$\theta_{\mathrm{acr}}$ distributions of the outlet and general directions.
Photons with $\theta_{\mathrm{acr}}<0.10$ rad are trapped into the outlet and
seldomly detected by any PMT.

###  5.3 Cosine distance for timing

In section 3.1, we ignore the timing information since PE is dominant in the
Jinping Prototype. We can expand the vector $\Lambda$ to include timing for
large detectors. The input is the average function $\lambda_{i}R_{i}(t)$ of
the inhomogeneous Poisson process on each PMT. For continuous form, the
$d_{\mathrm{cos}}$ is

$\displaystyle
d_{\mathrm{cos}}\left[\lambda_{i}(\mathbf{r}_{1})R_{i}(t;\mathbf{r}_{1})\mid\lambda_{i}(\mathbf{r}_{2})R_{i}(t;\mathbf{r}_{2})\right]$
(22) $\displaystyle=$ $\displaystyle
1-\frac{\sum_{i}\int\lambda_{i}(\mathbf{r}_{1})R_{i}(t;\mathbf{r}_{1})\lambda_{i}(\mathbf{r}_{2})R_{i}(t;\mathbf{r}_{2})\mathrm{d}t}{\sqrt{\sum_{i}\int[\lambda_{i}(\mathbf{r}_{1})R_{i}(t;\mathbf{r}_{1})]^{2}\mathrm{d}{t}}\sqrt{\sum_{i}\int[\lambda_{i}(\mathbf{r}_{2})R_{i}(t;\mathbf{r}_{2})]^{2}\mathrm{d}{t}}}.$

###  5.4 Upgrading Jinping Prototype

The criterion in section 3.3 is a guideline for future detectors. Due to the
space limits, PMT coverage over $40\text{\,}\%$ is difficult to achieve at a
small radius. Vertices in the central and the TR region make a huge difference
in energy reconstruction. If we focus on the most severe degeneracy at the
detector center and $600\text{\,}\mathrm{m}\mathrm{m}$ on the $z$-axis, the
main problem is that the gap between PMT1 and PMT3 is too large, as shown in
figure 8. Vertices at the bisector of the 1st and 2nd closest PMTs are in the
TR regions of both PMTs. Therefore the criterion could be updated to “the TR
regions of the neighboring PMTs do not overlap”. Compared to the criterion in
section 3.2 and figure 9 where the vertices are in the PMT-center direction,
the updated minimum number of the PMTs relaxes to 1/2 in a two-dimensional
slice and 1/4 in three dimensions. It shows that a minimum of 30 PMTs could
meet the demand, but should be uniformly distributed. In the current 30-PMT
arrangement, 0-to-$600\text{\,}\mathrm{m}\mathrm{m}$ large degeneracy is
caused by the absence of PMTs in the outlet direction, a big non-uniformity.
The Jinping Neutrino Experiment collaboration (JNE) plans to upgrade the
prototype to 60 PMTs and narrow the gaps between the PMTs around the outlet.

##  6 Conclusion

We obtain an accurate detector response model using regression, in which the
number of parameters is determined by validation. Compared to purely optics-
motivated models inevitably biased at the TR regions and data-driven ones
overwhelmed by the degrees of freedom, our model achieves an optimal balance.
It is shown to be unbiased in the TR regions while keeping the model
complexity under control. For $2\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$
$e^{-}$ Jinping prototype simulation, the vertex resolution is
$60\text{\,}\mathrm{m}\mathrm{m}$ at the center and
$20\text{\,}\mathrm{m}\mathrm{m}$ at the TR region except those around the
outlet. The energy resolution is $11\text{\,}\%$$/\sqrt{E\mathrm{[MeV]}}$. The
method is confirmed to work on Jinping prototype raw data by ^214Bi-^214Po
analysis. We believe our construction suits all the spherical detectors,
especially handles the optical complexity of those with the TR.

We investigate the reconstruction degeneracy at the Jinping prototype, and
confirm its origin to be the multimodality in likelihood functions. With a set
of carefully chosen approximations, we derive a necessary condition for a
detector to be free from reconstruction degeneracy: the expected PE ratio
between the 2 closest PMTs of any event should be less than 10. The criterion
is justified by comparing the Jinping prototype and an ideal detector setup.
That simple degeneracy criterion marks the first thorough and systematic study
on the mis-reconstruction of LS detectors known to us. It guides the upgrade
of the Jinping prototype and hopefully will serve as a valuable reference for
the design of PMT configurations at future LS detectors.

## Acknowledgments

We appreciate the development of the simulation tool of JSAP by Linyan Wan,
Ziyi Guo and Lei Guo. Many thanks to the efforts in the waveform
reconstruction by Aiqiang Zhang and Dacheng Xu. We are also grateful to Wentai
Luo and Xuewei Liu for discussions on reconstruction. The PMT arrangement is
inspired by Bohan Qi. The Corresponding author would like to thank the XMASS
collaboration for facilitating the idea of detector response modeling by
spherical harmonics. The early seed of this paper roots in exciting
discussions with Professors Shiro Ikeda, John Gregory Learned, Kai Uwe Martins
and Yoichiro Suzuki. We also thank the Jinping Neutrino Experiment
collaboration for sharing the data from the Jinping 1-ton prototype. This work
was supported in part by the National Natural Science Foundation of China (No.
12127808 and 12141503) and the Key Laboratory of Particle and Radiation
Imaging (Tsinghua University).

## References

  * [1] Borexino collaboration, _Science and technology of Borexino: a real-time detector for low energy solar neutrinos_ , _Astroparticle Physics_ 16 (2002) 205.
  * [2] KamLAND collaboration, _Precision Measurement of Neutrino Oscillation Parameters with KamLAND_ , _Phys. Rev. Lett._ 100 (2008) 221803.
  * [3] SNO+ collaboration, _Current Status and Future Prospects of the SNO+ Experiment_ , _Adv. High Energy Phys._ 2016 (2016) 6194250 [1508.05759].
  * [4] JUNO collaboration, _JUNO physics and detector_ , _Prog. Part. Nucl. Phys._ 123 (2022) 103927 [2104.02565].
  * [5] JUNO collaboration, _Neutrino Physics with JUNO_ , _Journal of Physics G: Nuclear and Particle Physics_ 43 (2016) 030401.
  * [6] JUNO collaboration, _JUNO physics and detector_ , _Progress in Particle and Nuclear Physics_ 123 (2022) 103927.
  * [7] KamLAND-Zen collaboration, _Search for majorana neutrinos near the inverted mass hierarchy region with kamland-zen_ , _Phys. Rev. Lett._ 117 (2016) 082503.
  * [8] KamLAND-Zen collaboration, _First Search for the Majorana Nature of Neutrinos in the Inverted Mass Ordering Region with KamLAND-Zen_ , arXiv:2203.02139.
  * [9] SNO+ collaboration, _The SNO+ experiment_ , _JINST_ 16 (2021) P08059 [2104.11687].
  * [10] J. F. Beacom, S. Chen, J. Cheng, S. N. Doustimotlagh, Y. Gao, S.-F. Ge et al., _Letter of Intent: Jinping Neutrino Experiment_ , _Chinese Phys. C_ 41 (2017) 023002.
  * [11] Q. Liu, M. He, X. Ding, W. Li and H. Peng, _A vertex reconstruction algorithm in the central detector of JUNO_ , _Journal of Instrumentation_ 13 (2018) T09005–T09005.
  * [12] H. S. Kim, _Finding an Event Vertex by Using a Weighting Method at RENO_ , _New Phys. Sae Mulli_ 62 (2012) 631.
  * [13] Z. Li, Y. Zhang, G. Cao, Z. Deng, G. Huang, W. Li et al., _Event vertex and time reconstruction in large-volume liquid scintillator detectors_ , _Nuclear Science and Techniques_ 32 (2021) 1.
  * [14] CHOOZ collaboration, _Search for neutrino oscillations on a long baseline at the CHOOZ nuclear power station_ , _Eur. Phys. J. C_ 27 (2003) 331 [hep-ex/0301017].
  * [15] Borexino collaboration, _Borexino calibrations: hardware, methods, and results_ , _Journal of Instrumentation_ 7 (2012) P10018–P10018.
  * [16] O. Tajima, _Measurement of electron anti-neutrino oscillation parameters with a large volume liquid scintillator detector, KamLAND_ , 2003.
  * [17] C. Galbiati and K. McCarty, _Time and space reconstruction in optical, non-imaging, scintillator-based particle detectors_ , _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ 568 (2006) 700.
  * [18] RENO collaboration, _RENO: An Experiment for Neutrino Oscillation Parameter $\theta_{13}$ Using Reactor Neutrinos at Yonggwang_, arXiv:1003.1391.
  * [19] G. Huang, Y. Wang, W. Luo, L. Wen, Z. Yu, W. Li et al., _Improving the energy uniformity for large liquid scintillator detectors_ , _Nucl. Instrum. Meth. A_ 1001 (2021) 165287 [2102.03736].
  * [20] Z. Qian, V. Belavin, V. Bokov, R. Brugnera, A. Compagnucci, A. Gavrikov et al., _Vertex and energy reconstruction in JUNO with machine learning methods_ , _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ (2021) 165527.
  * [21] T. Hastie, R. Tibshirani and J. Friedman, _The Elements of Statistical Learning: Data Mining, Inference, and Prediction_. Springer New York, 2013.
  * [22] J. A. Nelder and R. W. M. Wedderburn, _Generalized linear models_ , _Journal of the Royal Statistical Society. Series A (General)_ 135 (1972) 370.
  * [23] C. Davino, M. Furno and D. Vistocco, _Quantile Regression: Theory and Applications_. Wiley, 2013.
  * [24] R. J. Noll, _Zernike polynomials and atmospheric turbulence_ , _JOsA_ 66 (1976) 207.
  * [25] S. Agostinelli, J. Allison, K. Amako, J. Apostolakis, H. Araujo, P. Arce et al., _Geant4—a simulation toolkit_ , _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ 506 (2003) 250\.
  * [26] Z. Wang, Y. Wang, Z. Wang, S. Chen, X. Du, T. Zhang et al., _Design and analysis of a 1-ton prototype of the Jinping Neutrino Experiment_ , _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ 855 (2017) 81.
  * [27] Á. González, _Measurement of areas on a sphere using Fibonacci and latitude–longitude lattices_ , _Mathematical Geosciences_ 42 (2010) 49.
  * [28] Z. Guo, M. Yeh, R. Zhang, D.-W. Cao, M. Qi, Z. Wang et al., _Slow liquid scintillator candidates for MeV-scale neutrino experiments_ , _Astroparticle Physics_ 109 (2019) 33.
  * [29] JNE collaboration, _Measurement of muon-induced neutron production at China Jinping Underground Laboratory_ , _Chinese Physics C_ (2022) .
  * [30] D. Kraft, _A software package for sequential quadratic programming_. Tech Rep DFVLR-FB 88-28, 1988.
  * [31] JUNO collaboration, _TAO Conceptual Design Report: A Precision Measurement of the Reactor Antineutrino Spectrum with Sub-percent Energy Resolution_ , arXiv:2005.08745.
  * [32] XMASS collaboration, _XMASS detector_ , _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ 716 (2013) 78.
  * [33] XMASS collaboration, _Development of low radioactivity photomultiplier tubes for the XMASS-I detector_ , _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ 922 (2019) 171\.
  * [34] W. H. Richardson, _Bayesian-based iterative method of image restoration_ , _JoSA_ 62 (1972) 55.
  * [35] D. C. Xu, B. D. Xu, E. J. Bao, Y. Y. Wu, A. Q. Zhang, Y. Y. Wang et al., _Towards the ultimate PMT waveform analysis for neutrino and dark matter experiments_ , _Journal of Instrumentation_ 17 (2022) P06040.
  * [36] S. Kullback and R. A. Leibler, _On information and sufficiency_ , _The annals of mathematical statistics_ 22 (1951) 79.
  * [37] A. Basu, I. R. Harris, N. L. Hjort and M. C. Jones, _Robust and efficient estimation by minimising a density power divergence_ , _Biometrika_ 85 (1998) 549.
  * [38] Borexino collaboration, _Light propagation in a large volume liquid scintillator_ , _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ 440 (2000) 360\.