Prog. Theor. Exp. Phys. 2019, 053F01 (39 pages)
DOI: 10.1093/ptep/ptz015

M. Jiang1 , K. Abe2,3 , C. Bronner2 , Y. Hayato2,3 , M. Ikeda2 , K. Iyogi2 , J. Kameda2,3 , Y. Kato2 ,Y. Kishimoto2,3 , Ll. Marti2 , M. Miura2,3 , S. Moriyama2,3 , T. Mochizuki2 , M. Nakahata2,3 ,Y. Nakajima2,3 , Y. Nakano2 , S. Nakayama2,3 , T. Okada2 , K. Okamoto2 , A. Orii2 , G. Pronost2 ,H. Sekiya2,3 , M. Shiozawa2,3 , Y. Sonoda2 , A. Takeda2,3 , A. Takenaka2 , H. Tanaka2 , T. Yano2 ,R. Akutsu4 , T. Kajita3,4 , Y. Nishimura4 , K. Okumura3,4 , R. Wang4 , J. Xia4 , L. Labarga5 ,P. Fernandez5 , F. d. M. Blaszczyk6 , C. Kachulis6 , E. Kearns3,6 , J. L. Raaf6 , J. L. Stone3,6 ,S. Sussman6 , S. Berkman7 , J. Bian8 , N. J. Griskevich8 , W. R. Kropp8 , S. Locke8 , S. Mine8 ,P. Weatherly8 , M. B. Smy3,8 , H. W. Sobel3,8 , V. Takhistov8,9 , K. S. Ganezer10 , J. Hill10 ,J. Y. Kim11 , I. T. Lim11 , R. G. Park11 , B. Bodur12 , K. Scholberg3,12 , C. W. Walter3,12 ,M. Gonin13 , J. Imber13 , Th. A. Mueller13 , T. Ishizuka14 , T. Nakamura15 , J. S. Jang16 , K. Choi17 ,J. G. Learned17 , S. Matsuno17 , R. P. Litchﬁeld18 ,Y. Uchida18 , M. O. Wascko18 , N. F. Calabria19 ,M. G. Catanesi19 , R. A. Intonti19 , E. Radicioni19 , G. De Rosa20 , A. Ali21 , G. Collazuol21 ,F. Iacob21 , L. Ludovici22 , S. Cao23 , M. Friend23 , T. Hasegawa23 , T. Ishida23 , T. Kobayashi23 ,T. Nakadaira23 , K. Nakamura3,23 ,Y. Oyama23 , K. Sakashita23 , T. Sekiguchi23 , T. Tsukamoto23 ,KE. Abe24 , M. Hasegawa24 , Y. Isobe24 , H. Miyabe24 , T. Sugimoto24 , A. T. Suzuki24 ,Y. Takeuchi3,24 ,Y. Ashida1 , T. Hayashino1 , S. Hirota1 , T. Kikawa1 , M. Mori1 , KE. Nakamura1 ,T. Nakaya1,3 , R. A. Wendell1,3 , L. H. V. Anthony25 , N. McCauley25 , A. Pritchard25 ,K. M. Tsui25 , Y. Fukuda26 , Y. Itow27,28 , M. Murrase27 , P. Mijakowski29 , K. Frankiewicz29 ,C. K. Jung30 , X. Li30 , J. L. Palomino30 , G. Santucci30 , C. Vilela30 , M. J. Wilking30 ,C. Yanagisawa30,31 , D. Fukuda32 , K. Hagiwara32 , H. Ishino32 , S. Ito32 , Y. Koshio3,32 ,M. Sakuda32 ,Y. Takahira32 , C. Xu32 ,Y. Kuno33 , C. Simpson3,34 , D. Wark34,35 , F. Di Lodovico36 ,B. Richards36 , S. Molina Sedgwick36 , R. Tacik38,37 , S. B. Kim39 , M. Thiesse40 , L. Thompson40 ,H. Okazawa41 , Y. Choi42 , K. Nishijima43 , M. Koshiba44 , M. Yokoyama3,45 , A. Goldsack3,34 ,K. Martens3 , M. Murdoch3 , B. Quilain3 ,Y. Suzuki3 , M. R. Vagins3,8 , M. Kuze46 ,Y. Okajima46 ,T. Yoshida46 , M. Ishitsuka47 , J. F. Martin48 , C. M. Nantais48 , H. A. Tanaka48 , T. Towstego48 ,M. Hartz38 , A. Konaka38 , P. de Perio38 , S. Chen49 , L. Wan49 , A. Minamino50(The Super-Kamiokande Collaboration)
1

Department of Physics, Kyoto University, Kyoto, Kyoto 606-8502, Japan
Kamioka Observatory, Institute for Cosmic Ray Research, University ofTokyo, Kamioka, Gifu 506-1205, Japan
3
Kavli Institute for the Physics and Mathematics of the Universe (WPI), The University of Tokyo Institutes for
Advanced Study,University of Tokyo, Kashiwa, Chiba 277-8583, Japan
4
Research Center for Cosmic Neutrinos, Institute for Cosmic Ray Research, University of Tokyo, Kashiwa,
Chiba 277-8582, Japan
5
Department of Theoretical Physics, University Autonoma Madrid, 28049 Madrid, Spain
6
Department of Physics, Boston University, Boston, MA 02215, USA
7
Department of Physics and Astronomy, University of British Columbia, Vancouver, BC, V6T1Z4, Canada
8
Department of Physics and Astronomy, University of California, Irvine, Irvine, CA 92697-4575, USA
2© The Author(s) 2019. Published by Oxford University Press on behalf of the Physical Society of Japan.
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/),
which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.
Funded by SCOAP3

Downloaded from https://academic.oup.com/ptep/article-abstract/2019/5/053F01/5488936 by University of Tokyo Library user on 22 July 2019

Atmospheric neutrino oscillation analysis withimproved event reconstruction inSuper-Kamiokande IV

PTEP 2019, 053F01

M. Jiang et al.

9

Department of Physics and Astronomy, UCLA, CA 90095-1547, USA
Department of Physics, California State University, Dominguez Hills, Carson, CA 90747, USA
11
Department of Physics, Chonnam National University, Kwangju 500-757, Korea
12
Department of Physics, Duke University, Durham NC 27708, USA
13
Ecole Polytechnique, IN2P3-CNRS, Laboratoire Leprince-Ringuet, F-91120 Palaiseau, France
14
Junior College, Fukuoka Institute of Technology, Fukuoka, Fukuoka 811-0295, Japan
15
Department of Physics, Gifu University, Gifu, Gifu 501-1193, Japan
16
GIST College, Gwangju Institute of Science and Technology, Gwangju 500-712, Korea
17
Department of Physics and Astronomy, University of Hawaii, Honolulu, HI 96822, USA
18
Department of Physics, Imperial College London, London, SW7 2AZ, United Kingdom
19
Dipartimento Interuniversitario di Fisica, INFN Sezione di Bari and Università e Politecnico di Bari,
I-70125, Bari, Italy
20
Dipartimento di Fisica, INFN Sezione di Napoli and Università di Napoli, I-80126, Napoli, Italy
21
Dipartimento di Fisica, INFN Sezione di Padova and Università di Padova, I-35131, Padova, Italy
22
INFN Sezione di Roma and Università di Roma “La Sapienza”, I-00185, Roma, Italy
23
High Energy Accelerator Research Organization (KEK), Tsukuba, Ibaraki 305-0801, Japan
24
Department of Physics, Kobe University, Kobe, Hyogo 657-8501, Japan
25
Department of Physics, University of Liverpool, Liverpool, L69 7ZE, United Kingdom
26
Department of Physics, Miyagi University of Education, Sendai, Miyagi 980-0845, Japan
27
Institute for Space-Earth Environmental Research, Nagoya University, Nagoya, Aichi 464-8602, Japan
28
Kobayashi-Maskawa Institute for the Origin of Particles and the Universe, Nagoya University, Nagoya,
Aichi 464-8602, Japan
29
National Centre For Nuclear Research, 00-681 Warsaw, Poland
30
Department of Physics and Astronomy, State University of New York at Stony Brook, NY 11794-3800, USA
31
BMCC/CUNY, Science Department, NewYork, NewYork 10007, USA
32
Department of Physics, Okayama University, Okayama, Okayama 700-8530, Japan
33
Department of Physics, Osaka University, Toyonaka, Osaka 560-0043, Japan
34
Department of Physics, Oxford University, Oxford, OX1 3PU, United Kingdom
35
STFC, Rutherford Appleton Laboratory, Harwell Oxford, and Daresbury Laboratory, Warrington, OX11
0QX, United Kingdom
36
School of Physics and Astronomy, Queen Mary University of London, London, E1 4NS, United Kingdom
37
Department of Physics, University of Regina, 3737 Wascana Parkway, Regina, SK, S4SOA2, Canada
38
TRIUMF, 4004 Wesbrook Mall, Vancouver, BC, V6T2A3, Canada
39
Department of Physics, Seoul National University, Seoul 151-742, Korea
40
Department of Physics and Astronomy, University of Shefﬁeld, S3 7RH, Shefﬁeld, United Kingdom
41
Department of Informatics in Social Welfare, Shizuoka University of Welfare, Yaizu, Shizuoka, 425-8611,
Japan
42
Department of Physics, Sungkyunkwan University, Suwon 440-746, Korea
43
Department of Physics, Tokai University, Hiratsuka, Kanagawa 259-1292, Japan
44
The University of Tokyo, Bunkyo, Tokyo 113-0033, Japan
45
Department of Physics, University of Tokyo, Bunkyo, Tokyo 113-0033, Japan46
Department of Physics, Tokyo Institute of Technology, Meguro, Tokyo 152-8551, Japan47
Department of Physics, Faculty of Science and Technology, Tokyo University of Science, Noda, Chiba
278-8510, Japan48
Department of Physics, University of Toronto, ON, M5S 1A7, Canada
49
Department of Engineering Physics, Tsinghua University, Beijing, 100084, China
50
Faculty of Engineering, Yokohama National University, Yokohama, 240-8501, Japan∗
E-mail: mjiang@scphys.kyoto-u.ac.jp10...................................................................................................................
A new event reconstruction algorithm based on a maximum likelihood method has been developed for Super-Kamiokande. Its improved kinematic and particle identiﬁcation capabilities
enable the analysis of atmospheric neutrino data in a detector volume 32% larger than previous
analyses and increase the 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
2/39

Downloaded from https://academic.oup.com/ptep/article-abstract/2019/5/053F01/5488936 by University of Tokyo Library user on 22 July 2019

Received January 11, 2019; Revised February 15, 2019; Accepted February 15, 2019; Published May 14, 2019

PTEP 2019, 053F01

M. Jiang et al.

for the normal hierarchy, disfavoring the inverted hierarchy at 74% assuming oscillations at thebest ﬁt of the analysis.
...................................................................................................................
Subject Index

1.

C43, F20, F21, F23, H42

Introduction

2.

Neutrino oscillations

Neutrino oscillations are a result of the neutrino mass eigenstates differing from their weakinteraction (ﬂavor) eigenstates. For neutrinos in vacuum, their oscillation probability in the standard
three-ﬂavor paradigm can be described by six parameters, nearly all of which have been conﬁrmed to
3/39

Downloaded from https://academic.oup.com/ptep/article-abstract/2019/5/053F01/5488936 by University of Tokyo Library user on 22 July 2019

Neutrino oscillations have been conﬁrmed by a variety of experiments using both natural and artiﬁcial
sources. At present the data are well described assuming mixing among all three active neutrinos with
the Pontecorvo–Maki–Nakagawa–Sakata (PMNS) formalism [1,2]. Though most of its parameters
have been experimentally measured [3], the ordering of the mass states with the largest splitting(known as the mass hierarchy), the octant of the atmospheric mixing angle θ23 , and the value of its
CP-violating phase are unknown. These unresolved issues have been at the forefront of results fromthe T2K [4], NOvA [5], and Super-Kamiokande [6] (Super-K, SK) experiments and are the focus of
next-generation experiments planned in the USA [7], China [8], and Japan [9].
The atmospheric neutrino ﬂux 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 inﬂuence the oscillations of neutrinos passing through the Earth and produce a resonance enhancement of
the oscillation probability that depends on the values of θ13 , θ23 , and δCP . Importantly, this resonanceoccurs 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 ﬁnal states that complicate particle identiﬁcation
and where the ﬂux is orders of magnitude lower than its peak. Accordingly, recent Super-Kamiokande
results [6] 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 speciﬁc 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 identiﬁcation probability, and momentum resolution and doesso 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 kton·year
exposure with a 32% larger ﬁducial volume (FV) than previous Super-K analyses. Section 2 brieﬂyreviews atmospheric neutrino oscillations before the detector is introduced in Sect. 3. In Sect. 4 theatmospheric neutrino data reduction and event categorization are discussed. The new reconstruction
algorithm and its performance in comparison to the conventional reconstruction are detailed in Sect. 5.
Based on this performance improvement the expansion of the FV and corresponding systematic
errors are presented in Sects. 6 and 7, respectively. An analysis of the atmospheric neutrino data bythemselves is presented in Sect. 8.1 and is followed by an analysis in which θ13 is constrained by
the world data in Sect. 8.2, before concluding in Sect. 9.

PTEP 2019, 053F01

M. Jiang et al.

P(νμ → νμ ) ∼= 1 − 4 cos2 θ13 sin2 θ23 (1 − cos2 θ13 sin2 θ23 )

2 L
1.27m32
× sin2E
2 L
1.27m32
P(νμ ↔ νe ) ∼= sin2 θ23 sin2 2θ13 sin2
E(1)

using the assumption that m231 ≈ m232 since m221 is known to be small in comparison. In theseequations the neutrino travel length is represented by L (km) and its energy by E (GeV). The unit ofm2 is eV2 .
When neutrinos propagate in matter, the oscillation probabilities are modiﬁed due to interactions
with the electrons in the Earth. Indeed the forward scattering amplitude of νe differs from thoseof the other ﬂavors, producing an effective potential felt only by this species. This phenomenon is
known as the matter effect or the MSW effect of neutrino oscillations [10,11]. When neutrinos travel
through homogeneous matter, the oscillation parameters sin2 θ13 and m232 in Eq. 1 are replaced bytheir matter-equivalents [12]:sin2 2θ13,M =sin2 2θ13(cos 2θ13 − ACC /m232 )2 + sin2 2θ13

m232,M = m232 (cos 2θ13 − ACC /m232 )2 + sin2 2θ13 ,

(2)

(3)

√
where ACC = ±2 2GF Ne Eν and the electron density Ne is assumed to be constant. Here GF
represents the Fermi constant, and the sign of ACC is positive for neutrinos and negative forantineutrinos.
Resonant enhancement of these matter variables occurs when ACC /m232 = cos 2θ13 and isdependent upon the sign of m232 and whether the neutrino is a particle (ACC > 0) or an antiparticle (ACC < 0). Note that the enhancement only occurs for neutrinos if the hierarchy is normal(m232 > 0) whereas it only happens for antineutrinos if it is inverted (m232 < 0). Since atmospheric
neutrinos generally experience a varying matter proﬁle, and hence electron density Ne changes as
they travel through the Earth, they experience a variety of matter effects and similar enhancements
to the oscillation probability. The calculation of oscillation probability in this analysis takes suchvariation on matter density into consideration, with a simpliﬁed version of the preliminary reference
Earth model (PREM) (cf. Ref. [6]). Typically neutrinos with a few to ten GeV of energy that travel
through the Earth’s core experience matter effects the most strongly. From the equations above it canbe seen that atmospheric neutrinos are sensitive to sin2 2θ13 , sin2 θ23 , and the mass splitting m232
through the observation of the disappearance of νμ and appearance of νe .
4/39

Downloaded from https://academic.oup.com/ptep/article-abstract/2019/5/053F01/5488936 by University of Tokyo Library user on 22 July 2019

have non-zero values based on the results of reactor, atmospheric, solar, and long-baseline neutrino
experiments [3]. For brevity the details of the PMNS oscillation formalism will be omitted here;
more details about the treatment of neutrino oscillation in Super-K can be found in Ref. [6].
In vacuum the leading terms in the oscillations of atmospheric νe and νμ can be written as:
2 L
1.27m32
P(νe → νe ) ∼= 1 − sin2 2θ13 sin2
E

PTEP 2019, 053F01

M. Jiang et al.

3. The Super-Kamiokande detector

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 areseparated into three general categories, fully contained (FC), partially contained (PC), and upwardgoing muons (Up-μ).
For events classiﬁed as FC and PC, the neutrino interacts within the ﬁducial volume, deﬁned as
the region located more than 2 m from the ID wall in previous analyses. Events with no activity
in the outer detector are classiﬁed as FC. If energy deposition in the OD is observed, generally
from a high-energy muon exiting the ID, the event is classiﬁed 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-μ 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 speciﬁc to the three main sample categories are applied tothe data, details of which can be found in Ref. [23].
As an example, ﬁve steps of data reduction criteria are used for an FC sample. In the ﬁrst andsecond reduction steps, criteria on the total charge collected by ID and the number of hits on OD
5/39

Downloaded from https://academic.oup.com/ptep/article-abstract/2019/5/053F01/5488936 by University of Tokyo Library user on 22 July 2019

Super-Kamiokande is a cylindrical water Cherenkov detector with a total volume of 50 kilotons
and is 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 is 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 1885 8-inch PMTs coupled to wavelength-shifting plates are mounted on the inner surfaceof the OD, while its outer surface is covered with reﬂective 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 andended on 31 May 31 2018, when the detector began upgrade work for its next phase. Though the
topologies of the detector in the SK-III and -IV periods are the same, at the start of SK-IV the frontend electronics were upgraded to a system based on ASIC that uses a high-speed charge-to-time
converter [13]. After this upgrade Super-Kamiokande is able to collect all PMT hits above thresholdwithout incurring any dead time. More detailed descriptions of the detector and its electronics arepresented in Refs. [13–15].
Atmospheric neutrino interactions in the detector are simulated using the Honda et al. ﬂux calculation [16] and the NEUT [17] neutrino interaction software (version 5.4.0). Particles emerging
from the interactions are tracked through a simulation of the detector based on GEANT3 [18]. The
present work uses an updated version of NEUT relative to the previous atmospheric neutrino analysis(cf. Ref. [6]). In particular, charged-current quasi-elastic (CCQE) interactions in the new simulation
are on the local Fermi-gas model of Nieves [19,20], assuming an axial mass MA = 1.05 GeV/c2
and using the random phase approximation correction. Further, deep inelastic scattering (DIS) is
modeled using the GRV98 parton distribution function [21] and utilizing CKM matrix elements forthe calculation of structure functions. Corrections for low q2 scattering have been updated to thoseof Bodek and Yang [22], where q2 denotes the square of the transferred four-momentum.

PTEP 2019, 053F01

M. Jiang et al.

5.

FC event reconstruction

The Super-K event reconstruction algorithms determine an event’s physical properties such as theinteraction 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 SuperK, contributing to the discovery of atmospheric neutrino oscillation and has been used in both the K2K
and T2K experiments. APFit is a single-iteration ﬁtter based on the time and charge information of hit
PMTs [25]. The interaction vertex is reconstructed based on hit timing information after accounting
for the photon time of ﬂight. Then the direction of the ﬁrst ring found, 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 alikelihood 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
6/39

Downloaded from https://academic.oup.com/ptep/article-abstract/2019/5/053F01/5488936 by University of Tokyo Library user on 22 July 2019are 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 “ﬂasher” events
caused by internal corona discharge of a PMT. Failing tubes often produce multiple events withsimilar topologies as these discharges repeat. Typical ﬂasher events have a broad distribution in time
among hit PMTs and spatially similar hit distributions. Cuts based on the hit timing distribution andthe charge pattern correlation with other events that passed the third step reduction are applied toremove ﬂasher events. The ﬁfth reduction rejects the remaining cosmic ray muons, ﬂasher events,and electronic noise further using a more precise ﬁtter. Finally, the reconstructed vertex is required
to be within a speciﬁed ﬁducial volume to further reduce non-neutrino background events and toprevent 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
Sect. 6 in detail. An event’s visible energy (Evis ), which is deﬁned 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 arescanned by eye to determine the level of background contamination in the ﬁnal analysis sample andto estimate the uncertainty inherent in the reduction process.
After reduction, the FC data are subdivided 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 (μ-like), andsub-GeV (Evis < 1330.0 MeV) or multi-GeV (Evis > 1330.0 MeV) categories, as well as the numberof decay electrons (0, 1, 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-μ 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 preexisting algorithm and are divided into “stopping” and “through-going” subsamples based on theestimated muon stopping point for PC and Up-μ events. The “through-going” events in the Upμ sample are further divided into “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 the T2K experiment [24] and detects accelerator-generated neutrinos, they are excluded
based on event time information in this analysis. After all the selections there is a total of 13 FC
analysis samples, two PC samples and three Up-μ samples.

PTEP 2019, 053F01

M. Jiang et al.

5.1.

Likelihood function

An event topology hypothesis  (e.g., single-ring e-like) together with its associated kinematicparameters θ , 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 ﬁt. Based on the observed charge and hit time of each PMT, ﬁTQun constructs the following
likelihood function for a given hypothesis to estimate the kinematic variables:L(, θ ) =

unhit


Pj (unhit|, θ)

j

hit

{1 − Pi (unhit|, θ)}i× fq (qi |, θ)ft (ti |, θ).

(4)In this equation, the index j runs over all PMTs that did not register a hit (“unhit” PMTs) and for eachof these the probability that it does not register a hit given the ﬁtting hypothesis (, θ) is calculated
as Pj (unhit|, θ ). PMTs that did register a hit are indexed with i. For such PMTs the likelihooddensity for observing a charge qi under the ﬁtting hypothesis is represented by the charge likelihood
fq (qi |, θ ). The likelihood density of producing a hit at the observed time ti is deﬁned similarly asft (ti |, θ).
Since the processes of particle and optical photon propagation are decoupled from the responseof the PMT and the electronics, the likelihood can be rewritten in terms of the expected number ofphotoelectrons produced at the ith PMT given the hypothesis (the predicted charge), μi (, θ), asL(, θ ) =

unhit


Pj (unhit|μj )

j

hit

{1 − Pi (unhit|μi )}i

× fq (qi |μi )ft (ti |, θ).
7/39

(5)

Downloaded from https://academic.oup.com/ptep/article-abstract/2019/5/053F01/5488936 by University of Tokyo Library user on 22 July 2019ring pattern and opening angle. Rings from electrons tend to have rough edges produced by thelight from their electromagnetic showers, while rings from muons or charged pions predominantly
produce crisp edges. In the last step of the ﬁt, the momentum of each ring is evaluated based on
the observed charge inside a cone with a half angle of 70◦ drawn along the line connecting theinteraction 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 ﬁrst step to ﬁnd the initial vertex candidate in APFit.
The new algorithm, named ﬁTQun, 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 experiment [26], but has been developed from scratch for Super-K withadditional features such as multi-ring reconstruction for events with multiple ﬁnal-state particles.
Compared to APFit, ﬁTQun uses more information, including information from PMT hits outside ofthe expected Cherenkov cone and hit timing information, during the ﬁtting procedure. For a given
event ﬁTQun’s ﬁt procedure will run multiple times to determine the best kinematic parameters
for each possible particle conﬁguration hypothesis, while APFit only ﬁts those parameters once.
The remarkable evolution of computing power since the start of Super-Kamiokande has enabledﬁTQun to achieve higher reconstruction precision on similar timescales as APFit used to be. FiTQun
has already been used in the T2K analyses using an expanded ﬁducial volume due to its improved
resolution of reconstructed quantities and particle identiﬁcation [27].

PTEP 2019, 053F01

M. Jiang et al.

exp

ti

= t + smid /c + |RPMT− x − smid d|/cn ,i(6)where x and t are the vertex and creation time of the particle, respectively, d is the particle direction,is the position of the ith PMT, and smid represents half of the track length. Here c and cn
RPMTirepresent the light velocity in vacuum and in water, respectively. The time likelihood depends on
the predicted charge since a hit is recorded by the ﬁrst photon arriving at a PMT, which leads to anarrower distribution of hit time for higher numbers of incident photons and hence for more predictedcharge. The track length of a particle s, which is determined by the topology  and momentum p,
also affects the shape of the 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 thenmerged according to their relative intensities to obtain the ﬁnal 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 ﬁnal likelihood function assuming that the photons
expfrom a particle with earlier ti always arrive earlier than the photons from any other particles withexp
later ti values.
Once the likelihood function is deﬁned, the best set of kinematic parameters θ̂ for a given event
topology hypothesis is deﬁned as that which maximizes L(, θ). The best estimate for the particle
content of a given event is determined by comparison of L(, θ̂) among all hypotheses, .

5.2.

Fitting procedure

The ﬁTQun reconstruction process can be divided into four steps. The ﬁrst step, vertex pre-ﬁtting,roughly estimates the interaction vertex based on the PMT timing information. During the second
step, clusters of PMT hits in time are identiﬁed as candidates for particle activities. Thereafter,
the single-ring reconstruction, which performs ﬁts assuming event topologies with only a single
light-producting particle, and multi-ring reconstruction, which ﬁts using hypotheses with multiple
particles, are performed in sequence. During these ﬁts the negative log likelihood − ln L(, θ) is
minimized with respect to θ, using the MINUIT [28] package.
8/39

Downloaded from https://academic.oup.com/ptep/article-abstract/2019/5/053F01/5488936 by University of Tokyo Library user on 22 July 2019

Pj (unhit|μj ) and fq (qi |μ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 reﬂected (indirect light) are considered separately and summed to form the ﬁnal
μi . The predicted charge from direct light reaching a PMT is calculated by integrating the Cherenkov
emission proﬁle along the track while correcting for the distance from the track to the PMT, the lighttransmission 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 proﬁle and a scatteringfunction that has been generated in advance based on simulation and incorporates effects arisingfrom the relative position of the PMT and light source. For events with multiple Cherenkov rings thepredicted charge of each ring is ﬁrst calculated separately and then the sum from all the rings is used
to calculate the total expected μi . The ﬁnal charge likelihood fq (qi |μi ) is obtained by comparing theobserved charge in a PMT against the prediction assuming photoelectrons generated according toPoisson statistics.
exp
The time likelihood term can be expressed as ft (ti |ti , , p, μi ), where p is the momentum for theexp
topology  and ti is the expected hit time. The latter is deﬁned as the arrival time of unscatteredphotons emitted at the track midpoint and traveling directly to the PMT asPTEP 2019, 053F01

M. Jiang et al.

5.2.1. Vertex pre-ﬁtting
The vertex pre-ﬁtter is a fast algorithm that uses only the hit time information from PMTs around theprimary 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 x and time t that maximize
the goodness functione(−(Ti /σ ) /2) ,(7)− x|/cn
Tires = ti − t − |RPMTi(8)

G(x, t) =

res

2

i

where

is the residual hit time calculated on the assumption of a point-like light source and subtracting thephoton time of ﬂight. The value of parameter σ determines the precision of pre-ﬁtting. Here RPMTiis the position of the ith PMT. When the vertex and time get close to their true values by doinga grid search, the Tires distribute near zero, which results in a large goodness value. The pre-ﬁtter
is executed several times while gradually shrinking the grid size and σ to achieve high precision
efﬁciently. The ﬁtted vertex from this step, called the pre-ﬁt vertex, is just a rough estimation andwill be ﬁtted again with higher precision during minimization of L(, θ).
5.2.2. Hit clustering
Events in Super-K are deﬁned by detector activity in an O(10 μs) time window around an event
trigger but may contain multiple subevents representing clusters of PMT hits separated in time fromthe 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 subeventsfor further ﬁtting 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-ﬁnding algorithm. The vertex goodness G(x, t) from Eq. 7 is scanned in t withthe vertex position ﬁxed to the pre-ﬁt vertex to search for additional peaks from particle activity inthe detector. An example of the goodness distribution from a muon decay event is shown in Fig. 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 reﬂection processes within a clusterof hit activity, peaks are required to be above a minimum threshold F(t), which is deﬁned as:
F(t) := 0.25 argmax{i∈M

G(x, ti )} + η,(1 + ((t − ti )/γ )2(9)where M represents all local maxima of goodness function G(x, t). The value of the time constantγ is 25 ns when t < ti and 70 ns otherwise. An offset η = 9 is added to the threshold function to
suppress the effect of dark hits. The function F(t) is the blue curve in Fig. 1. The minimum goodness
between any two peaks must be lower than a second threshold, 0.6 × F(t), which is shown as a green
dashed curve. Under these criteria only the peaks labeled with a triangle are identiﬁed as candidates
for further ﬁtting.
During this procedure all vertex positions are assumed to lie close to the pre-ﬁt vertex when thepeak-ﬁnding algorithm runs. This assumption is broken when the primary particle travels a signiﬁcant
9/39

Downloaded from https://academic.oup.com/ptep/article-abstract/2019/5/053F01/5488936 by University of Tokyo Library user on 22 July 2019hit

PTEP 2019, 053F01

M. Jiang et al.

distance from the interaction vertex, as for high-momenta muons. Therefore, the vertex pre-ﬁtting
and peak-ﬁnding algorithm are rerun after masking the hits caused by the primary particle to improvedecay electron reconstruction efﬁciency. The vertex position used in the goodness function x is thenclose to the vertex of the secondary particle.
A time window deﬁned as −180 ns < Tires < 800 ns around each peak found is deﬁned to containits associated hits. Afterwards the vertex pre-ﬁtter and peak ﬁnder will be run once again in eachtime window using only the hits within each window. Peaks remaining after this step become theﬁnal candidates for full event reconstruction.

5.2.3. Single-ring reconstruction
The most basic reconstruction that is applied to the time windows deﬁned in the previous step is the
single-ring ﬁtter. This ﬁtter poses single-particle hypotheses for the likelihood function in Eq. 4. Threetypes of single-ring hypothesis are considered in ﬁTQun, 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 theobservation. Particle identiﬁcation (PID) is based on the best-ﬁt likelihood values for each of thesehypotheses. Electrons and muons, for instance, are separated by cutting on ln(Le /Lμ ), the logarithmof the likelihood ratio between the best-ﬁt electron and muon hypotheses. The distribution of thisvariable is shown for the FC atmospheric neutrino data and Monte Carlo (MC) simulation result with
sub-GeV and multi-GeV energies in Fig. 2. In both cases there is a clear separation of the likelihood
variable between electron-like (e-like) and muon-like (μ-like) events.
Figure 3 shows the vertex resolution as a function of visible energy for the FC single-ring chargedcurrent quasi-elastic (CCQE) event sample in the atmospheric neutrino MC whose true interaction
vertex lies within the conventional ﬁducial 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 ﬁTQun
for CCQE νe events is stable at 20.6 cm. For the CCQE νμ events, the vertex resolution of ﬁTQun,
10/39

Downloaded from https://academic.oup.com/ptep/article-abstract/2019/5/053F01/5488936 by University of Tokyo Library user on 22 July 2019

Fig. 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 ﬁxed to the pre-ﬁt 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.

PTEP 2019, 053F01(a)

M. Jiang et al.

(b)

(a)

(b)

Fig. 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 full triangles indicate the performanceof ﬁTQun and the open circles are for APFit. The true event vertex is required to be at least 200 cm away fromthe ID wall. The error bars represent the statistical uncertainty.

which changes from 29.2 cm to 15.9 cm, is better than that of APFit. In the same energy range,
the momentum resolution of ﬁTQun for CCQE νe events improves from 5.39% to 2.58% as thevisible energy increases, as shown in the top-left plot of Fig. 4. For CCQE νμ events the momentum
resolution is constant across the energy range, being lower than 2.5% for ﬁTQun and slightly worse
for APFit. Similarly, ﬁTQun shows an improved ability to discriminate between electrons and muons
with less than a 1% misidentiﬁcation rate for visible energies less than 1330 MeV (Fig. 5). Other
metrics used to measure ﬁt quality are generally stable and their typical values can be found in
Table 1. In general ﬁTQun performs as well as or better than APFit.
11/39

Downloaded from https://academic.oup.com/ptep/article-abstract/2019/5/053F01/5488936 by University of Tokyo Library user on 22 July 2019

Fig. 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 thenormal hierarchy assumed and the oscillation parameters taken to be m223 = 2.4 × 10−3 eV2 , sin2 θ23 = 0.5,
sin2 θ13 = 0.0210, and δCP = 0. The shaded histograms show charged-current νμ interactions. Error bars show
the statistical error. In this ﬁgure the reconstructed event vertex is required to be at least 200 cm away fromthe ID wall.

PTEP 2019, 053F01

M. Jiang et al.

(b)

(c)

(d)

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

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 ﬁnal states.
To save computing time, the ﬁTQun multi-ring ﬁtter is applied only to the time window around theprimary 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 froma previous ﬁt result. The hypothesis in the likelihood from Eq. 5 is updated to include a new ringand the likelihood is minimized again allowing the kinematic parameters for the rings to vary. Threehypotheses, including e-like, μ-like, and π + -like rings, are tested. The result of the updated ﬁt foreach of the particle hypotheses is compared to the original result to determine the validity of theadded ring. If the likelihood improves with the new ring it is accepted and the processes to search forfurther rings is repeated. The cycle of adding, ﬁtting, 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-ﬁt likelihoods from the two-ring and one-ring hypotheses. Based on an
optimization performed in MC, a cut at 9.35 (11.83) on this likelihood ratio is used to separate single12/39

Downloaded from https://academic.oup.com/ptep/article-abstract/2019/5/053F01/5488936 by University of Tokyo Library user on 22 July 2019(a)

PTEP 2019, 053F01(a)

M. Jiang et al.

(b)

Table 1. Summary of the basic performance of the APFit and ﬁTQun reconstruction algorithms on the fully
contained CCQE single-ring event sample with visible energy of 1 GeV.
Reconstruction
True CCQE νe sample
Vertex resolution
Direction resolution
Momentum bias
Momentum resolution
Mis-PID rate
True CCQE νμ sample
Vertex resolution
Direction resolution
Momentum bias
Momentum resolution
Mis-PID rate

ﬁTQun

APFit

20.6 cm
1.48◦
0.43%
2.90%
0.02%

24.9 cm
1.68◦
0.63%
3.56%
0.50%

15.8 cm
1.00◦−0.18%
2.26%
0.05%

17.3 cm
1.28◦
0.54%
2.60%
0.91%

and multi-ring events when the ﬁrst ring hypothesis is e-like (μ-like). However, ring candidates thathave an angular separation of less than 20◦ from the most energetic ring are discarded as spurious
since MC studies indicate that such candidates are typically due to particle scattering rather than a
new particle. When this occurs the two rings are merged and reﬁtted as one for all particle hypotheses
while keeping all other rings ﬁxed. This procedure is repeated for all rings in the event in descendingorder of their energy.
Table 2 shows the performance of the ring counting in both APFit and ﬁTQun using atmospheric
neutrino MC events. In this comparison only charged particles with energies more than 30 MeV above
the Cherenkov threshold in the ﬁnal states without any requirement on their angular separation withother particles are considered as true ring candidates. Compared to APFit, ﬁTQun 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 ﬁTQun tends to have more fake rings than APFit, theimprovements largely outweigh this slight downside.
13/39

Downloaded from https://academic.oup.com/ptep/article-abstract/2019/5/053F01/5488936 by University of Tokyo Library user on 22 July 2019

Fig. 5. Misidentiﬁcation 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 ﬁTQun. The reconstructed event vertex is required to be at least 200 cm away fromthe ID wall.

PTEP 2019, 053F01

M. Jiang et al.

(a)

(b)

Table 2. Summary table of ring counting performance of both APFit and ﬁTQun 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 energies more than 30 MeV above the Cherenkov
threshold in the ﬁnal states. The terms 1R, 2R, and 3R correspond to one ring, two rings, and three or morerings.
True number of rings

True 1R
True 2R
True ≥3R(a)ﬁTQun reconstruction
1R

2R

95.0%
27.8%
7.04%

4.64%
66.7%
25.5%

APFit reconstruction
≥ 3R
0.41%
5.56%
67.5%

1R

2R

95.9%
42.5%
20.2%

3.85%
52.8%
33.0%

≥ 3R
0.29%
4.63%
46.8%

(b)

Fig. 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 νμ charged-current interactions. Error bars are
statistical. The reconstructed event vertex is required to be at least 200 cm away from the ID wall.

14/39

Downloaded from https://academic.oup.com/ptep/article-abstract/2019/5/053F01/5488936 by University of Tokyo Library user on 22 July 2019

Fig. 6. The distribution of the likelihood ratio between the best-ﬁt 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 a single ring in the ﬁnal 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 atleast 200 cm away from the ID wall.

PTEP 2019, 053F01

M. Jiang et al.

(a)

(b)

Table 3. Summary table of performance of CC1π 0 events selection in three-ring events. The terms eee andμ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 π 0 . The purity column shows the purity of target
events in the passed invariant mass region.
ﬁTQun selection

Topology
eee events (target: νe CC1π )Sub-GeV
Multi-GeV
μee events (target: νμ CC1π 0 )Sub-GeV
Multi-GeV

APFit selection

Event rate

PurityEvent rate

Purity

278.5
112.1

52.1%
47.0%

175.9
32.8

42.9%
36.2%

384.5
143.5

54.6%
64.6%

201.2
51.6

38.0%
32.0%

0

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 μ-like is not as good as that for single-ring events.
The performance of the ﬁtter’s kinematic reconstruction on multi-ring events is checked by studying
the invariant mass of the π 0 in charged-current single-pion (CC1π 0 ) events. Figure 8 shows thereconstructed invariant mass calculated using the second and third most energetic rings for eventswith three e-like rings (eee events) for both reconstruction algorithms. Shaded histograms show trueCCνe 1π 0 events. Both ﬁtters have a peak in the invariant mass distribution near the π 0 mass, but
the ﬁTQun peak is larger, sharper, and has fewer backgrounds, indicating improved reconstruction
efﬁciency over APFit. The FWHM of the π 0 mass peak selected by ﬁTQun is 38 MeV, while this valueis 42 MeV for APFit. The event rate and the purity of the CCνe 1π 0 events with 85 < m < 185 MeV(indicated by the arrows in the ﬁgure) are 278.5 and 52.1% for ﬁTQun, which are much higher thanthe values for APFit: 175.9 and 42.9%. Similarly, CCνμ 1π 0 events from the three-ring sample whose
leading ring is μ-like and passing this selection are also studied and summarized in Table 3. The
ﬁTQun algorithm demonstrates a higher efﬁciency and purity for multi-ring event reconstruction.
15/39

Downloaded from https://academic.oup.com/ptep/article-abstract/2019/5/053F01/5488936 by University of Tokyo Library user on 22 July 2019

Fig. 8. Reconstructed invariant mass calculated using the second and third most energetic rings for eventswith three e-like rings by ﬁTQun (left) and APFit (right) from the sub-GeV atmospheric neutrino MC. Shaded
histograms show true CCνe 1π 0 events and arrows indicate the invariant mass range used to select events inwhich a π 0 is reconstructed properly. The reconstructed event vertex is required to be at least 200 cm awayfrom the ID wall.

PTEP 2019, 053F01

M. Jiang et al.

Fig. 10. Likelihood distribution used to separate multi-ring e-like events into neutrino-like and antineutrinolike subsamples. Error bars represent the statistical uncertainty of the data. Events with negative (positive)
likelihood values are designated ν̄e -like (νe -like).

Another important performance indicator for the multi-ring event reconstruction algorithm is theseparation of the neutrino and antineutrino components of the atmospheric neutrino sample, sinceat 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 ﬁrst stage of the separation is designed to extract CC νe + ν¯e interactions based on a likelihood
selection as in Ref. [29] using the APFit algorithm. Four variables, including the visible energy
fraction of the most energetic ring, the visible energy fraction of the most energetic charged pionlike ring, the number of decay electrons, and the distance to the farthest decay electron normalized
by the total visible energy, are used as the inputs to the likelihood function. However, in the present
study the inputs have been replaced by the equivalent variables from the ﬁTQun reconstruction.
Events that pass this selection are classiﬁed as “multi-ring e-like” while those that fail are termed“multi-ring other” as shown in Fig. 9. Both are used in the oscillation analysis discussed below.
16/39

Downloaded from https://academic.oup.com/ptep/article-abstract/2019/5/053F01/5488936 by University of Tokyo Library user on 22 July 2019

Fig. 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 aredesignated multi-ring e-like, while those with lower values are termed multi-ring other.

PTEP 2019, 053F01

M. Jiang et al.

Table 4. Summary of the multi-ring e-like event selection. The efﬁciency represents the fraction of events
passing the selection among all true charged-current νe events. Similarly, the purity is the fraction true CC νeout of all events passing the selection.
APFit

75.7%
77.8%

69.7%
69.5%

56.8%
58.8%

53.6%
52.6%

68.4%
30.0%

70.9%
25.9%

Neutrino and antineutrino interactions are separated from the multi-ring e-like sample during thesecond stage of the selection. The method is similar to the one discussed in Ref. [6] and uses thenumber of decay electrons, the number of reconstructed rings, and the event’s transverse momentum
as reconstructed by ﬁTQun. Figure 10 shows the ﬁnal likelihood distribution used in SK-IV. The
efﬁciency and purity for selecting e-like events and identifying true CC ν̄e (νe ) events as ν̄e -like(νe -like) are summarized in Table 4. FiTQun shows improved performance in each of these metrics
that will translate to better sensitivity to sin2 θ13 and the mass hierarchy.

6.

Fiducial volume expansion

In previous analyses using the APFit algorithm the ﬁducial volume (FV) was deﬁned as the region
located 200 cm from the wall of the inner detector. With this deﬁnition approximately 30% of theSuper-K inner detector mass is not used for atmospheric neutrino analysis. On the other hand, thesensitivity to several open questions is currently limited by a lack of statistics. The analysis sensitivity
can be increased by taking advantage of ﬁTQun’s improved reconstruction performance to expand
the ﬁducial volume while keeping signal purities high and backgrounds low.
The event selection and categorization in this study are similar to those one used in Ref. [6], but
all cut variables used to deﬁne the FC samples are reconstructed by ﬁTQun. Partially containedand Up-μ events, on the other hand, are deﬁned using APFit in the same manner as in the previous
publications.

6.1.

MC study for ﬁducial volume expansion

Though a larger ﬁducial volume provides higher statistics, the purity of signal interactions may
decrease if the reconstruction performance deteriorates for events close to the ID wall. Events thatare 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 tothe wall, the possibility of events with interaction vertices outside of the ID wall being reconstructedwithin the FV increases. Although the OD is designed to veto such events, as introduced in Sect. 4,
events that 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
17/39

Downloaded from https://academic.oup.com/ptep/article-abstract/2019/5/053F01/5488936 by University of Tokyo Library user on 22 July 2019

First stage
Multi-ring e-like events
Efﬁciency
PuritySecond stage
Multi-ring true CCνe events
Efﬁciency
PurityMulti-ring true CCν̄e events
Efﬁciency
PurityﬁTQun

PTEP 2019, 053F01

M. Jiang et al.

analysis sample. Such “entering events” are a type of misreconstructed background in the oscillation
analysis with potentially large systematic errors.
Figure 11 shows the purity of several components of the FC sub-GeV μ-like 0 decay-e sample as afunction 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, νμ CC and ν̄μ CC interactions (red) for this sample,
remains stable as the cut position is decreased to smaller values until it reaches 50 cm, where thebackground from entering ν events increases dramatically. All 13 FC event categories exhibit similar
behavior, indicating that a FV cut at 50 cm is an acceptable limit due to entering ν events. The purityof the 13 FC samples deﬁned using this cut and the conventional FV deﬁnition are summarized inTable 5 (cf. Table II of Ref. [6]).
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 fromthe standard SK analysis, which uses a cut at 200 cm. The rise in sensitivity as the FV is expanded
indicates that an FV cut at 50 cm can be expected to provide improved hierarchy sensitivity as longas 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.

6.2.

Background estimation

In addition to entering backgrounds and events with incorrectly assigned PID, non-neutrino (non-ν)
backgrounds from cosmic ray muons and ﬂashing 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 rejectedby cuts on activity in the OD. Any remaining events are rejected by the FV cut, since their true vertex
position, i.e., the entrance point of the muon, is expected to be on the wall of the 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.
18/39

Downloaded from https://academic.oup.com/ptep/article-abstract/2019/5/053F01/5488936 by University of Tokyo Library user on 22 July 2019

Fig. 11. The fraction of each component of the FC sub-GeV μ-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 nearestwall.

PTEP 2019, 053F01

M. Jiang et al.

Table 5. Sample purity broken down by neutrino ﬂavor assuming neutrino oscillations with m232 = 2.52 ×
10−3 eV2 and sin2 θ23 = 0.51 with an FV cut of 200 cm. The number in the parentheses is the result within the
new region, where the distance to the nearest wall is between 50 cm and 200 cm.
SampleCC νeFully contained (FC) multi-GeV
Single-ringνe -like
0.726(0.748)ν̄e -like
0.553(0.566)
μ-like0.000(0.000)Multi-ringνe -like
0.588(0.609)ν̄e -like
0.526(0.541)μ-like0.010(0.016)Other
0.283(0.302)CC νμ + ν̄μ

CC ντ

NC

0.242
(0.227)0.020(0.015)0.001(0.025)0.033(0.208)0.000(0.001)0.001(0.001)0.028(0.045)0.040(0.063)0.004(0.011)0.000(0.000)0.000(0.000)0.795(0.805)0.974(0.968)0.984(0.980)0.001(0.001)0.000(0.000)0.000(0.000)0.191(0.150)0.026(0.031)0.016(0.019)0.016(0.036)0.011(0.018)0.000(0.000)0.922(0.837)0.077(0.066)0.379(0.371)0.000(0.000)0.058(0.064)0.003(0.003)0.996(0.995)0.027(0.016)0.008(0.007)0.003(0.004)0.113(0.105)0.056(0.053)0.001(0.001)

0.117(0.112)0.300(0.301)0.001(0.002)0.026(0.032)0.054(0.059)0.021(0.023)0.959(0.946)0.342(0.342)0.036(0.032)0.020(0.016)0.004(0.005)0.053(0.051)0.204(0.188)0.134(0.118)0.026(0.031)
0.295(0.274)Flasher events are caused by the internal electrical discharge of a PMT and represent another kindof non-neutrino background. The reconstructed vertex for such events is also expected to be on thewall of the detector, so, as was the case for cosmic ray muons, more ﬂasher 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 modeledin the atmospheric neutrino MC, the distribution of the distance from the reconstructed vertex to thenearest wall of FC events is compared to MC in Fig. 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 forthe presence of a large unmodeled component to the data is observed.
19/39

Downloaded from https://academic.oup.com/ptep/article-abstract/2019/5/053F01/5488936 by University of Tokyo Library user on 22 July 2019

Fully contained (FC) sub-GeV
e-like, single-ring
0 decay-e
0.728(0.702)
1 decay-e
0.907(0.712)μ-like, single-ring
0 decay-e
0.010(0.034)
1 decay-e
0.000(0.001)
2 decay-e
0.000(0.000)π 0 -likeTwo-ring
0.051(0.109)CC ν̄ePTEP 2019, 053F01

M. Jiang et al.

Fig. 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 ﬁgure shows the ratio of the MC to the data.

To evaluate the background rate more precisely all FC events within the ﬁducial volume were
scanned by eye using a graphical event display program. In total, one ﬂasher event and 24 cosmic
ray muon events were found in the FC sample with dwall > 50 cm. Most of the cosmic ray muon
backgrounds were reconstructed as multi-ring downward-going μ-like events with more than 5 GeV/cof momentum. Eye scans of these events indicated that their vertices were incorrectly reconstructedwithin the new ﬁducial volume due to the presence of an erroneous second ring found by the ﬁtter.
Both the rate and type of such misreconstructions are consistent with those found in independent
20/39

Downloaded from https://academic.oup.com/ptep/article-abstract/2019/5/053F01/5488936 by University of Tokyo Library user on 22 July 2019

Fig. 12. Sensitivity to rejecting 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m223 = 2.4 × 10−3 eV2 , sin2 θ23 = 0.5, sin2 θ13 = 0.0210, and δCP = 0. Systematic errors are assumed to be
the same as those for the analysis with an FV cut at 200 cm.

PTEP 2019, 053F01

M. Jiang et al.

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

cosmic ray muon samples. Using the dwall distribution of the cosmic ray muon events, the expectedcosmic ray background in the ﬁnal analysis sample is statistically removed. This is done under theassumption that the FC event reduction process does not depend strongly on the muon kinematics.
By comparison of the event classiﬁcation of this control sample with that observed in the eyescanned background, this assumption has been validated for all particle directions and momenta.
As an example, Fig. 14 shows the dwall distribution for multi-ring μ-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 shapes of these two samples
are in good agreement. Due to their limited contamination in the ﬁnal FC sample, the effect of ﬂasher
PMTs is treated with a systematic error on their rate in the oscillation analysis presented below.
21/39

Downloaded from https://academic.oup.com/ptep/article-abstract/2019/5/053F01/5488936 by University of Tokyo Library user on 22 July 2019

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

PTEP 2019, 053F01

M. Jiang et al.

Fig. 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 partiallycontained (upward-going muon) event rates. The number in the parentheses for the fully contained event shows
the ﬁducial volume cut value. Orange (red) circles denote the fully contained event rate within the conventional
FV (new region).

22/39

Downloaded from https://academic.oup.com/ptep/article-abstract/2019/5/053F01/5488936 by University of Tokyo Library user on 22 July 2019

Fig. 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 ﬁfth column) and other samples are shown as reconstructed momentum
distributions (ﬁrst column). Cyan lines denote the best-ﬁt 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.

PTEP 2019, 053F01

M. Jiang et al.

Expanding the FV opens the possibility that events interacting within the ID but close to its
wall may produce particles that escape into the OD. Poor modeling of the response of the OD canthereby potentially introduce biases and relative inefﬁciencies in the FC sample. Figure 15 shows
the distribution of OD hits used to deﬁne the FC sample after the reduction processes for data withdwall > 50 cm. The detection inefﬁciency due to the cut on OD hits for both data and MC are
conﬁrmed to be consistent, around 0.2%. Based on this result in conjunction with the stability of thereconstruction algorithm, its systematic errors, and sample purities across the detector, the ﬁducial
volume deﬁnition is expanded from its conventional value to dwall > 50 cm (expanded FV) in theanalysis presented below.
Zenith angle distributions for each analysis sample using the expanded FV are shown in Fig. 16.
Their event rates since the start of SK-IV have been stable at 8.3 (2.2) FC events per day in theconventional FV (new region), 0.73 PC events per day, and 1.49 Up-μ events per day, as shown in
Fig. 17.

7.

Systematic error

While systematic errors related to the atmospheric neutrino ﬂux and cross section model are thesame as those used in Ref. [6], systematic errors on the event selection have been updated for FC
events using the expanded ﬁducial volume outlined above. Since ﬁTQun is not used to reconstruct
PC or Up-μ 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 acrossthree calibration samples spanning momenta up to 10 GeV/c: electrons from cosmic ray muon decay,
atmospheric neutrino events producing single π 0 s from neutral current interactions, and stopping
cosmic ray muons. The difference between data and MC for the calibration samples reconstructedwithin the expanded FV is shown in Fig. 18. The calculation of the total systematic error follows thatof Ref. [6] and is the sum in quadrature of the absolute energy scale error, which is deﬁned as thelargest data–MC difference across all samples, and momenta with the average of the time variationof these samples throughout SK-IV. The uncertainty from the up/down asymmetry of the detector is
23/39

Downloaded from https://academic.oup.com/ptep/article-abstract/2019/5/053F01/5488936 by University of Tokyo Library user on 22 July 2019

Fig. 18. Absolute energy scale error measurement for SK-IV with ﬁTQun 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.

PTEP 2019, 053F01

M. Jiang et al.

Table 6. Summary of uncertainty from the energy scale. The energy scale error is the quadratic sum of theabsolute energy scale error and the time variation. The absolute energy scale error value is obtained from themost discrepant sample, i.e., Michel electrons from cosmic ray muons in this study.
ﬁTQun>50 cmﬁTQun> 200 cm

APFit [6]> 200 cm

Energy calibration error
Absolute energy scale error
Time variation

2.02%
1.92%
0.62%

2.17%
2.09%
0.59%

2.1%

Up/down asymmetry energy calibration error

0.67%

0.58%

0.4%

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 deviateddirection is taken as the uncertainty from the detector asymmetry. As summarized in Table 6, theuncertainty of the energy scale on the conventional (expanded) FV is 2.17% (2.02%) using ﬁTQunreconstruction, which is similar to the value with APFit: 2.1% [6].
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 identiﬁcation (PID),
estimated number of rings (ring counting), and the two-stage separation of the multi-GeV multi-ringe-like event samples described above are evaluated using the atmospheric neutrino data by comparing
the likelihood distributions of the atmospheric neutrino data and MC. The MC distributions are shifted
and smeared during a ﬁt to data distribution. The fractional change in the MC after the cut used in
the likelihood selection before and after ﬁtting is deﬁned as the systematic error.
Systematic errors are estimated for each of the FC analysis samples for both the conventional
ﬁducial 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 isassumed between the systematic error categories. The sizes of the systematic errors evaluated using
ﬁTQun and APFit are consistent within the conventional ﬁducial volume, with improvements seen
in ﬁTQun for many samples. For many error categories the ﬁTQun reconstruction shows consistenterrors between the conventional ﬁducial 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 beneﬁt to theoscillation sensitivity from including events in the new region in the analyses. A quantitative estimation of the sensitivity is presented in the next section. The numbers for the different FV regions
are merged to form the ﬁnal response function used in the analysis, which will be described in thenext section. Systematic errors and their sizes at the best-ﬁt point of the analysis are presented inTable A1.

8. Atmospheric neutrino oscillation analysis
The Super-K atmospheric neutrino data are separately ﬁtted against both the normal and inverted
hierarchy hypotheses. Since θ13 has been measured and is well constrained by reactor experiments,
two types of analyses are done: one with θ13 as a free ﬁtting parameter and the other with it constrained
24/39

Downloaded from https://academic.oup.com/ptep/article-abstract/2019/5/053F01/5488936 by University of Tokyo Library user on 22 July 2019

Reconstruction
Dwall range

PTEP 2019, 053F01

M. Jiang et al.

Table 7. Summary of systematic errors related to the ring counting, particle identiﬁcation, and multi-GeV
multi-ring e-like separation. The sign of the number represents the correlation among samples, with negative
numbers denoting fully anti-correlated samples.
Reconstruction
Dwall range

Particle ID of single-ring events
Sub-GeV
e-likeμ-likeMulti-GeV
e-likeμ-likeﬁTQun> 200 cm

APFit
> 200 cm

1.94%
0.59%
1.08%
1.25%

1.20%
0.48%
0.42%
1.21%

1.6%
1.0%
−3.0%
0.6%

−3.58%−2.32%−2.39%−1.66%

−1.9%
2.3%

8.61%
2.65%

1.21%
−2.33%

1.0%−1.2%−4.21%−3.07%−0.62%
0.72%

−0.9%
2.4%

0.99%
−0.89%

0.36%−0.37%−0.28%
0.22%

0.23%−0.21%

0.06%−0.06%−0.35%
0.35%

Particle ID of the brightest ring in multi-ring events
Sub-GeV
e-like−3.19%
μ-like
1.31%
Multi-GeV
e-like1.94%
μ-like−1.06%−0.72%
0.31%

4.19%−1.56%

1.10%−0.66%

3.33%−1.56%

Multi-GeV, multi-ring e-like—other separatione-like−0.88%
other
0.50%−0.67%
0.53%

3.0%
−3.4%

Multi-GeV, multi-ring νe -like—ν̄e -like separation
νe -like−3.64%
ν¯e -like
4.51%−2.33%
2.10%

6.82%−6.04%

to 0.0210±0.0011 [3]. Though the constraint on this parameter is much stronger than can be expectedfrom Super-K alone, the appearance of upward-going electron events characteristic of the mass
hierarchy signal is driven by θ13 and therefore measurements with atmospheric neutrinos represent
an important test of the analysis’ hierarchy preference. The solar mixing parameters, m221 andsin2 θ12 , on the other hand, have little impact on the oscillations of neutrinos in the energy range
considered here and will be ﬁxed in the analysis as described below.
25/39

Downloaded from https://academic.oup.com/ptep/article-abstract/2019/5/053F01/5488936 by University of Tokyo Library user on 22 July 2019

Ring counting
Sub-GeV, single-ringe-like, p < 400 MeV
e-like, p > 400 MeV
μ-like, p < 400 MeV
μ-like, p > 400 MeV
Sub-GeV, multi-ringe-likeμ-like, Evis > 600 MeV
Multi-GeV, single-ringe-likeμ-likeMulti-GeV, multi-ringe-likeμ-likeﬁTQun
50–200 cm

PTEP 2019, 053F01

M. Jiang et al.
Table 8. Parameter grid used in the ﬁt.

ParameterRange|m232,31 |sin2 θ13
sin2 θ23−3

Scan points

1.5 × 10 –3.5 × 10
0.0–0.06
0.3–0.725
0.0–2πδCP−340 points
25 points
35 points
19 points

Parameter

Value

m221
sin2 θ12
sin2 θ13(7.53 ± 0.18) × 10−5 eV20.307 ± 0.013
0.0210 ± 0.0011

Data are ﬁtted to the MC expectation with a binned χ 2 method assuming Poisson statistics. The
effects of systematic errors are regarded as scaling factors on the MC in each bin [30], with different
error sources assumed to be independent. The deﬁnition of χ 2 used in the ﬁts is   2

On
i
2
χ =2
En − On + On ln
+
,
(10)
En
σi
n
i

where
En = En,0 (1 +

fni i ).

(11)i

In this equation On is the observed number of events in the nth analysis bin. Similarly, En,0 represents
the nominal MC expectation in that bin and En is the expectation after incorporating the effect of
the systematic errors. The coefﬁcient fni describes the fractional change in the bin’s MC under a 1σivariation of the ith systematic error source. The sum over ( i /σi )2 penalizes the χ 2 for adjusting thesystematic errors when bringing the MC into agreement with the data. For each set of oscillation
parameters tested, the MC are ﬁtted to the data by minimizing χ 2 over the error parameters i . Allﬁts are performed over 515 analysis bins in SK-IV.

8.1. Analysis with unconstrained θ13
The ﬁt with θ13 unconstrained by reactor measurements is performed over four parameters, |m232,31 |,sin2 θ13 , sin2 θ23 , and δCP with a total of 71 systematic errors. Fits to the normal hierarchy use m232
while those for the inverted hierarchy use m231 . The agreement between the data and MC is evaluatedusing Eq. 10 at each point in the grid of points shown in Table 8, while keeping the solar mixing
parameters, m212 and sin2 θ12 , ﬁxed to the values in Table 9. Systematic errors representing theuncertainty in these parameters are included as Eq. 10 in the analysis. All systematic errors havebeen proﬁled over in the ﬁgures presented below. The best-ﬁt parameters for each hierarchy aredeﬁned as the parameter set returning the smallest χ 2 value. Between the two hierarchy ﬁts the one
with smallest minimum χ 2 value is taken as the global best ﬁt and is taken to indicate the mass
hierarchy preference.
26/39

Downloaded from https://academic.oup.com/ptep/article-abstract/2019/5/053F01/5488936 by University of Tokyo Library user on 22 July 2019

Table 9. Values of oscillation parameters ﬁxed in the analysis and their systematic errors. Note that sin2 θ13 is
only ﬁxed in the “θ13 constrained”analyses described in Sect. 8.2.

PTEP 2019, 053F01

M. Jiang et al.

Table 10. νμ + ν̄μ CC backgrounds in a multi-GeV e-like sample.
Reconstruction
Dwall range
Single-ringνe -likeν̄e -like
Multi-ringνe -likeν̄e -like

ﬁTQun
50–200 cm

ﬁTQun> 200 cm

APFit
> 200 cm

0.064
0.003

0.058
0.003

0.109
0.009

0.059
0.023

0.054
0.021

0.121
0.040

Fig. 20. Constraints on neutrino oscillation parameters from the Super-K atmospheric neutrino data with noassumed constraint on θ13 . The solid blue and dashed orange lines denote the normal and inverted hierarchy ﬁtresults, respectively. The latter has been offset from the former by the difference in their minimum χ 2 values.

Results and discussion
Figure 19 shows the expected lower limit of sin2 θ13 at 1σ as a function of the true value of sin2 θ23
for the ﬁTQun-based analysis samples with both the conventional and expanded FV as well as ananalysis with APFit-reconstructed samples with the conventional FV. The sensitivities for the twoalgorithms using the same conventional FV are similar since the fraction of events common between
them is larger than 97%. The expected numbers of atmospheric neutrino events selected by APFit
27/39

Downloaded from https://academic.oup.com/ptep/article-abstract/2019/5/053F01/5488936 by University of Tokyo Library user on 22 July 2019

Fig. 19. Expected lower limit of sin2 θ13 at 68.3% conﬁdence level assuming that its true value is 0.0210 asa function of the true value of sin2 θ23 for a livetime of 3118.5 days. Gray and blue lines show the sensitivity
for samples reconstructed with ﬁTQun and using the conventional FV and the expanded FV, respectively. Theorange line denotes the sensitivity when events are reconstructed with APFit in the conventional FV.

PTEP 2019, 053F01

M. Jiang et al.

Table 11. Parameter estimates for each analysis and hierarchy hypothesis considered. Here NH (IH) refers to
the normal (inverted) hierarchy ﬁt. The terms “free” and “constrained” refer to ﬁts without and with a constraint
on sin2 θ13 , respectively. For sin2 θ23 parameter ranges are shown for both octants, with the best-ﬁtting octant
enclosed in a box. The expected absolute χ 2 value for the constrained ﬁt is 523.6 and the probability for
obtaining a larger value is 0.244 for NH.
θ13 freeθ13 constrained

NH

IH

NH

IH

χ
sin2 θ13
sin2 θ23 (1st oct.)sin2 θ23 (2nd oct.)|m232,31 | [×10−3 eV2 ]
δCP

576.3
0.008+0.025−0.0050.425+0.051−0.034
0.588+0.030−0.062
2.63+0.10−0.21
3.84+2.00
−2.14

578.1
0.008+0.015−0.0070.425+0.075−0.0270.575+0.034−0.0752.53+0.14−0.08
4.19+2.09
−4.19

576.5
–
0.425+0.046−0.0370.600+0.013−0.0302.53+0.22−0.12
3.14+2.67
−1.35

579.0
–
0.425+0.055−0.036
0.588+0.022−0.0372.53+0.14−0.31
4.89+1.51−3.46

2

and ﬁTQun are 24 188.7 and 24 184.3, respectively, while the number of the events passing both cutsof the two algorithms is 23 520.2. The 3% difference is mainly due to the vertex resolution near the
FV boundary.
However, since the contamination of νμ events in the e-like samples most sensitive to θ13 is lower forthe ﬁTQun-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 isclear.
One dimensional allowed regions for θ13 , |m232,31 |, sin2 θ23 , and δCP from the ﬁt to the data areshown in Fig. 20. The normal hierarchy hypothesis yields a better data–MC agreement than the22
inverted hierarchy hypothesis with χNH,min− χIH,min= −1.81. The 1σ allowed region for sin2 θ13 isfrom 0.003 to 0.033 (from 0.001 to 0.023) for the normal (inverted) hierarchy ﬁt, which is consistent
with the globally preferred value. The point at sin2 θ13 = 0.0 is disfavored at approximately 1.8σ
(1.2 σ ) for the normal (inverted) ﬁt. A summary of the best-ﬁt information and parameter constraints
is presented in Table 11.
As discussed in Sect. 2, the determination of sin2 θ13 and the mass hierarchy using atmospheric
neutrinos is achieved using the 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-likesamples, where the asymmetry is deﬁned as the ratio of the difference of the number of upward-goingand downward-going events relative to their sum. Here upward-going (downward-going) events are
events whose zenith angle satisfy cos θz < −0.4 (cos θz > 0.4). Excesses between a few and ten
GeV in the multi-GeV e-like νe and the multi-ring e-like νe samples drive the normal hierarchy
preference.
The best-ﬁt atmospheric mixing parameters from the normal hierarchy ﬁt are m232 = 2.63+0.10−0.21 ×+0.05110−3 eV2 and sin2 θ23 = 0.588+0.030
forthesecondoctant(best-ﬁt)
and
0.425
for
the
ﬁrst
−0.062−0.0342octant. Maximal mixing (sin θ23 = 0.5) is weakly disfavored at around 1σ signiﬁcance. Besides thedata excesses in the upward-going regions of the single-ring e-like νe sample and multi-ring othersample, the data deﬁcits in the multi-GeV μ-like sample also contribute to this preference. It should
be noted that the preference for sin2 θ23 is coupled to the θ13 measurement, since both parameters
feature in the νμ → νe oscillation probability. The jagged nature of the |m232 | curve is due to thelimited statistics at low energies in both data and MC. In this region, the oscillation probability is
28/39

Downloaded from https://academic.oup.com/ptep/article-abstract/2019/5/053F01/5488936 by University of Tokyo Library user on 22 July 2019

Hierarchy

PTEP 2019, 053F01

M. Jiang et al.

roughly proportional to sin2 (1.27m232 L/E) and changes wildly for small changes in the neutrinoenergy, E, or the neutrino travel length, L. Since there is no sufﬁcient statistics to fully sample andaverage over these oscillations, a rapid change in the χ 2 is observed in the ﬁgure for small changesin m232 .
The best-ﬁt value of δCP is 3.84 radians for both the normal and inverted hierarchy ﬁts. Comparing
with the least preferred parameter value of 0.8 radians, more electron neutrino appearances in the
sub-GeV e-like samples are expected to be observed due to νμ → νe oscillations. In the multi-GeV
region, δCP similarly modulates the amount of electron neutrino appearance but it is subdominant to
the effects of θ13 . In the inverted hierarchy ﬁt, the data are less sensitive to θ13 and δCP since only
antineutrino events experience matter effects in the Earth while there are more neutrino events thanantineutrino events.

8.2. Analysis with constrained θ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 theneutrino 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 hierarchyand δCP . Indeed, Eq. 2 shows that the size of the MSW resonance relies on the value of sin2 θ13 directly.
Since reactor neutrino experiments constrain this parameter more precisely than the analysis in the
29/39

Downloaded from https://academic.oup.com/ptep/article-abstract/2019/5/053F01/5488936 by University of Tokyo Library user on 22 July 2019

Fig. 21. Ratio of upward- (cos θ < −0.4) to downward-going (cos θ > 0.4) events of multi-GeV e-like events
as a function of energy for the mass hierarchy-sensitive analysis samples. The error bars are statistical. The
momentum estimator for multi-ring samples is the sum of the momenta of all rings. The orange line shows
the best-ﬁt result assuming the inverted hierarchy hypothesis and the cyan that from the normal hierarchy
hypothesis.

PTEP 2019, 053F01

M. Jiang et al.

Fig. 23. Expected sensitivity to rejecting the wrong octant as a function of the true value of sin2 θ23 for thetrue normal mass hierarchy. The other conditions are the same as in Fig. 22.

previous section, the analysis presented here is performed with the value of sin2 θ13 constrained
to 0.0210 ± 0.0011 [3]. During the ﬁt 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 ﬁt uses the same data
samples, binning, and parameter grid for all other oscillation parameters as the analysis presented
above.

Results and discussion
The expected sensitivities to the mass hierarchy and θ23 octant for the true normal hierarchy are illustrated in Figs. 22 and 23 for the ﬁTQun-reconstructed analysis samples using both the conventionaland expanded FV and for the APFit-based analysis with the conventional FV. Improved sensitivity
is found with the ﬁTQun-based analyses regardless of the ﬁducial volume used for the same reasons
as produced the improved sensitivity to θ13 in the previous section.
30/39

Downloaded from https://academic.oup.com/ptep/article-abstract/2019/5/053F01/5488936 by University of Tokyo Library user on 22 July 2019

Fig. 22. Expected sensitivity to the normal mass hierarchy (left) and inverted hierarchy (right) as a function of
the true value of sin2 θ23 . Here sin2 θ13 = 0.0210 ± 0.0011 and the assumed livetime is 3118.5 days. Gray andblue bands show the sensitivity of the analysis with event samples reconstructed with ﬁTQun 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 widths of the bands correspond to the uncertainty from δCP .

PTEP 2019, 053F01

M. Jiang et al.

Fig. 25. Constraints on the atmospheric mixing parameters using the SK-IV atmospheric neutrino and the
expanded FV. The solid blue (dashed orange) line shows the 90% CL for the normal (inverted) hierarchy. The
star denotes the best-ﬁt value, which is at the same point for the normal and inverted hierarchies, as shown in
Table 11. In each contour sin2 θ13 is constrained to be 0.0210 ± 0.0011. The contours have both been drawn
relative to the global best ﬁt.

Figure 24 shows the χ 2 value as a function of the atmospheric neutrino mixing parameters andδCP in the θ13 -constrained ﬁt. As in the unconstrained ﬁt the data prefer the normal mass hierarchy
with χ 2 = −2.45. The 1σ allowed region for |m232 | is from 2.41–2.75 ×10−3 eV2 (2.36–2.67
×10−3 eV2 ) for the normal (inverted) hierarchy ﬁt, which is consistent with the result of the unconstrained ﬁt. The contour for the allowed region with 90% CL is shown in Fig. 25. The result of theθ13 -constrained ﬁt shows a weak preference for the ﬁrst octant of θ23 , while it is in good agreementwith other measurements, as shown in Fig. 26. Adding the constraint on sin2 θ13 , the preferred valueof θ23 shifts from the second to the ﬁrst octant, though both are allowed at 1σ . Indeed, the difference2in χ 2 between the minimum value of each octant is χ1st−2nd= −0.73 (0.13) for the constrained(unconstrained) ﬁt. This result can be explained by degeneracies in the appearance probability thatarise for certain combinations of sin2 2θ13 and sin2 θ23 as shown in Eq. 1; an increase in one can be
31/39

Downloaded from https://academic.oup.com/ptep/article-abstract/2019/5/053F01/5488936 by University of Tokyo Library user on 22 July 2019

Fig. 24. Constraints on neutrino oscillation parameters from SK-IV atmospheric neutrino data using the
expanded FV and assuming sin2 θ13 = 0.0210 ± 0.0011. The solid blue and dashed orange lines denotethe normal and inverted hierarchy ﬁt results, respectively. The latter has been offset from the former by thedifference in their minimum χ 2 values.

PTEP 2019, 053F01

M. Jiang et al.

Fig. 27. Distributions of the difference in best-ﬁt χ 2 values between ﬁrst octant and second octant ﬁts to pseudodata sets used in the generation of the CLOs value for the SK θ13 constrained analysis. In the cyan (orange)
histogram the pseudo data have been generated assuming sin2 θ23 = 0.4 (sin2 θ23 = 0.6). Shaded portions of
the histograms denote the fraction of pseudo data sets with more extreme values than those observed in the2= −0.73.
data, χdata

compensated by a decrease in the other to describe the data. The same upward-going event excesses(deﬁcits) in the e-like (μ-like) samples that drive the hierarchy, θ13 , and octant preferences in theunconstrained analysis are not strong enough to support the increased excesses (deﬁcits) expected forhigher values of θ13 and the second octant. Since the θ13 is now ﬁxed, the data can be accommodated
by reducing the expected number of e-like events and increasing the number of μ-like events bymoving θ23 to the ﬁrst octant.
The best-ﬁt value of δCP is 3.14 (4.89) for the normal (inverted) hierarchy hypothesis, with a tighter
constraint compared with the unconstrained ﬁt. Parameter values and their 1σ errors are summarizedin Table 11.
The CLs method [34] is used to address the signiﬁcance of the octant and mass hierarchy preferences observed in the data. As in the previous Super-K study [6] 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
32/39

Downloaded from https://academic.oup.com/ptep/article-abstract/2019/5/053F01/5488936 by University of Tokyo Library user on 22 July 2019

Fig. 26. Constraints on the atmospheric mixing parameters at the 90% CL. The ﬁTQun result (solid red) is takenfrom the analysis with the normal mass hierarchy hypothesis and sin2 θ13 constrained to be 0.0210 ± 0.0011.
Constraints from the Super-K I–IV with APFit (dashed cyan) [6], T2K (dotted yellow) [4], NOvA (dashedgreen) [31], IceCube (dashed black) [32], and MINOS+ (dashed blue) [33] experiments are also shown.

PTEP 2019, 053F01

M. Jiang et al.

Table 12. Summary of the probability of observing a χ 2 preference for the NH more extreme than that observedin the data assuming an IH, p0 (IH), and CLHs values for a range of assumed parameters. The ﬁrst row shows
the true θ23 used to generate MC ensembles used in the calculations. Other oscillation parameters are takenfrom the analysis best ﬁt. The value of θ13 is constrained to 0.0210 ± 0.0011.
True sin2 θ23

0.4

0.425

0.5

0.6

p0 (IH)
CLHs

0.025
0.308

0.033
0.260

0.065
0.229

0.072
0.143

the current detector exposure and systematic errors have been varied as Gaussian parameters withwidths speciﬁed by their uncertainties. The CLs parameter for the octant study, CLO
s , is deﬁned as:CLOs =p0 (O2 ),1 − p0 (O1 )(12)where p0 (Oi ) represents the p-value with respect to obtaining a difference in the best-ﬁt χ 2 values
between both octant hypotheses smaller (larger) than the value of data assuming that the true valueof sin2 θ23 is O2 = 0.6 (O1 = 0.4). These pseudo MC distributions of these hypotheses are shownas the orange and cyan shaded histograms in Fig. 27, respectively. The CLO
s value is found to be just
0.098 and, accordingly, the data show no strong preference for the octant.
Similarly, the parameter for the mass hierarchy is deﬁned asCLHs =p0 (IH),1 − p0 (NH)(13)

where p0 (IH) and p0 (NH) are p-values for obtaining a difference in the χ 2 values of the best-ﬁt mass
hierarchies more extreme than those of the data assuming a true IH and NH, respectively.
Figure 28 shows the distribution for the mass hierarchy determination. Due to the large uncertaintyon θ23 , MC ensembles have been generated with different assumed values θ23 and with the otheroscillation parameters ﬁxed to their best-ﬁt values. The p-values and CLH
s values for the hierarchy
2
test are summarized in Table 12. As sin θ23 ranges from 0.4 to 0.6 the observed CLH
s values decrease
from a 30.8% CL preference for the inverted mass hierarchy to 14.3%. At the analysis best-ﬁt value
33/39

Downloaded from https://academic.oup.com/ptep/article-abstract/2019/5/053F01/5488936 by University of Tokyo Library user on 22 July 2019

Fig. 28. Distributions of the difference in best-ﬁt χ 2 values between normal and inverted hierarchy ﬁts topseudo data sets used in the generation of the CLHs value for the SK θ13 constrained analysis. In the cyan(orange) histogram the pseudo data have been generated assuming the normal (inverted) hierarchy at theanalysis best ﬁt shown in Table 11. Shaded portions of the histograms denote the fraction of pseudo data sets2
with more extreme values than those observed in the data, χdata= −2.45.

PTEP 2019, 053F01

M. Jiang et al.

the preference is 26.0%. A smaller value of sin2 θ23 predicts fewer electron neutrino appearances,which results in smaller p0 (IH) and larger CLH
s values due to the multi-GeV e-like sample event
excesses seen in the data.

9.

Conclusion

A new reconstruction algorithm based on a maximum likelihood method has been developed forSuper-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 identiﬁcation. The new algorithm has further demonstrated reliable performance overa larger volume of the detector than the previous volume. Accordingly, in the present analysis theﬁducial volume deﬁnition 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 ﬁducial volume deﬁnition 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 a normal mass hierarchy, the constraints on the atmospheric mixing parameters are+0.013+0.22
2
sin2 θ23 = 0.425+0.046−0.037 (0.600−0.030 ) for the ﬁrst (second) octant and m32 = 2.53−0.12 , with
+2.67
δCP = 3.14−1.35 . The data show a weak preference for the normal mass hierarchy, disfavoring theinverted mass hierarchy at 74.0% assuming oscillation parameters at the analysis best-ﬁt point. No
strong preference for the θ23 octant is observed.
Figure 29 shows the expected sensitivity to mass hierarchy as a function of livetime for bothreconstruction algorithms and their ﬁducial volumes. The expected improvement in sensitivity with
the new algorithm becomes more apparent as data are 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 2800 days of livetime.
Acknowledgements
We gratefully acknowledge the cooperation of the Kamioka Mining and Smelting Company. The SuperKamiokande experiment has been built and operated from funding by the Japanese Ministry of Education,
34/39

Downloaded from https://academic.oup.com/ptep/article-abstract/2019/5/053F01/5488936 by University of Tokyo Library user on 22 July 2019

Fig. 29. Expected sensitivity to the normal mass hierarchy as a function of livetime assuming sin2 θ23 = 0.5.
Gray and blue lines show the sensitivity of the atmospheric neutrino sample reconstructed with ﬁTQun in
the conventional FV and expanded FV, respectively. Orange lines denote the sensitivity when events arereconstructed with APFit in the conventional FV.

PTEP 2019, 053F01

M. Jiang et al.

Funding
Open Access funding: SCOAP3 .

Appendix A. Systematic uncertainties
The systematic errors uses in the atmospheric neutrino analysis within the expanded FV and with aconstrained θ13 are summarized in Table A1.
Table A1. Systematic error used in this analysis. The second column shows the best-ﬁt value of the systematic
error parameter in % and the third column shows the estimated 1σ error size in %. The ﬁt result is within the
expanded FV and with θ13 constrained to 0.0210 ± 0.0011. The systematic errors in ring counting, particle
identiﬁcation, and multi-ring likelihood selection within the conventional ﬁducial volume (dwall > 200 cm)
and within the new region (50 cm < dwall < 200 cm) have been merged.
Systematic error
Flux normalization(νμ + ν̄μ )/(νe + ν̄e )ν̄e /νeν̄μ /νμUp/down ratioEν < 1 GeV
Eν > 1 GeV2
Eν < 1 GeV
1 < Eν < 10 GeV
Eν > 10 GeV3
Eν < 1 GeV
1 < Eν < 10 GeV
Eν > 10 GeV4
Eν < 1 GeV
1 < Eν < 10 GeV
Eν > 10 GeV5< 400 MeV
1

> 400 MeV

Multi-GeV
Multi-ring sub-GeV
Multi-ring multi-GeV

Horizontal/vertical ratioPC< 400 MeV

e-likeμ-like0-decay μ-likee-likeμ-like0-decay μ-likee-likeμ-likee-likeμ-likee-likeμ-likee-likeμ-like0-decay μ-likee-likeμ-like0-decay μ-like> 400 MeV

Fit value (%)

σ (%)

7.3
9.6
1.3
2.7
6.6
−1.1
−2.7
−0.34
0.36
3.24
9.2
0.079
0.24
0.87
0.63
0.40
1.34
0.55
0.16
0.32
0.16
0.24
0.16
0.16
−0.023−0.023−0.069
−0.32
−0.44
−0.32

25
15
2
3
5
5
5
8
2
6
15
0.1
0.3
1.1
0.8
0.5
1.7
0.7
0.2
0.4
0.2
0.3
0.2
0.2
0.1
0.1
0.3
1.4
1.9
1.4

Continued on next page

35/39

Downloaded from https://academic.oup.com/ptep/article-abstract/2019/5/053F01/5488936 by University of Tokyo Library user on 22 July 2019

Culture, Sports, Science and Technology, the US Department of Energy, and the US 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 ResearchCouncil (NSERC) of Canada, the Scinet and Westgrid consortia of Compute Canada, and the National ScienceCentre of Poland (2015/17/N/ST2/04064, 2015/18/E/ST2/00758).

PTEP 2019, 053F01

M. Jiang et al.
Table A1. Continued from previous page

Systematic error
Multi-GeV
Multi-ring sub-GeV

PC
K/π ratio in ﬂux calculation6
Neutrino path length
Sample-by-sample

FC multi-GeV
PC + stopping UP-μ

Matter effects
Solar activity
MA in Quasi-elastic interaction
Single π production, axial coupling
Single π production, CA5
Single π production, Background scale factor error
CCQE cross section7
CCQE ν̄/ν ratio7
CCQE μ/e ratio7
DIS cross section
DIS model comparisons8
DIS Q2 distribution (high W )9DIS Q2 distribution (low W )9
Coherent π production
NC/CCντ cross section
Single π production, π 0 /π ±Single π production, ν̄i /νi (i = e, μ)10
NC fraction from hadron simulation
π + decay uncertainty sub-GeV 1-ringe-like 0-decay
μ-like 0-decay
e-like 1-decay
μ-like 1-decay
μ-like 2-decay
Final state and secondary interactions11
1-ringMulti-ring
Meson exchange current12
m221 [3]sin2 (θ12 ) [3]sin2 (θ13 ) [3]
FC reduction
PC reduction
FC/PC separation
PC stopping/through-going separation (bottom)
PC stopping/through-going separation (barrel)
PC stopping/through-going separation (top)Non-ν background
Sub-GeV μ-like
Multi-GeV μ-like
Sub-GeV 1-ring 0-decay μ-like
PC

σ (%)

−0.74
−0.53
−0.32
−0.30
−0.64
−0.35
−0.39
−5.3
−2.3
−0.22
2.46
0.06
0.31
−4.79
−7.03
9.33
−3.03
7.37
4.81
−2.74
5.72
4.84
−5.70
−9.21
−9.41
−0.10
9.08
1.75
2.30
−5.07

3.2
2.3
1.4
1.3
2.8
1.5
1.7
10
10
5
5
6.8
7
10
10
10
10
10
10
10
10
10
10
10
100
20
25
40
10
10

−0.050−0.066−0.34−0.025
−0.47

0.6
0.8
4.1
0.3
5.7

−6.83
6.93
1.08
0.0
0.05
0.00
0.27
−0.26
0.000
−1.75
1.13
−29.0
−0.001−0.002−0.001−0.0151010100.0002
1.3
0.11
1.3
1
0.02
6.8
8.5
40
0.02
0.07
0.02
0.49

Continued on next page

36/39

Downloaded from https://academic.oup.com/ptep/article-abstract/2019/5/053F01/5488936 by University of Tokyo Library user on 22 July 2019

Multi-ring multi-GeV

e-likeμ-likee-likeμ-likee-likeμ-likeFit value (%)

PTEP 2019, 053F01

M. Jiang et al.
Table A1. Continued from previous pageσ (%)

e-likeμ-likee-likeμ-likee-likeμ-likee-likeμ-likee-likeμ-like−0.000−0.001
−0.40
−0.44
0.43
0.19
1.03
1.32
−0.74
−0.42
0.42
−0.66
0.59

0.03
0.07
2.00
−0.86
0.84
0.37
2
2.57
−1.44
−0.82
0.82
−1.28
1.14

e-likeμ-likee-likeμ-like

−0.11
0.10
−0.02
0.02

0.52
−0.49
0.09
−0.09

e-likeμ-likee-likeμ-likeνe ,ν̄e

0.29
−0.12
−0.27
0.16
−0.42
0.32
−0.24
−0.53
0.048
0.029
0.018
0.30
−0.33
3.26

−1.33
0.53
1.25
−0.73
−0.71
0.52
2.02
0.67
0.5
0.3
0.6
1.7
1.5
3

−4.9
2.3
−6.55
1.89
−0.039
0.29
89.2

11
17
24
2.57
1.03
0.7
67

Sub-GeV e-like (ﬂasher event)Multi-GeV e-like (ﬂasher event)
Fiducial volume
Ring separation< 400 MeV
> 400 MeV
Multi-GeV
Multi-ring sub-GeV
Multi-ring multi-GeV

Particle identiﬁcation (1 ring)
Sub-GeV
Multi-GeV
Particle identiﬁcation (multi-ring)
Sub-GeV
Multi-GeV
Multi-ring likelihood selection

Multi-ring e-like
Multi-ring other

Energy calibration
Up/down asymmetry energy calibration
UP-μ reduction
Stopping
Through-going
UP-μ stopping/through-going separation
Energy cut for stopping UP-μ
Path length cut for through-going UP-μThrough-going UP-μ showering separation
Background subtraction for UP-μStopping13
Non-showering13
Showering13
νe /ν̄e separation
Sub-GeV 2-ring π 0
Decay-e tagging
FC stopping muon subtraction
1

Uncertainty decreases linearly with log Eν from 25% (0.1 GeV) to 7% (1 GeV).
Uncertainty is 7% up to 10 GeV, increases linearly with log Eν from 7% (10 GeV) to 12% (100 GeV) and then to 20%(1 TeV).
3
Uncertainty increases linearly with log Eν from 5% (30 GeV) to 30% (1 TeV).
4
Uncertainty increases linearly with log Eν from 8% (100 GeV) to 20% (1 TeV).
5
Uncertainty increases linearly with log Eν from 6% (50 GeV) to 40% (1 TeV).
6
Uncertainty increases linearly from 5% to 20% between 100 GeV and 1 TeV.
7
Difference from the Nieves [35] model is set to 1.0.
8
Difference from CKMT [36] parametrization is set to 1.0.
9
Difference from GRV98 [21] is set to 1.0.
10
Difference from the Hernandez [37] model is set to 1.0.
11
See Ref. [38].
12
Difference from NEUT without the model from Ref. [35] is set to 1.0.
13
The uncertainties in background subtraction for upward-going muons are only for the most horizontal bin, −0.1 <cos θ < 0.
2

37/39

Downloaded from https://academic.oup.com/ptep/article-abstract/2019/5/053F01/5488936 by University of Tokyo Library user on 22 July 2019

Fit value (%)

Systematic error

PTEP 2019, 053F01

M. Jiang et al.

References

38/39

Downloaded from https://academic.oup.com/ptep/article-abstract/2019/5/053F01/5488936 by University of Tokyo Library user on 22 July 2019[1] B. Pontecorvo, Sov. Phys. JETP 26, 984 (1968).
[2] Z. Maki, M. Nakagawa, and S. Sakata, Prog. Theor. Phys. 28, 870 (1962).
[3] C. Patrignani et al. [Particle Data Group], Chin. Phys. C 40, 100001 (2016).
[4] K. Abe et al. [T2K Collaboration], Phys. Rev. Lett. 118, 151801 (2017) [arXiv:1701.00432 [hep-ex]][Search INSPIRE].
[5] P. Adamson et al. [NOvA Collaboration], Phys. Rev. Lett. 118, 231801 (2017) [arXiv:1703.03328[hep-ex]] [Search INSPIRE].
[6] K. Abe et al. [Super-Kamiokande Collaboration], Phys. Rev. D 97, 072001 (2018) [arXiv:1710.09126[hep-ex]] [Search INSPIRE].
[7] R. Acciarri et al., arXiv:1601.05471 [physics.ins-det] [Search INSPIRE].
[8] F. An et al., J. Phys. G: Nucl. Part. Phys. 43, 030401 (2016) [arXiv:1507.05613 [physics.ins-det]][Search INSPIRE].
[9] K. Abe et al. [Hyper-Kamiokande Proto-Collaboration], arXiv:1805.04163 [physics.ins-det] [Search
INSPIRE].
[10] L. Wolfenstein, Phys. Rev. D 17, 2369 (1978).
[11] S. P. Mikheev and A. Yu. Smirnov, Sov. J. Nucl. Phys. 42, 913 (1985).
[12] C. Giunti, C.W. Kim, and M. Monteno, Nucl. Phys. B 521, 3 (1998) [arXiv:hep-ph/9709439] [Search
INSPIRE].
[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. A 610, 710 (2009) [arXiv:0911.0986 [physics.ins-det]][Search INSPIRE].
[14] S. Fukuda et al., Nucl. Instrum. Meth. A 501, 418 (2003).
[15] K. Abe et al., Nucl. Instrum. Meth. A 737, 253 (2014) [arXiv:1307.0162 [physics.ins-det]] [Search
INSPIRE].
[16] M. Honda, T. Kajita, K. Kasahara, and S. Midorikawa, Phys. Rev. D 83, 123001 (2011)[arXiv:1102.2688 [astro-ph.HE]] [Search INSPIRE].
[17] Y. Hayato, Nucl. Phys. Proc. B. Suppl. 112, 171 (2002).
[18] R. Brun, F. Bruyant, F. Carminati, S. Giani, M. Maire, A. McPherson, G. Patrick, and L. Urban, GEANT
Detector Description and Simulation Tool. CERN-W5013, CERN-W-5013, W5013, W-5013 (1994).
[19] R. Gran, J. Nieves, F. Sanchez, and M. J. Vicente Vacas, Phys. Rev. D 88, 113007 (2013)[arXiv:1307.8105 [hep-ph]] [Search INSPIRE].
[20] J. Nieves, I. Ruiz Simo, and M. J. Vicente Vacas, Phys. Rev. C 83, 045501 (2011) [arXiv:1102.2777[hep-ph]] [Search INSPIRE].
[21] M. Glück, E. Reya, and A. Vogt, Eur. Phys. J. C 5, 461 (1998) [arXiv:hep-ph/9806404] [Search INSPIRE].
[22] A. Bodek and U.-k. Yang, AIP Conf. Proc. 792, 257 (2005) [arXiv:hep-ph/0508007] [Search INSPIRE].
[23] Y. Ashie et al. [Super-Kamiokande Collaboration], Phys. Rev. D 71, 112005 (2005)[arXiv:hep-ex/0501064] [Search INSPIRE].
[24] K. Abe et al. [T2K Collaboration], Phys. Rev. Lett. 112, 061802 (2014) [arXiv:1311.4750 [hep-ex]][Search INSPIRE].
[25] M. Shiozawa, Nucl. Instrum. Meth. A 433, 240 (1999).
[26] R. B. Patterson, E. M. Laird, Y. Liu, P. D. Meyers, I. Stancu, and H. A. Tanaka, Nucl. Instrum. Meth. A
608, 206 (2009) [arXiv:0902.2222 [hep-ex]] [Search INSPIRE].
[27] K. Abe et al. [T2K Collaboration], Phys. Rev. Lett. 121, 171802 (2018) [arXiv:1807.07891 [hep-ex]][Search INSPIRE].
[28] F. James, MINUIT Function Minimization and Error Analysis: Reference Manual Version 94.1,
CERN-D-506, CERN-D506 (1994).
[29] R. Wendell et al. [The Super-Kamiokande Collaboration], Phys. Rev. D 81, 092004 (2010)[arXiv:1002.3471 [hep-ex]] [Search INSPIRE].
[30] G. L. Fogli, E. Lisi, A. Marrone, D. Montanino, and A. Palazzo, Phys. Rev. D 66, 053010 (2002)[arXiv:hep-ph/0206162] [Search INSPIRE].
[31] P. Adamson et al. [NOvA Collaboration], Phys. Rev. Lett. 118, 151802 (2017) [arXiv:1701.05891[hep-ex]] [Search INSPIRE].
[32] M. G. Aartsen et al., Nucl. Phys. B 908, 161 (2016).
[33] P. Adamson et al., Phys. Rev. Lett. 112, 191801 (2014) [arXiv:1403.0867 [hep-ex]] [Search INSPIRE].

PTEP 2019, 053F01

M. Jiang et al.

39/39

Downloaded from https://academic.oup.com/ptep/article-abstract/2019/5/053F01/5488936 by University of Tokyo Library user on 22 July 2019[34] A. L. Read, J. Phys. G: Nucl. Part. Phys. 28, 2693 (2002).
[35] J. Nieves, J. E. Amaro, and M. Valverde, Phys. Rev. C 70, 055503 (2004); 72, 019902 (2005) [erratum][arXiv:nucl-th/0408005] [Search INSPIRE].
[36] A. Capella, A. Kaidalov, C. Merino, and J. Tran Thanh Van, Phys. Lett. B 337, 358 (1994)[arXiv:hep-ph/9405338] [Search INSPIRE].
[37] E. Hernández, J. Nieves, and M. Valverde, Phys. Rev. D 76, 033005 (2007) [arXiv:hep-ph/0701149][Search INSPIRE].
[38] E. Pinzon, Measurement of pion-carbon cross sections at DUET and measurement of neutrino oscillation
parameters at the T2K experiment, Ph.D. Thesis, York University Toronto, Canada (2017).