# Atmospheric Neutrino Oscillation Analysis With Improved Event Reconstruction
in Super-Kamiokande IV

M. Jiang18  K. Abe1  C. Bronner1  Y. Hayato1  M. Ikeda1  K. Iyogi1  J. Kameda1
Y. Kato1  Y. Kishimoto1  Ll. Marti1  M. Miura1  S. Moriyama1  T. Mochizuki1
M. Nakahata1  Y. Nakajima1  Y. Nakano1  S. Nakayama1  T. Okada1  K. Okamoto1
A. Orii1  G. Pronost1  H. Sekiya1  M. Shiozawa1  Y. Sonoda1  A. Takeda1  A.
Takenaka1  H. Tanaka1  T. Yano1  R. Akutsu2  T. Kajita2  Y. Nishimura2  K.
Okumura2  R. Wang2  J. Xia2  L. Labarga4  P. Fernandez4  F. d. M. Blaszczyk6
C. Kachulis6  E. Kearns6  J. L. Raaf6  J. L. Stone6  S. Sussman6  S. Berkman5
J. Bian7  N. J. Griskevich7  W. R. Kropp7  S. Locke7  S. Mine7  P. Weatherly7
M. B. Smy7  H. W. Sobel7  V. Takhistov7 also at Department of Physics and
Astronomy, UCLA, CA 90095-1547, USA. K. S. Ganezer8  J. Hill8  J. Y. Kim9  I.
T. Lim9  R. G. Park9  B. Bodur10  K. Scholberg10  C. W. Walter10  M. Gonin44
J. Imber44  Th. A. Mueller44  T. Ishizuka11  T. Nakamura12  J. S. Jang13  K.
Choi14  J. G. Learned14  S. Matsuno14  R. P. Litchfield15  Y. Uchida15  M. O.
Wascko15  N. F. Calabria45  M. G. Catanesi45  R. A. Intonti45  E. Radicioni45
G. De Rosa46  A. Ali48  G. Collazuol48  F. Iacob48  L. Ludovici47  S. Cao16
M. Friend16  T. Hasegawa16  T. Ishida16  T. Kobayashi16  T. Nakadaira16  K.
Nakamura16  Y. Oyama16  K. Sakashita16  T. Sekiguchi16  T. Tsukamoto16  KE.
Abe17  M. Hasegawa17  Y. Isobe17  H. Miyabe17  T. Sugimoto17  A. T. Suzuki17
Y. Takeuchi17  Y. Ashida18  T. Hayashino18  S. Hirota18  T. Kikawa18  M.
Mori18  KE. Nakamura18  T. Nakaya18  R. A. Wendell18  L. H. V. Anthony19  N.
McCauley19  A. Pritchard19  K. M. Tsui19  Y. Fukuda20  Y. Itow21  M. Murrase21
P. Mijakowski23  K. Frankiewicz23  C. K. Jung24  X. Li24  J. L. Palomino24  G.
Santucci24  C. Vilela24  M. J. Wilking24  C. Yanagisawa24 also at BMCC/CUNY,
Science Department, New York, New York, USA. D. Fukuda25  K. Hagiwara25  H.
Ishino25  S. Ito25  Y. Koshio25  M. Sakuda25  Y. Takahira25  C. Xu25  Y.
Kuno26  C. Simpson27  D. Wark27  F. Di Lodovico28  B. Richards28  S. Molina
Sedgwick28  R. Tacik29  S. B. Kim30  M. Thiesse31  L. Thompson31  H. Okazawa32
Y. Choi34  K. Nishijima41  M. Koshiba35  M. Yokoyama36  A. Goldsack3  K.
Martens3  M. Murdoch3  B. Quilain3  Y. Suzuki3  M. R. Vagins3  M. Kuze37  Y.
Okajima37  T. Yoshida37  M. Ishitsuka38  J. F. Martin39  C. M. Nantais39  H.
A. Tanaka39  T. Towstego39  M. Hartz40  A. Konaka40  P. de Perio40  S. Chen42
L. Wan42  A. Minamino43
(The Super-Kamiokande Collaboration)  1  2  3  4  5  6  7  8  9  10  11  12
13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32
33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48
<EMAIL_ADDRESS>Kamioka Observatory, Institute for Cosmic Ray
Research, University of Tokyo, Kamioka, Gifu 506-1205, Japan Research Center
for Cosmic Neutrinos, Institute for Cosmic Ray Research, University of Tokyo,
Kashiwa, Chiba 277-8582, Japan Kavli Institute for the Physics andMathematics
of the Universe (WPI), The University of Tokyo Institutes for Advanced
Study,University of Tokyo, Kashiwa, Chiba 277-8583, Japan Department of
Theoretical Physics, University Autonoma Madrid, 28049 Madrid, Spain
Department of Physics and Astronomy, University of British Columbia,
Vancouver, BC, V6T1Z4, Canada Department of Physics, Boston University,
Boston, MA 02215, USA Department of Physics and Astronomy, University of
California, Irvine, Irvine, CA 92697-4575, USA Department of Physics,
California State University, Dominguez Hills, Carson, CA 90747, USA Department
of Physics, Chonnam National University, Kwangju 500-757, Korea Department of
Physics, Duke University, Durham NC 27708, USA Junior College, Fukuoka
Institute of Technology, Fukuoka, Fukuoka 811-0295, Japan Department of
Physics, Gifu University, Gifu, Gifu 501-1193, Japan GIST College, Gwangju
Institute of Science and Technology, Gwangju 500-712, Korea Department of
Physics and Astronomy, University of Hawaii, Honolulu, HI 96822, USA
Department of Physics, Imperial College London , London, SW7 2AZ, United
Kingdom High Energy Accelerator Research Organization (KEK), Tsukuba, Ibaraki
305-0801, Japan Department of Physics, Kobe University, Kobe, Hyogo 657-8501,
Japan Department of Physics, Kyoto University, Kyoto, Kyoto 606-8502, Japan
Department of Physics, University of Liverpool, Liverpool, L69 7ZE, United
Kingdom Department of Physics, Miyagi University of Education, Sendai, Miyagi
980-0845, Japan Institute for Space-Earth Environmental Research, Nagoya
University, Nagoya, Aichi 464-8602, Japan Kobayashi-Maskawa Institute for the
Origin of Particles and the Universe, Nagoya University, Nagoya, Aichi
464-8602, Japan National Centre For Nuclear Research, 00-681 Warsaw, Poland
Department of Physics and Astronomy, State University of New York at Stony
Brook, NY 11794-3800, USA Department of Physics, Okayama University, Okayama,
Okayama 700-8530, Japan Department of Physics, Osaka University, Toyonaka,
Osaka 560-0043, Japan Department of Physics, Oxford University, Oxford, OX1
3PU, United Kingdom School of Physics and Astronomy, Queen Mary University of
London, London, E1 4NS, United Kingdom Department of Physics, University of
Regina, 3737 Wascana Parkway, Regina, SK, S4SOA2, Canada Department of
Physics, Seoul National University, Seoul 151-742, Korea Department of Physics
and Astronomy, University of Sheffield, S3 7RH, Sheffield, United Kingdom
Department of Informatics inSocial Welfare, Shizuoka University of Welfare,
Yaizu, Shizuoka, 425-8611, Japan STFC, Rutherford Appleton Laboratory, Harwell
Oxford, and Daresbury Laboratory, Warrington, OX11 0QX, United Kingdom
Department of Physics, Sungkyunkwan University, Suwon 440-746, Korea The
University of Tokyo, Bunkyo, Tokyo 113-0033, Japan Department of Physics,
University of Tokyo, Bunkyo, Tokyo 113-0033, Japan Department of Physics,Tokyo
Institute of Technology, Meguro, Tokyo 152-8551, Japan Department of Physics,
Faculty of Science and Technology, Tokyo University of Science, Noda, Chiba
278-8510, Japan Department of Physics, University of Toronto, ON, M5S 1A7,
Canada TRIUMF, 4004 Wesbrook Mall, Vancouver, BC, V6T2A3, Canada Department of
Physics, Tokai University, Hiratsuka, Kanagawa 259-1292, Japan Department of
Engineering Physics, Tsinghua University, Beijing, 100084, China Faculty of
Engineering, Yokohama National University, Yokohama, 240-8501, Japan Ecole
Polytechnique, IN2P3-CNRS, Laboratoire Leprince-Ringuet, F-91120 Palaiseau,
France Dipartimento Interuniversitario di Fisica, INFN Sezione di Bari and
Università e Politecnico di Bari, I-70125, Bari, Italy Dipartimento di Fisica,
INFN Sezione di Napoli and Università di Napoli, I-80126, Napoli, Italy INFN
Sezione di Roma and Università di Roma “La Sapienza”, I-00185, Roma, Italy
Dipartimento di Fisica, INFN Sezione di Padova and Università di Padova,
I-35131, Padova, Italy

###### Abstract

A new event reconstruction algorithm based on a maximum likelihood method has
been developed for Super-Kamiokande. Its improved kinematic and particle
identification capabilities enable the analysis of atmospheric neutrino data
in a detector volume 32% larger than previous analyses and increases
sensitivity to the neutrino mass hierarchy. Analysis of a 253.9 kton-year
exposure of the Super-Kamiokande IV atmospheric neutrino data has yielded a
weak preference for the normal hierarchy, disfavoring the inverted hierarchy
at 74% assuming oscillations at the best fit of the analysis.

##  1 Introduction

Neutrino oscillations have been confirmed by a variety of experiments using
both natural and artificial sources. At present the data are well described
assuming mixing among all three active neutrinos with the Pontecorvo-Maki-
Nakagawa-Sakata (PMNS) formalism Pontecorvo:1967fh ; Maki:1962mu . Though most
of its parameters have been experimentally measuredOlive:2016xmw , the
ordering of the mass states with the largest splitting (known as the mass
hierarchy), the octant of the atmospheric mixing angle $\theta_{23}$, and the
value of its CP-violating phase are unknown. These unresolved issues have been
at the forefront of results from the T2K Abe:2017uxa , NOvA Adamson:2017gxd
and Super-Kamiokande Abe:2017osc  (Super-K, SK) experiments and are the focus
of next-generation experiments planned in the U.S. Acciarri:2016crz , China
An:2015jdp , and Japan Abe:2018uyc .

The atmospheric neutrino flux provides a source of both neutrinos and
antineutrinos with a wide variety of energies and path lengths suitable to
probe each of these open questions. Matter effects influence the oscillations
of neutrinos passing through the Earth and produce a resonance enhancement of
the oscillation probability that depends on the value of $\theta_{13}$,
$\theta_{23}$, and $\delta_{CP}$. Importantly this resonance occurs only for
neutrinos if the hierarchy is normal (largest splitting is between the two
heaviest states) and only for antineutrinos if it is inverted. The hierarchy
signal manifests most strongly as an upward-going excess of electron neutrino
or antineutrino events with multi-GeV energies, a region where interactions
frequently have multi-particle final states that complicate particle
identification and where the flux is orders of magnitude lower than its peak.
Accordingly, recent Super-Kamiokande results Abe:2017osc  have been limited by
both a lack of statistics and event miscategorization in the signal region.

A new event reconstruction algorithm has been developed based on a maximum
likelihood method designed to extract specific event topologies and determine
the best set of kinematic parameters. The new algorithm demonstrates improved
reconstruction performance across a variety of metrics including vertex
resolution, particle identification probability, and momentum resolution and
does so across a larger volume of the Super-K detector than its predecessor.
As a result, it is expected to improve Super-K’s sensitivity to the remaining
parameters of the PMNS mixing paradigm.

This paper presents an analysis of 3118.5 days of SK-IV atmospheric data for a
253.9 kiloton$\cdot$year exposure with a 32% larger fiducial volume (FV) that
previous Super-K analyses. Section 2 briefly reviews atmospheric neutrino
oscillations before the detector is introduced in Section 3. In Section 4 the
atmospheric neutrino data reduction and event categorization are discussed.
The new reconstruction algorithm and its performance in comparison to the
conventional reconstruction are detailed in Section 5. Based on this
performance improvement the expansion of the FV and corresponding systematic
errors are presented in Sections 6 and 7, respectively. An analysis of the
atmospheric neutrino data by themselves is presented in Section 8.1 and is
followed by an analysis in which $\theta_{13}$ is constrained by the world
data in Section 8.2, before concluding in Section 9.

##  2 Neutrino Oscillations

Neutrino oscillations are a result of the neutrino mass eigenstates differing
from their weak-interaction (flavor) eigenstates. For neutrinos in vacuum,
their oscillation probability in the standard three-flavor paradigm can be
described by six parameters, nearly all of which have been confirmed to have
non-zero values based on the results of reactor, atmospheric, solar, and long-
baseline neutrino experiments Olive:2016xmw . For brevity the details of the
PMNS oscillation formalism will be omitted here, more details about the
treatment of neutrino oscillation in Super-K are at Abe:2017osc .

In vacuum the leading terms in the oscillations of atmospheric $\nu_{e}$ and
$\nu_{\mu}$ in can be written as:

$\displaystyle P(\nu_{e}\rightarrow\nu_{e})$ $\displaystyle\cong$
$\displaystyle 1-\sin^{2}2\theta_{13}\sin^{2}\left(\frac{1.27\Delta
m_{32}^{2}L}{E}\right)$ $\displaystyle P(\nu_{\mu}\rightarrow\nu_{\mu})$
$\displaystyle\cong$ $\displaystyle
1-4\cos^{2}\theta_{13}\sin^{2}\theta_{23}(1-\cos^{2}\theta_{13}\sin^{2}\theta_{23})$
$\displaystyle\times\sin^{2}\left(\frac{1.27\Delta m_{32}^{2}L}{E}\right)$
$\displaystyle P(\nu_{\mu}\leftrightarrow\nu_{e})$ $\displaystyle\cong$
$\displaystyle\sin^{2}\theta_{23}\sin^{2}2\theta_{13}\sin^{2}\left(\frac{1.27\Delta
m_{32}^{2}L}{E}\right)$ (1)

using the assumption that $\Delta m_{31}^{2}\approx\Delta m_{32}^{2}$ since
$\Delta m_{21}^{2}$ is known to be small in comparison. In these equations the
neutrino travel length is represented by $L$ (km) and its energy by $E$ (GeV).
The unit of $\Delta m^{2}$ is eV2. Atmospheric neutrinos are sensitive to
$\sin^{2}2\theta_{13}$ and the mass splitting $\Delta m_{32}^{2}$ via
$\nu_{\mu}\leftrightarrow\nu_{e}$ oscillations, which depend on
$\sin^{2}\theta_{23}$.

When neutrinos propagate in matter, the oscillation probabilities are modified
due to interactions with the electrons in the earth. Indeed the forward
scattering amplitude of $\nu_{e}$ differs from those of the other flavors,
producing an effective potential felt only by this species. This phenomenon is
known as the matter effect or the MSW effect of neutrino oscillations
Wolfenstein:1978msw ; Mikheev:1985msw . When neutrinos travel through
homogeneous matter, the oscillation parameters $\sin^{2}\theta_{13}$ and
$\Delta m_{32}^{2}$ in Equation 1 are replaced by their matter-equivalents
Giunti:1997MSW ,

$\sin^{2}2\theta_{13,M}=\frac{\sin^{2}2\theta_{13}}{(\cos
2\theta_{13}-A_{CC}/\Delta m^{2}_{32})^{2}+\sin^{2}2\theta_{13}}$ (2) $\Delta
m_{32,M}^{2}=\Delta m_{32}^{2}\sqrt{(\cos 2\theta_{13}-A_{CC}/\Delta
m^{2}_{32})^{2}+\sin^{2}2\theta_{13}}$ (3)

where $A_{CC}=\pm 2\sqrt{2}G_{F}N_{e}E_{\nu}$ and the electron density $N_{e}$
is assumed to be constant. Here $G_{F}$ represents the Fermi constant.

Resonant enhancement of these matter variables occurs when $A_{CC}/\Delta
m^{2}_{32}=\cos 2\theta_{13}$ and is dependent upon the sign of $\Delta
m_{32}^{2}$ and whether the neutrino is a particle ($A_{CC}>0$) or an
antiparticle ($A_{CC}<0$). Note that the enhancement only occurs for neutrinos
if the hierarchy is normal ($\Delta m_{32}^{2}>0$) whereas it only happens for
antineutrinos if it is inverted ($\Delta m_{32}^{2}<0$). Though atmospheric
neutrinos generally experience a varying matter profile, and hence electron
density $N_{e}$ changes as they travel through the earth, this phenomena is
nonetheless present. The calculation of oscillation probability in this
analysis takes such variation on matter density into consideration, with a
simplified version of the preliminary reference Earth model (PREM) (c.f.
Abe:2017osc ). Typically neutrinos with a few to ten GeV of energy that travel
through the Earth’s core experience matter effects the most strongly.

##  3 The Super-Kamiokande Detector

Super-Kamiokande is a cylindrical water Cherenkov detector with a total volume
of 50 kilotons and located inside the Kamioka mine in Gifu Prefecture, Japan.
It is optically separated into two regions, an inner detector (ID), which
forms the primary neutrino target and has a total volume of 32 kilotons and an
outer detector (OD), which is a two-meter cylindrical shell surrounding the ID
and used primarily as a veto. The ID is instrumented with more than 11,000
inward-facing 20-inch photomultiplier tubes (PMTs), representing 40%
photocathode coverage of the target volume. A total of 1,885 8-inch PMTs
coupled to wavelength-shifting plates are mounted on the inner surface of the
OD, while its outer surface is covered with reflective sheets to increase
light collection.

Since the start of operations in 1996, Super-Kamiokande has gone through four
data-taking periods, SK-I, -II, -III, and -IV. The present work focuses on the
latter, which started in September 2008 and ended on May 31${}^{\text{st}}$
2018 when the detector began upgrade work for its next phase. Though the
topology of the detector in the SK-III and -IV periods are the same, at the
start of SK-IV the front end electronics were upgraded to a system based on
ASIC that uses a high-speed charge-to-time converter Nishino:2009zu . After
this upgrade Super-Kamiokande is able to collect all PMT hits above threshold
without incurring any dead time. More detailed descriptions of the detector,
its electronics are presented in Fukuda:2002uc ; Nishino:2009zu ; Abe:2013gga
.

Atmospheric neutrino interactions in the detector are simulated using the
Honda et. al flux calculation Honda:2011nf  and the NEUT Hayato:2002sd
neutrino interaction software (version 5.4.0). Particles emerging from the
interactions are tracked through a simulation of the detector based on GEANT3
Brun:1994aa . The present work uses an updated version of NEUT relative to the
previous atmospheric neutrino analysis (c.f. Abe:2017osc ). In particular,
charged-current quasi-elastic (CCQE) interactions in the new simulation are on
the local Fermi-Gas model of Nieves Nieves:2013wx ; Nieves:2011wx , assuming
an axial mass $M_{A}=1.05~{}\mbox{GeV/c}^{2}$ and using the random phase
approximation correction. Further, deep inelastic scattering (DIS) is modeled
using the GRV98 parton distribution function Gluck:1998xa  and utilizing CKM
matrix elements for the calculation of structure functions. Corrections for
low $q^{2}$ scattering have been updated to those of Bodek and Yang
Bodek:2005de , where $q^{2}$ denotes the square of the transferred four-
momentum.

##  4 Event sample and reduction

The current analysis utilizes atmospheric neutrino data collected during the
SK-IV period with a total livetime of 3118.45 days. As in previous Super-K
analyses, the atmospheric neutrino data are separated into three general
categories, fully contained (FC), partially contained (PC), and upward-going
muons (Up-$\mu$).

For events classified as FC and PC, the neutrino interacts within the fiducial
volume, defined as the region located more than 2 m from the ID wall in
previous analyses. Events with no activity in the outer detector are
classified as FC. If energy deposition in the OD is observed, generally from a
high energy muon exiting the ID, the event is classified as PC. Muons created
by neutrino interactions in the rock around SK or in the OD water and
traveling upward through the detector form the Up-$\mu$ sample.

At trigger level the Super-Kamiokande event sample consists mainly of
downward-going cosmic ray muons and low energy radioactivity from contaminants
in the water such as radon. To remove these backgrounds, reduction processes
specific to the three main sample categories are applied to the data, whose
detail can be found in Ashie:2005osc .

As an example, five steps of data reduction criteria are used for FC sample.
In the first and second reduction steps, criteria on the total charge
collected by ID and the number of hits on OD are applied to remove most cosmic
ray muons. In the third reduction step, a dedicated cosmic ray muon
reconstruction algorithm is used to apply a cut on OD hits near the entrance
or exit point of a reconstructed track. The fourth step is designed to remove
so-called “flasher” events caused by internal corona discharge of a PMT.
Failing tubes often produce multiple events with similar topologies as these
discharges repeat. Typical flasher events have a broad distribution in time
among hit PMTs and spatially similar hit distributions. Cuts based on the hit
timing distribution and the charge pattern correlation with other events that
passed the third step reduction are applied to remove flasher events. The
fifth reduction rejects the remaining cosmic ray muons, flasher events and
electronic noise further using a more precise fitter. Finally, the
reconstructed vertex is required to be within a specified fiducial volume to
further reduce non-neutrino background events and to prevent the
reconstruction performance from deteriorating near the detector wall. These
steps are particularly relevant to the new reconstruction algorithm considered
here and will be discussed in Section 6 in detail. An event’s visible energy
($E_{\mbox{vis}}$), which is defined as the energy of an electromagnetic
shower producing the same amount of Cherenkov light as observed in the event,
is also required to be larger than 30 MeV to remove low energy backgrounds.
All events passing the reduction are scanned by eye to determine the level of
background contamination in the final analysis sample and to estimate the
uncertainty inherent to the reduction process.

After reduction, the FC data are sub-divided based upon the number of observed
Cherenkov rings, the particle ID (PID) of the most energetic ring, and the
number of observed electrons from muon decays into combinations of single- or
multi-ring, electron-like ($e$-like) or muon-like ($\mu$-like), and sub-GeV
($E_{\mbox{vis}}<1330.0$ MeV) or multi-GeV ($E_{\mbox{vis}}>1330.0$ MeV)
categories, as well as the number of decay electrons (0, 1, or 2 or more). In
the present study, all FC events have been reconstructed using the new
algorithm, which will be introduced in the next section.

The PC and Up-$\mu$ samples also have their own reduction processes, which are
optimized for the topology of OD activity of neutrino events. These two
samples are reconstructed using the pre-existing algorithm and are divided
into “stopping” and “through-going” subsamples based on the estimated muon
stopping point for PC and Up-$\mu$ events. The “through-going” events in the
Up-$\mu$ sample are further divided in to “showering” and “non-showering”
based on whether the event induces an electromagnetic shower while traversing
the ID. Although Super-K also works as the far detector for T2K experiment
Abe:2013nka  and detects accelerator-generated neutrinos, they are excluded
based on event time information in this analysis. After all selections there
are a total of 13 FC analysis samples, two PC samples and three Up-$\mu$
samples.

##  5 FC Event Reconstruction

The Super-K event reconstruction algorithms determine an event’s physical
properties such as the interaction vertex, number of particles, the particle
types and momenta based on PMT hit information.

The conventional event reconstruction algorithm, APFit, was introduced at the
beginning of Super-K, contributing to the discovery of atmospheric neutrino
oscillation and has been used in both the K2K and T2K experiments. APFit is a
single iteration fitter based on the time and charge information of hit PMTs
Shiozawa:1999ap . The interaction vertex is reconstructed based on hit timing
information after accounting for the photon time of flight. Then the direction
of the first found ring, usually also the most energetic ring, is determined
based on the observed charge distribution with respect to the interaction
vertex. Additional ring candidates are found using a Hough-transform based
method and selected by a likelihood function optimized to reject spurious ring
candidates. After determining the number of Cherenkov rings in the event, the
particle type of each ring is determined based upon the Cherenkov ring pattern
and opening angle. Rings from electrons tend to have rough edges produced by
the light from their electromagnetic showers, while rings from muons or
charged pions predominately produce crisp edges. In the last step of the fit,
the momentum of each ring is evaluated based on the observed charge inside
cone with a half angle of $70^{\circ}$ drawn along the line connecting the
interaction vertex to each ring’s center. Corrections and adjustments are made
to account for charge sharing between overlapping rings. It is worth noting
that the hit time information is only used during the first step to find the
initial vertex candidate in APFit.

The new algorithm, named fiTQun, employs a maximum likelihood method to
reconstruct particle types and determine kinematics in the detector
simultaneously. The algorithm is based on methods developed for the MiniBooNE
experimentPat:2009mini , but has been developed from scratch for Super-K with
additional features such as multi-ring reconstruction for events with multiple
final-state particles. Compared to APFit, fiTQun uses more information,
including information from PMT hits outside of the expected Cherenkov cone and
hit timing information, during the fitting procedure. For a given event
fiTQun’s fit procedure will run multiple times to determine the best kinematic
parameters for each possible particle configuration hypothesis, while APFit
only fits those parameters once. The remarkable evolution of computing power
since the start of Super-Kamiokande has enabled fiTQun to achieve higher
reconstruction precision on similar time scales as APFit used to be. FiTQun
has already been used in the T2K analyses using an expanded fiducial volume
due its improved resolution of reconstructed quantities and particle
identification Abe:2018wpn .

###  5.1 Likelihood Function

An event topology hypothesis $\Gamma$ (e.g. single-ring $e$-like) together
with its associated kinematic parameters $\theta$, which include the vertex
position, particle creation times, the azimuthal and zenith angles of the
particle directions, as well as their momenta are considered in the likelihood
function during a fit. Based on the observed charge and hit time of each PMT,
fiTQun constructs the following likelihood function for a given the hypothesis
to estimate the kinematic variables:

$\displaystyle L(\Gamma,\theta)$ $\displaystyle=$
$\displaystyle\prod^{\text{unhit}}_{j}P_{j}(\text{unhit}|\Gamma,\theta)\prod^{\text{hit}}_{i}\\{1-P_{i}(\text{unhit}|\Gamma,\theta)\\}$
(4) $\displaystyle\times
f_{q}(q_{i}|\Gamma,\theta)f_{t}(t_{i}|\Gamma,\theta).$

In this equation, the index $j$ runs over all PMTs which did not register a
hit (“unhit” PMTs) and for each of these the probability that it does not
register a hit given the fitting hypothesis ($\Gamma,\theta$) is calculated as
$P_{j}(\text{unhit}|\Gamma,\theta)$. PMTs which did register a hit are indexed
with $i$. For such PMTs the likelihood density for observing a charge $q_{i}$
under the fitting hypothesis is represented by charge likelihood
$f_{q}(q_{i}|\Gamma,\theta)$. The likelihood density of producing a hit at the
observed time $t_{i}$ is defined similarly as $f_{t}(t_{i}|\Gamma,\theta)$.

Since the processes of particle and optical photon propagation are decoupled
from the response of the PMT and the electronics, the charge likelihood can be
rewritten in terms of the expected number of photoelectrons produced at the
$i$-th PMT given the hypothesis (the predicted charge),
$\mu_{i}(\Gamma,\theta)$, as

$\displaystyle L(\Gamma,\theta)$ $\displaystyle=$
$\displaystyle\prod^{\text{unhit}}_{j}P_{j}(\text{unhit}|\mu_{j})\prod^{\text{hit}}_{i}\\{1-P_{i}(\text{unhit}|\mu_{i})\\}$
(5) $\displaystyle\times f_{q}(q_{i}|\mu_{i})f_{t}(t_{i}|\Gamma,\theta)$

$P_{j}(\text{unhit}|\mu_{j})$ and $f_{q}(q_{i}|\mu_{i})$ are properties of the
PMT and the electronics, and therefore do not explicitly depend on the process
of Cherenkov photon emission and propagation in water.

In the calculation of the predicted charge, the contributions from direct
light and light that has scattered or been reflected (indirect light) are
considered separately and summed to form the final $\mu_{i}$. The predicted
charge from direct light reaching a PMT is calculated by integrating the
Cherenkov emission profile along the track while correcting for the distance
from the track to the PMT, the light transmission in water, and the PMT
angular acceptance. Charge produced by indirect light reaching the PMT is
predicted by integrating the product of the direct light emission profile and
a scattering function that has been generated in advance based on simulation
and incorporates effects arising from the relative position of the PMT and
light source. For events with multiple Cherenkov rings the predicted charge of
each ring is first calculated separately and then the sum from all the rings
is used to calculate the total expected $\mu_{i}$. The final charge likelihood
$f_{q}(q_{i}|\mu_{i})$ is obtained by comparing the observed charge in a PMT
against the prediction assuming photoelectrons generated according to Poisson
statistics.

The time likelihood term can be expressed as
$f_{t}(t_{i}|t^{\text{exp}}_{i},\Gamma,p,\mu_{i})$, where $p$ is the momentum
for the topology $\Gamma$ and $t^{\text{exp}}_{i}$ is the expected hit time.
The latter is defined as the arrival time of unscattered photons emitted at
the track midpoint and traveling directly to the PMT as

$\displaystyle
t^{\text{exp}}_{i}=t+s_{\text{mid}}/c+|\bm{R}_{i}^{\text{PMT}}-\bm{x}-s_{\text{mid}}\bm{d}|/c_{n},$
(6)

where $\bm{x}$ and $t$ are the vertex and creation time of the particle,
respectively, $\bm{d}$ is the particle direction, $\bm{R}_{i}^{\text{PMT}}$ is
the position of the $i$-th PMT and $s_{\text{mid}}$ represents half of the
track length. Here $c$ and $c_{n}$ represent the light velocity in vacuum and
in water, respectively. The time likelihood depends on the predicted charge
since a hit is recorded by the first photon arriving at a PMT, which leads a
narrower distribution of hit time for higher numbers of incident photons and
hence, more predicted charge. The track length of a particle $s$, which is
determined by the topology $\Gamma$ and momentum $p$, also affects the shape
of time likelihood since not all photons are generated at the track midpoint.
Contributions to the likelihood from direct and indirect photon hits are
calculated separately and then merged according to their relative intensities
to obtain the final time likelihood. The time likelihoods are determined using
particle gun simulations. For multi-particle hypotheses, the time likelihood
is calculated ring-by-ring and then merged to a final likelihood function
assuming the photons from a particle with earlier $t^{\text{exp}}_{i}$ always
arrive earlier than the photons from any other particles with later
$t^{\text{exp}}_{i}$ values.

Once the likelihood function is defined, the best set of kinematic parameters
$\hat{\theta}$ for a given event topology hypothesis is defined as that which
maximizes $L(\Gamma,\theta)$. The best estimate for the particle content of a
given event is determined by comparison of $L(\Gamma,\hat{\theta})$ among all
hypotheses, $\Gamma$.

###  5.2 Fitting procedure

The fiTQun reconstruction process can be divided into four steps. The first
step, vertex pre-fitting, roughly estimates the interaction vertex based on
the PMT timing information. During the second step, clusters of PMT hits in
time are identified as candidate for particle activity. Thereafter the single-
ring reconstruction, which performs fits assuming event topologies with only a
single light-producting particle, and multi-ring reconstruction, which fits
using hypotheses with multiple particles, are performed in sequence. During
these fits the negative log likelihood -ln $L(\Gamma,\theta)$ is minimized
with respect to $\theta$, using the MINUITJames:1994vla  package.

####  5.2.1 Vertex pre-fitting

The vertex pre-fitter is a fast algorithm which uses only the hit time
information from PMTs around the primary event trigger to estimate an initial
vertex position assuming all observed light was emitted from a common point.
This is done by searching for the vertex position $\bm{x}$ and time $t$ which
maximize the goodness function

$\displaystyle G(\bm{x},t)$ $\displaystyle=$
$\displaystyle\sum^{\text{hit}}_{i}e^{(-(T_{i}^{\text{res}}/\sigma)^{2}/2)},$
(7)

where

$\displaystyle T_{i}^{\text{res}}$ $\displaystyle=$ $\displaystyle
t_{i}-t-|\bm{R}_{i}^{\text{PMT}}-\bm{x}|/c_{n}$ (8)

is the residual hit time calculated on the assumption of a point-like light
source and subtracting the photon time of flight. The value of parameter
$\sigma$ determines the precision of pre-fitting. Here
$\bm{R}_{i}^{\text{PMT}}$ is the position of the $i^{th}$ PMT. When the vertex
and time get close to their true values by doing a grid search, the
$T_{i}^{\text{res}}$ distribute near zero, which results in a large value of
the goodness. The pre-fitter is executed several times with gradually
shrinking the grid size and $\sigma$ to achieve high precision efficiently.
The fitted vertex from this step, called pre-fit vertex, is just a rough
estimation and will be fitted again with higher precision during minimization
of $L(\Gamma,\theta)$.

####  5.2.2 Hit clustering

Events in Super-K are defined by detector activity in an $O(10~{}\mu$s) time
window around an event trigger but may contain multiple subevents representing
clusters of PMT hits separated in time from the primary trigger. A common
example is muon decay in which the muon produces the primary event trigger and
its decay produces additional delayed detector hits that form a subevent. A
hit clustering algorithm is used to search for activity around the primary
trigger and locates any additional subevents for further fitting with the more
precise reconstruction methods discussed below.

The hit clustering algorithm starts by searching in time for subevent activity
around the event trigger using a peak-finding algorithm. The vertex goodness
$G(\bm{x},t)$ from Equation 7 is scanned in $t$ with the vertex position fixed
to the pre-fit vertex to search for additional peaks from particle activity in
the detector. An example of the goodness distribution from a muon decay event
is shown in Figure 1. The two dominant peaks in the distribution correspond to
the parent muon and its Michel electron, respectively. To avoid counting peaks
created by scattering or reflection processes within a cluster of hit
activity, peaks are required to be above a minimum threshold $F(t)$, which is
defined as:

$F(t):=0.25~{}{\operatorname*{argmax}_{i\in
M}}\\{\frac{G(\bm{x},t_{i})}{(1+((t-t_{i})/\gamma)^{2}}\\}+\eta$ (9)

where $M$ represents all local maxima of goodness function $G(\bm{x},t)$. The
value of time constant $\gamma$ is 25 ns when $t<t_{i}$ and 70 ns otherwise.
An offset $\eta$ = 9 is added to the threshold function to suppress the effect
of dark hits. The function $F(t)$ is the blue curve in Figure 1. The minimum
goodness between any two peaks must be lower than a second threshold,
$0.6\times F(t)$, which is shown as green dashed curve. Under these criteria
only the peaks labelled with a triangle are identified as candidates for
further fitting.

Figure 1:  Goodness as a function of hit time for an event with a muon and a
Michel electron. The black line shows the goodness as a function of hit time
with the vertex fixed to the pre-fit vertex. Blue (dashed) and green (dotted)
curves denote thresholds for identifying candidate peaks. Red triangles denote
the time of candidates and gray vertical lines indicate time windows around
those candidates as determined by the algorithm. If such windows overlap, the
candidates are merged into one.

During this procedure all vertex positions are assumed to lie close to the
pre-fit vertex when the peak-finding algorithm runs. This assumption is broken
when the primary particle travels a significant distance from the interaction
vertex, as for high momenta muons. Therefore, the vertex pre-fitting and peak-
finding algorithm are rerun after masking the hits caused by the primary
particle to improve decay electron reconstruction efficiency. The vertex
position used in the goodness function $\bm{x}$ is then close to the vertex of
the secondary particle.

A time window defined as -180 ns $<T_{i}^{\text{res}}<$ 800 ns around each
found peak is defined to contain its associated hits. Afterwards the vertex
pre-fitter and peak finder will be run once again in each time window using
only the hits within each window. Peaks remaining after this step become the
final candidates for full event reconstruction.

####  5.2.3 Single-ring reconstruction

The most basic reconstruction which is applied to the time windows defined in
the previous step is the single-ring fitter. This fitter poses single-particle
hypotheses for the likelihood function in Equation 4. Three types of single-
ring hypothesis are considered in fiTQun, that of an electron, a muon, and a
charged-pion. For each hypothesis the kinematic parameters of the event,
including the interaction vertex, the particle momentum, and its direction are
varied to maximize the likelihood function against the observation. Particle
identification (PID) is based on the best-fit likelihood values for each of
these hypotheses. Electrons and muons, for example, are separated by cutting
on ln($L_{e}/L_{\mu}$), the logarithm of the likelihood ratio between the
best-fit electron and muon hypotheses. The distribution of this variable is
shown for FC atmospheric neutrino data and Monte Carlo (MC) simulation result
with sub-GeV and multi-GeV energies in Figure 2. In both cases there is a
clear separation of the likelihood variable between electron-like ($e$-like)
and muon-like ($\mu$-like) events.

(a) Sub-GeV events (b) Multi-GeV events

Figure 2:  PID likelihood distribution for FC single-ring sub-GeV (left) and
multi-GeV (right) comparing data (points) and atmospheric neutrino MC
(histograms). Neutrino oscillations are taken into account with the normal
hierarchy assumed and the oscillation parameters taken to be $\Delta
m_{23}^{2}=2.4\times 10^{-3}\mbox{eV}^{2}$, $\mbox{sin}^{2}\theta_{23}=0.5$,
$\mbox{sin}^{2}\theta_{13}=0.0210$, and $\delta_{CP}=0$. The shaded histograms
show charged-current $\nu_{\mu}$ interactions. Error bars show the statistical
error. In this figure the reconstructed event vertex is required to be at
least 200 cm away from the ID wall.

(a) True electron event with single ring (b) True muon event with single ring

Figure 3:  Vertex resolution of single-ring electron (left) and muon (right)
events in the FC CCQE event sample in the atmospheric neutrino MC, plotted as
a function of visible energy. The red full triangle indicate the performance
of fiTQun and the black open circle are for APFit. The true event vertex is
required to be at least 200 cm away from the ID wall.

(a) Momentum resolution of true single-electron events. (b) Momentum
resolution of true single-muon events.

(c) Momentum bias of true single-electron events. (d) Momentum bias of true
single-muon events

Figure 4:  Momentum resolution (top) and bias (bottom) of single-ring electron
(left) and muon (right) events in the FC CCQE event sample in the atmospheric
neutrino MC, plotted as a function of visible energy. The red full triangle
indicate the performance of fiTQun and the black open circle are for APFit.
The reconstructed event vertex is required to be at least 200 cm away from the
ID wall.

(a) True single-electron events. (b) True single-muon events.

Figure 5:  Mis-identification rate of single-ring electron (left) and muon
(right) events from FC CCQE atmospheric neutrino MC events plotted as a
function of visible energy. The open circles indicate the performance of APFit
and the triangles are for fiTQun. The reconstructed event vertex is required
to be at least 200 cm away from the ID wall.

Figure 3 shows the vertex resolution as the function of visible energy for the
FC single-ring charged-current quasi-elastic (CCQE) event sample in the
atmospheric neutrino MC whose true interaction vertex lies within the
conventional fiducial volume, the region located more than 2 m from the ID
wall. As the visible energy increases from 100 MeV to 1330 MeV the vertex
resolution of fiTQun for CCQE $\nu_{e}$ events is stable at 20.64 cm , while
the resolution of APFit varies from 34.63 cm to 25.26 cm. For the CCQE
$\nu_{\mu}$ events, the vertex resolution of fiTQun, which changes from 29.24
cm to 15.89 cm, is better than that of APFit, which varies from 34.21 cm to
18.16 cm. In the same energy range, the momentum resolution of fiTQun for CCQE
$\nu_{e}$ events improves from 5.39% to 2.58% as the visible energy increases,
while the resolution of APFit changes from 7.04% to 3.32%, as shown in the top
left plot of Figure 4. For CCQE numu events the momentum resolution is
constant across the energy range, being lower than $2.5\%$ for fiTQun and
slightly worse for APFit. Similarly fiTQun shows an improved ability to
discriminate between electrons and muons with less than a 1% mis-
identification rate for visible energies less than 1330 MeV (Figure 5). Other
metrics used to measure fit quality are generally stable and their typical
values can be found in Table 1. In general fiTQun performs as well or better
than APFit.

Table 1:  Summary of the basic performance of the APFit and fiTQun reconstruction algorithms on the fully contained CCQE single ring event sample with visible energy of 1 GeV.  Reconstruction | fiTQun | APFit
---|---|---
True CCQE $\nu_{e}$ sample
Vertex Resolution | 20.64 cm | 24.86 cm
Direction Resolution | 1.48∘ | 1.68∘
Momentum Bias | 0.43% | 0.63%
Momentum Resolution | 2.90% | 3.56%
Mis-PID rate | 0.02% | 0.50%
True CCQE $\nu_{\mu}$ sample
Vertex Resolution | 15.83 cm | 17.32 cm
Direction Resolution | 1.00∘ | 1.28∘
Momentum Bias | -0.18% | 0.54%
Momentum Resolution | 2.26% | 2.60%
Mis-PID rate | 0.05% | 0.91%

####  5.2.4 Multi-ring reconstruction

In atmospheric neutrino analyses it is essential to reconstruct events with
multiple light-producing particles since a large fraction of events with
multi-GeV energies have multi-particle final states. To save computing time,
the fiTQun multi-ring fitter is applied only to the time window around the
primary event trigger and not to any additional subevents. The process of
reconstructing a multi-ring event starts by performing an iterative search for
an additional ring on top of any existing rings from a previous fit result.
The hypothesis in the likelihood from Equation 5 is updated to include a new
ring and minimized again allowing the kinematic parameters for the rings to
vary. Three hypotheses, including $e$-like, $\mu$-like and $\pi^{+}$-like
rings are tested. The result of the updated fit for each of the particle
hypotheses is compared to the original result to determine the validity of the
added ring. If the likelihood improves with the new ring it is accepted and
the processes to search for further rings is repeated. The cycle of adding,
fitting and examining new rings iterates until either a newly added ring fails
the likelihood criterion or six rings are found in the event. Figure 6 shows
the difference in the best-fit likelihoods from the 2-ring and a 1-ring
hypotheses. Based on an optimization performed in MC, a cut at 9.35 (11.83) on
this likelihood ratio is used to separate single- and multi-ring event when
the first ring hypothesis is $e$-like ($\mu$-like). However, ring candidates
that have an angular separation of less than $20^{\circ}$ from the most
energetic ring are discarded as spurious since MC studies indicate such
candidates are typically due to the particle scattering rather than a new
particle. When this occurs the two rings are merged and refit as one for all
particle hypotheses while keeping all other rings fixed. This procedure is
repeated for all rings in the event in descending order of their energy.

Table 2 shows the performance of the ring counting in both APFit and fiTQun
using atmospheric neutrino MC events. In this comparison only charged
particles with energy more than 30 MeV above Cherenkov threshold in the final
states without any requirement on their angular separation with other
particles are considered as true ring candidates. Compared to APFit, fiTQun
shows a greater ability to reconstruct multi-ring events, while the ability to
correctly identify single-ring topologies is the same in both algorithms.
Although fiTQun tends to have more fake rings than that of APFit, the
improvements largely outweigh this slight downside.

Table 2:  Summary table of ring counting performance of both APFit and fiTQun on fully contained atmospheric neutrino events. Columns denote the number of reconstructed rings and rows the number of true rings. The true number of rings is counted using only particles with energy more than 30 MeV above Cherenkov threshold in the final states. The terms 1R, 2R, and 3R correspond to one ring, two rings, and three or more rings. True Number | fiTQun Reconstruction | APFit Reconstruction
---|---|---
of Rings | 1R | 2R |  $\geq$ 3R | 1R | 2R |  $\geq$ 3R
True 1R | 95.0% | 4.64% | 0.41% | 95.9% | 3.85% | 0.29%
True 2R | 27.8% | 66.7% | 5.56% | 42.5% | 52.8% | 4.63%
True $\geq$3R | 7.04% | 25.5% | 67.5% | 20.2% | 33.0% | 46.8%

Figure 7 shows the PID likelihood variable distribution of the most energetic
ring in fully contained multi-ring events. Due to the overlap of Cherenkov
photons from multiple particles, the separation between $e$-like and
$\mu$-like is not as good as that for single-ring events.

(a) Sub-GeV events (b) Multi-GeV events

Figure 6:  The distribution of the likelihood ratio between the best-fit
single-ring hypothesis and multi-ring hypothesis for sub-GeV (left) and multi-
GeV (right) $e$-like events. The points show fully contained atmospheric
neutrino data and the histograms show the MC prediction including neutrino
oscillations. The shaded histograms show events with one a single-ring in the
final state. Error bars show the statistical error. The terms 1R and MR
represent single-ring and multi-ring, respectively. The reconstructed event
vertex is required to be at least 200 cm away from the ID wall.

(a) Sub-GeV events (b) Multi-GeV events

Figure 7:  PID likelihood distributions for the most energetic ring in sub-GeV
(left) and multi-GeV (right) fully contained multi-ring events. Points show
the data and the histograms show the atmospheric neutrino MC including
neutrino oscillations. The shaded histograms show $\nu_{\mu}$ charged-current
interactions. Error bars are statistical. The reconstructed event vertex is
required to be at least 200 cm away from the ID wall.

(a) Event selected by fiTQun (b) Event selected by APFit

Figure 8:  Reconstructed invariant mass calculated using the second and the
third most energetic rings for events with three e-like rings by fiTQun (left)
and APFit (right) from the sub-GeV atmospheric neutrino MC. Shaded histograms
show true CC$\nu_{e}1\pi^{0}$ events and arrows indicate the invariant mass
range used to select events in which a $\pi^{0}$ is reconstructed properly.
The reconstructed event vertex is required to be at least 200 cm away from the
ID wall.

The performance of the fitter’s kinematic reconstruction on multi-ring events
is checked by studying the invariant mass of $\pi^{0}$s in charged-current
single pion (CC1$\pi^{0}$) events. Figure 8 shows the reconstructed invariant
mass calculated using the second and the third most energetic rings for events
with three $e$-like rings ($eee$ events) for both reconstruction algorithms.
Shaded histograms show true CC$\nu_{e}1\pi^{0}$ events. Both fitters have a
peak in the invariant mass distribution near the $\pi^{0}$ mass, but the
fiTQun peak is larger, sharper, and has fewer backgrounds indicating improved
reconstruction efficiency over APFit. The FWHM of the $\pi^{0}$ mass peak
selected by fiTQun is 38 MeV, while this value is 42 MeV for APFit. The event
rate and the purity of the CC$\nu_{e}1\pi^{0}$ events with $85<m<185$ MeV
(indicated by the arrows in the figure) are 278.5 and 52.1% for fiTQun, which
are much higher than those values of APFit: 175.9 and $42.9\%$. Similarly,
CC$\nu_{\mu}1\pi^{0}$ events from the three ring sample whose leading ring is
$\mu$-like and passing this selection are also studied and summarized in Table
3. The fiTQun algorithm demonstrates a higher efficiency and purity for multi-
ring event reconstruction.

Table 3:  Summary table of performance of CC1$\pi^{0}$ events selection in three ring events. The terms $eee$ and $\mu ee$ represent the reconstructed PID of the three rings in order of momentum. The event rate column shows the number of events passing the invariant mass cut to select $\pi^{0}$. The purity column shows the purity of target events in the passed invariant mass region.  Topology | fiTQun selection | APFit selection
---|---|---
Event rate | Purity | Event rate | Purity
$\bm{eee}$ events (target: $\nu_{e}$CC1$\pi^{0}$)
Sub-GeV | 278.5 | 52.1% | 175.9 | 42.9%
Multi-GeV | 112.1 | 47.0% | 32.8 | 36.2%
$\bm{\mu ee}$ events (target: $\nu_{\mu}$CC1$\pi^{0}$)
Sub-GeV | 384.5 | 54.6% | 201.2 | 38.0%
Multi-GeV | 143.5 | 64.6% | 51.6 | 32.0%

Another important performance indicator for multi-ring event reconstruction
algorithm is the separation of neutrino and antineutrino components of the
atmospheric neutrino sample, since at multi-GeV energies they have the most
sensitivity to the mass hierarchy. A two-stage likelihood method has been
developed to purify the neutrino and antineutrino components for multi-ring
events. The first stage of the separation is designed to extract CC
$\nu_{e}+\bar{\nu_{e}}$ interactions based on a likelihood selection as in
Wendell:2010md  using the APFit algorithm. However, in the present study the
inputs have been replaced by the equivalent variables from the fiTQun
reconstruction. Events that pass this selection are classified as “multi-ring
$e$-like” while those that fail are termed “multi-ring other” as shown in
Figure 9. Both are used in the oscillation analysis discussed below.

Table 4:  Summary of the Multi-ring e-like event selection. The efficiency represents the fraction of events passing the selection among all true charged-current $\nu_{e}$ events. Similarly, the purity is the fraction true CC $\nu_{e}$ out of all events passing the selection.  | fiTQun | APFit
---|---|---
First stage
Multi-ring $e$-like events |  |
Efficiency | 75.7% | 69.7%
Purity | 77.8% | 69.5%
Second stage
Multi-ring True CC$\nu_{e}$ events |  |
Efficiency | 56.8% | 53.6%
Purity | 58.8% | 52.6%
Multi-ring True CC$\bar{\nu}_{e}$ events |  |
Efficiency | 68.4% | 70.9%
Purity | 30.0% | 25.9%

Neutrino and antineutrino interactions are separated from the multi-ring
$e$-like sample during the second stage of the selection. The method is
similar to the one discussed in Abe:2017osc  and uses the number of decay
electrons, the number of reconstructed rings, electrons, and the event’s
transverse momentum as reconstructed by fiTQun. Figure 10 shows the final
likelihood distribution used in SK-IV. The efficiency and purity for selecting
$e$-like events and identifying true CC $\bar{\nu}_{e}$ ($\nu_{e}$) events as
$\bar{\nu}_{e}$-like ($\nu_{e}$-like) is summarized in Table 4. FiTQun shows
improved performance in each of these metrics which will translate to better
sensitivity to sin${}^{2}\theta_{13}$ and the mass hierarchy.

Figure 9:  Likelihood distribution used to separate SK-IV multi-ring events
whose most energetic ring is e-like. Error bars represent the statistical
uncertainty of the data. Events with likelihood values larger than -0.25 are
designated multi-ring e-like, while those with lower values are termed multi-
ring other.  Figure 10:  Likelihood distribution used to separate multi-ring
$e$-like events into neutrino-like and antineutrino-like subsamples. Error
bars represent the statistical uncertainty of the data. Events with negative
(positive) likelihood values are designated $\bar{\nu}_{e}$-like
($\nu_{e}$-like).

##  6 Fiducial Volume expansion

In previous analyses using the APFit algorithm the fiducial volume (FV) was
defined as the region located 200 cm from the wall of the inner detector. With
this definition approximately $30\%$ of the Super-K inner detector mass is not
used for atmospheric neutrino analysis. On the other hand, the sensitivity to
several open questions is currently limited by a lack of statistics. The
analysis sensitivity can be increased by taking advantage of fiTQun’s improved
reconstruction performance to expand the fiducial volume while keeping signal
purities high and backgrounds low.

The event selection and categorization in this study is similar to the one
used in Abe:2017osc , but all cut variables used to define the FC samples are
reconstructed by fiTQun. Partially contained and Up-$\mu$ events on the other
hand, are defined using APFit in the same manner as in the previous
publications.

###  6.1 MC study for fiducial volume expansion

Though a larger fiducial volume provides higher statistics, the purity of
signal interactions may decrease if the reconstruction performance
deteriorates for events close to the ID wall. Events that are too close to the
wall and producing particles that travel toward it will only illuminate a few
PMTs and as a result will have poorly imaged Cherenkov rings. Further, as the
FV cut is moved closer to the wall, the possibility of events with interaction
vertices outside of the ID wall being reconstructed within the FV increases.
Although the OD is designed to veto such events, as introduced in Section 4,
events which originate in the 55 cm dead region between the ID and OD optical
barriers or interact very near the OD wall may not produce enough light to
trigger the OD cut and may remain in the analysis sample. Such “entering
events” are a type of mis-reconstructed background in the oscillation analysis
with potentially large systematic errors.

Figure 11 shows the purity of several components of the FC sub-GeV $\mu$-like
0 decay-e sample as a function of the cut on the reconstructed distance from
the event vertex to the nearest ID wall (dwall) as an example. The fraction of
signal events, $\nu_{\mu}$CC and $\bar{\nu}_{\mu}$ CC interactions (red) for
this sample, remains stable as the cut position is decreased to smaller values
until it reaches 50 cm, where the background from entering $\nu$ events
increases dramatically. All of 13 FC event categories exhibit similar
behavior, indicating that a FV cut at 50 cm is an acceptable limit due to
entering $\nu$ events. The purity of the 13 FC samples defined using this cut
and the conventional FV definition are summarized in Table 5 (c.f. Table II of
Ref Abe:2017osc ).

Figure 11:  The fraction of each component of the FC sub-GeV $\mu$-like 0 decay-e event sample in the atmospheric neutrino MC simulation as a function of the reconstructed distance from the interaction vertex to the nearest wall.  Table 5: Sample purity broken down by neutrino flavor assuming neutrino oscillations with $\Delta m^{2}_{32}=2.52\times 10^{-3}\mbox{eV}^{2}$ and $\sin^{2}\theta_{23}=0.51$ with a FV cut of 200 cm. The number in the parentheses is the result within the new region, where the distance to nearest wall is between 50 cm and 200 cm.  Sample | CC $\nu_{e}$ | CC $\bar{\nu}_{e}$ | CC $\nu_{\mu}+\bar{\nu}_{\mu}$ | CC $\nu_{\tau}$ | NC
---|---|---|---|---|---
Fully Contained (FC) Sub-GeV
e-like, Single-ring |  |  |  |  |
0 decay-e | 0.728 | 0.242 | 0.001 | 0.000 | 0.028
| (0.702) | (0.227) | (0.025) | (0.001) | (0.045)
1 decay-e | 0.907 | 0.020 | 0.033 | 0.001 | 0.040
| (0.712) | (0.015) | (0.208) | (0.001) | (0.063)
$\mu$-like, Single-ring |  |  |  |  |
0 decay-e | 0.010 | 0.004 | 0.795 | 0.001 | 0.191
| (0.034) | (0.011) | (0.805) | (0.001) | (0.150)
1 decay-e | 0.000 | 0.000 | 0.974 | 0.000 | 0.026
| (0.001) | (0.000) | (0.968) | (0.000) | (0.031)
2 decay-e | 0.000 | 0.000 | 0.984 | 0.000 | 0.016
| (0.000) | (0.000) | (0.980) | (0.000) | (0.019)
$\pi^{0}$-like |  |  |  |  |
Two-ring | 0.051 | 0.016 | 0.011 | 0.000 | 0.922
| (0.109) | (0.036) | (0.018) | (0.000) | (0.837)
Fully Contained (FC) Multi-GeV
Single-ring |  |  |  |  |
$\nu_{e}$-like | 0.726 | 0.077 | 0.058 | 0.027 | 0.113
| (0.748) | (0.066) | (0.064) | (0.016) | (0.105)
$\bar{\nu}_{e}$-like | 0.553 | 0.379 | 0.003 | 0.008 | 0.056
| (0.566) | (0.371) | (0.003) | (0.007) | (0.053)
$\mu$-like | 0.000 | 0.000 | 0.996 | 0.003 | 0.001
| (0.000) | (0.000) | (0.995) | (0.004) | (0.001)
Multi-ring |  |  |  |  |
$\nu_{e}$-like | 0.588 | 0.117 | 0.054 | 0.036 | 0.204
| (0.609) | (0.112) | (0.059) | (0.032) | (0.188)
$\bar{\nu}_{e}$-like | 0.526 | 0.300 | 0.021 | 0.020 | 0.134
| (0.541) | (0.301) | (0.023) | (0.016) | (0.118)
$\mu$-like | 0.010 | 0.001 | 0.959 | 0.004 | 0.026
| (0.016) | (0.002) | (0.946) | (0.005) | (0.031)
Other | 0.283 | 0.026 | 0.342 | 0.053 | 0.295
| (0.302) | (0.032) | (0.342) | (0.051) | (0.274)

An MC study was performed to investigate the effect of decreases in purity
against improved statistical power of an expanded FV. Figure 12 shows the
sensitivity to the normal mass hierarchy assuming different FV cut values.
Here systematic errors are assumed to be the same as those from the standard
SK analysis, which uses a cut at 200 cm. The rise in sensitivity as the FV is
expanded indicates that a FV cut at 50 cm can be expected to provide improved
hierarchy sensitivity as long as systematic errors in the detector region
outside of the conventional cut boundary are stable. A more precise
sensitivity study with updated systematic errors for this region is presented
below.

Figure 12:  Sensitivity to reject the wrong mass hierarchy as a function of
the FV cut for a 3118.5 day exposure. The normal hierarchy is assumed to be
the true hierarchy with oscillation parameters taken to be $\Delta
m_{23}^{2}=2.4\times 10^{-3}\mbox{eV}^{2}$, $\mbox{sin}^{2}\theta_{23}=0.5$,
$\mbox{sin}^{2}\theta_{13}=0.0210$, and $\delta_{CP}=0$. Systematic error are
assumed to be the same as those for the analysis with a FV cut at 200 cm.

###  6.2 Background estimation

In addition to entering backgrounds and events with incorrectly assigned PID,
non-neutrino (non-$\nu$) backgrounds from cosmic ray muons and flashing PMTs
are potential background sources when the FV is expanded.

Cosmic ray muons must traverse the OD before reaching the ID and are therefore
mostly rejected by cuts on activity in the OD. Any remaining events are
rejected by the FV cut, since their true vertex position, that is the entrance
point of the muon, is expected to be on the wall of tank. However, when the FV
cut is near the wall cosmic ray background events might be introduced into the
FC sample due to the vertex resolution of the reconstruction algorithm.

Flasher events are caused by the internal electrical discharge of a PMT and
represent another kind of non-neutrino background. The reconstructed vertex
for such events is also expected to be on the wall of the detector, so as was
the case for cosmic ray muons, more flasher backgrounds might exist in the FC
sample when the FV is enlarged.

In order to observe possible excesses in the data due to such backgrounds,
which are not modeled in the atmospheric neutrino MC, the distribution of the
distance from the reconstructed vertex to the nearest wall of FC events is
compared to MC in Figure 13. Though not shown here, the distributions of the
events’ momentum, direction, and particle ID have been checked for both the
conventional FV region (dwall $>$ 200 cm) and new region (50 cm $<$ dwall $<$
200 cm), respectively. No evidence for the presence of a large unmodeled
component to the data is observed.

Figure 13:  The distribution of the distance from the reconstructed event
vertex to the nearest ID wall of FC data (points) and MC (solid line). The
lower figure shows the ratio of the MC to the data.

To evaluate the background rate more precisely all FC events within the
fiducial volume were scanned by eye using a graphical event display program.
In total, one flasher event and 24 cosmic ray muon events were found in the
observed in the FC sample with dwall $>$ 50 cm. Most of the cosmic ray muon
backgrounds were reconstructed as multi-ring downward-going $\mu$-like events
with more than 5 GeV/c of momentum. Eye scans of these events indicated that
their vertices were incorrectly reconstructed within the new fiducial volume
due to the presence of an erroneous second ring found by the fitter. Both the
rate and type of such mis-reconstructions is consistent with that found in
independent cosmic ray muon samples. Using the dwall distribution of the
cosmic ray muon events, the expected cosmic ray background in the final
analysis sample is statistically removed. This is done under the assumption
that the FC event reduction process does not depend strongly on the muon
kinematics. By comparison of the event classification of this control sample
with that observed in the eye-scanned background, this assumption has been
validated for all particle directions and momenta. As an example, Figure 14
shows the dwall distribution for multi-ring $\mu$-like events. More than 20,
000 cosmic ray muon events (solid line) are shown here and their dwall
distribution has been normalized to the eye-scanned background events (points)
in the FC sample. The shape of these two samples are in good agreement. Due to
their limited contamination in the final FC sample, the effect of flasher PMTs
is treated with a systematic error on their rate in the oscillation analysis
presented below.

Figure 14:  Distribution of the distance to the nearest ID wall from the
reconstructed vertex of events passing the multi-ring $\mu$-like sample
selection. The solid line is from a cosmic ray muon control sample and the
points show the distribution of events identified as backgrounds in the FC
data. The histograms have been normalized by the total number of events seen
in the FC sample.

Expanding the FV opens the possibility that events interacting within the ID
but close to its wall may produce particles which escape into the OD. Poor
modeling of the response of the OD can thereby potentially introduce biases
and relative inefficiencies in the FC sample. Figure 15 shows the distribution
of OD hits used to define the FC sample after the reduction processes for data
with dwall $>50$ cm. The detection inefficiency due to the cut on OD hits for
both data and MC are confirmed to be consistent, around 0.2%. Based on this
result in conjunction with the stability of the reconstruction algorithm, its
systematic errors, and sample purities across the detector, the fiducial
volume definition is expanded from its conventional value to dwall $>50$ cm
(expanded FV) in the analysis presented below.

Figure 15:  Distribution number of OD hits in cluster for events within new
region (50 cm $<$ dwall $<$ 200 cm) of data (points) and atmospheric neutrino
MC (solid line). The red dashed line shows the threshold to select FC events.

Zenith angle distributions for each analysis sample using the expanded FV are
shown in Figure 16. Their event rates since the start of SK-IV of the
experiment have been stable at 8.3 (2.2) FC events per day in the conventional
FV (new region), 0.73 PC events per day, and 1.49 Up-$\mu$ events per day, as
shown in Figure 17.

Figure 16:  Data and MC comparisons for the SK-IV data divided into 18
analysis samples. The expanded FV, where dwall $>$ 50 cm, is shown here.
Samples with more than one zenith angle bin are shown as zenith angle
distributions (second through fifth column) and other samples are shown as
reconstructed momentum distributions (first column). Cyan lines denote the
best-fit MC assuming the normal hierarchy. Narrow panels below each
distribution show the ratio relative to the normal hierarchy MC. In all panels
the error bars represent the statistical uncertainty.  Figure 17:  Final event
rates as a function of time since the start of SK-IV operations. The error
bars are statistical. Circles denote the fully contained event rate and
upward-facing (downward-facing) triangles show the partially contained
(upward-going muon) event rates. The number in the parentheses for fully
contained event shows the fiducial volume cut value.

##  7 Systematic error

In order to use the fiTQun reconstruction algorithm for physics analysis the
corresponding systematic errors must be evaluated. While systematic errors
related to the atmospheric neutrino flux and cross section model are the same
as those used in Abe:2017osc , systematic errors on the event selection have
been updated for FC events using the expanded fiducial volume outlined above.
Since fiTQun is not used to reconstruct PC or Up-$\mu$ events in the present
work, the errors for those samples are the same as in the previous
publications.

The uncertainty on the absolute energy scale is estimated by comparison of
data and MC across three control samples spanning momenta up to 10 GeV/c:
electrons from cosmic ray muon decay, atmospheric neutrino events producing
single $\pi^{0}$s from neutral current interactions, and stopping cosmic ray
muons. The difference between data and MC for the calibration samples
reconstructed within the expanded FV is shown in Figure 18. The total
systematic error follows that of Abe:2017osc  and is the sum in quadrature of
the absolute energy scale, which is defined as the largest data-MC difference
across all samples, and momenta with the average time variation of these
samples throughout SK-IV. The uncertainty from the up/down asymmetry of the
detector is measured with Michel electrons from cosmic ray muons by comparison
of their momenta for data and MC as a function of zenith angle. The difference
between data and MC at the most deviated direction is taken as the uncertainty
from the detector asymmetry. As summarized in Table 6, the uncertainty of
energy scale on the conventional (expanded) FV is $2.17\%$ (2.02%) using
fiTQun reconstruction, which is similar to the value with APFit: 2.1%
Abe:2017osc .

Table 6: Summary of uncertainty from the energy scale. The energy scale error is the quadratic-sum of absolute energy scale and the time variation. The absolute energy scale value is obtained from the most discrepant sample, i.e. Michel electrons from cosmic ray muons in this study.  Reconstruction | fiTQun | fiTQun | APFitAbe:2017osc
---|---|---|---
Dwall range |  $>$50 cm |  $>$ 200 cm |  $>$ 200 cm
Energy calibration error | 2.02% | 2.17% | 2.1%
Absolute energy scale | 1.92% | 2.09% |
Time variation | 0.62% | 0.59% |
Up/down asymmetry energy calibration error | 0.67% | 0.58% | 0.4%
Figure 18:  Absolute energy scale measurement for SK-IV with fiTQun and the
expanded FV (dwall $>$ 50 cm). Vertical error bars denote the statistical
uncertainty and horizontal error bars shows the momentum range spanned by each
control sample.

Lacking other control samples that span the same energies and event topologies
as the atmospheric neutrino data, systematic errors on the event selection
including the particle identification (PID), estimated number of rings (ring
counting), and the two-stage separation of the multi-GeV multi-ring e-like
event samples described above are evaluated using the atmospheric neutrino
data.

Systematic errors are estimated for each of the FC analysis samples for both
the conventional fiducial volume and the additional region formed between the
conventional and expanded FV boundaries. The results are summarized in Table
7. Errors for the individual samples subject to a particular uncertainty, such
as the particle ID for single ring-events, are assumed to be fully correlated
(or fully anti-correlated). Anti-correlations are indicated by negative
numbers in the table. No correlation is assumed between the systematic error
categories. The sizes of the systematic errors evaluated using fiTQun and
APFit are consistent within the conventional fiducial volume, with
improvements seen in fiTQun for many samples. For many error categories the
fiTQun reconstruction shows consistent errors between the conventional
fiducial volume and the new region (50 cm $<$ dwall $<$ 200 cm). In some cases
however, such as the ring counting uncertainty for multi-GeV single-ring
events, larger systematic errors are observed. At present these are not large
enough to offset the benefit to the oscillation sensitivity from including
events in the new region in the analyses. A quantitative estimation of the
sensitivity is presented in next the section. The numbers for the different FV
regions are merged to form the final response function used in the analysis,
which will be described in next section. Systematic errors and their sizes at
the best fit point of the analysis are presented in Tables LABEL:tab:sysa.

Table 7: Summary of systematic errors related to the ring counting, particle identification and multi-GeV multi-ring $e$-like separation. The sign of the number represents the correlation between these systematic errors, with negative numbers denoting fully anti-correlated samples.  Reconstruction | fiTQun | fiTQun | APFit
---|---|---|---
Dwall range | 50 cm$\sim$200 cm |  $>$ 200 cm |  $>$ 200 cm
Ring counting
Sub-GeV , Single-ring |  |  |
$e$-like, p $<$ 400 MeV | 1.94% | 1.20% | 1.6%
$e$-like, p $>$ 400 MeV | 0.59% | 0.48% | 1.0%
$\mu$-like, p $<$ 400 MeV | 1.08% | 0.42% | -3.0%
$\mu$-like, p $>$ 400 MeV | 1.25% | 1.21% | 0.6%
Sub-GeV, Multi-ring |  |  |
$e$-like | -3.58% | -2.39% | -1.9%
$\mu$-like, $E_{\text{vis}}>$ 600MeV | -2.32% | -1.66% | 2.3%
Multi-GeV Single-ring |  |  |
$e$-like | 8.61% | 1.21% | 1.0%
$\mu$-like | 2.65% | -2.33% | -1.2%
Multi-GeV Multi-ring |  |  |
$e$-like | -4.21% | -0.62% | -0.9%
$\mu$-like | -3.07% | 0.72% | 2.4%
Particle ID of single-ring events
Sub-GeV |  |  |
$e$-like | 0.99% | 0.36% | -0.28%
$\mu$-like | -0.89% | -0.37% | 0.22%
Multi-GeV |  |  |
$e$-like | 0.23% | 0.06% | -0.35%
$\mu$-like | -0.21% | -0.06% | 0.35%
Particle ID of the brightest ring in multi-ring events
Sub-GeV |  |  |
$e$-like | -3.19% | -0.72% | 4.19%
$\mu$-like | 1.31% | 0.31% | -1.56%
Multi-GeV |  |  |
$e$-like | 1.94% | 1.10% | 3.33%
$\mu$-like | -1.06% | -0.66% | -1.56%
Multi-GeV Multi-ring $e$-like - other separation
$e$-like | -0.88% | -0.67% | 3.0%
other | 0.50% | 0.53% | -3.4%
Multi-GeV Multi-ring $\nu_{e}$-like - $\bar{\nu}_{e}$-like separation
$\nu_{e}$-like | -3.64% | -2.33% | 6.82%
$\bar{\nu_{e}}$-like | 4.51% | 2.10% | -6.04%

##  8 Atmospheric Neutrino Oscillation Analysis

The Super-K atmospheric neutrino data are separately fit against both the
normal and inverted hierarchy hypotheses. Since $\theta_{13}$ has been
measured and is well constrained by reactor experiments, two types of analyses
are done: one with $\theta_{13}$ as a free fitting parameter and the other
with it constrained to $0.0210\pm 0.0011$ Olive:2016xmw . Though the
constraint on this parameter is much stronger than can be expected from
Super-K alone, the appearance of upward-going electron events characteristic
of the mass hierarchy signal is driven by $\theta_{13}$ and therefore
measurements with atmospheric neutrinos represent an important test of the
analysis’ hierarchy preference. The solar mixing parameters, $\Delta
m^{2}_{12}$ and $\sin^{2}\theta_{12}$, on the other hand, have little impact
on the oscillations of neutrinos in the energy range considered here and will
be fixed in the analysis as described below.

Data are fit to the MC expectation with a binned $\chi^{2}$ method assuming
Poisson statistics. The effect of systematic errors are regarded as scaling
factors on the MC in each bin Fogli:2002pt , with different error sources
assumed to be independent. The definition of $\chi^{2}$ used in the fits is

$\displaystyle\chi^{2}=2\displaystyle\sum_{n}\left(E_{n}-\mathcal{O}_{n}+\mathcal{O}_{n}\ln\frac{\mathcal{O}_{n}}{E_{n}}\right)+\displaystyle\sum_{i}\left(\frac{\epsilon_{i}}{\sigma_{i}}\right)^{2},$
(10)

where,

$\displaystyle E_{n}$ $\displaystyle=$ $\displaystyle\displaystyle
E_{n,0}(1+\displaystyle\sum_{i}f^{i}_{n}\epsilon_{i})$ (11)

In this equation $\mathcal{O}_{n}$ is the observed number of events in the
$n^{th}$ analysis bin. Similarly, $E_{n,0}$ represents the nominal MC
expectation in that bin and $E_{n}$ is the expectation after incorporating the
effect of the systematic errors. The coefficient $f^{i}_{n}$ describes the
fractional change in the bin’s MC under a $1{\sigma_{i}}$ variation of the
$i^{th}$ systematic error source. The sum over $(\epsilon_{i}/\sigma_{i})^{2}$
penalizes the $\chi^{2}$ for adjusting the systematic errors when bringing the
MC into agreement with the data. For each set of oscillation parameters
tested, the MC fit to the data by minimizing $\chi^{2}$ over the error
parameters $\epsilon_{i}$. All fits are performed over 515 analysis bins in
SK-IV.

###  8.1 Analysis with unconstrained $\theta_{13}$

Figure 19:  Expected sensitivity to reject $\sin^{2}\theta_{13}$ = 0 assuming
its true value is $0.0210$ as a function of the true value of
$\sin^{2}\theta_{23}$ for a livetime of 3118.5 days. Grey and blue lines show
the sensitivity for samples reconstructed with fiTQun and using the
conventional FV and the expanded FV, respectively. The orange line denotes the
sensitivity when events are reconstructed with APFit in the conventional FV.

Figure 20:  Constraints on neutrino oscillation parameters from the Super-K atmospheric neutrino data with no assumed constraint on $\theta_{13}$. The solid blue and dashed orange lines denote the normal and inverted hierarchy fit results, respectively. The latter has been offset from the former by the difference in their minimum $\chi^{2}$ values. Table 8: Parameter grid used in the fit. Parameter | Range | Scan points
---|---|---
$|\Delta m^{2}_{32,31}|$ | $1.5\times 10^{-3}\sim 3.5\times 10^{-3}$ | 40 points
$\sin^{2}\theta_{13}$ | $0.0\sim 0.06$ | 25 points
$\sin^{2}\theta_{23}$ | $0.3\sim 0.725$ | 35 points
$\delta_{CP}$ | $0.0\sim 2\pi$ | 19 points
Table 9: Values of oscillation parameters fixed in the analysis and their systematic errors. Note that $\sin^{2}\theta_{13}$ is only fixed in the “$\theta_{13}$ constrained”analyses described in Section 8.2. Parameter | Value
---|---
$\Delta m^{2}_{21}$ | $(7.53\pm 0.18)\times 10^{-5}\mbox{eV}^{2}$
$\sin^{2}\theta_{12}$ | $0.307\pm 0.013$
$\sin^{2}\theta_{13}$ | $0.0210\pm 0.0011$

The fit with $\theta_{13}$ unconstrained by reactor measurements is performed
over four parameters, $|\Delta m^{2}_{32,31}|$,
$\sin^{2}\theta_{13}$,$\sin^{2}\theta_{23}$ and $\delta_{CP}$ with a total of
71 systematic errors. Fits to the normal hierarchy use $\Delta m^{2}_{32}$
while those for the inverted hierarchy use $\Delta m^{2}_{31}$. The agreement
between the data and MC is evaluated using Equation 10 at each point in the
grid of points shown in Table 8, while keeping the solar mixing parameters,
$\Delta m^{2}_{12}$ and $\sin^{2}\theta_{12}$, fixed to the values in Table 9.
Systematic errors representing the uncertainty in these parameters are
included in the analysis. The best-fit parameters for each hierarchy are
defined as the parameter set returning the smallest $\chi^{2}$ value. Between
the two hierarchy fits the one with smallest minimum $\chi^{2}$ value is taken
as the global best-fit and is taken to indicate the mass hierarchy preference.

### Results and Discussion

Figure 19 shows the expected sensitivity to reject a zero value of
$\sin^{2}\theta_{13}$ as a function of the true value of
sin${}^{2}\theta_{23}$ for the fiTQun-based analysis samples with both the
conventional and expanded FV as well as an analysis with APFit-reconstructed
samples with the conventional FV. The sensitivities for the two algorithms
using the same conventional FV are similar since the fraction of events common
between them is larger than 97%. The expected number of atmospheric neutrino
events selected by APFit and fiTQun are 24188.7 and 24184.3, respectively,
while the number of the events passed both the cut of the two algorithms is
23520.2. The 3% difference is mainly due to the vertex resolution near FV
boundary.

However, since the contamination of $\nu_{\mu}$ events in the e-like samples
most sensitive to $\theta_{13}$ is lower for the fiTQun-selected sample as
shown in Table 10, there is a slight boost in its sensitivity. Expanding the
FV naturally incorporates more signal events in the analysis and the
sensitivity improvement is clear.

Table 10:  $\nu_{\mu}+\bar{\nu}_{\mu}$ CC backgrounds in multi-GeV e-like sample.  Reconstruction | fiTQun | fiTQun | APFit
---|---|---|---
Dwall range | 50 $\sim$ 200 cm |  $>200$ cm |  $>$ 200 cm
Single-ring |  |  |
$\nu_{e}$-like | 0.064 | 0.058 | 0.109
$\bar{\nu}_{e}$-like | 0.003 | 0.003 | 0.009
Multi-ring |  |  |
$\nu_{e}$-like | 0.059 | 0.054 | 0.121
$\bar{\nu}_{e}$-like | 0.023 | 0.021 | 0.040

One-dimensional allowed regions for $\theta_{13}$, $|\Delta m^{2}_{32,31}|$,
$\sin^{2}\theta_{23}$ and $\delta_{CP}$ from the fit to the data are shown in
Figure 20. The normal hierarchy hypothesis yields better data-MC agreement
than the inverted hierarchy hypothesis with
$\chi^{2}_{\text{NH,min}}-\chi^{2}_{\text{IH,min}}=-1.81$. The 1$\sigma$
allowed region for $\sin^{2}\theta_{13}$ is from 0.003 to 0.033 (from 0.001 to
0.023) for the normal (inverted) hierarchy fit, which is consistent with the
globally preferred value. The point at $\sin^{2}\theta_{13}=0.0$ is disfavored
at approximately $1.8{\sigma}$ (1.2 $\sigma$) for normal (inverted) fit. A
summary of the best-fit information and parameter constraints is presented in
Table 11.

As discussed in Section 2, the determination of $\sin^{2}\theta_{13}$ and the
mass hierarchy using atmospheric neutrinos is achieved using upward-going
electron neutrino appearance at multi-GeV energies. Figure 21 shows the up-
down asymmetry of the multi-GeV single- and multi-ring electron-like samples,
where the asymmetry is defined as the ratio of the difference of the number of
downward-going and downward-going events relative to their sum. Here downward-
going (upward-going) events are events whose zenith angle satisfy
$\mbox{cos}\theta_{z}<-0.4$ ($\mbox{cos}\theta_{z}>0.4$). Excesses between a
few and ten GeV in the multi-GeV e-like $\nu_{e}$ and the Multi-Ring e-like
$\nu_{e}$ samples, drive the normal hierarchy preference.

The best-fit atmospheric mixing parameters from the normal hierarchy fit are
$\Delta m^{2}_{32}=2.63^{+0.10}_{-0.21}\times 10^{-3}\mbox{eV}^{2}$ and
$\sin^{2}\theta_{23}=0.588^{+0.030}_{-0.062}$ for second octant (best-fit) and
$0.425^{+0.051}_{-0.034}$ for first octant. Maximal mixing
($\sin^{2}\theta_{23}=0.5$) is weakly disfavored at around $1{\sigma}$
significance. Besides the data excesses in the upward-going regions of the
single-ring e-like $\nu_{e}$ sample and multi-ring other sample, the data
deficits in the multi-GeV $\mu$-like sample also contribute to this
preference. It should be noted that the preference for $\sin^{2}\theta_{23}$
is coupled to the $\theta_{13}$ measurement, since both parameters feature in
the $\nu_{\mu}\rightarrow\nu_{e}$ oscillation probability.

The best-fit value of $\delta_{CP}$ is 3.84 radians for both the normal and
inverted hierarchy fits. Comparing with the least preferred parameter value of
0.8 radians, more electron neutrino appearance in the sub-GeV e-like samples
are expected to be observed due to $\nu_{\mu}\rightarrow\nu_{e}$ oscillations.
In the multi-GeV region, $\delta_{CP}$ similarly modulates the amount of
electron neutrino appearance but it is subdominant to the effects of
$\theta_{13}$. In the inverted hierarchy fit, the value of $\theta_{13}$ can
be more easily adjusted to bring the MC prediction into better agreement with
data since only antineutrino events experience matter effects in the earth.
Since there are more neutrino events than antineutrino events in the
atmospheric sample, the freedom to move $\theta_{13}$ ultimately leads to a
weaker constraint on $\delta_{CP}$.

Figure 21:  Ratio of upward- ($\mbox{cos}\theta<-0.4$) to downward-going
($\mbox{cos}\theta>0.4$) events as a function of energy for the mass
hierarchy-sensitive analysis samples. energy. The error bars are statistical.
For the single-ring samples the energy is taken to be the visible energy
assuming the observed event was an electron. The energy estimator for multi-
ring samples is the total observed energy summed over each reconstructed ring
after adjusting for each ring’s PID. The orange line shows the best-fit result
assuming the inverted hierarchy hypothesis and the cyan that from the normal
hierarchy hypothesis.

###  8.2 Analysis with constrained $\theta_{13}$

Since atmospheric neutrinos come from all directions, carry a wide range of
energies, and often produce particles that are invisible to Super-K it is not
possible in general to fully reconstruct the neutrino kinematics. As a result
atmospheric neutrinos themselves do not typically have oscillation parameter
sensitivity at the same level as that in long-baseline or reactor experiments,
such that the introduction of constraints from those measurements can improve
sensitivity to the mass hierarchy and and $\delta_{CP}$. Indeed, Equation 2
shows that the size of MSW resonance relies on the value of
$\mbox{sin}^{2}\theta_{13}$ directly. Since reactor neutrino experiments
constrain this parameter more precisely than the analysis in the previous
section, the analysis presented here constrains the value of
$\sin^{2}\theta_{13}$ to $0.0210\pm 0.0011$ Olive:2016xmw . During the fit the
parameter is held at its central value and a systematic error that imparts the
impact of its uncertainty is added to the analysis. The fit uses the same data
samples, binning, and parameter grid for all other oscillation parameters as
the analysis presented above.

### Results and Discussion

Figure 22:  Expected sensitivity to the normal mass hierarchy (left) and
inverted hierarchy (right) as a function of the true value of
$\sin^{2}\theta_{23}$. Here $\sin^{2}\theta_{13}=0.0210\pm 0.0011$ and the
assumed livetime is 3118.5 days. Grey and blue bands show the sensitivity of
the analysis with event samples reconstructed with fiTQun in the conventional
FV and expanded FV, respectively. Orange lines denote the sensitivity when
events are reconstructed using the APFit algorithm with the conventional FV.
The width of the bands corresponds to the uncertainty from $\delta_{CP}$.
Figure 23:  Expected sensitivity to reject wrong octant as a function of the
true value of $\sin^{2}\theta_{23}$ for true normal mass hierarchy. The other
conditions are the same as Figure 22.

The expected sensitivity to the mass hierarchy and $\theta_{23}$ octant for
true normal hierarchy are illustrated in Figure 22 and Figure 23 for the
fiTQun-reconstructed analysis samples using both the conventional and expanded
FV and for the APFit-based analysis with the conventional FV. Improved
sensitivity is found with the fiTQun-based analyses regardless of the fiducial
volume used for the same reasons as produced the improved sensitivity to
$\theta_{13}$ in the previous section.

Figure 24 shows the $\chi^{2}$ value as a function of the atmospheric neutrino
mixing parameters and $\delta_{CP}$ in the $\theta_{13}$-constrained fit. As
in the unconstrained fit the data prefer the normal mass hierarchy with
$\Delta\chi^{2}=-2.45$. The 1$\sigma$ allowed region for $|\Delta m^{2}_{32}|$
is from 2.41 to 2.75 $\times 10^{-3}\mbox{eV}^{2}$ (2.36 to 2.67 $\times
10^{-3}\mbox{eV}^{2}$) for the normal (inverted) hierarchy fit, which is
consistent with the result of unconstrained fit. The contour for the allowed
region with 90% CL is shown in Figure 25. The preference for the second octant
of $\theta_{23}$ in the unconstrained fit is changed. Adding the constraint on
$\sin^{2}\theta_{13}$ the preferred value of $\theta_{23}$ shifts from the
second to first octant, though both are allowed at 1$\sigma$. Indeed the
difference in $\chi^{2}$ between the minimum value of each octant is
$\Delta\chi^{2}_{1^{st}-2^{nd}}$ = -0.73 (0.13) for the constrained
(unconstrained) fit. This result can be explained by degeneracies in the
appearance probability that arise for certain combinations of
$\sin^{2}2\theta_{13}$ and $\sin^{2}\theta_{23}$ as shown in Equation 1; an
increase in one can be compensated by a decrease in the other to describe the
data. The same upward-going event excesses (deficits) in the $e$-like
($\mu$-like) samples that drive the hierarchy, $\theta_{13}$, and octant
preferences in the unconstrained analysis are not strong enough to support the
increased excess (deficits) expected for higher values of $\theta_{13}$ and
the second octant. Since the $\theta_{13}$ is now fixed, the data can be
accommodated by reducing the expected number of $e$-like events and increasing
the number of $\mu$-like events by moving $\theta_{23}$ to the first octant.

The best-fit value of $\delta_{CP}$ is 3.14 (4.89) for the normal (inverted)
hierarchy hypothesis, with a tighter constraint compared with the
unconstrained fit. Parameter values and their $1{\sigma}$ errors are
summarized in Table 11.

Figure 24:  Constraints on neutrino oscillation parameters from SK-IV atmospheric neutrino data using the expanded FV and assuming $\mbox{sin}^{2}\theta_{13}=0.0210\pm 0.0011$. The solid blue and dashed orange lines denote the normal and inverted hierarchy fit results, respectively. The latter has been offset from the former by the difference in their minimum $\chi^{2}$ values. Figure 25:  Constraints on the atmospheric mixing parameters using SK-IV atmospheric neutrino and the expanded FV. The solid blue (dashed orange) line shows 90% C.L. for the normal (inverted) hierarchy. The star denotes the best-fit value, which is at the same point for normal and inverted hierarchy, as shown in Table 11. In each contour $\mbox{sin}^{2}\theta_{13}$ is constrained to be $0.0210\pm 0.0011$. The contours have both been drawn relative to the global best-fit. Table 11: Parameter estimates for each analysis and hierarchy hypothesis considered. Here NH (IH) refers to the normal (inverted) hierarchy fit. The terms “Unconstrained” and “Constrained” refer to fits without and with a constraint on $\mbox{sin}^{2}\theta_{13}$, respectively. For $\sin^{2}\theta_{23}$ parameter ranges are shown for both octants, with the best-fit octant enclosed in a box. The expected absolute $\chi^{2}$ value for the constrained fit is 523.6 and the probability for obtaining larger value is 0.244 for the NH.  |  $\theta_{13}$ Free |  $\theta_{13}$ Constrained
---|---|---
Hierarchy | NH | IH | NH | IH
$\chi^{2}$ | 576.3 | 578.1 | 576.5 | 579.0
$\mbox{sin}^{2}\theta_{13}$ | $0.008^{+0.025}_{-0.005}$ | $0.008^{+0.015}_{-0.007}$ | – | –
$\mbox{sin}^{2}\theta_{23}$ ($1^{st}$ oct.) | $0.425^{+0.051}_{-0.034}$ | $0.425^{+0.075}_{-0.027}$ |  $0.425^{+0.046}_{-0.037}$ |  $0.425^{+0.055}_{-0.036}$
$\mbox{sin}^{2}\theta_{23}$ ($2^{nd}$ oct.) |  $0.588^{+0.030}_{-0.062}$ |  $0.575^{+0.034}_{-0.075}$ | $0.600^{+0.013}_{-0.030}$ | $0.588^{+0.022}_{-0.037}$
$|\Delta m^{2}_{32,31}|$ [$\times 10^{-3}$ eV2] | $2.63^{+0.10}_{-0.21}$ | $2.53^{+0.14}_{-0.08}$ | $2.53^{+0.22}_{-0.12}$ | $2.53^{+0.14}_{-0.31}$
$\delta_{CP}$ | $3.84^{+2.00}_{-2.14}$ | $4.19^{+2.09}_{-4.19}$ | $3.14^{+2.67}_{-1.35}$ | $4.89^{+1.51}_{-3.46}$

The $\mbox{CL}_{s}$ method Read:2002hq  is used to address the significance of
the octant and mass hierarchy preferences observed in the data. As in the
previous Super-K study Abe:2017osc  MC ensembles assuming different parameter
combinations of the octant and hierarchy were generated with statistical and
systematic error variations applied to each pseudo data set. Statistical
variations have been applied assuming the current detector exposure and
systematic errors have been varied as Gaussian parameters with widths
specified by their uncertainties. The $\mbox{CL}_{s}$ parameter for the octant
study, $\mbox{CL}^{\mbox{O}}_{s}$, is defined as:

$\mbox{CL}^{\mbox{O}}_{s}=\frac{p_{0}(O_{2})}{1-p_{0}(O_{1})},$ (12)

where $p_{0}(O_{i})$ represents the $p$-value with respect to obtaining a
difference in the best-fit $\chi^{2}$ values between both octant hypotheses
smaller (larger) than the value of data assuming the true value of
$\mbox{sin}^{2}\theta_{23}$ is $O_{2}=0.6$ ($O_{1}=0.4$). These pseudo MC
distributions of these hypotheses are shown as the orange and cyan shaded
histograms in Figure 26, respectively. The $\mbox{CL}^{\mbox{O}}_{s}$ value is
found to be just 0.098 and accordingly, the data show no strong preference for
the octant.

Figure 26:  Distributions of the difference in best-fit $\chi^{2}$ values
between first octant and second octant fits to pseudo data sets used in the
generation of the $\mbox{CL}^{\mbox{O}}_{s}$ value for the SK $\theta_{13}$
constrained analysis. In the cyan (orange) histogram the pseudo data have been
generated assuming $\sin^{2}\theta_{23}=0.4$ ($\sin^{2}\theta_{23}=0.6$).
Shaded portions of the histograms denote the fraction of pseudo data sets with
more extreme values than that observed in the data,
$\Delta\chi^{2}_{data}=-0.73$.

Similarly the parameter for the mass hierarchy is defined as

$\mbox{CL}^{\mbox{H}}_{s}=\frac{p_{0}(\mbox{IH})}{1-p_{0}(\mbox{NH})},$ (13)

where $p_{0}(\mbox{IH})$ and $p_{0}(\mbox{NH})$ are p-values for obtaining a
difference in the $\chi^{2}$ of the best-fit mass hierarchies more extreme
than that of the data assuming a true IH and NH, respectively.

Figure 27:  Distributions of the difference in best-fit $\chi^{2}$ values
between normal- and inverted-hierarchy fits to pseudo data sets used in the
generation of the $\mbox{CL}^{\mbox{H}}_{s}$ value for the SK $\theta_{13}$
constrained analysis. In the cyan (orange) histogram the pseudo data have been
generated assuming the normal (inverted) hierarchy at the analysis best-fit
shown in Table 11. Shaded portions of the histograms denote the fraction of
pseudo data sets with more extreme values than that observed in the data,
$\Delta\chi^{2}_{data}=-2.45$.

Figure 27 shows the distribution for the mass hierarchy determination. Due to
the large uncertainty on $\theta_{23}$, MC ensembles have been generated with
different assumed values $\theta_{23}$ and with the other oscillation
parameters fixed to their best-fit values. The $p$-values and
$\mbox{CL}^{\mbox{H}}_{s}$ values for hierarchy test are summarized in Table
12. As sin${}^{2}\theta_{23}$ ranges from 0.4 to 0.6 the observed
$\mbox{CL}^{\mbox{H}}_{s}$ values decrease from a 30.8% C.L. preference for
the inverted mass hierarchy to 14.3%. At the analysis best-fit value the
preference is 26.0%. A smaller value of $\mbox{sin}^{2}\theta_{23}$ predicts
less electron neutrino appearance, which results in a smaller
$p_{0}(\mbox{IH})$ values and larger $\mbox{CL}^{\mbox{H}}_{s}$ due to the
multi-GeV $e$-like samples’ event excesses seen in data.

Table 12: Normal hierarchy significance summarized in terms of the probability of observing a $\chi^{2}$ preference for the NH more extreme than that observed in the data assuming an IH, $p_{0}(\mbox{IH})$ , and $\mbox{CL}^{\mbox{H}}_{s}$ values for a range of assumed parameters. The first row shows the true $\theta_{23}$ used to generate MC ensembles used in the calculations. Other oscillation parameters are taken from the analysis’ best-fit. The value of $\theta_{13}$ is constrained to 0.0210 $\pm$ 0.0011. True sin${}^{2}\theta_{23}$ | 0.4 | 0.425 | 0.5 | 0.6
---|---|---|---|---
$p_{0}(\mbox{IH})$ | 0.025 | 0.033 | 0.065 | 0.072
$\mbox{CL}^{\mbox{H}}_{s}$ | 0.308 | 0.260 | 0.229 | 0.143

##  9 Conclusion

Figure 28:  Expected sensitivity to the normal mass hierarchy as a function of
livetime assuming $\sin^{2}\theta_{23}$=0.5. Grey and blue bands show the
sensitivity of the atmospheric neutrino sample reconstructed with fiTQun in
the conventional FV and expanded FV, respectively. Orange lines denote the
sensitivity when events are reconstructed with APFit in the conventional FV.

A new reconstruction algorithm based on a maximum likelihood method has been
developed for Super-K. Compared to the conventional reconstruction algorithm,
the new algorithm shows improved performance in a variety of metrics including
event vertex resolution, particle momentum resolution, and particle
identification. The new algorithm has further demonstrated reliable
performance over a larger volume of the detector than the previous volume.
Accordingly, in the present analysis the fiducial volume definition has been
expanded to include all events reconstructed more than 50 cm from any ID wall.
This represents a 32% increase in the number of events available for analysis
relative to the 200 cm threshold used in previous SK analyses.

Using the new algorithm with its expanded fiducial volume definition an
analysis of a 253.9 kton-year exposure of the SK-IV atmospheric data has
yielded oscillation parameter estimates consistent with both previous Super-K
measurements and results from other experiments. Assuming the normal mass
hierarchy constraints on the atmospheric mixing parameters are
$\mbox{sin}^{2}\theta_{23}=0.425^{+0.046}_{-0.037}$
($0.600^{+0.013}_{-0.030}$) for first (second) octant and $\Delta
m^{2}_{32}=2.53^{+0.22}_{-0.12}$, with $\delta_{CP}=3.14^{+2.67}_{-1.35}$. The
data show a weak preference for the normal mass hierarchy, disfavoring the
inverted mass hierarchy at 74.0% assuming oscillation parameters at the
analysis best-fit point. No strong preference for the $\theta_{23}$ octant is
observed.

Figure 28 shows the expected sensitivity to mass hierarchy as the function of
livetime for both reconstruction algorithms and their fiducial volumes. The
expected improvement in sensitivity with the new algorithm becomes more
apparent as data is accumulated, even assuming the conventional FV. While the
new reconstruction has only been applied to the 3118.5 day SK-IV fully
contained data set, future efforts will expand this study to include other
Super-K samples and run periods, which constitute an additional 2,800 days of
livetime.

##  10 Acknowledgments

We gratefully acknowledge the cooperation of the Kamioka Mining and Smelting
Company. The Super-Kamiokande experiment has been built and operated from
funding by the Japanese Ministry of Education, Culture, Sports, Science and
Technology, the U.S. Department of Energy, and the U.S. National Science
Foundation. Some of us have been supported by funds from the National Research
Foundation of Korea NRF-2009-0083526 (KNRC) funded by the Ministry of Science,
ICT, and Future Planning, the European Union H2020 RISE-GA641540-SKPLUS, the
Japan Society for the Promotion of Science, the National Natural Science
Foundation of China under Grants No. 11235006, the National Science and
Engineering Research Council (NSERC) of Canada, the Scinet and Westgrid
consortia of Compute Canada, and the National Science Centre of Poland
(2015/17/N/ST2/04064, 2015/18/E/ST2/00758).

## References

  * (1) B. Pontecorvo, Sov. Phys. JETP, 26, 984–988 (1968).
  * (2) Z. Maki, M. Nakagawa, and S. Sakata, Prog. Theor. Phys., 28, 870 (1962).
  * (3) C. Patrignani et al., Chin. Phys., C40(10), 100001 (2016).
  * (4) K. Abe et al., Phys. Rev. Lett., 118(15), 151801 (2017), arXiv:1701.00432.
  * (5) P. Adamson et al., Phys. Rev. Lett., 118(23), 231801 (2017), arXiv:1703.03328.
  * (6) K. Abe et al., Phys.Rev., D97, 072001 (2018), arXiv:1710.09126.
  * (7) R. Acciarri et al. (2016), arXiv:1601.05471.
  * (8) Fengpeng An et al., J. Phys., G43(3), 030401 (2016), arXiv:1507.05613.
  * (9) K. Abe et al. (2018), arXiv:1805.04163.
  * (10) L. Wolfenstein, Phys. Rev., D17, 2369–2374 (1978).
  * (11) S. P. Mikheev and A. Yu. Smirnov, Sov. J. Nucl. Phys., 42, 913–917 (1985).
  * (12) C. Giunti, C.W. Kim, and M. Monteno, Nucl.Phys, B521, 3–36 (1998), arXiv:9709439.
  * (13) H. Nishino, K. Awai, Y. Hayato, S. Nakayama, K. Okumura, M. Shiozawa, A. Takeda, K. Ishikawa, A. Minegishi, and Y. Arai, Nucl. Instrum. Meth., A610, 710–717 (2009), arXiv:0911.0986.
  * (14) Y. Fukuda et al., Nucl. Instrum. Meth., A501, 418–462 (2003).
  * (15) K. Abe, Y. Hayato, T. Iida, K. Iyogi, J. Kameda, et al., Nucl.Instrum.Meth., A737, 253–272 (2014), arXiv:1307.0162.
  * (16) M. Honda, T. Kajita, K. Kasahara, and S. Midorikawa, Phys.Rev., D83, 123001 (2011), arXiv:1102.2688.
  * (17) Y. Hayato, Nucl. Phys. Proc. Suppl., 112, 171–176, [,171(2002)] (2002).
  * (18) Rene Brun, F. Bruyant, Federico Carminati, Simone Giani, M. Maire, A. McPherson, G. Patrick, and L. Urban (1994).
  * (19) R. Gran, J. Nieves, F. Sanchez, and M. J. Vicente Vacas, Phys. Rev., D88, 113007 (2013), arXiv:hep-ph/1307.8105.
  * (20) J. Nieves, I. Ruiz Simo, and M. J. Vicente Vacas, Phys. Rev., C83, 045501 (2011), arXiv:hep-ph/1102.2777.
  * (21) M. Gluck, E. Reya, and A. Vogt, Eur.Phys.J., C5, 461–470 (1998), arXiv:hep-ph/9806404.
  * (22) Arie Bodek and Un-ki Yang, AIP Conf. Proc., 792, 257–260, [,257(2005)] (2005), arXiv:hep-ph/0508007.
  * (23) Y. Ashie. et al., Phys.Rev., D71, 112005 (2005), arXiv:0501064.
  * (24) K. Abe et al. (2013), arXiv:1311.4750.
  * (25) M. Shiozawa et al., Nucl.Instrum.Meth., A433, 240–246 (1999).
  * (26) R. B. Patterson, E. M. Laird, Y. Liu, P. D. Meyers, I. Stancu, and H. A. Tanaka, Nucl.Instrum.Meth, A608, 206–224 (2009), arXiv:0902.2222.
  * (27) K. Abe et al., Phys. Rev. Lett., 121(17), 171802 (2018), arXiv:1807.07891.
  * (28) F. James (1994).
  * (29) R. Wendell et al., Phys.Rev., D81, 092004 (2010), arXiv:1002.3471.
  * (30) G. L. Fogli, E. Lisi, A. Marrone, D. Montanino, and A. Palazzo, Phys.Rev., D66, 053010 (2002), arXiv:hep-ph/0206162.
  * (31) Alexander L. Read, J. Phys., G28, 2693–2704, [11(2002)] (2002).
  * (32) J. Nieves, Jose Enrique Amaro, and M. Valverde, Phys. Rev., C70, 055503, [Erratum: Phys. Rev.C72,019902(2005)] (2004), arXiv:nucl-th/0408005.
  * (33) A. Capella, A. Kaidalov, C. Merino, and J. Tran Thanh Van, Phys. Lett., B337, 358 (1994), hep-ph/9405338.
  * (34) E. Hernandez, J. Nieves, and M. Valverde, Phys. Rev., D76, 033005 (2007), hep-ph/0701149.
  * (35) E. Pinzon,  Measurement of Pion-Carbon Cross Sections at DUET and Measurement of Neutrino Oscillation Parameters at the T2K Experiment,  PhD thesis, York U. (2017).

##  11 Systematic Uncertainties

Table 13: Systematic error used in this analysis. The second column shows the best-fit value of the systematic error parameter in percent and the third column shows the estimated $1{\sigma}$ error size in percent. The fit result within expanded FV and $\theta_{13}$ constrained to $0.0210\pm 0.0011$. The same applies below. The systematic errors about ring counting, particle identification and Multi-ring likelihood selection within conventional fiducial volume (dwall $>$ 200 cm) and within new region (200 cm $>$ dwall $>$ 50 cm) have been merged. Systematic Error | Fit Value (%) |  $\sigma$ (%)
---|---|---
Flux normalization |  $E_{\nu}<1\;\mathrm{GeV}$111Uncertainty decreases linearly with $\log E_{\nu}$ from 25 %(0.1 GeV) to 7 %(1 GeV). |  | $7.3$ | $25$
|  $E_{\nu}>1\;\mathrm{GeV}$222Uncertainty is 7 % up to 10 GeV, linearly increases with $\log E_{\nu}$ from 7 %(10 GeV) to 12 %(100 GeV) and then to 20 %(1 TeV) |  | $9.6$ | $15$
$(\nu_{\mu}+\bar{\nu}_{\mu})/(\nu_{e}+\bar{\nu}_{e})$ | $E_{\nu}<1\;\mathrm{GeV}$ | $1.3$ | $2$
| $1<E_{\nu}<10\;\mathrm{GeV}$ |  | $2.7$ | $3$
|  $E_{\nu}>10\;\mathrm{GeV}$333Uncertainty linearly increases with $\log E_{\nu}$ from 5 %(30 GeV) to 30 %(1 TeV). |  | $6.6$ | $5$
$\bar{\nu}_{e}/\nu_{e}$ | $E_{\nu}<1\;\mathrm{GeV}$ | $-1.1$ | $5$
| $1<E_{\nu}<10\;\mathrm{GeV}$ |  | $-2.7$ | $5$
|  $E_{\nu}>10\;\mathrm{GeV}$444Uncertainty linearly increases with $\log E_{\nu}$ from 8 %(100 GeV) to 20 %(1 TeV). |  | $-0.34$ | $8$
$\bar{\nu}_{\mu}/\nu_{\mu}$ | $E_{\nu}<1\;\mathrm{GeV}$ | $0.36$ | $2$
| $1<E_{\nu}<10\;\mathrm{GeV}$ |  | $3.24$ | $6$
|  $E_{\nu}>10\;\mathrm{GeV}$555Uncertainty linearly increases with $\log E_{\nu}$ from 6 %(50 GeV) to 40 %(1 TeV). |  | $9.2$ | $15$
Up/down ratio | $<400\;\mathrm{MeV}$ |  $e$-like | $0.079$ | $0.1$
|  |  $\mu$-like | $0.24$ | $0.3$
|  | 0-decay $\mu$-like | $0.87$ | $1.1$
| $>400\;\mathrm{MeV}$ |  $e$-like | $0.63$ | $0.8$
|  |  $\mu$-like | $0.40$ | $0.5$
|  | 0-decay $\mu$-like | $1.34$ | $1.7$
| Multi-GeV |  $e$-like | $0.55$ | $0.7$
|  |  $\mu$-like | $0.16$ | $0.2$
| Multi-ring Sub-GeV |  $e$-like | $0.32$ | $0.4$
|  |  $\mu$-like | $0.16$ | $0.2$
| Multi-ring Multi-GeV |  $e$-like | $0.24$ | $0.3$
|  |  $\mu$-like | $0.16$ | $0.2$
| PC |  | $0.16$ | $0.2$
Horizontal/vertical ratio | $<400\;\mathrm{MeV}$ |  $e$-like | $-0.023$ | $0.1$
|  |  $\mu$-like | $-0.023$ | $0.1$
|  | 0-decay $\mu$-like | $-0.069$ | $0.3$
| $>400\;\mathrm{MeV}$ |  $e$-like | $-0.32$ | $1.4$
|  |  $\mu$-like | $-0.44$ | $1.9$
|  | 0-decay $\mu$-like | $-0.32$ | $1.4$
| Multi-GeV |  $e$-like | $-0.74$ | $3.2$
|  |  $\mu$-like | $-0.53$ | $2.3$
| Multi-ring Sub-GeV |  $e$-like | $-0.32$ | $1.4$
|  |  $\mu$-like | $-0.30$ | $1.3$
| Multi-ring Multi-GeV |  $e$-like | $-0.64$ | $2.8$
|  |  $\mu$-like | $-0.35$ | $1.5$
| PC |  | $-0.39$ | $1.7$
K/$\pi$ ratio in flux calculation666Uncertainty increases linearly from 5$\%$ to 20$\%$ between 100GeV and 1TeV. | $-5.3$ | $10$
Neutrino path length | $-2.3$ | $10$
Sample-by-sample | FC Multi-GeV | $-0.22$ | $5$
| PC + Stopping UP-$\mu$ |  | $2.46$ | $5$
Matter effects | $0.06$ | $6.8$
Solar Activity | $0.31$ | $7$
$M_{A}$ in QE |  |  | $-4.79$ | $10$
Single $\pi$ Production, Axial Coupling | $-7.03$ | $10$
Single $\pi$ Production, $C_{A5}$ | $9.33$ | $10$
Single $\pi$ Production, BKG | $-3.03$ | $10$
CCQE cross section777Difference from the Nieves Nieves:2004wx  model is set to 1.0 | $7.37$ | $10$
CCQE $\bar{\nu}/\nu$ ratio777Difference from the Nieves Nieves:2004wx  model is set to 1.0 | $4.81$ | $10$
CCQE $\mu/e$ ratio777Difference from the Nieves Nieves:2004wx  model is set to 1.0 | $-2.74$ | $10$
DIS cross section | $5.72$ | $10$
DIS model comparisons888Difference from CKMT CKMT94  parametrization is set to 1.0 | $4.84$ | $10$
DIS $Q^{2}$ distribution (high W)999Difference from GRV98 Gluck:1998xa  is set to 1.0 | $-5.70$ | $10$
DIS $Q^{2}$ distribution (low W)999Difference from GRV98 Gluck:1998xa  is set to 1.0 | $-9.21$ | $10$
Coherent $\pi$ production | $-9.41$ | $100$
NC/CC | $-0.10$ | $20$
$\nu_{\tau}$ cross section | $9.08$ | $25$
Single $\pi$ production, $\pi^{0}/\pi^{\pm}$ | $1.75$ | $40$
Single $\pi$ production, $\bar{\nu}_{i}/\nu_{i}$ (i=$e,\mu$)101010Difference from the HernandezHernandez07  model is set to 1.0 | $2.30$ | $10$
NC fraction from hadron simulation | $-5.07$ | $10$
$\pi^{+}$ decay uncertainty Sub-GeV 1-ring |  |
|  $e$-like 0-decay | $-0.050$ | $0.6$
|  $\mu$-like 0-decay |  | $-0.066$ | $0.8$
|  $e$-like 1-decay |  | $-0.34$ | $4.1$
|  $\mu$-like 1-decay |  | $-0.025$ | $0.3$
|  $\mu$-like 2-decay |  | $-0.47$ | $5.7$
Final state and seconday interactions111111Ref. Pinzon:2017fsi  |  |
| 1-ring | $-6.83$ | $10$
| Multi-ring |  | $6.93$ | $10$
Meson exchange current121212Difference from NEUT without model from Nieves:2004wx  is set to 1.0 | $1.08$ | $10$
$\Delta m^{2}_{21}$ Olive:2016xmw  | $0.0$ | $0.0002$
$\sin^{2}(\theta_{12})$ Olive:2016xmw  | $0.05$ | $1.3$
$\sin^{2}(\theta_{13})$ Olive:2016xmw  | $0.00$ | $0.11$
FC reduction |  |  | $0.27$ | $1.3$
PC reduction | $-0.26$ | $1$
FC/PC separation | $0.000$ | $0.02$
PC stopping/through-going separation (bottom) | $-1.75$ | $6.8$
PC stopping/through-going separation (barrel) | $1.13$ | $8.5$
PC stopping/through-going separation (top) | $-29.0$ | $40$
Non-$\nu$ background | Sub-GeV $\mu$-like | $-0.001$ | $0.02$
| Multi-GeV $\mu$-like | $-0.002$ | $0.07$
| Sub-GeV 1-ring 0-decay $\mu$-like | $-0.001$ | $0.02$
| PC | $-0.015$ | $0.49$
| Sub-GeV $e$-like (flasher event) | $-0.000$ | $0.03$
| Multi-GeV $e$-like (flasher event) | $-0.001$ | $0.07$
Fiducial Volume | $-0.40$ | $2.00$
Ring separation | $<400\;\mathrm{MeV}$ |  $e$-like | $-0.44$ | $-0.86$
|  |  $\mu$-like | $0.43$ | $0.84$
| $>400\;\mathrm{MeV}$ |  $e$-like | $0.19$ | $0.37$
|  |  $\mu$-like | $1.03$ | $2$
| Multi-GeV |  $e$-like | $1.32$ | $2.57$
|  |  $\mu$-like | $-0.74$ | $-1.44$
| Multi-ring Sub-GeV |  $e$-like | $-0.42$ | $-0.82$
|  |  $\mu$-like | $0.42$ | $0.82$
| Multi-ring Multi-GeV |  $e$-like | $-0.66$ | $-1.28$
|  |  $\mu$-like | $0.59$ | $1.14$
Particle identification (1 ring) |  |
| Sub-GeV |  $e$-like | $-0.11$ | $0.52$
|  |  $\mu$-like | $0.10$ | $-0.49$
| Multi-GeV |  $e$-like | $-0.02$ | $0.09$
|  |  $\mu$-like | $0.02$ | $-0.09$
Particle identification (multi-ring) |  |
| Sub-GeV |  $e$-like | $0.29$ | $-1.33$
|  |  $\mu$-like | $-0.12$ | $0.53$
| Multi-GeV |  $e$-like | $-0.27$ | $1.25$
|  |  $\mu$-like | $0.16$ | $-0.73$
Multi-ring likelihood selection | Multi-ring $e$-like |  $\nu_{e}$,$\bar{\nu}_{e}$ | $-0.42$ | $-0.71$
| Multi-ring Other |  | $0.32$ | $0.52$
Energy calibration | $-0.24$ | $2.02$
Up/down asymmetry energy calibration | $-0.53$ | $0.67$
UP-$\mu$ reduction | Stopping | $0.048$ | $0.5$
| Through-going |  | $0.029$ | $0.3$
UP-$\mu$ stopping/through-going separation | $0.018$ | $0.6$
Energy cut for stopping UP-$\mu$ | $0.30$ | $1.7$
Path length cut for through-going UP-$\mu$ | $-0.33$ | $1.5$
Through-going UP-$\mu$ showering separation | $3.26$ | $3$
Background subtraction for UP-$\mu$ |  |
| Stopping131313The uncertainties in BG subtraction for upward-going muons are only for the most horizontal bin, $-0.1<\cos\theta<0$. | $-4.9$ | $11$
| Non-showering131313The uncertainties in BG subtraction for upward-going muons are only for the most horizontal bin, $-0.1<\cos\theta<0$. |  | $2.3$ | $17$
| Showering131313The uncertainties in BG subtraction for upward-going muons are only for the most horizontal bin, $-0.1<\cos\theta<0$. |  | $-6.55$ | $24$
$\nu_{e}/\bar{\nu}_{e}$ Separation | $1.89$ | $2.57$
Sub-GeV 2-ring $\pi^{0}$ | $-0.039$ | $1.03$
Decay-e tagging | $0.29$ | $0.7$