学位論文

Search for Proton Decay Using an Improved Event
Reconstruction Algorithm in Super-Kamiokande（スーパーカミオカンデにおける改良型
事象再構成アルゴリズムを用いた陽子崩壊探索）

平成２９年５月博士（理学）申請

東京大学大学院理学系研究科
物理学専攻
須田

祐介

Abstract
Most of grand unified theories predict proton decay as a result of baryon number
violation, which has not been found so far. The proton decay mode of p → e+ π 0
is predicted as the most dominant decay mode in many non-super-symmetric
grand unified theories. In this thesis, we discuss the result of a search for
p → e+ π 0 in Super-Kamiokande. We collect 483 live days data from 4973
live days to 5456 live days since the latest published paper [1]. Almost halfof all data is collected by a new electronics which enables us to reduce thebackground rate by approximately half using a neutron tagging technique. A
newly developed event reconstruction algorithm is applied to the data with the
new electronics. This new algorithm is able to extend the fiducial volume by
10% from 22.5 kton to 24.7 kton and reduce the background rate by about
30% compared to the conventional analysis while keeping the same level of thesignal efficiency. We analyze the 179.5 kton·years data (from September 2008 to
August 2016) by using the new algorithm and find no candidate for this decay
mode. By combining with zero-observation results of a conventional algorithm
in the 172.5 kton·years data (from the beginning of the experiment to August2008) [1], the lifetime limit for p → e+ π 0 is set to be τ /B > 1.88 × 1034 years at
90 % confidence level. This limit is 18% longer than that of the latest publishedpaper and the current most stringent constraint on non-super-symmetric grandunified theories in the world.

Acknowledgments
I am deeply grateful to my supervisor, Masashi Yokoyama for his enormous
support and patience. He gave me a lot of appropriate advice every important
point. I strongly remember the encouragement he gave me when I was about
to have a breakdown after I opened the data. My thesis would never completed
without his guidance.
I would like to thank all of previous and current Super-Kamiokande collaborators, especially the following people: Masayuki Nakahata, Masato Shiozawa,
Shigetaka Moriyama, Yoshinari Hayato, Roger Wendell, Makoto Miura, Jun
Kameda, Shoei Nakayama, Hidekazu Tanaka, Ed Kearns, Chris Walter andShunichi Mine. Shiozawa-san’s insightful comments were always helpful for myanalysis. I learned a lot of things related to our proton decay analysis from
Miura-san. I was trained by Mine-san the presentation skill.
I also thank Kimihiro Okumura, Yasuhiro Nishimura and Euan Richard forhelping my works on the FC reduction during the expert period.
A crucial element of this thesis is the new event reconstruction algorithm,
fiTQun. I would like to express my gratitude to the members of the fiTQun
working group: Sophie Berkman, Patrick de Perio, Miao Jiang, Akira Konaka,
Andy Missert, Hirohisa Tanaka, Shimpei Tobayama, Mike Wilking and Eric
Zimmerman. My special thanks go to Tobayama-san who introduced me fiTQun
and gave me many important suggestions. His passionate attitude to research
was highly impressive. I appreciate the support from Mine-san, Tobayama-san
and Jiang-san on many fiTQun studies, especially on its energy scale study.
I spend most of time during my master and doctoral programs in Kamioka.
I thank the following people made my Kamioka life happy: Tomonobu Tomura, Hiroyoshi Uchida, Kunxian Huang, Seiko Hirota, Shota Takahashi, Keishi
Hosokawa, Yuki Nakano, Osamu Takachio, Hiroki Takiya, Naoya Oka, Katsuma Nakagawa, Yuto Haga, Isao Kametani, Keiji Tateishi, Yosuke Ohnishi,
Masatoshi Kobayashi, Yuji Okajima, Asato Orii, Ryosuke Akutsu, Daisuke
Fukuda, Fumimasa Muto, Yutaro Sonoda and all of the Kamioka Observatory
member. Talking with Tomura-san and my friends was always fun and essential
to survive in Kamioka. Life in Kamioka made me stronger as a person.
I also thank my college friends: Shion Chen, Chikuma Kato, Chihiro Kozakai, Tatsuya Mori, Miki Nishimura, Takanao Tsuyuki and Harumichi Yokoyama.
I would like to show my appreciation to the members in my lab: Hiroaki
Aihara, Yoshiyuki Onuki, Denis Epifanov, Ian Watson, Konosuke Iwamoto, Sogoi

Figure 1: Word by Yoji Totsuka. It is pronounced as “Yohshi-yo Hohkai-seyo”
which literally means “Decay, protons”.

Mineo, Clement Ng, Yasutaka Kanazawa, Nobuhiro Shimizu, Taichiro Koga,
Junya Sasaki, Yifan Jin, Naruhiro Chikuma, Fuminao Hosomi, Hiroko Niikura,
Tianmeng Lou, Alejandro Mora Grimaldo Johnny, Wan Kun, Tatsuki Ose, Riku
Tamura, Ryo Fujita and Kuniko Kono.
My deep appreciation goes to my family and friends Kohei Kato, Takeru
Kageyama, Kaito Sasaki and Kanako Katayama.
Finally, I want to thank the word written by Yoji Totsuka on February 17,
2003, as shown in Figure 1. He was the first spokesperson of Super-Kamiokande
and hoped that proton decay is observed. This word strongly motivated me tosearch for proton decay and always gave me a hope. I hope that protons hear
our wish someday in the future.

ii

Contents
1 Introduction

12 Physics motivation
2.1 Standard Model and its Problems . . . . . . . . . . . . . . . . . .
2.2 Grand Unified Theories . . . . . . . . . . . . . . . . . . . . . . .
2.2.1 SU (5) Model . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.2 SUSY SU (5) Model . . . . . . . . . . . . . . . . . . . . .
2.2.3 SO(10) Model . . . . . . . . . . . . . . . . . . . . . . . .
2.2.4 Proton Decay Search . . . . . . . . . . . . . . . . . . . . .

4
4
5
5
6
7
8

3 Super-Kamiokande Experiment10
3.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
3.2 Principle of the Detector . . . . . . . . . . . . . . . . . . . . . . . 11
3.3 Water and Air . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.4 Photosensors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.5 Data Acquisition . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
3.5.1 System for SK-I to SK-III . . . . . . . . . . . . . . . . . . 17
3.5.2 System for SK-IV . . . . . . . . . . . . . . . . . . . . . . 19
3.6 Detector Calibration . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.6.1 PMT Gain . . . . . . . . . . . . . . . . . . . . . . . . . . 20
3.6.2 PMT Timing . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.6.3 Water Transparency . . . . . . . . . . . . . . . . . . . . . 24
4 Monte Carlo Simulations28
4.1 Signal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
4.2 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
4.2.1 Neutrino Flux . . . . . . . . . . . . . . . . . . . . . . . . . 31
4.2.2 Neutrino Interactions . . . . . . . . . . . . . . . . . . . . 32
4.3 Detector . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
5 Data Set
37
5.1 First Reduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
5.2 Second Reduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
5.3 Third Reduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

iii

CONTENTS

5.4
5.5
5.6

Fourth Reduction . . . . . . . . . . . . . . . . . . . . . . . . . . .
Fifth Reduction . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Final FC Data Set . . . . . . . . . . . . . . . . . . . . . . . . . .

41
41
45

6 Event Reconstruction
47
6.1 Overview of Two Algorithms . . . . . . . . . . . . . . . . . . . . 47
6.2 APfit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
6.2.1 Vertex Fitting . . . . . . . . . . . . . . . . . . . . . . . . 48
6.2.2 Ring Counting . . . . . . . . . . . . . . . . . . . . . . . . 49
6.2.3 Particle Identification . . . . . . . . . . . . . . . . . . . . 51
6.2.4 Momentum Determination . . . . . . . . . . . . . . . . . . 53
6.2.5 Michel Electron Tagging . . . . . . . . . . . . . . . . . . . 55
6.3 fiTQun . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
6.3.1 Likelihood Function . . . . . . . . . . . . . . . . . . . . . 55
6.3.2 Vertex Pre-fitter . . . . . . . . . . . . . . . . . . . . . . . 69
6.3.3 Hit Clustering . . . . . . . . . . . . . . . . . . . . . . . . 69
6.3.4 Single-Ring Fitter . . . . . . . . . . . . . . . . . . . . . . 70
6.3.5 Multi-Ring Fitter . . . . . . . . . . . . . . . . . . . . . . . 74
6.4 Performance in MC Simulation: APfit vs. fiTQun . . . . . . . . 77
6.5 Fiducial Volume . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
6.6 Energy Scale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
6.6.1 High Energy Stopping Muon Sample . . . . . . . . . . . . 84
6.6.2 Low Energy Stopping Muon Sample . . . . . . . . . . . . 86
6.6.3 Michel Electron Energy Spectrum . . . . . . . . . . . . . 88
6.6.4 Neutral Pion Mass . . . . . . . . . . . . . . . . . . . . . . 88
6.6.5 Energy Scale Error . . . . . . . . . . . . . . . . . . . . . . 90
6.6.6 Uniformity of Energy Scale Error . . . . . . . . . . . . . . 91
7 Data Analysis
94
7.1 Selection Criteria . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
7.2 Systematic Errors . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
7.3 MC-Data Agreement in the Outside the Signal Region . . . . . . 118
7.4 Result . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121
7.5 Lifetime Limit . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
8 Discussion
129
8.1 Comparison with GUT models . . . . . . . . . . . . . . . . . . . 129
8.2 Further Improvements for fiTQun . . . . . . . . . . . . . . . . . . 130
8.3 Discovery Potential in Future . . . . . . . . . . . . . . . . . . . . 130
9 Conclusion

138

A Special Handling for the Blasting in the Mine140

B Time Variation of PMT Hit Rate

143

C Development of FiTQun for SK-I to SK-III

146

iv

CONTENTS

D Mass Production of Atmospheric Neutrino MCv

149

Chapter 1

Introduction
Where do we come from? What are we? Where are we going? – Paul Gauguin.
These might be philosophical questions, but we can paraphrase them in a physics
point of view; How does our universe begin and end? What are fundamental
particles and interactions between them?
As for the first question, it is considered that our universe at birth was atiny space with an extremely high temperature and has been expanding at anaccelerating rate. The end of the universe is not clear, but dark energy and dark
matters which are unknown thing may play a crucial role. Such cosmologicaldiscussions rely on the general relativity proposed by A. Einstein in 1916. In
its 100th anniversary year, the LIGO and Virgo collaborations announced thatthey observed gravitational waves from a binary black hole merger [2]. By thisobservation, the last remaining prediction of the general relativity is confirmed.
About the second question, we have a beautiful theory of elementary particles, so-called the standard model (SM). The SM contains three fundamental
forces except for gravity; strong, weak and electromagnetic forces. Weak andelectromagnetic forces are unified at a high energy. Most of elementary particle
phenomena can be explained by the SM very precisely. The last missing pieceof the SM prediction was a Higgs boson, but it was recently discovered by the
ATLAS [3] and the CMS [4] experiments in 2012.
However, there are still a lot of things which imply the existence of physics
beyond the SM; the unification of three forces, neutrino masses, electric charge
quantization and so on. Neutrinos are the one of the fundamental particle in
the SM which assumes neutrinos are massless. This assumption is completely
excluded by the observation of neutrino oscillations at the Super-Kamiokande
experiment in 1998 [5]. The theory beyond the SM must exist. The most attractive model is a grand unified theory (GUT) which merges three fundamental
forces at very high energy. Most of GUTs predict that a nucleon spontaneously
decays into lighter particles, “nucleon decay” as shorthand, while protons are
considered as stable in the SM. Since the temperature of the early universe isextremely high, the verification of GUT is one of the important issues to answer
the questions above. Moreover, a proton’s fate might be related to the end of1

CHAPTER 1. INTRODUCTION

universe. Therefore, we search for proton decay particular in p → e+ π 0 whichis the most dominant decay mode in many GUTs.
Since protons live a long time, we need to monitor a large number of protons.
Our experimental equipment is the Super-Kamiokande (SK) detector which is
the largest volume detector for nucleon decay search. The most sensitive proton
decay mode in the SK detector is p → e+ π 0 . Figure 1.1 shows a schematic
diagram of p → e+ π 0 and its typical Monte Carlo simulation event display in SK.
There are four distinct observation periods; SK-I, -II, -III and -IV. The current
period is SK-IV which has the most stable operation and the best detector
performance. As of today, the amount of data taken in the SK-IV period is
almost half of the whole period. We have been searching for many nucleon
decay modes not only the p → e+ π 0 mode since the experiment started in 1996.
No significant data excess was found. However, the nucleon lifetime predictedby many GUTs is still in the SK search region. It can be said, nucleon decay
could happen at anytime.
We have to improve our analysis method, not just waiting to observe nucleon
decay. We developed a new event reconstruction algorithm, called fiTQun.
FiTQun has better performance than a conventional algorithm, APfit, which
has been used since the experiment started. At present, fiTQun works only for
the SK-IV data. So, we search for p → e+ π 0 with the data of SK-I to SK-III
by APfit and that of SK-IV by fiTQun in this thesis.
Based on the huge efforts of the SK collaboration over the past 20 years,
this thesis exists. The author highly contributed to the followings related to thecontents of this thesis:• Quality control of a data set for proton decay and atmospheric neutrino
analyses as an expert (1 year)• Validation of fiTQun with some collaborators through a proper calibration• Evaluation of a fiTQun’s energy scale error• Generation of a 3000 years atmospheric neutrino Monte Carlo simulation
as the background of p → e+ π 0• First application of fiTQun to the proton decay analysis (p → e+ π 0 )through detailed checks and a tuning of fiTQun• First proton decay search (p → e+ π 0 ) including an unexplored region of
the SK detector (10% expansion of a conventional fiducial volume)
In Chapter 2, we discuss the physics motivation of this thesis, including abrief history and the importance of proton decay search. Our experimental set
up and its calibration are explained in Chapter 3. Chapter 4 describes Monte
Carlo simulations for both the signal and background events. The data set
for this analysis is prepared through a proper data reduction process describedin Chapter 5. Chapter 6 explains the event reconstruction algorithms. The
analysis and result are shown in Chapter 7. Chapter 8 discusses about theresult and the future prospects. Finally, we conclude this thesis in Chapter 9.
2

CHAPTER 1. INTRODUCTION

Figure 1.1: Schematic diagram of p → e+ π 0 (bottom left) and its typicalMonte Carlo simulation event display in SK-IV (center). A free proton decayedat around the tank center. It emitted a positron on the left and a neutral pionon the right. As a result of the pion decay, three electromagnetic shower rings
are observed. Each dot corresponds to a PMT and the color shows the amountof the observed charge.

3

Chapter 2

Physics motivation
2.1

Standard Model and its Problems

The Standard Model (SM) is a quantum field theory with the gauge group
SU (3)C × SU (2)L × U (1)Y , three generations of fermions and a scalar Higgs
boson. SU (3)C and SU (2)L ×U (1)Y describe the strong and electro-weak forces,
respectively, which are carried by the gauge bosons G1∼8 , W1∼3 and B. The
SM fermions are listed in Table 2.1 with their representations of each groupand hypercharge Y = 2(Q − T3 ), where Q is the electric charge and T3 is theSU (2) isospin. The Higgs boson makes some fermions and gauge bosons massive
through the spontaneous symmetry breaking of SU (2)L × U (1)Y .
Most of high energy phenomena can be very precisely explained by the SM. In
summer 2012, the Higgs boson as the last missing piece of the SM was discoveredby the ATLAS [3] and the CMS [4] experiments. The SM was completed,however, it is obviously not an ultimate theory because there are no gravity,
dark matter and neutrino masses. Furthermore, there are many theoreticalquestions such as the hierarchy problem, triangle anomalies, the number of free

Table 2.1: Fermion fields in the standard model and their representations of
non-Abelian groups and hypercharge. Bold numbers indicate the dimension ofrepresentations.

u, d
uc
dc
νe , eec

SU (3)CSU (2)LU (1)Y

3
3̄
3̄
1
1

2
1
1
2
1

1
3
− 43
2
3

4

−1
2

CHAPTER 2. PHYSICS MOTIVATION

parameters and so on. In particular, the electric charge quantization of fermionshighly motivates us to unify quarks and leptons. Since the electric charges are
the arbitral coupling constants of U (1)Y , there is no reason that the electriccharge of down quarks is the third of that of electrons. But, the magnitudeof electric charge of protons is measured as the same of that of electrons withan accuracy of 1 × 10−21 [6]. In addition to that, it can be said that thebaryon/lepton number symmetry in the SM is just an accidental symmetry(not a local gauge symmetry), which might be violated in a high energy limit.
Also, the three running coupling constants in the SM run across each other ina very high energy region > 1014 GeV. These circumstantial evidences imply
that there is a larger gauge group which is approximated to the SM group in a
low energy limit. Such models are called Grand Unified Theories (GUTs). We
discuss several GUTs in the following.

2.2

Grand Unified Theories

There are many GUT models and they predict nucleon decay in general. In thissection, we focus on SU (5), SUSY SU (5) and SO(10) models. Comprehensive
discussions on GUTs (flipped SU (5), E6 , etc.) can be found in the document [7].

2.2.1

SU (5) Model

The first grand unified theory was proposed by H. Georgi and S. Glashow in1974 [8]. The rank of a gauge group represents the number of independent
quantum numbers; the rank of U (1) is one and that of SU (n) is n − 1. So, asimple gauge group which contains the SM groups as subgroups should have arank more than three. In addition to that, complex representations are necessary
in that group. They found out that SU (5) is the minimal answer.
Each family of fermions and gauge bosons in SU (5) are represented by 10+ 5̄
and 24, respectively, as follows;
 c 


d1
0
uc3 −uc2 u1 d1
−uc3
 dc2 0
uc1
u2 d2 
 c 
 c

c



0
u3 d3 5 :  d3 
and10 :  u2 −u1
(2.1)
,
 e 
−u1 −u2 −u3
0
ec 
−νe
−d1 −d2 −d3 −ec 0
G11 − √2B
30
 G2

1

G31
24 : 

 X1

Y1


G12
2
G2 − √2B
30
G32

G13
G23
G33 − √2B
30

X2

X3

Y2

Y3

X1
X2
X3
3
3B
W
√ + √
2
30
−

W

Y1
Y2
Y3
W+
3
3B
√ + √
−W
2
30





 , (2.2)

where the suffixes of quarks show the SU (3)C color, Gs correspond to gluons,
X and Y are 12 new gauge bosons. Since TrQ = Tr(T3 + Y /2) = 0, this model
5

CHAPTER 2. PHYSICS MOTIVATION

u
X

u
d

e+

u

d

d

d

u

e+
Y

u
u

Figure 2.1: Feynman diagrams for p → e+ π 0 via X and Y bosons.

can explain the electric charge quantization as follow,
5X

Qa = 3Qdc + Qe + Qνe = 0

a=1

→

1
Qd = −Qdc = − .
3

(2.3)

In the gauge sector, the upper left 3 × 3 matrix is for SU (3)C rotations and thelower right 2 × 2 matrix is for SU (2)L rotations whereas the other off-diagonal
elements (Xs and Y s) connect quarks and leptons. Thus, the baryon/lepton
number is no longer a conserved quantity, and nucleons spontaneously decay.
Because we, made of matters, exist stably, both of X and Y should be massive.
Through a spontaneous symmetry breaking of SU (5) → SU (3)C × SU (2)L ×
U (1)Y occurred at a GUT energy scale (> 1014 GeV), they earn masses.
In the minimal SU (5), the most dominant decay mode is p → e+ π 0 whose
predicted proton lifetime is about 1031±1 years [9]. Figure 2.1 shows the Feynman diagrams for p → e+ π 0 . However, this was already excluded by the previous experiments as seen in Section 2.2.4. Furthermore, SU (5) predicts the
Weinberg angle (θW ) as sin2 θW = 3/8, but this prediction was also excluded
by experimental results of about 0.23 [10]. The minimal SU (5) GUT model iscompletely ruled out.

2.2.2

SUSY SU (5) Model

One attractive idea to save the SU (5) model is to introduce a symmetry between
fermions and bosons, so-called “super-symmetry (SUSY).” t t SUSY makes thenumber of both fermions and bosons double, however, three running coupling
constants can completely match at around 2 × 1016 GeV [11–13] as shown in
Figure 2.2.
Since we do not observe any of light leptons with spin-zero, SUSY must be
broken at a low energy, typically > 1 TeV. In order to observe SUSY particles,
direct searches have been conducted by the ATLAS and the CMS experiments,
but no new particles were found so far. The GUT energy scale increases by
SUSY, and the proton lifetime is predicted to be longer than that of non-SUSY
SU (5).
In the minimal SUSY SU (5) model, the most dominant decay mode of protons is p → ν̄K + since SUSY partners of Higgs (Higgsinos) which are mediators
6

CHAPTER 2. PHYSICS MOTIVATION

Figure 2.2: Running coupling constants in the SM (left) and the minimal SUSY
SM (right) powered by SOFTSUSY [16]. Taken from [7].

prefer to couple with heavier particles [14]. Its predicted upper limit of thepartial lifetime of p → ν̄K + is < 2.9 × 1030 years whereas that of p → e+ π 0
is > 4.1 × 1033 years [14]. Looking at Super-Kamiokande results, it seems thatthe minimal SUSY SU (5) model is also excluded as well as the minimal SU (5)model [15].

2.2.3

SO(10) Model

The minimal SU (5) model can not explain why neutrinos have mass. The seesaw
mechanism can solve this question by introducing right-handed neutrino [17–19]. A right-handed neutrino is naturally implemented in the next larger GUT,SU (10) which has rank five [20, 21]. All fermions including the right-handed
neutrino can be embedded in a single 16 representation (= 10 + 5 + 1).
There are many patterns of SO(10) symmetry breaking as follow.
SU (5) × U (1),SU (4) × SU (2)L × SU (2)R ,SO(10) →SU(3) × SU (2)L × SU (2)R × U (1),

(2.4)


SU (4) × SU (2)L × U (1)→ SU (3)C × SU (2)L × U (1)YThe second term SU (4) × SU (2)L × SU (2)R is known as the Pati-Salam modelin which SU (2)L × SU (2)R expresses a left-right symmetry (Parity is recoveredat very high energy). Since there are more than one symmetry breaking energy
scales in SO(10), it is possible to choose them so that three running coupling
constants unify.
The most dominant decay mode of SO(10) is p → e+ π 0 as well as the SU (5)model and its lifetime is predicted in a large range 1032 − 1039 years [22, 23].
7

CHAPTER 2. PHYSICS MOTIVATION

2.2.4

Proton Decay Search

The baryon number conservation law (Weyl [24] (1929), Stueckelberg [25] (1938),
Wigner [26] (1949)) had been tested by various proton decay searches before theminimal SU (5) model was proposed in 1974. The first experimental search for
proton decay was conducted by M. Goldhaber [27] in 1954. He searched forfission products of Th232 induced by proton decay, and set the limit of > 1020
years [27]. This kind of inclusive indirect search was done for several materials
by using geochemical and radiochemical methods which are described in [28];Th232 (Flerov et al [29]), Te130 → Xe129 (Evans and Steinberg [30]), K39 → Ar37
(Fireman [31]) and Mica (Bennett [32]).
At the same time, the direct search by using liquid scintillator was in fashion[27, 33–40]. The first one is Reines et al [27], they measured charged particles
with kinetic energy of above 100 MeV from proton decay by 300 liters of liquid
scintillator in the depth of 200 m.w.e.
After 1974, many GUT motivated experiments began to be proposed. The
first partial lifetime limit for p → e+ π 0 is > 6 × 1029 years measured by Learned
et al using liquid scintillators located in deep underground [40]. Thereafter,
four experiments were conducted by using fine grained iron calorimeter which
consists of alternative iron layers (NUSEX [41], KGF [42], Soudan [43] and Frejus experiment [44]). Although such calorimeter has excellent energy resolution
and particle identification, it is not feasible to scale up the detector volume due
to cost. Compared to that technique, it is much easier to construct large volume
of water Cherenkov detector discussed in Chaper 3 with lower cost. Thus, thecurrent main stream in proton decay search is to use water Cherenkov technique.
The previous water Cherenkov experiments are Homestake (150 ton in fiducial mass) [45], IMB (3.3 kton in fiducial mass) [46] and KAMIOKANDE (1.0kton in fiducial mass) [47]. They searched for proton decay in many decay
modes, but no evidences were found. As a successor of the KAMIOKANDE
experiment, the Super-Kamiokande (SK) experiment has been leading proton
decay search with its fiducial mass of 22.5 kton [48, 49]. Unfortunately, we
found no evidence and thus set the partial lifetime limit to p → e+ π 0 [1] andp → νK + [50] are,
τ /B(p → e+ π 0 ) > 1.6 × 1034 years(2.5)τ /B(p → νK + ) > 5.9 × 1033 years,

(2.6)

and

which ruled out the minimal SU (5) and minimal SUSY SUSY(5) model, respectively.
The summary plot of previous proton decay searches is shown in Figure 2.3.
SK has led this field and given the stringent constraints to GUTs for almost
20 years. This superiority of SK is expected to continue after about 10 years
from today before starting a next generation of proton decay search experimentdiscussed in Chapter 8. Since the sensitivity of SK is within the range of the8Lifetime Limit (years)

CHAPTER 2. PHYSICS MOTIVATION

1037

Minimal SU(5), p→ e + π0
Minimal SO(10), p → e + π0
Goldhaber
Flerov et al
Evans & Steinberg

33

Fireman

10

Bennett
Reines et al (1954)
Reines et al (1958)
Backenstoss et al
Giamati & Reines
Kropp & Reines

29

10

Gurr et al
Bergamasco & Picchi
Reines & Crouch
Learned et alBattistoni et al (NUSEX)

25

10Krishnaswamy et al (KGF)
Bartelt et al (Soudan)
Berger et al (Frejus)
Cherry et al (Homestake)
Hirata et al (Kamiokande)
McGrew et al (IMB)

21

10

Shiozawa et al (Super-K I)
Nishino et al (Super-K I-II)
Abe et al (Super-K I-IV)

1940

1960

1980

2000

2020

2040
Year

Figure 2.3: Summary of proton lifetime limits by the previous experiments.
The black (red) region shows the predicted lifetime range in the minimal SU (5)(SO(10)) model. The red, green, black and blue points correspond to the geochemical and radiochemical method, liquid scintillator, fine grained calorimeter
and water Cherenkov experiments. The results of Learned et al and KGF,
Soudan, Frejus, Homestake, KAMIOKANDE, IMB and Super-Kamiokande experiments are for p → e+ π 0 whereas the others are the other modes or inclusive
results.

minimal SO(10) model, there is still a chance to discover proton decay. However,the contamination of background events goes on increasing year after year. In
this thesis, we search for p → e+ π 0 which is the most dominant decay modein the minimal SO(10) model and the most sensitive decay mode in SK. We
expand the fiducial volume and reduce the background rate by using a new event
reconstruction algorithm described in Chaper 6. This analysis is the cleanest
search for p → e+ π 0 with the highest statistic of data in the world.

9

Chapter 3

Super-Kamiokande
Experiment3.1

Overview

Super-Kamiokande (SK) is a large water Cherenkov detector proposed in 1987,
whose main physics target is to search for proton decay. The SK detector
mainly consists of a cylindrical stainless steel tank (41.4 m in hight and 39.3
m in diameter), 50 kton ultra-pure light water (H2 O) and the photomultipliertubes (PMTs), which is shown in Figure 3.1 [51]. In order to avoid cosmic
ray muon background, the detector is located in the Kamioka mine 1,000 m
underground of Mt. Ikenoyama, in Kamioka town, Gifu Prefecture, Japan. At
that point (2,700 m.w.e.), the cosmic ray muon flux is 6 × 10−8 cm−2 s−1 sr−1
which is five orders of magnitude smaller than that on the surface of the earth,and its rate observed in SK is about 2 kHz.
There are two concentric cylinders inside the tank; the inner detector (ID)and the outer detector (OD). The ID is a main detector which contains 32 ktonof water with 36.2 m in height and 33.8 m in diameter, whereas the OD is acosmic ray muon veto detector which covers the region surrounding the ID withthe thickness of about 2 m. Both detectors are optically separated by opaquematerials. The OD is also used to prevent backgrounds from the surroundingrock such as neutrons and gamma rays.
On the ID surface, 20-inch PMTs are uniformly mounted facing inward viewing the ID volume while 8-inch PMTs are uniformly mounted on the OD surface
facing outward viewing the OD, as shown in Figure 3.2. In order to prevent thegeomagnetic field effect to PMTs, 26 sets of Helmholtz coils are equipped on
the wall of the tank and reduce the magnetic field from 450 mG to 50 mG.
The observation started on April 1996 with 11,146 ID PMTs and 1,885
OD PMTs and stopped on July 2001 for a maintenance work to replace bad
PMTs. This phase is called SK-I. While refilling the tank with the water after
the replacement work, an accident happened on November 12th, 2001 at 11:01

10

CHAPTER 3. SUPER-KAMIOKANDE EXPERIMENT

Figure 3.1: Schematic view of the Super-Kamiokande detector. Taken from [51].

JST. One of the ID PMTs at the bottom collapsed, a vacuum space suddenly
emerged and was immediately crushed, and then shock wave was generated anddestroyed PMTs. As a result of chain reaction of shock wave, about 60% of all
PMTs were broken. In order to avoid same incident, fiber reinforced plastics(FRPs) and acrylic cases come into use to protect ID PMTs hereafter. The next
phase of SK, SK-II, started from October 2002 after redistributing the surviving
ID PMTs and installing new OD PMTs and stopped on October 2005. After
reconstructing the SK, the third phase (SK-III) started from June 2006 with
11,129 ID PMTs. On September 2008, the read-out electronics and the dataacquisition system were upgraded. The SK-IV phase has begun at that timeand continues until this day. The photo-coverage of SK is 40%, except SK-II(19%). Table 3.1 shows the summary of the SK phases.

3.2

Principle of the Detector

Super-Kamiokande observes charged particles by detecting the Cherenkov radiation using PMTs. In a medium with a refractive index n, a charged particle
emits electromagnetic shock wave when the velocity of the particle (v = βc) exceeds the phase velocity c/n of light, as shown in Figure 3.3. This phenomenon is
known as the Cherenkov radiation. Along with the particle track, the Cherenkov
light is emitted as a cone with its opening angle θC which satisfies,cos θC =1.
nβSince n is about 1.34 in water, θC is 42◦ for the particle with β ' 1.
11

(3.1)

CHAPTER 3. SUPER-KAMIOKANDE EXPERIMENT

Figure 3.2: Supporting frames of PMTs. Taken from [51].

Table 3.1: Summary of the four Super-Kamiokande observation phases. The
live time of SK-IV corresponds to the data collecting from Sep. 2008 to Aug.
2016.

Observation
Live time (days)
Number of PMTs
Photo Coverage
FRP&Acrylic Case
Electronics
Trigger

Start
End
ID
OD

ID
OD

SK-I

SK-II

SK-III

SK-IV

Apr. 1996
Jul. 2001
1489.2
11,146
1,885
40%
No
ATM
OD QTC
HardwareOct. 2002
Oct. 2005
798.6
5,182
1,885
19%
Yes
ATM
OD QTC
HardwareJun. 2006
Aug. 2008
518.1
11,129
1,885
40%
Yes
ATM
OD QTC
HardwareSep. 2008(running)
2650.4
11,129
1,885
40%
Yes
QBEE
QBEE
Software12

CHAPTER 3. SUPER-KAMIOKANDE EXPERIMENT

Cherenkov light

Charged particle

Figure 3.3: Schematic view of the Cherenkov radiation.

The momentum threshold for the Cherenkov radiation is derived from thecondition of β > 1/n. For example, the Cherenkov threshold for electrons inwater is 0.57 MeV/c. Table 3.2 summarizes the Cherenkov threshold for various
charged particles in water.
The number of emitted photons per wavelength λ per unit travel distance xof a charged particle is given as,


d2 N
2πα
1
= 2 1− 2 2 ,
(3.2)
dxdλ
λ
n β
where α is the fine structure constant. For the wavelength between 300 nm and
600 nm, 340 photons are emitted per unit cm by the particle with β ' 1 inwater (n = 1.34).
Figure 3.4 shows the Cherenkov light patterns for two different types of theparticle observed in SK. The difference between those two is the sharpness of
the Cherenkov ring edge. Electrons or high energy gamma-rays create fuzzy
rings (“shower”) as shown in the top display since they produce electromagnetic shower and are deflected by the multiple scattering, whereas muons or
charged pions create sharper rings (“non-shower”) as shown in the bottom display since they do not produce electromagnetic shower. Not only the particle
type, the kinematics of particle such as the vertex, direction and momentum aredetermined by using the Cherenkov light arrival time and its amount received
at each PMT as described in Chapter 6.

13

CHAPTER 3. SUPER-KAMIOKANDE EXPERIMENT

Figure 3.4: Shower (top) and non-shower (bottom) events in the fully-contained
data set discussed in Chap 5. Each dot shows a PMT which detects photons,
and its color indicates the amount of the charge. According to a reconstruction
algorithm (fiTQun discussed in Chap 6), the momentum of the particle on the
top is 468 MeV/c assuming a single electron ring and that on the bottom is
506 MeV/c assuming a single muon ring.
14

CHAPTER 3. SUPER-KAMIOKANDE EXPERIMENT

Table 3.2: Momentum thresholds of the Cherenkov radiation in water for various
charged particles assuming n = 1.34.

Momentum (MeV/c)

3.3

e

µ

π+

p

0.57

118

156

1052

Water and Air

Since Cherenkov light travels in a large distance before reaching the PMTs, ahigh transparency of the water is necessary. As for the water in SK, spring waterin the mine is used. The water is continuously purified and circulated at a rateof about 30 tons/hour by a complex system [51]. The water purification system
mainly consists of several filters, a UV sterilizer and two degasifiers in orderto remove particles larger than 0.2 µm, bacteria and radon dissolved in water,
respectively. The water temperature is controlled at around 13 ◦ C by heatexchangers in order to reduce PMT dark noise and suppress bacteria growth.
In the mine, the air is radon-rich due to the leaking from the rock and canbe dissolved in the water since there is a gap between the top of the SK tankand the water surface. Radiation from radon decay becomes the low energy
background. In order to reduce such background, the rock of the experimental
hall was coated with a polyurethane material, and radon-free air is continuously
supplied from outside the mine into the experimental area.

3.4

Photosensors

For the Cherenkov light detection in the ID, the 20-inch PMTs manufactured byHamamatsu Photonics K.K. (R3600) are used [51]. Figure 3.5 shows a picture of
the 20-inch PMT. The bialkali (Sb-K-Cs) photocathode has a quantum efficiencyof 22% at a peak where wavelength is 360 nm, as shown in Figure 3.6. The
dynode structure is a 11-stage Venetian blind type. The operation of the ID
PMTs is conducted applying high voltage ranging in 1700-2000V between thephotocathode and the last dynode. The collection efficiency of the first dynode,gain, transit time spread and dark noise rate are over 70%, 107 , 2.2 ns (1σ) fora single photoelectron and 3 kHz at 0.25 photoelectrons threshold, respectively.
Although the acrylic covers are attached to the ID PMTs used in SK-II to SKIV, the effect on the photon detection is small because the transparency of thecover for a normal incident photon in water is over 96% at 350 nm.
The Hamamatsu 8-inch PMTs (R1408) are used for the OD [52]. In order
to increase the photon collection efficiency, a 1.3 cm thickness 60 cm square
wavelength-shifting plate (ultra-violet light → blue-green light) is equipped toeach OD PMT as shown in Figure 3.7. Due to the reemission process of theplate, the timing resolution of the OD PMT gets worse from 13 ns to 15 ns

15CHAPTER 3. SUPER-KAMIOKANDE EXPERIMENT

Figure 3.5: Photos of the 20-inch PMT (top) and the FRP and acrylic case(bottom).
16

CHAPTER 3. SUPER-KAMIOKANDE EXPERIMENT

Figure 3.6: Quantum efficiency as a function of wavelength for the 20-inch PMT.
Taken from [51].

in FWHM. However, the collection efficiency of the OD PMT is increased by
60%. Furthermore, the reflective material (Tyvek) sheets which have a 90%
reflectivity at 400 nm are put on the tank wall and the spaces between the OD
PMTs in order to enhance the photon collection efficiency. From SK-III, theTyvek sheets are also used to separate the barrel and the top and bottom region
in the OD for better rejection of “corner-clipping” cosmic-ray muon events.

3.5

Data Acquisition

A PMT signal is fed into front-end electronics via a 70 m RG58 coaxial cable, and
it is considered as “hit” when the pulse height exceeds a 0.25 photoelectrons
level. An event trigger is issued if the number of the hit PMTs surpasses a
trigger threshold. The event consists of the time and charge of the hit PMTs
within a given time window around the trigger timing. There are several triggers
depending on the visible energy; the low energy (LE), high energy (HE), outerdetector (OD) trigger and so on. In this analysis, we use the HE triggered
events. Since the data acquisition (DAQ) system which was used in SK-I to
SK-III was upgraded for SK-IV, the following sentences discuss the two different
DAQ systems, separately.

3.5.1

System for SK-I to SK-III

The front-end electronics of the ID PMTs for SK-I to SK-III is the analog
timing module (ATM) [53]. One ATM has 12 channels of the charge-to-analog
17

CHAPTER 3. SUPER-KAMIOKANDE EXPERIMENT

Figure 3.7: Schematic drawing of the OD PMT with the wavelength-shifting
plate. Taken from [52].

converters (QACs) and time-to-analog converters (TACs), and 960 ATMs wereused for the ID PMTs in total.
For each PMT, a 200 ns rectangular pulse with a 15 mV pulse height isgenerated when the signal is a hit. The ATM outputs the analog sum of all
the rectangular pulse (HITSUM). Then, the event trigger is issued if the sum of
HITSUMs from all ATMs exceeds a trigger threshold. Once the event trigger
happened, each ATM starts to digitize and record the time whose the hit PMT
signal surpasses the threshold and the charge integrated in a 400 ns time window
around the PMT trigger timing. In order to avoid missing successive events suchas Michel electrons, two switching pairs of both the QAC and TAC are preparedfor each channel because the process time for one channel is 5.5 µs. The time
window length of events is 1.3 µs around the event trigger [−400 ns, +900 ns].
The signal arrived from 400 ns to 900 ns after the hit is neglected due to a cablereflection pulse caused by a impedance mismatch between an ATM and a PMT.
Similar to the ATM system, signals from OD PMTs are processed by chargeto-time converter modules. When an event trigger happened, the time andcharge information are converted to digital by multi-hit time-to-digital converter
modules (LeCroy 1877).
There are three types of the event trigger in this system; the super low energy(SLE), LE, HE and OD triggers. Table 3.3 shows the trigger thresholds of SLE,LE and HE for each SK period. The SLE threshold corresponds to 4.6 MeV

18

CHAPTER 3. SUPER-KAMIOKANDE EXPERIMENT

Table 3.3: Summary of the hardware triggers used in SK-I to SK-III.

SLE
LE
HE

SK-I

SK-II

SK-III

186 mV
320 mV
340 mV

110 mV152 mV
180 mV

186 mV
302 mV
320 mV

electrons. The OD threshold corresponds to 19 hits.

3.5.2

System for SK-IV

Since the ATM has an electronics dead time and a cable reflection issue, anew electronics QBEE [54] was developed and installed for both the ID and
OD PMTs. One QBEE has 24 channels, and each channel consists of a custom charge-to-time converter (QTC) and a multi-hit time-to-digital converter(TDC). The timing and charge resolutions of the QTC are 0.3 ns for 2 pC andabout 0.2 pC for below 50 pC, respectively. The charge dynamic range of the
QTC is 0.2 - 2500 pC, and the charge non-linearity is better than 1% for theoverall range.
The QTC records all the hits and immediately integrates the charge fromthe signals in a 400 ns time window for each hit, and TDC converts the QTC
output into the digitized time and charge. The digitized time and charge arerecorded when a software trigger is happened. The software trigger is issued
when the number of hit PMTs in a 200 ns sliding time window (N200 ) exceeds
a threshold. In this system, the SLE, LE and HE triggers are defined as well as
the previous system and these event time windows become larger. In addition
to these triggers, special high energy (SHE) and “after trigger (AFT)” triggers
were introduced for a neutron tagging algorithm which is used in this analysis.
The AFT trigger is issued only for an event triggered by a SHE trigger without
an OD trigger and save an additional 500 µs of data. The software triggers
employed in this system are summarized in Table 3.4.

3.6

Detector Calibration

It is necessary for any physics analysis to keep the charge and time responseof each PMT are the same for the same intensity and timing of the incident
light. Thus, the detector calibration is important. Moreover, calibration results
discussed below are used as inputs to the detector simulation and the event
reconstruction algorithms.
In this section, we focus on measurements in the ID. The PMT gain calibration is discussed at first since the hit timing depends on the amount of thehit charge. The relative timing calibration and the water transparency mea-

19

CHAPTER 3. SUPER-KAMIOKANDE EXPERIMENT

Table 3.4: Summary of the software triggers used in SK-IV.
Trigger Type

N200 Threshold

Time Window (µs)

SLE
LE
HE
SHE
AFT
OD

34 → 3147
50
70 → 58
SHE without OD
22−0.5 ∼ +1.0
−5 ∼ +35−5 ∼ +35−5 ∼ +35+35 ∼ +535
−5 ∼ +35

surement and its time variation are then discussed. Details of the calibration
methods and other calibration items (e.g. relative QE correction, PMT&black
sheet reflection, OD calibration) can be found in [51, 55].

3.6.1

PMT Gain

The first step of the PMT gain calibration is to set the high-voltage (HV) for
each PMT so that the output charge of each PMT is the same for the same lightintensity. As for the calibration source, a xenon lamp is used. Its flash light ispassed through an ultra-violet filter and injected into a scintillator ball which isinstalled at the tank center via an optical fiber, and then isotropic diffused lightis emitted. Since the SK detector is cylindrical, the arrived light intensity to each
PMT depends on the PMT location. In order to compensate this geometricaleffect, 420 “standard PMTs” whose HVs were adjusted to output the samecharge for the same light intensity were mounted in the tank. The HVs of other
PMTs in the tank are adjusted in order to match the output charge to that of
the nearest standard PMT. The reproducibility of the observed charge for all
ID PMTs with respect to their reference values is 1.3%.
The next step is to determine the relative gain difference among PMTs. The
calibration source for this relative gain calibration is a nitrogen laser (USHO
KEC-100, a pulse width of 0.4 ns in FWHM at a wavelength of 337 nm). The
output light is monitored by a faster response 2-inch PMT (Hamamatsu H243150, rise time: 0.7 ns). The wavelength of this light is shifted to 398 nm with apulse width of 0.2 ns by a dye, and the light is passed through neutral filters andinjected into a diffuser ball placed in near the tank center via an optical fiber.
This shifted wavelength corresponds to almost maximum point of the effective
hit taking into account of the Cherenkov spectrum, light absorption and QE
spectrum. By using this laser system, two measurements are conducted. The
first one uses high-intensity flashes (Is ) for every PMT and records the averageobserved charge Qops (i) for each PMT i which is given as;
Qobs (i) ∝ Is × a(i) × QE (i) × G(i),(3.3)

where a(i), QE (i) and G(i) are the acceptance, QE and gain of the i-th PMT.
20

CHAPTER 3. SUPER-KAMIOKANDE EXPERIMENTThe second one uses low-intensity flashes (Iw ) which make a single photoelectron(p.e.) hit for a few PMTs, and counts the number of hits Nobs (i) for each PMT
which is given as;
Nobs (i) ∝ Iw × a(i) × QE (i).

(3.4)

Then, the gain of each PMT is obtained by taking the ratio of the above two
equations namely,
G(i) ∝

Qobs (i).
Nobs (i)(3.5)

The relative gain for each PMT is calculated by normalizing G(i) to the averagegain over all PMTs, and is used as a PMT-by-PMT correction factor in conversion from the output charge to the number of p.e.s. The RMS of the relative
gain distribution is 5.9%.
Finally, the absolute gain for all PMTs is determined. The calibration source
is a nickel-californium ball located at the tank center. From this ball, 6.19.0 MeV low energy gamma-rays are uniformly emitted by capturing neutrons
from the decay of 252 Cf by a nickel nucleus. After applying the relative gain
factor to each PMT, the single p.e. charge distribution for all PMTs is obtained
as shown in Figure 3.8. The averaged single p.e. charge value is determinedfrom the single p.e. distribution after subtracting a similar distribution taken
in an off-timing region of the calibration. PMT gain conversion factors are
determined as 2.055, 2.297, 2.243 and 2.658 pC per single p.e. for SK-I, II,
III and IV, respectively. This single p.e. distribution is used in the detector
simulation.
The time variation of the PMT gain is measured by using the cosmic-ray
muons going through the detector, so-called “through-going” muons. Sincemuons deposit almost constant energy of about 2 MeV/cm independent of theirenergy, the muon track can be regarded as a constant light source. Whileassuming the observed light is the direct light, the observed charge of the i-th
PMT is given as follows,


f (θ)
li
Q(i) = Q0,iexp − 0 ,
(3.6)
li
L
where Q0,i is a constant value, f (θ) is the PMT acceptance as a function of thelight incident angle θ to the PMT, li is the photon travel length and L0 is themeasured attenuation length. Here, the vertical downward muons are used andtheir tracks are reconstructed by connecting the entering and exiting points on
the ID wall. From this equation, the effective observed charge (ln(Q(i)li /f (θ)))
is calculated as a function of l. The intercept of that function represents the
PMT gain. Figure 3.9 shows the time variation of the averaged PMT gain foreach SK. The time dependence of the averaged PMT gain is about 2% per year,
however, the reason of the increase is unknown. In the analysis in SK-I to SKIII, we take into account of the PMT gain correction for these variations in theevent reconstruction step.
21

CHAPTER 3. SUPER-KAMIOKANDE EXPERIMENT

106
105
104
103
102
10
1

0

5

10 15 20 25 30 35 40 45 50
pC

Figure 3.8: Single photoelectron distribution for all PMTs measured by thecalibration using the nickel-californium ball. Taken from [55].

For SK-IV, another method to evaluate the PMT gain time variation isavailable. By taking the early off timing hit distribution for each PMT in thethrough-going muon sample, one can estimate the PMT absolute gain by looking
at the single photoelectron peak. Recently, it revealed that the PMT gain time
variation depends on their production year as shown in Figure 3.10. In theanalysis in SK-IV, the observed charge of each PMT is corrected by using these
PMT gain variations (five production year sets are used).

3.6.2

PMT Timing

The PMT hit timing of each PMT varies due to the variation of the cable length,
electronics processing time and timing response dependence of the charge which
is known as a time-walk effect. In order to correct the hit timing variation, a
two-dimensional time versus charge (TQ) plot is taken for each PMT by using
the nitrogen laser system described above. Time-of-flight (TOF) is subtracted
from the travel time between the diffuser ball to the respective PMT and thenthe hits within a ±50 ns time window of the monitor PMT hit timing arechosen. Figure 3.11 shows a TQ plot. For each charge bin, the TOF subtracted
timing distribution is fitted by an asymmetric Gaussian which is used in thedetector simulation as shown in Figure 3.12 (in the simulation for SK-I to SKIII, a symmetric Gaussian is used instead). A polynomial function as a function

22

PMT Gain (a.u.)

CHAPTER 3. SUPER-KAMIOKANDE EXPERIMENT

1.2

1.15

1.1

1.05

1

0.95
0

1

2

3

4

5

6

7
8
Elapsed Years

Figure 3.9: Time variation of the averaged PMT gain measured by using the
PMTs in the bottom of the tank. The horizontal axis is the elapsed years fromthe beginning of each SK. The colors show the detector periods; black, red,green and blue are for SK-I, -II, -III and -IV, respectively.

Dark noise charge peak w.r.t. April 2009

1.2

1.15

1992-1995 PMT
1996-1997 PMT
2003 PMT
2004 PMT2005 PMT

1.1

1.05

1

0.95
01/01/09

01/01/10

01/01/11

01/01/12

01/01/13

01/01/14

01/01/15

01/01/16

date

Figure 3.10: PMT production year dependence of the PMT gain time variationin SK-IV. Five different production year groups are displayed.

23

CHAPTER 3. SUPER-KAMIOKANDE EXPERIMENT

Q (p.C.)

tqmap of cable # 00010
0 1 2 3 4 5 6 7 8 9 10

1.0

1230

0.9
1220

0.8
0.7

1210

T (ns)

0.6
1200

0.5
0.4

1190

0.3

Entries 561679
1180

0.2
0.1

1170

0

20

40

60

80

100

120

140

160

180

0.0

QBin

Figure 3.11: Example of the TQ plot. The vertical and horizontal axes represent to the TOF subtracted hit timing and the hit charge, respectively. Taken
from [55].

of charge is then obtained by fitting those fitted peak timings and standard
deviations, so-called “TQ map.” The TQ map of each PMT is used to correct
the hit timing.

3.6.3

Water Transparency

The light absorption and scattering in the water are measured by comparing
calibration data results with the detector simulation. Detailed discussion can beseen in the document [51, 55]. As for the data, a nitrogen laser light is injected
into the tank via an optical fiber and the hit timing distribution of each PMT
are recorded for each five different wavelength of the laser; 337, 375, 405, 445
and 473 nm. When the light of wavelength λ travels a distance l, its intensity
is simulated to decrease as exp(−l/L(λ)) in the detector simulation. Here, L(λ)
is the attenuation length in the water which is defined as;
L(λ) = (αsys (λ) + αasy (λ) + αabs (λ))−1 ,(3.7)

where αsys (λ), αasy (λ) and αabs (λ) are empirical coefficients of the light “symmetric” scattering, “asymmetric” scattering and absorption, respectively. Rayleighscattering and the symmetric component of Mie scattering are considered as thesymmetric scattering, and the asymmetric component of Mie scattering is forthe asymmetric scattering. These coefficients are given as the following empirical functions, and their nine parameters P0−8 are then obtained by fitting to24

CHAPTER 3. SUPER-KAMIOKANDE EXPERIMENT

Figure 3.12: Timing distribution of the charge bin 14 (single photoelectronlevel) for all ID PMTs. The red points show the data and the blue line shows
the fitted asymmetric Gaussian function. Taken from [55].

the data;
αabs (λ) =αsym (λ) =αasy (λ) = P6

P 0 P1+ C(λ),
λ4

P4
λ4


1.0 +

P5
λ2

(3.8)


,


P7
2
1.0 + 4 (λ − P8 ) ,
λ

(3.9)

(3.10)

where C(λ) is the factor derived by the measurement result for λ ≥ 464 nm [56],whereas the following equation is employed for λ ≤ 464 nm,
C(λ) = P0 P2

λ
500

P 3
.

(3.11)

Figure 3.13 shows the fitted coefficients which are used in the detector simulation.
The time variation of the water transparency is measured by using thethrough-going muon sample. From Equation 3.6, the water attenuation length
L can be evaluated. Figure 3.14 shows the time variation of the measured water25

CHAPTER 3. SUPER-KAMIOKANDE EXPERIMENT

Inverse water transparency (1/m)

10-1
Total

Absorption

Symmetric

Asymmetric

10-2

10-3

10-4
300

325

350

375

400

425

450

475

500

Wavelength (nm)Figure 3.13: Fitted coefficients of the water transparency as a function of the
wavelength of the nitrogen laser. The dotted lines are the fitted functions andthe solid line is the sum of them. The points are the data measured in April
2009. Taken from [55].

attenuation length in the detector, which is used in the event reconstruction
algorithm for the water attenuation correction. Thanks to the sophisticatedwater control system and closing the tank since the beginning of SK-III, thevariation of the water quality in SK-IV is the smallest.

26

Water Attenuation Length (m)CHAPTER 3. SUPER-KAMIOKANDE EXPERIMENT

130
120
110
100
90
80
70
60
50
1996 1998 2000 2002 2004 2006 2008 2010 2012 2014 2016 2018
Date

Figure 3.14: Time variation of the water attenuation length from the beginningof SK to Aug. 31st, 2016. The colors show the detector periods; black, red,green and blue are for SK-I, -II, -III and -IV, respectively. Since the detector
conditions such as the water circulation, OD segmentation Tyvek sheet, electronics and so on have been changed, it is hard to compare the absolute water
attenuation length between each SK period.

27

Chapter 4Monte Carlo Simulations
In this chapter, we describe how the signal and background events for p → e+ π 0
are simulated in this analysis, including the detector behavior. We employ aMonte Carlo (MC) method. As for the background, the atmospheric neutrinos
are only considered since other noise events are negligible and can be highly
rejected by the data reduction described in Chapter 5.

4.1

SignalIn this analysis, we assume that proton decay equally happened in each protonof H2 O which contains two and eight protons in hydrogen and oxygen nuclei,
respectively. Protons in hydrogen which are called “free protons” are stationary
and do not interact with other nucleons. Free protons decay into a positron anda neutral pion back-to-back with a monochromatic momentum of 459.43 MeV/c.
Then, neutral pions immediately decay into two gammas (the branching ratioof this mode is 98.8% and the rest is for Dalitz decay). On the other hand,
protons in oxygen which are called “bound protons” suffer from several effects
as described below.
For the decay of bound protons, we first consider the Fermi motion and thenuclear binding energy in 16 O. We simulate the nucleon momentum in 16 O
based on an electron-12 C scattering experiment [57]. Figure 4.1 (left) shows
the simulated proton momentum distributions for both the s-state and p-stateof 16 O. The binding energy is simulated as a Gaussian distribution, and thiseffect is considered by subtracting the binding energy from the proton mass in
the rest. The mean and σ of the binding energy distributions are (39.0 MeV/c,
10.2 MeV/c) for the s-state and (15.5 MeV/c, 3.82 MeV/c) for the p-state. The
simulated proton invariant mass in 16 O is shown in Figure 4.1 (right).
In addition to the Fermi motion, we consider an effect from a spectator
nucleon in 16 O. When a proton is going to decay and its wave function isoverlapping the other one, this proton decay becomes a two-body system decay
and the momentum of the decaying proton is changed. This decay is called the28Number of Events

Number of Events

CHAPTER 4. MONTE CARLO SIMULATIONS

104
103

104
103

102

102

10

10

1

1
0

0

100 200 300 400 500 600 700 800 900 1000
True Proton Momentum (MeV/c)

200

400
600
800
1000
1200
True Proton Invariant Mass (MeV/c2)

Figure 4.1: Simulated proton momentum (left) and invariant mass (right) distributions. The cyan, black, green, red and purple histograms correspond to
free, bound, s-state, p-state and correlated decay protons, respectively.
correlated decay. According to T. Yamazaki and Y. Akaishi [58], the predictedprobability of the correlated decay is about 10%. Figure 4.1 shows the proton
momentum and invariant mass distributions in 16 O.
While a positron emerged in 16 O can immediately escape without any interaction, a neutral pion often interacts hadronically with nucleons via inelastic scattering, absorption or charge exchange which converts a neutral pion toa charged pion. These pion final state interactions (π-FSI) are simulated in
NEUT [59], and the cross sections for each interaction are calculated from acascade model [60] which was tuned by various π-(p, n) and π-N scattering
data [61].
Figure 4.2 shows the comparison between the tuned π ± -12 C scattering crosssection in NEUT and the experimental data. The cumulative fractions of theπ 0 -FSI are shown in Figure 4.3. The momentum and angular distributions of
the scattered pions are determined from the results of the phase shift analysis ofthe π-N scattering experiments [62]. Because of the Pauli exclusion principle,
the FSI nucleon momentum should be larger than the Fermi surface momentum
given as follows;

 31
3 2
pF (r) =π ρ(r),(4.1)
2
where ρ(r) is the nuclear density as a function of the distance r from the nuclear
center. We use the Woods-Saxon nuclear density model [63] for the calculationof ρ as follows.
ρ(r) =

ρ(0),1 + exp( r−a
b )

(4.2)

where ρ(0) = 0.48m3π , mπ is the invariant mass of pions, a = 2.69 fm is the
maximum nuclear radius of 16 O and b = 0.41 fm is the surface thickness of 16 O.
Figure 4.4 shows both ρ and pF as a function of r.
29

CHAPTER 4. MONTE CARLO SIMULATIONS

Figure 4.2: Cross sections for π + -12 C (left) and π − -12 C (right) as a function ofpion initial momentum. The dots are for the experimental data, the solid lines
are for the calculation by NEUT and the dotted lines are for the old FSI model.
Taken from [61].

1Charge exchange

π production

0.8

Fraction

Scattering

0.6
Absorption

0.4

0.2

0
300

No interaction

350

400

450

500

550

π0 Initial Momentum (MeV/c)

600

Figure 4.3: Cumulative fractions of the π 0 -FSI as a function of π 0 initial momentum. Taken from [1].

30

0.08

220

0.07

200
180

0.06

160

0.05

140
120

0.04

100
0.03

80

0.02

60
40

0.01
0
0

20
1

2

Fermi Surface Momentum (MeV/c)-3

Nuclear Density (fm )CHAPTER 4. MONTE CARLO SIMULATIONS

0
3
4
5
6
7
8
9
10
Radial Distance from the Center of Oxygen (fm)Figure 4.4: Woods-Saxon nuclear density (black) and the Fermi surface momentum (red) as a function of the radial distance from the center of 16 O.
After proton decay in 16 O, the remaining nuclei can be in an excited stateand it de-excites through gamma-rays or nucleon emission. We simulate thoseemissions based on H. Ejiri[64].

4.2

BackgroundAs for the background of proton decay, atmospheric neutrinos are only considered since other noise events are negligible. In this section, we discuss how
many atmospheric neutrinos are expected to come in the SK detector and theirinteractions in water.

4.2.1

Neutrino Flux

The atmospheric neutrino flux at the Super-Kamiokande site is calculated by
M. Honda et al. [65, 66], called the Honda flux. The Honda flux is estimatedbased on the primary cosmic ray flux measurements by the AMS [67, 68] andBESS [69, 70] experiments, and takes into account the effects of the solar activity
and the geomagnetic field to the primary cosmic ray flux. The interaction of
the primary cosmic ray with the atmosphere is calculated based on JAM [71](primary cosmic ray energy for below 32 GeV) and DPMJET-III [72] (for above
32 GeV). Neutrino oscillations are not considered in the atmospheric neutrinoMC, but those effects are taken into account by reweighing each MC event in
the analysis step.
The expected zenith angle distributions of the atmospheric neutrino flux at
SK by the Honda flux are shown in Figure 4.5. The peaks at the horizontal
31

CHAPTER 4. MONTE CARLO SIMULATIONS

4x10

3

3x10

2

0.32 GeV

3.2 GeV

1.0 GeV

νµ
10

3

νe

νµ

φν (m−2s−1sr−1GeV −1)φν (m−2s−1sr−1GeV−1)φν (m−2s−1sr−1GeV−1)

10

νµ

10

νµ

2

νe

1

νµ
νµ

5x10

0

νe

νe

νe
2x10

νe
4x10

2

−1.0

−0.50

0

0.50

1.0

3x10

0

1

−1.0

−0.50

0

cos θ

cos θ

0.50

1.0

−1.0

−0.50

0

0.50

1.0

cos θ

Figure 4.5: Expected zenith angle distributions of the atmospheric neutrino flux
at SK by the Honda flux for the neutrino energies of 0.32 GeV (left), 1.0 GeV(middle) and 3.2 GeV (right). Taken from [66].

direction are due to a longer travel length of cosmic ray particles coming horizontally. Since there is the geomagnetic field, a large up-down asymmetry canbe seen in the lowest energy plot. The direction-averaged atmospheric neutrino
flux at SK by the Honda flux is shown in Figure 4.6 with other flux models
such as the FLUKA [73] and the Bartol [74]. For this proton decay analysis,
the neutrinos with the order 1GeV energy can be the dominant background.

4.2.2

Neutrino Interactions

The interactions of neutrinos with the water are simulated by NEUT [59]. In
the atmospheric neutrino MC, the interactions of neutrinos with hydrogen andoxygen nuclei are only simulated while neglecting the interaction of neutrinos
with electrons since the cross section for neutrino-electron interactions is threeorders of magnitude smaller than that for neutrino-nuclei interactions.
As described in the following sentences, there are roughly three categories forboth Charged Current (CC) and Neutral Current (NC) interactions; (Quasi-)
Elastic (QE) scattering, resonant (RES) productions and Deep Inelastic Scattering (DIS). Figure 4.7 shows the distributions of the CC νµ and ν¯µ cross sections
per nucleon as a function of neutrino energy in NEUT with the experimental
data [75].
The particle type of neutrinos generated in the atmosphere are only νe andνµ . However, ντ can be arrived in SK through neutrino oscillations. The decayof produced τ is complex and is simulated by TAUOLA [76].

32CHAPTER 4. MONTE CARLO SIMULATIONS

10

1.5xν µνe

Neutrino Flux Ratiosφν x E ν2 (m−2sec−1sr−1GeV)

5
2

0.75xν e101νµThis WorkHKKMS06
Bartol
Fluka

0

10

−1

10

0

10
Eν (GeV)

10This WorkHKKMS06
BartolFluka

νe + νeνµ2

ν µ ( x1.5)
νe

1

1

νµ + νµ

−1

10

0

10
Eν (GeV)

10

νe

1

Figure 4.6: Direction-averaged atmospheric neutrino flux at SK by the Honda
flux. The absolute flux is on the left and the the flux ratio is on the right. Taken
from [66].

Figure 4.7: Cross sections per nucleon divided by neutrino energy as a functionof neutrino energy for CC νµ (left) and ν¯µ (right). Taken from [77].

33CHAPTER 4. MONTE CARLO SIMULATIONSElastic and Quasi-Elastic Scattering
In NC elastic scattering, a neutrino just scatters off a nucleon without creating
any particles. On the other hand, a neutrino is converted into a corresponding
charged lepton in CC Quasi-Elastic (CCQE) scattering.
ν + N → l + N0

(4.3)

Here, N and N 0 are the initial and final nucleons. This interaction is the dominant process up to around 1 GeV. For free protons, the CCQE cross sections
are simulated based on the Llewellyn-Smith model [78]. The modification model
provided by Smith and Moniz [79] is used for a nucleon in 16 O (a relativistic
Fermi gas model). The differential cross section of the model is a function of
the square of the four-momentum transfer. The axial vector mass contained ina form factor of this model is set at 1.21 GeV based on the measurements of
K2K [80] and MiniBooNE [81].
In the simple CCQE interaction described above, we assume that a neutrino
interacts only with a single nucleon. However, the nucleons which are not thetarget nucleon may affect the interaction process, and such effect is suggested by
MiniBooNE result [81]. This is known as the Meson Exchange Current (MEC)
model [82]. Recently, NEUT has been modified to include the MEC effect which
treats the interaction of a neutrino with a pair of two nucleons only for CCQE
interactions as follows.
ν + N N 0 → l + N 00 N 000(4.4)

Resonant Productions
By the delta resonance through neutrino interactions, a single meson or photon
is produced as follows,ν + N → l + N 0 + (π, γ, K, η).

(4.5)

In the previous version of NEUT used in the previous p → e+ π 0 analysis [1],this is simulated based on the Rein and Sehgal model [83]. This model has
been tuned using form factors from Graczyk et al. [84] in NEUT used in thisthesis. The parameters in form factors of the modified model are obtained fromthe Argonne National Laboratory and Brookhaven National Laboratory bubble
chamber experiments [85]. In this p → e+ π 0 analysis, νe CC events with 1 π 0(νe CC1π 0 ) become the dominant background as described in Chapter 7. By
adopting this better model, the momentum transfer of W bosons is decreased
and the outgoing lepton momentum becomes larger compared to the previous
model. Thus, the total momentum of the νe CC1π 0 events goes higher. Figure 4.8 shows the true outgoing electron momentum and true total momentumof νe CC1π 0 events within a neutrino energy from 1 GeV to 3 GeV for the previous and current MC. As a result, the background for p → e+ π 0 is reduced.

34

CHAPTER 4. MONTE CARLO SIMULATIONS

MC11d=1147.4,MC16a=1097.5,Error=0.045
Number of Events

Number of Events

MC11d=361.2,MC16a=327.1,Error=0.104
70
60
50

240
220
200
180
160
140

40

120
100

30

80
20

60
40

10

20
0
0

100

0
0

200
300
400
500
600
700
800
True Lepton Momentum (MeV/c), Ptot<500MeV/c

500

1000

1500

2000
2500
3000
True Total Momentum (MeV/c)

Figure 4.8: True outgoing electron momentum (left) and true total momentum(right) of νe CC1π 0 events within a neutrino energy from 1 GeV to 3 GeV for
the previous (black) and current MC (red). In the left figure, the events withthe total momentum of less than 500 MeV/c are displayed.
Coherent pion production is a process in which a neutrino interacts with theentire oxygen nucleus and produces a forward-going pion as follows,ν +16 O → l +16 O + π.

(4.6)

The Rein and Sehgal model [86] is used in NEUT for the coherent pion production.

Deep Inelastic Scattering
In a higher energy region above around 10 GeV, DIS becomes the dominant
process. Such a high energy neutrino can interact with a constituent quark in
the target nucleon and often produce multiple hadrons as follows,ν + N → l + N 0 + hadrons.

(4.7)

This process is simulated based on the GRV98 parton distribution functions [87]and the corrections by Bodek and Yang [88]. The produced hadrons are only
considered as pions in the region of W < 2 GeV, and its multiplicity is determined by the bubble chamber experiments [89, 90]. In the region of W > 2 GeV,
heavier mesons (K and η) are also considered as the hadron production, which
are simulated by PYHIA/JETSET [91].

4.3

Detector

Particles escaped from the nucleus in concern are simulated by a GEANT3 [92]based detector simulation, called SKDETSIM. Particle interactions with water,
particle decays and Cherenkov photon production are simulated by GEANT. For
the hadronic interaction, GCALOR [93] and the NEUT pion cascade model areused for the interaction with the pions above and below 500 MeV/c, respectively.
35

CHAPTER 4. MONTE CARLO SIMULATIONSThe Cherenkov photon propagation in water is simulated by a custom code
which includes Rayleigh/Mie scattering and the absorption of photons by H2 O
and is tuned by the calibration results. For photons with short wavelength λ(≤ 450 nm), Rayleigh scattering is the dominant effect caused by small particles(r  λ, where r is the radius of a particle), which has the λ−4 dependence andscatters photons symmetrically in the forward and backward directions. The
absorption of photons by H2 O is dominant for long wavelength (≥ 450 nm).
Mie scattering is caused by large particles (r  λ) and makes a sharp peakin the forward region. The reflectivity of both the PMTs and black sheets ismodeled and tuned by using the calibration results.
The PMT responses in both the time and charge and the electronics are also
simulated based on the calibration measurements.
The output data structure of SKDETSIM is the same as that of the observed
data. Thus, we can analyze both the data and MC events in the same way. In
this thesis, 30,000 events of the p → e+ π 0 MC generated randomly in the ID
and the 3000 (500) years atmospheric neutrino MC for SK-IV (SK-I to -III) are
used.

36

Chapter 5

Data Set
In this analysis, we use fully contained (FC) events whose vertices and all visible particles are contained inside the ID, in order to reconstruct the invariant
mass and momentum of protons precisely. Most of FC events consist of atmospheric neutrino events. There are two other types in atmospheric neutrino
events; partially contained (PC) and upward-going muon (UPMU). PC events
are the events whose vertices are contained inside the ID, but any of the charged
particles go outside. UPMU events are the high energy muons generated by
upward-going neutrinos interacted with the rock surrounded the detector.
In order to prepare the FC data set, we need to remove background events
caused by cosmic ray muons, radioactivities, solar and atmospheric neutrinos
and so on. Through the following five reduction steps, we can reject background
events except FC atmospheric neutrinos. Figure 5.1 shows the schematic drawing for various event types in SK. The reduction criteria are common for thedata, atmospheric neutrino MC and proton decay MC.

5.1

First Reduction

The event trigger is fired at a rate of 106 times a day, most of which are low
energy events coming from radioactive background and cosmic ray muon events.
Those low energy background events can be largely rejected by applying thefollowing criteria.
1. PE300 ≥ 200 (≥ 100 for SK-II)
2. NHITA800 ≤ 50 (≤ 55 for SK-IV) or OD trigger is off
where PE300 is the maximum number of total photoelectrons (p.e.s) observedwith ID in a sliding 300 ns time window and NHITA800 is the number of hit
OD PMTs within a time window [-400 ns, +400 ns]. In the first cut, thecut threshold roughly corresponds to 22 MeV/c electrons. Therefore, most ofradioactive background, low energy Michel electron from cosmic ray muons andsolar neutrinos are rejected. Since the PMT density of SK-II is half that of37

CHAPTER 5. DATA SET

FC neutrino
Cable hole μ
Invisible μ
VETO

Elec.

Stopping μ
Through-going μRadioactive BG
PC neutrino

Michel e
FlasherProton decay
Michel e

Rock
UPMU neutrino

Figure 5.1: Schematic drawing for various event types in SK. Blue shows theemission of Cherenkov light.

38

CHAPTER 5. DATA SET

the other detector periods, the half threshold is applied for SK-II. The second
cut rejects not only cosmic ray muons, but also PC and UPMU events. The
threshold in the second cut for SK-IV is adjusted to correct the changes in the
OD electronics boards and the PMT gains. After the first reduction, the event
rate is decreased to the order of 103 events/day.

5.2

Second Reduction

A more stringent rejection is performed in the second reduction compared to
the first one. The events passed through the following criteria proceed to thethird reduction.
1. PEmax /PE300 ≤ 0.5
2. NHITA800 ≤ 25 (≤ 30 for SK-IV)
or PEtot ≥ 100, 000 (≥ 50, 000 for SK-II)
where PEmax is the highest observed p.e.s among all hit ID PMTs and PEtot is
the sum of the observed p.e.s in the ID inside the time window. The first cut
can reject electrical noise events caused by a single PMT which has a large hit.
More rejection in OD activated events is done by the second cut while accepting
very high energy FC events which have light leakage into the OD side. Afterthe second reduction, the event rate decreases to the order of 100 events/day.

5.3

Third Reduction

The third reduction removes specific background events such as; hard muons,through-going muons, stopping-muons, cable hole muons, coincidence muons,
flasher and low energy events. After the third reduction, the event rate decreases
to around 40 events/day.

Hard Muon Cut
Cosmic ray muons which have an energy of ≥ 1 TeV, “hard” muons, would be
rejected by the following cut;
• NTHITA500 ≥ 40
where NTHITA500 is the number of OD observed p.e.s in a sliding 500 ns time
window.

Through-Going Muon Cut
Through-going muons which satisfy the following all criteria are rejected.
• PEmax > 231• Number of hit ID PMTs ≥ 1000
39

CHAPTER 5. DATA SET• NHITAin > 10 or NHITAout > 10
• Goodness of the through-going muon fit > 0.75
where NHITAin (NHITAout ) is the number of hit OD PMTs in the muon entrance (exit) hit cluster. The goodness of a through-going muon event is evaluated by using the observed hit time, its resolution and the expected hit time
calculated by the muon entering time and its track information.

Stopping Muon Cut
Similar to through-going muon cut, stopping muon events are rejected by thefollowing cut;
• NHITAin ≥ 10 when the goodness of the stopping-muon fit ≥ 0
and only for SK-I, NHITAin ≥ 5 when the goodness of the stopping-muon
fit ≥ 0.5

Cable Hole Muon Cut
Twelve cable holes for PMT’s signal and high voltage cables are located on the
top of the tank. No OD PMTs were installed under four holes out of twelve.
These holes admit cosmic ray muons to enter the detector without OD activity.
To reject those “cable hole muons,” four veto counters (2.0 m × 2.5 m plastic
scintillation counters) were installed on each hole. The reduction criteria are as
follows;• One veto counter hit
and a distance from the reconstructed vertex to the cable holes in boththe x and y direction is ≤ 4 m

Coincidence Muon Cut
Cosmic ray muons sometimes enter the detector just after a SLE trigger bycoincidence. The following cuts reject these events.
• Number of hit OD PMTs in a fixed time window between +300 ns and+800 ns ≥ 20
and number of ID observed p.e.s in a fixed time window between +300 ns
and +800 ns > 5000 p.e.s (2500 p.e.s for SK-II)

Flasher Event Cut
Light emission from a PMT itself occurs sometimes due to an electrical discharge, which may look a FC-like event. This is called “flasher” event. Usually
the hit timing distribution of flasher events is longer than that of particle induced events. Flasher events which emit light for a longer time can be rejectedby the following cuts;40CHAPTER 5. DATA SET• NMIN100 ≥ 20 (15 for SK-I)or, only for SK-I, NMIN100 ≥ 10 if the number of hit ID PMTs ≤ 800• Goodness of the time-of-flight fitter ≤ 0.4
where NMIN100 is the minimum number of hit ID PMTs in a sliding 100 ns time
window between +200 ns and +700 ns. Since the TOF vertex fit for a flasher
event will not work well, the second term is applied for the reduction.

Low Energy Event Cut
The remaining low energy events such as the radioactive background and electrical noise are subtracted by;• NHIT50 < 50 (25 for SK-II)
where NHIT50 is the number of time-of-flight corrected ID hits in a sliding 50 ns
time window. This threshold corresponds to the visible energy of about 9 MeV.

5.4

Fourth Reduction

Flasher events are further rejected by the fourth reduction using a likelihood
method. It is empirically known that flasher PMTs repeatedly emit light in aparticular way and leave similar hit patterns. A typical flasher event is displayedin Figure 5.2. The flasher likelihood consists of three parameters evaluatedbetween two events; a correlation factor r using the charge patterns, a parameterdKS comparing the time distribution of the charge by a Kolmogorov-Smirnov
test and a parameter which indicates whether the PMT with the highest p.e. isconsistent between two events or not. Dividing the ID wall into 1450 patches of
2 m × 2 m square, the correlation factor r is calculated by;Nr=

1 X (QA,i − hQA i) (QB,i − hQB i)
N i
σA σB

(5.1)

where N is the number of patches, QA(B),i is the total observed p.e.s in thei-th patches for the event A (B), QA(B) and σA(B) are the averaged observed
p.e.s and their RMS. The likelihood cut value is determined by comparing theatmospheric neutrino MC and a hand selected real flasher sample. Both two
events are rejected as flasher when their likelihood value is higher than the cut
value. The fourth reduction removes a few events/day.

5.5

Fifth Reduction

The fifth reduction rejects the remaining cosmic ray muons and flasher events,
electronic noise and calibration related events. The fifth reduction removes a
few events/day.
41

CHAPTER 5. DATA SET

Figure 5.2: Typical flasher event display in SK-IV. The color indicates the hit
timing (red: early hit, blue: late hit). A PMT located in the center of the tanksuddenly flashed light. This type of the hit pattern hardly occurs by proton
decay and atmospheric neutrino events.

42

CHAPTER 5. DATA SET

Stopping Muon Cut
The improved reduction criterion is used for stopping muon events.
• NHITA0in ≥ 10or,
• Goodness of the stopping muon fit ≥ 0.5 and NHITA0in ≥ 5
where NHITA0in is the number of hit OD PMTs in the muon entrance (exit) hit
cluster without an OD timing correction.

Cable Hole Muon Cut
Considering the case of no veto counter hit, but cosmic ray muons enter thedetector through the cable holes, tighter cuts based on the reconstruction parameters are necessary. The following criteria are used only for SK-IV.
• Goodness of stopping muon fit ≥ 0.4
and PEtot > 4000
and cos(z) of muon’s direction < −0.6
and the distance between the muon entering point and the nearest cablehole < 250 cm• NHITAAPfit ≥ 4
where NHITAAPfit is the number of hit OD PMTs within 8 m from the muon entering point fitted by an event reconstruction algorithm, APfit (see Chapter 6),
in a sliding 200 ns time window from −400 ns to +400 ns.

Coincidence Muon Cut
In addition to the coincidence muon cut in the third reduction, a tighter cut isapplied as follows.
• PE500 < 300 p.e.s (150 p.e.s for SK-II) and PElate ≥ 20 p.e.s
where PE500 is the total number of observed ID p.e.s in a fixed time window
between −100 ns and +400 ns and PElate is the maximum number of hit OD
PMTs in a sliding 200 ns time window from +400 ns to +1600 ns.

Invisible Muon Cut
Cosmic ray muons below the Cherenkov threshold, “invisible” muons, can enter
the ID without PMT hits and decay into Michel electrons. These Michel electrons can emit light inside the ID, then it looks like FC event. The invisible
muon events which satisfy the following criteria are rejected;• PEtot < 1000 p.e.s (500 p.e.s for SK-II)
43

CHAPTER 5. DATA SET• NHITACearly ≥ 5• NHITACearly + NHITAC500 ≥ 10 (if DISTclust < 500 cm)
where NHITACearly is the maximum number of hit OD PMTs in a sliding 200 ns
time window from −8900 ns to −100 ns, NHITAC500 is the number of hit OD
PMTs in a fixed time window between −100 ns and +400 ns and DISTclust
is a distance between two OD hit clusters used for the NHITACearly and the
NHITAC500 .

Long-Tail Flasher Event
A tighter cut for the flasher events than the third reduction is conducted asfollows;• NMIN100 ≥ 6 and the goodness of the TOF fitter < 0.4
or,• NMIN100 < 6
and the goodness of the TOF fitter < 0.3 for SK-II to SK-IV only

Electronics Related Event
The electronic noise from the high voltage systems or electronics boards cancreate a lot of fake hits. The following criterion removes the electronic noise
events.
• N0 ≥ 250 (125 for SK-II) and N0 − N1 ≥ 100 (50 for SK-II)
where N0 (N1 ) is the number of hit ID PMTs with a less (greater) than a single
p.e.

Dead Electronics Hut
Due to some trouble, it could happen the electronics huts turn to be off. In orderto check whether the hut is active or not, the following criterion is necessary.
• All electronics huts have at least one OD hits.

Wrong OD Timing
In the phase of SK-I to SK-III, the OD timing was occasionally out of thesynchronization with the ID timing. This problem caused large early timing
clusters in the OD or the abnormal number of sub-events. Thus, the following
criteria are applied for SK-I to SK-III only.
• PEtot ≥ 1000 (500 for SK-II)
and NHITACearly ≥ 20 and the number of sub-events is 7.

44

CHAPTER 5. DATA SET

Calibration Related Event Cut
Some events for the calibration by a laser or a xenon lamp can be contaminatedin the data set due to some trouble. But, using the information from themonitored PMTs for the calibration sources, it is possible to remove those events
as follows.
• Number of reference PMT hits for the laser and the xenon lamp is zero.

5.6

Final FC Data Set

Finally, the FC data set for this analysis is obtained after applying the following
criteria.
• Dwall > 1.5 m for fiTQun (2.0 m for APfit)• NHITAC < 16 (< 10 for SK-II)• Evis > 30 MeV
where Dwall is the distance from the wall and the reconstructed vertex, NHITAC
is the number of hit OD PMTs in the largest charge cluster and Evis is the visibleenergy which is defined as the total energy of all reconstruction rings assuming
all of them are originated by electrons. FiTQun is a newly developed event
reconstruction algorithm and has better performances than APfit as describedin Chapter 6.
In addition to the automated FC reduction, a real time check in the data
quality is conducted by human eyes for the events after the third reduction.
The scanners check the hit charge and timing distribution and the hit pattern
for both the ID and OD using an event display tool, judge whether the reduction algorithm works properly or not and be aware of something unusual (forinstance, Appendix A).
The final event rate is about 8 events/day. The reduction efficiencies of eachreduction step for the data, atmospheric neutrino MC and p → e+ π 0 MC in
SK-IV are shown in Table 5.1.

45

CHAPTER 5. DATA SET

Table 5.1: Reduction efficiencies of each reduction step for the data (from Sep.
2008 to Aug. 2016, 2650.4 live days), atmospheric neutrino MC and p → e+ π 0
MC in SK-IV. The true-Dwall > 2 m cut is used to the calculation of the rates
for MC. Here, the fiducial volume cut in the final reduction is Dwall > 2 m.
Reduction
1st
2nd
3rd
4th
5th
Final

Data
1771.31 day−1
274.40 day−1
39.55 day−1
37.19 day−1
34.46 day−1
8.08 day−1(4694674 events)(727266 events)(104830 events)(98575 events)(91336 events)(21410 events)

46

Atm.-ν MC

p → e+ π 0 MC

99.99%
99.99%
99.83%
98.87%
98.83%
98.07%

99.98%
99.93%
99.72%
98.19%
97.94%
97.94%

Chapter 6

Event Reconstruction
The event reconstruction algorithm (fitter) works to determine the event properties such as the vertex, number of Cherenkov rings, particle type and momentum. There are two different fitters; one is an existing fitter (APfit) used since
the beginning of SK, the other one is a newly developed fitter based on the
maximum likelihood method (fiTQun).
In this Chapter, we first describe the reconstruction scheme of both fittersand then discuss their performance. After that, the fiducial volume and energy
scale for fiTQun are discussed. The same algorithm is applied to both the
fully-contained data and MC. FiTQun is applied only for SK-IV while APfit isapplied for SK-I to SK-IV.

6.1

Overview of Two Algorithms

The conventional event reconstruction algorithm, APfit, was developed in late1990’s and contributed to the discovery of atmospheric neutrino oscillation and
the K2K and T2K experiments. APfit is a step-by-step fitter by using the timeand charge information of hit PMTs. Thanks to the remarkable evolution of
CPU power, we can reconstruct events more precisely. The improved event
reconstruction algorithm, fiTQun, is a maximum likelihood fitter which uses
not only hit T&Q information but also the information of non-hit (“unhit”)
PMTs. Using simultaneous fitting method and more information, the precisionof reconstruction of fiTQun is much better than that of APfit. FiTQun is able tofit the vertex ring-by-ring, number of rings up to six and particle identification
under the three particle hypotheses (e, µ, charged-π). On the other hand, APfit
can fit the single-vertex, number of rings up to five and particle identification
under the two particle hypotheses (e, µ).
Based on the reconstruction method used in the MiniBooNE experiment [94],fiTQun was initially developed by the T2K collaboration for the SK-IV detector
as a far detector of T2K. For SK analyses, we have improved fiTQun not only for
SK-IV but also the other detector periods. We have validated fiTQun for SK-IV

47

CHAPTER 6. EVENT RECONSTRUCTION

Table 6.1: Characteristics of fiTQun and APfit.

Type
VertexRing Counting
Particle Identification
Availability (2016)

fiTQunAPfit

Maximum Likelihood Fitter
Multi-Vertices
6-rings (max.)e, µ, charged-π
SK-IV

Step-by-Step Fitter
Single-Vertex
5-rings (max.)e, µSK-I to SK-IV

analysis, but we need more tuning of fiTQun including the detector simulation
for the other SK phases (see Appendix C). Therefore, we use fiTQun for thedata analysis in SK-IV. Table 6.1 shows the characteristics of fiTQun and APfit.
The process speed of the fiTQun (APfit) multi-ring fit in our computer system
is about 8.8 min (1.3 min) per events for the atmospheric neutrino MC events
preselected as the background events of p → e+ π 0 described in Appendix D.

6.2

APfit

In this section, we review the reconstruction scheme of APfit. The detail information can be found in [95]. The event reconstruction of APfit is done asfollowing steps.
1. Vertex fitting
2. Ring counting
3. Particle identification
4. Momentum determination
5. Michel electron search
For single-ring events determined in the second step, a more precise vertex fit
is done by using the Cherenkov ring pattern depending on the particle type.
However, it is not available for multi-ring events and thus not used in thisanalysis.

6.2.1

Vertex Fitting

The first step of the reconstruction is to define the vertex position. There arethree processes: first, “point-fit” roughly finds the vertex, then a search for theedge of the most energetic Cherenkov ring and its direction is performed usingthat vertex. Finally, “TDC-fit” redetermins the vertex taking into account thefinite particle track length and the scattered light.

48

CHAPTER 6. EVENT RECONSTRUCTION

Point-fit is performed assuming that Cherenkov light is emitted from a single point. For a test vertex x, the time-of-flight (TOF) subtracted hit timing
distribution is calculated. That distribution becomes narrower when x comesclose to the true vertex. By maximizing the following goodness function, thepoint-fit vertex can be derived.

1 X(ti (x) − t0 )2
Goodness =
exp −
(6.1)
N i
2(1.5σ)2
where N is the number of hit PMTs, ti is the TOF subtracted timing of thei-th PMT, t0 is the interaction time maximizing the goodness, the factor 1.5 ischosen to optimize the fitting performance and σ is the PMT timing resolutionof 2.5 ns.
The next step is to search the ring edge and the direction of the most energetic Cherenkov ring. The observed charge distribution depending on the
Cherenkov angle θ, PE(θ), can be calculated using the point-fit vertex and assuming a test direction, which is shown in Figure 6.1. Then, the ring edge θedge
under that assumption can be found where the following conditions are satisfied.
d2 PE(θ)
=0
d2 θ

(6.2)

By maximizing the following estimator Q(θedge ), the edge and the direction of
the most energetic Cherenkov ring are determined.
!2
R θedgePE(θ)dθ(θedge − θexp )2dPE(θ)0exp −Q(θedge ) =(6.3)
sin θedgedθ
2σθ2θ=θedgewhere θexp and σθ are the expected Cherenkov angle using the total charge
inside the cone and its resolution, respectively.
Finally, TDC-fit adjusts the vertex more precisely taking into account thefinite particle track length and the scattered light. The initial track lengthis estimated from the total charge and an assumption of 3 MeV/cm energy
deposits. The timing residual for PMTs inside the cone is calculated assuming
the Cherenkov light is emitted along the track with the same angle. For PMTs
outside the cone, the timing residual is calculated assuming the Cherenkov lightis emitted at the vertex considering the scattering parameters. The final vertexof the most energetic ring is determined by maximizing an estimator.

6.2.2

Ring Counting

The next step is to search for additional Cherenkov rings. Candidate rings
are found by a shape extraction algorithm known as the Hough transform.
Figure 6.2 shows an illustration of this idea. For each hit PMT, virtual circles(dashed line) centered on each PMT position are drawn with 42◦ half angle.
The direction of the intersection of those virtual rings is used as the directionof the additional ring. Then, the candidate ring are judged whether it is a fake49

CHAPTER 6. EVENT RECONSTRUCTION

number of p.e.s

(i)
150
100
50
0

second derivative

10

0

20

40

60

80

(ii)

0

-10

detected edge
0

20

40

60

80

θ(degree)Figure 6.1: An observed charge distribution depending on the Cherenkov angle(top) and its second derivative (bottom). Taken from [96].

50CHAPTER 6. EVENT RECONSTRUCTION
42 deg. ringhit PMT(possible center)

centerCherenkov ring

(most probable)

Figure 6.2: Illustration of the ring finding principle. Taken from [96].

ring or not by a log likelihood method. The likelihood consists of the observed
and expected charge. This ring counting continues up to five rings by default.
Figure 6.3 shows the ring counting likelihood distribution for the FC sub-GeV(visible energy < 1.33 GeV) sample of the data and the atmospheric neutrino
MC.

6.2.3

Particle Identification

After the ring counting, the particle type of each ring is estimated by the particle
identification (PID) process. Electrons or gamma-rays create fuzzy Cherenkov
rings by electromagnetic showering (e-like rings). On the other hand, Cherenkov
rings which have sharp ring edges are created by muons or charged-pions (µ-likerings). Moreover, the Cherenkov angle of e-like rings are almost constant at
42◦ , but that of µ-like rings can vary depending on their momenta and energy
losses. PID can classify the particle type into e-like or µ-like by using thesedifferences in the ring patterns and the angles.
PID likelihoods are built with two particle type hypotheses for each ring andused for PID. The likelihood consists of the observed and expected charge. The
expected charge is evaluated from the SK MC simulations and given as follows.

qexp,i (e) = αe Qexp (pe , θi )

R
ri

3/2

51

e−ri /L f (Θi ) + qscat,i(6.4)

Number of Events

CHAPTER 6. EVENT RECONSTRUCTION

SK-IV Data
ν eCC w/o π
ν eCC w/ 0π0oth.
ν eCC w/ π0
ν µCC w/o πν µCC w/ 0π0 oth.
ν µCC w/ π0
NC

500
400
300
200
100
0
-15

-10

-5

0

5
10
15
Ring counting likelihood

Figure 6.3: Ring counting likelihood distribution for the FC sub-GeV sample ofthe data (black points) and the atmospheric neutrino MC (stacked histograms)in SK-IV. The events whose reconstructed vertices are within the region 2 m
inside from the nearest wall without tagged Michel electrons are selected. The
events in the positive (negative) side are identified as multi- (single-) ring events.

52

CHAPTER 6. EVENT RECONSTRUCTION




2
α
sin
θ
µ
x
i
 + qknock,i  e−ri /L f (Θi )qexp,i (µ) =   dθ 
ri sin θxi + ri dx
x=x

(6.5)

i

+ qscat,iwhere
αe , αµ : normalization factor
Qexp (pe , θi ): expected charge distribution depending on the electron momentum and the angle θi for the i-th PMT
R: virtual sphere radius (16.9 m)
ri : distance from the vertex to the i-th PMT
L: light attenuation length in waterf (Θi ): angular acceptance as a function of the photon injection angle Θi for
the i-th PMT
qscat,i : expected charge due to scattering light for the i-th PMT
θxi : Cherenkov angle of a muon track length at x
qknock,i : expected charge due to knock-on electrons for the i-th PMT
Figure 6.4 shows the PID likelihood distributions for the FC sub-GeV single-ringand multi-ring sample of the data and the atmospheric neutrino MC.

6.2.4

Momentum Determination

The momentum of each particle is evaluated by the observed charge inside theCherenkov cone with a half angle of 70◦ . The observed charge of the i-th PMT
for the n-th ring is separated based on the expected charge distributions fromeach ring, described as follows;
qexp,i,n.
0n0 qexp,i,nqobs,i,n = qobs,i P(6.6)

By correcting the light attenuation and angular acceptance of the PMTs, thetotal charge for each ring (RTOT) is calculated as follows;




X
X
cos ΘiGMC qobs,i,n eri /Lα−Si (6.7)
RT OT =GDataf (Θi )◦◦
θi,n <70 ,−50 ns<ti <250 ns

θi <70where
GMC , GData : relative PMT gain parameter for the data and the MC simulation
α: normalization factor
θi,n : angle between the i-th PMT and the n-th ring direction
ti : TOF subtracted hit timing of the i-th PMT
Si : expected charge of the i-th PMT due to scattering light.
RTOT is converted to the corresponding momentum depends on the PID of thering (e-like or µ-like) using conversion tables obtained from MC.

53

Number of Events

CHAPTER 6. EVENT RECONSTRUCTION

SK-IV Data
ν eCC w/o π
ν eCC w/ 0π0oth.
ν eCC w/ π0
ν µCC w/o πν µCC w/ 0π0 oth.
ν µCC w/ π0
NC

400
350
300
250
200
150
100
50

Number of Events

0
-10

-8

-6

-4

-2

0

2

4

6
8
10
PID likelihood

SK-IV Data
ν eCC w/o π
ν eCC w/ 0π0oth.
ν eCC w/ π0
ν µCC w/o πν µCC w/ 0π0 oth.
ν µCC w/ π0
NC

350
300
250
200
150
100
50
0
-300

-200

-100

0
100
200
300
PID likelihood for all rings

Figure 6.4: PID likelihood distributions for the FC sub-GeV single-ring (top)and multi-ring (bottom) sample of the data (black points) and the atmospheric
neutrino MC (stacked histograms) in SK-IV. Here, the events whose reconstructed vertices are within the region 2 m inside from the nearest wall are
selected. The events in the positive (negative) side are identified as muon(electron-) like events.

54

CHAPTER 6. EVENT RECONSTRUCTION

6.2.5

Michel Electron Tagging

Michel electrons are an indicator of the existence of a muon. The tagging ofMichel electrons is conducted based on the number of hits in a cluster, the total
charge and the goodness of the vertex fit. Due to an impedance mismatch seenin the ATM boards, the time interval between 800 ns and 1,200 ns from theprimary event timing is excluded from the tagging for SK-I to SK-III. But, thattime interval is available for SK-IV since the mismatch was solved in the QBEE
boards. Therefore, the tagging efficiency for SK-IV is higher than the otherdetector periods. The tagging efficiency of Michel electrons is 80% for µ+ and
63% for µ− for SK-I to SK-III, 96% for µ+ and 80% for µ− for SK-IV.

6.3

fiTQunFirst of all, we discuss the likelihood function since fiTQun is kind of maximum
likelihood fitter. Then, the reconstruction procedure of fiTQun is shown as thefollowing steps.
1. Vertex pre-fitter
2. Hit clustering
3. Single-ring fitter
4. Multi-ring fitter
The document [97] can provide much more information about fiTQun.

6.3.1

Likelihood Function

The event reconstruction is conducted by maximizing the following likelihood
function,L(x) =unhit
Y

Pj (unhit|x)

j

hit
Y
{1 − Pi (unhit|x)}fq (qi |x)ft (ti |x).

(6.8)

i

where x is a track parameter set of { vertex position (x, y, z), time (t), zenith
angle (θ), azimuth angle (φ), momentum (p)}, P (unhit|x) (1 − P (unhit|x))
indicates the unhit (hit) probability, fq (qi |x) and ft (qi |x) are the probability
density functions (PDFs) of the observed charge and time, respectively. The hit(unhit) probability function is multiplied for all hit (unhit) PMTs. This likelihood function can be simplified by introducing a predicted charge (µ) insteadof x.
L(x) =unhitYj

Pj (unhit|µj )hitY{1 − Pi (unhit|µi )}fq (qi |µi )ft (ti |x).
i

55

(6.9)

CHAPTER 6. EVENT RECONSTRUCTION

Then, the unhit probability function and charge PDF become independent of theCherenkov photon emission profile, these are determined by the characteristicsof PMTs and its readout electronics. The time likelihood also uses the predictedcharge as described below. For the likelihood calculation, the predicted charge
for each PMT for a given track parameter set is evaluated at first, and thelikelihood is then evaluated based on the predicted charge and the observed
charge and time information.

Predicted Charge
For a given single particle hypothesis x, the predicted charge of each PMT isestimated for the direct (µdir ) and indirect light (µsct ), separately. In the case ofa multi-particle hypothesis, the predicted charge is independently calculated forthe n-th particle hypothesis and summed for all particle hypotheses as follows.
X
sct
dark
(µdir
µi =
(6.10)
i,n + µi,n ) + µnwhere µdark indicates the contribution from the PMT dark rate taken fromSKDETSIM (5.7195 kHz for SK-IV) for a time window which is described in
Section 6.3.3. For a single particle hypothesis, µdir and µsct are given as anintegral along the particle track length s from its initial position;Zµdir = Φ(p) ds g(p, s, cos θ)Ω(R)T (R)(η)(6.11)
and
sct

µ

Z
= Φ(p)ds1
ρ(p, s)Ω(R)T (R)(η)A(s)4π(6.12)

where the variables θ, R, η indicate the relative orientation of the particle and
the PMT which are all function of s as shown in Figure 6.5. Φ(p) is a normalization factor which is proportional to the average total number of photons forthe particle with initial momentum p. In the following sentences, we describe
the following factors: the Cherenkov emission profile g, the PMT solid angle
factor Ω, the light transmission factor T , the PMT angular acceptance , thephoto emission fraction ρ and the light scattering factor A.
First of all, we need to know how much Cherenkov photons are emitted along
the particle track s. The Cherenkov emission profile g(p, s, cos θ) is the numberof emitted photons per unit track length per unit solid angle at angle θ withrespect to the particle direction for the particle with initial momentum p. The
normalization of g is given as;Zg(p, s, cos θ) dsdΩ = 1.
(6.13)
Figure 6.6 and 6.7 show the Cherenkov profiles for electrons and muons which are
generated by injecting a single particle with various momenta in SKDETSIM.
56

CHAPTER 6. EVENT RECONSTRUCTION

s

R

PMT
Figure 6.5: Schematic diagram of the relation between the variables for thepredicted charge. The vertex position of the particle is shown as the white dot.
Taken from [97].

The photon emission angles for electrons are almost constant and have the broad
distributions, whereas that for muons are depending on the momentum and havethe sharper distribution. These differences enable us to distinguish the particle
type.
Next, we need to know how Cherenkov photons are received by the PMT.
For that purpose, the solid angle Ω(R) for the PMT is roughly evaluated asfollows.
Ω(R) =

πa2
R2 + a2(6.14)

where a is the radius of the ID PMT (0.254 m). This formula is valid at distanceR > 1 m. While traveling the distance from the emitting point to the PMT,
photons can be absorbed or scattered in water. Such effects are considered by
the light transmission factor T (R).
T (R) = exp(−R/Latt )(6.15)

where Latt is the light attenuation length in water used in SKDETSIM (74.9646 m
for SK-IV). The PMT angular acceptance (η) should be taken into account.
Figure 6.8 shows the angular acceptance as a function of the angle η between the
PMT normal and the direction of the particle position viewed from the PMT,
which is obtained from SKDETSIM. The fitted polynomial function is used in
57

0.8

400

300

0.6

300

200

0.4

200

100

0.2

100

0

0

0

0.6

0.8

e 50MeV/cs (cm)s (cm)

400

1
cosθ

400

0.25
0.2

300

0.15

200

0.5
0.4
0.3
0.2
0.1
0.6

0.8

e 200MeV/cs (cm)s (cm)

CHAPTER 6. EVENT RECONSTRUCTION

1
cosθ

400
300

0.1

200

0.1
100
0

0.05

0.05
0.6

e 500MeV/c

0.8

1
cosθ

0

0

100
0

0.6

e 1000MeV/c

0.8

1
cosθ

0

Figure 6.6: Cherenkov emission profiles for electrons at various initial momenta.
In all the figures, the photon emission angles have board distributions and arepeaked at around cos θ = 0.75 which corresponds to the Cherenkov angle for aparticle with β = 1 in water. Taken from [97].

58

400

s (cm)s (cm)

CHAPTER 6. EVENT RECONSTRUCTION

4

300

400
0.8

300

0.6

3
200

0.6

0.8

s (cm)µ 160MeV/c

1
cosθ

400

1

100

0

0

0.3
0.25

300

0.4
0.2
0.6

0.8

µ 300MeV/cs (cm)

100
0

200

2

1
cosθ

400
0.15
300

0.2
200

0.15
0.1

100

0

0.1

200

0.05

100

0.05
0

0.6

µ 600MeV/c

0.8

1
cosθ

0

0

0.6

µ 1000MeV/c

0.8

1
cosθ

0

Figure 6.7: Cherenkov emission profiles for muons at various initial momenta.
In all the figures, the photon emission angles have shaper distributions andstronger dependence on momentum compared to electrons. Taken from [97].

59CHAPTER 6. EVENT RECONSTRUCTION

2500
2000
1500
1000
500
0
0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
cosη

Figure 6.8: PMT angular acceptance as a function of cos η obtained fromSKDETSIM. The vertical axis is in arbitrary unit. The solid line shows a fitted
polynomial function, and it is normalized to be 1 at η = 0. Taken from [97].

the predicted charge calculation and its normalization condition is (η = 0) = 1.
Finally, the indirect light caused by light scattering and reflection by detector
components should be considered. This contribution is taken into account inA(s) as,A(s) = A(xPMT , zVTX , RVTX , ϕ, θ, φ) ≡dµsct.
dµiso,dir(6.16)µiso,dir is the predicted charge from isotropic light source along with the same
particle track and total light intensity as µdir , which is obtained by removing
A(s) from Equation 6.12 and defining ρ(p, s) as follows.
Zρ(p, s) ≡ g(p, s, cos θ) dΩ(6.17)In the indirect light evaluation, we assume that the Cherenkov angle does notdepend on momentum and is fixed at the value of β = 1. By taking the ratio
in Equation 6.16, the momentum dependence is factored out and then A(s)becomes purely a function of geometrical variables between the particle track
and the PMT in concern. Since the detector is cylindrical, A(s) can be describedby six-parameters as shown in Figure 6.9. In order to obtain A(s), we randomly
generate 3 MeV/c electrons within the detector while multiple scattering is60

CHAPTER 6. EVENT RECONSTRUCTION

Rvtx
φ
zvtx
(θ,ϕ)

zPMT

Figure 6.9: Schematic diagram of the relation between the variables for the
scattering table in case of the PMTs on the barrel. For the PMTs on the top or
bottom wall, xPMT in Equation 6.16 is replaced to the PMT’s distance RPMT
from the vertical axis at the center of the tank. The vertex position of theparticle is shown as the white dot. Taken from [97].

turned off. Such electrons can be regarded as point light sources with β = 1.
Then, six-dimensional histograms for both µdir and µiso,dir are made by counting
the direct and indirect light arriving at the PMTs, respectively. At each binin the six dimensional parameter space, the numerator dµsct is obtained fromthe histogram for the indirect light, while the denominator dµiso,dir is calculated
by averaging the histogram for the direct light over the particle direction. The
ratio is linearly interpolated in terms of the six parameters. The scattering table
A(s) is thus obtained.
The integral along the particle track for the predicted charge takes too muchtime to calculate it. In order to avoid such a problem, we employ an approximation of the predicted charge. The photon acceptance factor J(s) is defined
as follows.
J(s) ≡ Ω(R)T (R)(η) ≈ j0 + j1 s + j2 s2(6.18)As shown in Figure 6.10, J(s) is well approximated as the parabola function.
The coefficients j0 , j1 and j2 are determined by evaluating J(s) at three points
along the particle track; the initial particle position, the position at which 90 %of Cherenkov light is emitted, and their mid-position. Then, the formulation of
the predicted charge for the direct light (Equation 6.11) becomes simpler after

61

CHAPTER 6. EVENT RECONSTRUCTION

the approximation.
dir

µ

Z
ds g(s)J(s) ≈ Φ(p)(I0 j0 + I1 j1 + I2 j2 )= Φ(p)(6.19)

where
Z
In ≡

ds g(s)sn .

(6.20)In is obtained as a function of three parameters; p, R and θ. After integrating
Equation 6.20 at discrete those three parameter bins, In is fitted as a functionof p at each R, θ bin for smoothing the likelihood surface. Finally, In is linearlyinterpolated in terms of R and θ. In a similar manner to the predicted chargeof the direct light, an approximation is introduced for that of the indirect light
as follows.
J(s)A(s) ≈ k0 + k1 s + k2 s2 .

(6.21)The predicted charge for the indirect light in Equation 6.12 reduces as follows.
Z
1
1
sct
ρ(s)J(s)A(s) = Φ(p) (k0 + K1 k1 + K2 k2 ),µ = Φ(p) ds(6.22)
4π
4π
where
Z
Kn ≡

ds ρ(s)sn .

(6.23)

Because Kn is a function of p only, it is fitted as a function of p for smoothing
the likelihood surface. The coefficients kn are determined in the same way as
jn .

Unhit Probability
Photoelectron counting by PMTs obeys a Poisson distribution with mean µcalculated as the predicted charge. Thus, the probability of no photoelectronis given as e−µ . Due to the threshold of the PMT signal, there are possibilities
that a photoelectron does not surpass the threshold and make a hit. Correcting
up to the third order of µ for such threshold effect, the unhit probability is
estimated as;
P (unhit|µ) = (1 + a1 µ + a2 µ2 + a3 µ3 )e−µ(6.24)

where the coefficients an are obtained from SKDETSIM. Figure 6.11 shows agood agreement between Equation 6.24 and the true values taken from SKDETSIM.

62

CHAPTER 6. EVENT RECONSTRUCTION

J(s)0.0012Full

0.00118

Parabolic0.001160.001140.001120.00110.001080.00106

0

20

40

60

80

100

120

140

160
s (cm)

Figure 6.10: Photon acceptance factors as a function of s for the initial conditions R = 500 cm, θ = 90◦ , η = 0◦ . The black and red lines correspond to J(s)and the parabola function (j0 + j1 s + j2 s2 ), respectively. The coefficients jn aredetermined by evaluating J(s) at the three points indicated by the circles. J(s)can be well-reproduced by the parabola function. Taken from [97].

63P(unhit| µ )CHAPTER 6. EVENT RECONSTRUCTION1
10-1

w/o corr.
w/ 3rd order corr.
MC true

10-2
10-3
10-4
10-5
0

2

4

6

8

10

12
µ (p.e.)Figure 6.11: Unhit probabilities as a function of the predicted charge. The dots
indicate the values taken from SKDETSIM. The red and blue lines are for withand without the 3rd order correction of the PMT threshold effect, respectively.
The corrected function reproduces the MC true values. Taken from [97].

Charge Likelihood
The charge likelihood fq (q|µ) is the probability to observe the charge q of the
PMT when its predicted charge is µ. In order to obtain fq (q|µ), normalized
observed charge distributions are prepared by generating photoelectrons obeying
a Poisson distribution with discrete values of µ at the PMTs in SKDETSIM asshown in Figure 6.12. Then, those distributions are fitted by a polynomial as afunction of µ at each fixed value of q. The final polynomial is obtained after alinear interpolation of the fit parameters in terms of the observed charge q.

Time Likelihood
The time likelihood ft (ti |x) depends on the event hypothesis x and the positionof the PMT. However, it is impractical to completely consider such complex dependencies. In order to approximate ft (ti |x), we assume all Cherenkov photons
are emitted at the midpoint of the particle track, and then the time likelihood
can be expressed by the residual hit time given as,itres
i = ti − t − smid /c − |RPMT − x − smid d|/(c/n)(6.25)iwhere ti , RPMT
are the raw hit time and position of the i-th PMT, x, t are thevertex position and time, d is the particle direction and smid is the half of the64CHAPTER 6. EVENT RECONSTRUCTION

0.7

0.7

0.6

µ = 0.30 p.e.

0.5

0.6

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0
0

5

10

15

0
0

20
25
q (p.e.)0.7µ = 1.50 p.e.

0.5

10

0.6

15

20
25
q (p.e.)

0.4

0.3

0.3

0.2

0.2

0.1

0.1
5

10

15

0
0

20
25
q (p.e.)0.7µ = 3.00 p.e.

0.5

0.4

5

10

15

20
25
q (p.e.)0.7

0.6µ = 6.00 p.e.

0.5

0.6

0.4

0.3

0.3

0.2

0.2

0.1

0.1
5

10

15

µ =12.00 p.e.

0.5

0.4

0
0

5

0.7

0.6

0
0

µ = 0.70 p.e.

0.5

0
0

20
25
q (p.e.)

5

10

15

20
25
q (p.e.)Figure 6.12: Normalized charge likelihoods as a function of observed charge
at various predicted charges. The dots are taken from SKDETSIM and are
well-fitted by the fitted function indicated by the solid line. Taken from [97].

65CHAPTER 6. EVENT RECONSTRUCTION

particle track length. By using tres
i , time likelihoods for the direct and indirect
light are calculated, and the final likelihood is obtained by merging those two
likelihoods as described below.
For the direct light, we impose the assumption that the time likelihood
ftdir (tres ) follows a Gaussian and its parameters depend only on µdir and p. The
width of the residual time distribution decreases when arriving photons increasesince a hit is recorded by the first photon arriving at a PMT. That width
becomes wider for a long track event since the assumption in Equation 6.25
becomes less accurate. In order to obtain the tres distribution, particles are
generated at various fixed momenta in SKDETSIM. At each momentum, thetwo dimensional histogram of tres − log µdir is filled for every hit caused by adirect photon. Then, the mean and σ of the tres distribution are obtained byfitting a Gaussian at each log µdir . After smoothing those two parameters interms of log µdir and p, the final ftdir (tres ) is obtained. Figure 6.13 and 6.14
show the ftdir (tres ) distributions for electrons and muons, respectively.
In contrast to the direct light time likelihood, we do not assume any dependence of µsct and p for the indirect light time likelihood, and it is defined as anempirical function given by,r πexp(τ 2 /2σ 2 )(τ < 0)sct resσ + 2γ ×, (6.26)ft (ti ) = 1/
(τ /γ + 1) exp(−τ /γ) (τ > 0)2where τ = tresi − 5 ns, σ = 8 ns and γ = 25 ns. In order to reproduce the effectof reflected light, this function has a long tail on the right side.
In order to obtain the total time likelihood for a single particle hypothesis,
the direct and indirect light time likelihoods should be merged assuming thata hit is produced by an indirect photon if there is no direct photon, which is
given as,dir res
sct resft (tresi ) = ωft (ti ) + (1 − ω)ft (ti ),(6.27)

wheredir

1 − e−µi(6.28)sct .
dir
1 − e−µi e−µi
For a multi-particle hypothesis, it is possible to extend Equation 6.27. At
first, tres
i ’s of each particle hypothesis are calculated separately. Then, theseparticles are ordered by tres
i . The total time likelihood for a multi-particle
hypothesis is defined as follows.
ω≡ft (tres
i )=

1
1−
×

sct
−(µdiri,j +µi,j )
j=1 e

Qall

all
X

j−1
Y

"

j=1

+

all
Y
k=1

!
−µdir
i,k

e

dir

dir res(1 − e−µi,j )ft,j(ti )(6.29)

k=1

!
−µdir
i,k

e

j−1
Y

!
−µsct
i,k

e

k=1

66

#
−µsct
i,j

(1 − esct res
)ft,j
(ti )

CHAPTER 6. EVENT RECONSTRUCTION

0.15

ep=300MeV/cµ = 0.26 p.e.

0.1

0.05

-20

0

20

0

40
t [ns]
res

0.2

e-

0.15

p=300MeV/cµ = 3.80 p.e.

-20

0

20

40
tres [ns]ep=300MeV/cµ =15.14 p.e.

0.3
0.2

0.1

0.1

0.05
0

p=300MeV/cµ = 1.05 p.e.

0.1

0.05

0

e-

-20

0

20

0

40
t [ns]
res

-20

0

20

40
tres [ns]Figure 6.13: Normalized direct light time likelihoods for 300 MeV/c electrons
as a function of the residual time at various predicted charges for the direct
light. The dots are taken from SKDETSIM. The red curves are the Gaussianfunction parameterized by µdir and p. Taken from [97].

67CHAPTER 6. EVENT RECONSTRUCTION

0.15µp=450MeV/cµ = 0.26 p.e.

0.1

0.05

-20

0

20

0

40
t [ns]
res

0.2

µ-

0.15

p=450MeV/cµ = 3.80 p.e.

-20

0

20

40
tres [ns]µp=450MeV/cµ =15.14 p.e.

0.3
0.2

0.1

0.1

0.05
0

p=450MeV/cµ = 1.05 p.e.

0.1

0.05

0

µ-

-20

0

20

0

40
t [ns]
res

-20

0

20

40
tres [ns]Figure 6.14: Normalized direct light time likelihoods for 450 MeV/c muons as
a function of the residual time at various predicted charges for the direct light.
The dots are taken from SKDETSIM. The red curves are the Gaussian function
parameterized by µdir and p. Taken from [97].

68CHAPTER 6. EVENT RECONSTRUCTION

6.3.2

Vertex Pre-fitter

The minimization of − ln L(x) starts with a single vertex position as an initial
seed in order to avoid stack at a local minimum. For that purpose, the vertex
pre-fitter is firstly performed in the likelihood calculation by assuming a point
light source, which is estimated by maximizing the following goodness function.
hitXiexp(−(Tres/σ)2 /2),(6.30)ii
Tres= ti − t − |RPMT− x|/(c/n).

(6.31)G(x, t) =

i

where

i
Tres
is the residual hit time of the i-th PMT. A coarse grid scan in the space of{x, t} is iteratively done for the goodness while shrinking the grid size and σ.
After the grid scan, the final seed vertex is obtained by minimizing of −G(x, t)with σ = 4 ns.

6.3.3

Hit Clustering

Since the time window for an event is order 10 µs around an event trigger, there
might be multiple hit clusters. For example, a Michel electron produced by amuon as a primary event are counted as “subevent” since their time interval
is enough long, typically 2 µs. The hit clustering algorithm is thus performed
before going to the maximum likelihood reconstruction.
First, the goodness G(x, t) discussed in Equation (6.30) is evaluated every
8 ns step with an optimized σ (=6.3 ns) while fixing the vertex position givenby the vertex pre-fitter. The evaluation is started at ∼300 ns before the trigger
using all hits in the time window between −200 ns and 15,000 ns around thetrigger. An example of the goodness distribution as a function of t for an eventwith a parent muon and a Michel electron shown in Figure 6.15. The goodness
peak is identified by using a goodness threshold curve function. Two peaks in
Figure 6.15 are consistent with true particle times.
After the peaks are searched, a time window is set around each peak, containing the hits belong to the peak. The hits outside the time window are considered
as unhit in the likelihood calculation. The residual time is re-evaluated for theipeak time, and then the time window is defined as −180 ns < Tres
< 800 ns. Iftime windows are overlapped, they are merged into a single window.
Finally, the vertex pre-fitter and the peak search are conducted again for a
precise peak finding. The final time windows are fixed for each remaining peak,
and the maximum likelihood reconstruction is done in each of them. Peakdetection efficiencies of Michel electrons from sub-GeV stopping muons for thedata and MC are 88.41 ± 0.08% and 87.81 ± 0.15%, respectively.

69

CHAPTER 6. EVENT RECONSTRUCTION

800
700
600
500
400
300
200
100
0

700

800

900 1000 1100 1200 1300 1400
t (ns)

Figure 6.15: Goodness distribution as a function of hit time for an event witha parent muon and a Michel electron. The dots show the scanned goodness
points. The blue and green curves are the goodness threshold functions which
are used to determine the event time indicated by the red lines. The green
dotted lines are for the true time for each particle. Taken from [97].

6.3.4

Single-Ring FitterFor each subevent described in Section 6.3.3, the single-ring fitter is performedby maximizing the likelihood under a single particle hypothesis. There are threeparticle hypotheses which fiTQun can handle: e, µ and π + . The single-ring fitterfor e and µ hypotheses reconstructs seven parameters {x, t, θ, φ, p}. In addition
to these parameters, the energy loss parameter Eloss is fitted for π + hypothesis.
For the subevent which has multiple event peaks in the time window described above, the single-ring fitters are performed on each peak using all thehits related to each peak without overlapping. The hits for the peak after theiiprimary one are selected in a range of −30 ns < Tres
< 60 ns where Tresisthe residual time calculated for the peak. The single-ring fitters are performed
using those hits while the other hits are regarded as unhits. For the primarypeak, the hits which are not selected by the other peaks are used to the fit.
First of all, the single electron fitter is conducted. As for the initial values ofthe vertex {x, t} to the fit, the results of the vertex pre-fitter are used. While
fixing the vertex and the momentum estimated by the total observed charge, theinitial value of the direction {θ, φ} is set by scanning the likelihood at 400 pointsof the direction which are equally spaced on the unit sphere. The initial value ofthe momentum p is determined by another likelihood scanning fixing the otherparameters. Using these initial seven values and electron hypothesis, the single
electron track parameters are then simultaneously fitted by minimizing − ln L.

70

CHAPTER 6. EVENT RECONSTRUCTION

Next, the single muon fitter is performed. By using the fitted track parameters of the single electron fitter except the momentum, the initial value ofthe momentum is obtained by a likelihood scan. After that, the simultaneous
minimization of − ln L is done assuming muon hypothesis.
The last single-ring fitter is for the π + hypothesis. Since the mass of π + s issimilar to that of muons, π + s usually have similar Cherenkov emission profilesof muons. However, a π + can interact hadronically with the nuclei in water,
and then change its direction or be absorbed. In such a case, a π + produces
two rings; the ring before the hadron interaction is called the “upstream-track”
ring and the other ring after that is called the “downstream-track” ring. The
single π + fitter reconstructs this upstream-track ring using the single muon fit
result as a seed. The predicted charge of the upstream-track is obtained bysubtracting the predicted charge of the downstream-track from the predictedcharge of an imaginary full track without any hadron interaction. Here, theadditional fit parameter Eloss is introduced which is defined as the kinetic energy
deposit in the upstream-track. Using Eloss , the vertex of the downstream-track
is calculated by assuming constant kinetic energy loss per unit track length. For
the time likelihood, the residual time in Equation 6.25 is estimated using the
midpoint of the upstream-track. All eight parameters of the upstream-track aresimultaneously fitted by minimizing of − ln L.
By taking the ratio of resulting log likelihoods, the particle type of thesing-ring can be identified. Figure 6.16 shows ln(Le /Lµ ) as a function of thereconstructed momentum for the single electron fit in the FC true-fiducial νeand νµ CCQE in the atmospheric neutrino MC. Both ln(Le /Lµ ) for νe and νµhave the opposite dependence on the momentum. Therefore, electron events
are separated by the following cut which is decided by taking into account themisidentification rate for both electrons and muons.
ln(Le /Lµ ) > −10.

(6.32)

Although we do not distinguish π + and µ in this analysis, it is worth to look atln(Lπ+ /Lµ ) as a function of Eloss in the T2K accelerator neutrino MC as shownin Figure 6.17. The separation is not clear compared to ln(Le /Lµ ), but somefraction of π + events can be distinguished from muons by making a thresholdline. In the left figure, ln(Lπ+ /Lµ ) is distributed near zero, which means thatmuons are fitted almost equally well by the µ and π + hypotheses. Since both µand π + have similar Cherenkov emission profiles, muon rings are fitted well by
the upstream-track π + hypothesis while increasing Eloss to the point where themomentum at the assumed upstream-track point drops below the Cherenkov
threshold.
Examples of the data-MC agreement related to the PID likelihood for the FC
sub-GeV single-ring and multi-ring sample in SK-IV are shown in Figure 6.18.
Our atmospheric neutrino MC agrees well with the data regarding PID performance. The mis-PID rate in single-ring events is evaluated by using the
stopping muon sample described in Section 6.6. As shown in Figure 6.19, themis-PID rates for both muons and electrons in the data (MC) are 0.10 ± 0.03%(0.15 ± 0.03%) and 0.61 ± 0.08% (0.48 ± 0.04%), respectively.
71

6

2000

ln(Le/Lµ )ln(Le/Lµ )CHAPTER 6. EVENT RECONSTRUCTION

4

0

30

2000

20

0

10

2
-2000
0

-2000
200 400 600 800 1000
Rec. e momentum (MeV/c)00

200 400 600 800 1000
Rec. e momentum (MeV/c)

0

Figure 6.16: Log likelihood ratio ln(Le /Lµ ) as a function of the reconstructed
momentum of the sing-ring electron fit for the FC true-fiducial νe (νµ ) CCQE
in the atmospheric neutrino MC on the left (right) figure. The magenta lines
are the cut criteria for the e/µ separation. Taken from [97].

Figure 6.17: Log likelihood ratio ln(Lπ+ /Lµ ) as a function of the energy loss
for the FC true-fiducial νµ CCQE (NC π + ) in the T2K accelerator neutrino
MC on the left (right) figure. The black lines are the cut criteria for the µ/π +separation. Taken from [97].

72

Number of Events

CHAPTER 6. EVENT RECONSTRUCTION

800

SK-IV Data
ν eCC w/o π
ν eCC w/ 0π0oth.
ν eCC w/ π0
ν µCC w/o πν µCC w/ 0π0 oth.
ν µCC w/ π0
NC

700
600
500
400
300
200
100

Number of Events

0
-3000

-2000

-1000

0

1000

2000
3000
ln(Le/Lµ)

SK-IV Data
ν eCC w/o π
ν eCC w/ 0π0oth.
ν eCC w/ π0
ν µCC w/o πν µCC w/ 0π0 oth.
ν µCC w/ π0
NC

250
200
150
100
50
0
-50 -40 -30 -20 -10

0

10

20

30 40 50ln(Lall-e/ Lothers)

Figure 6.18: PID likelihood distributions for the FC sub-GeV single-ring (top)and multi-ring (bottom) sample of the data (black points) and the atmospheric
neutrino MC (stacked histograms) in SK-IV. Here, the events whose reconstructed vertices are within the region 1.5 m inside from the nearest wall are
selected. In the top figure, the events in the above (below) -10 are identified aselectron- (muon-) like events. The horizontal axis in the bottom figure is for the
likelihood ratio between an all e-like hypothesis and the best fit among otherhypotheses for two or three ring events.

73

300
250

Number of Events

Number of Events

CHAPTER 6. EVENT RECONSTRUCTION

Mis-PID rate
Data: 0.10±0.03%
MC: 0.15±0.03%

200
150

600

Mis-PID rate
Data: 0.61±0.08%
MC: 0.48±0.04%

500
400
300

100

200

50

100

0
-25000 -20000 -15000 -10000 -5000

0

0
-600

5000 10000ln(Le/Lµ) + 10

-400

-200

0

200

400
600
ln(Le/Lµ) + 10

Figure 6.19: Distributions of the log likelihood ratio ln(Le /Lµ ) + 10 for bothmuons (left) and electrons (right) in the stopping muon sample of the data
(black points) and the MC (red histogram).

6.3.5

Multi-Ring Fitter

The reconstruction of multi-particle events is crucial for proton decay and atmospheric neutrino analyses. FiTQun searches up to six rings only in the firstsubevent since subsequent subevents usually do not contain multi-particles.
The ring counting of fiTQun is conducted in a brute force way as shownin Figure 6.20. Electrons and gamma-rays are reconstructed as the electron
ring, and muons and π + s are reconstructed as the upstream-track π + ring since
they can be fitted by the π + hypothesis. In this figure, the single electron
fit result is used as the first ring and the brute force ring search is done, but
this search is also repeated for the single upstream-track π + fit. First of all,the direction of an additional second ring is roughly evaluated by scanning the
likelihood at 400 points of the direction, assuming the same vertex of the singlering fit result and the second ring momentum is assumed to be 50 MeV/c.
Then, the second ring momentum is fitted by minimizing − ln L while all theother parameters are fixed. After that, the direction and the momentum of
the second ring are simultaneously fitted fixing all the other parameters. When
there are additional rings in an event, the momentum of the single ring fit tends
to be overestimated due to the existence of the charge from the other rings.
Therefore, the simultaneous fit of the momenta of the two rings is performed
while all the other parameters are fixed. Finally, all the parameters of the tworings including the common vertex are then fitted simultaneously. In the caseof the π + ring, Eloss is also fitted simultaneously with the other parameters.
In order to check whether the new ring is true ring or fake ring, the log likelihood ratio of the single-ring fit to the better of the two-ring fit, ln(L2Rex /L1Re )is evaluated. Figure 6.21 shows the square root of ln(L2Rex /L1Re ) distribution
for FC sub-GeV atmospheric neutrino MC events with no Michel electrons. True
single electron events such as νe CC events are concentrated to a lower region
while multi-ring events such as NC events have higher values. The new ringwhich passed the following criterion is identified as the true ring.
74

CHAPTER 6. EVENT RECONSTRUCTION

e

1st

ee

2nd

eee

3rd
4th

e

eeee

ee
eee

ee e

e e
ee

e ee

e e

e
e

e

e

Figure 6.20: Schematic diagram of the ring finding in the multi-ring fitter. The
fifth and the sixth ring are assumed to be electrons. In this diagram, the firstring is assuming an electron ring. The same procedure is repeated for the caseof assuming π + as the first ring.

• ln(L2R /L1R ) > 87.5, if the second ring is the electron-ring.
or,• ln(L2R /L1R ) > 140, if the second ring is the π + -ring.
For the (n + 1)-th ring, the criterion becomes as,• ln(L(n+1)R /LnR ) > 70.
Following the second ring search, additional ring search is repeated until
either a fake ring is found or up to six rings in a similar manner. But, in orderto save computation time, the fifth and the sixth ring are fitted by the electron
hypothesis only, and the simultaneous fit of all rings is conducted only for the
momentum of them after the two-ring fit. Figure 6.21 also shows the squareroot of ln(L3Reex /L2Ree ) distribution for FC sub-GeV atmospheric neutrino MC
events with no Michel electrons. It can be seen that both νe CC with π 0 eventsand NC events are dominant interaction modes above the threshold. The multiring hypothesis with the smallest − ln L is used in a further precise fit as the“seeding multi-ring hypothesis”.
Once the seeding multi-ring hypothesis is obtained, the rings in that hypothesis are reordered by the visible energy. For the most energetic ring, if theangle between the ring and a lower energy ring is smaller than 20◦ , they are
merged into a single ring by adding the visible energy. This angle thresholdis decided to reduce an over ring-counting seen in p → e+ π 0 analysis at thedevelopment stage as shown in Figure 6.22. Figure 6.23 shows the number ofrings distributions for free proton events in the p → e+ π 0 analysis before andafter the ring correction change. After the ring merge procedure, the most energetic ring is fitted including the vertex by three particle hypotheses (e, µ, π + )while the other rings are fixed. Then, the ring type of the most energetic ring is75

Number of Events

CHAPTER 6. EVENT RECONSTRUCTION

700

SK-IV Data
ν eCC w/o π
ν eCC w/ 0π0oth.
ν eCC w/ π0
ν µCC w/o πν µCC w/ 0π0 oth.
ν µCC w/ π0
NC

600
500
400
300
200
100

Number of Events

0

0

10

20

250

30
40
50
sgn × ln(L2Rex/ L1Re)

SK-IV Data
ν eCC w/o π
ν eCC w/ 0π0oth.
ν eCC w/ π0
ν µCC w/o πν µCC w/ 0π0 oth.
ν µCC w/ π0
NC

200
150
100
50
0

0

10

20

30
40
50
sgn × ln(L3Reex / L2Ree)

Figure 6.21: Square root of ln(L2Rex /L1Re ) (top) and ln(L3Reex /L2Ree ) (bottom) distributions for the FC sub-GeV sample without tagged Michel electrons.
The dots show the SK-IV data. Stacked histograms are the atmospheric neutrino MC, which are categorized by final states. The events whose reconstructedvertices are within the region 1.5 m inside from the nearest wall are selected.

76Number of Events

CHAPTER 6. EVENT RECONSTRUCTION

25

20

15

10

5

0
0

20

40

60

80

100

120
140
160
180
Minimum angle (deg)Figure 6.22: Minimum opening angle between each ring only for free proton
events in the p → e+ π 0 MC. The black, red, green and blue histograms correspond to 2, 3, 4 and 5-ring events, respectively. By introducing the opening
angle cut at 20◦ , we can merge the fake rings.

determined whether it is a showering (e and γ) or a non-showering (µ and π + )ring. If the following criterion is satisfied for the ring, then it is marked as ashower ring and the original hypothesis is replaced by the selected hypothesis.
ln(Le /Lπ+ ) > −10.

(6.33)The simultaneous fit and particle identification procedures as above are sequentially performed for all the remaining rings in the descending order of the visibleenergy.
In the multi-ring fitter, the e/π + separation is only used. However, in thecase of the ring is judged as a non-showering ring, the ring is always assumed
to be a muon if it is the most energetic ring, whereas it is always assumed to bea π + if it is a lower energy ring. Because νµ CC events usually have a muon as
the most energetic ring and NC events rarely have a π + as the most energetic
ring.
The final fit result obtained as above is used in this analysis.

6.4

Performance in MC Simulation: APfit vs.
fiTQunIn this section, the SK-IV MC based performances of both APfit and fiTQun
are shown. Here, we use the FC true-fiducial CCQE single electron and muon
77

Number of Events

CHAPTER 6. EVENT RECONSTRUCTION

1200
Fraction of (2R,3R,4R) events (%):
APfit

1000

: 40.3, 58.4, 1.3

fiTQun(before): 16.7, 68.9, 13.7
fiTQun(after) : 20.1, 76.5, 3.3

800
600
400
200
0
1

2

3

4

5
6
7
Number of Rings

Figure 6.23:
Number of rings distributions for free proton events in the
p → e+ π 0 MC. The black histogram shows the APfit result and the blue (red)histogram shows the fiTQun result before (after) the ring correction change.
events in the atmospheric neutrino MC for Figure 6.24-6.29. The fiducial volume
for APfit is defined as the region 2 m inside from the wall while we show thefiTQun performance in two different fiducial volume definitions; 1.5 m and 2 m.
Figure 6.24 shows the distributions of the distance between the reconstructed
and true vertices. The vertex resolutions defined at 1σ deviation position for
single electron (muon) events are; 20.1 (19.3) cm for the best single ring fit
of fiTQun and 29.4 (22.7) cm for APfit. The momentum dependencies of thevertex resolution are shown in Figure 6.25.
The angle between the reconstructed and true directions, and their momentum dependencies are also shown in Figure 6.26 and 6.27. For single electron
and muon events, the direction resolutions defined at an 1σ deviation angle are
2.81◦ (2.91◦ ) and 1.73◦ (1.75◦ ) for the best single ring fit of fiTQun (APfit).
In order to evaluate the momentum resolution, the fractional difference between the reconstructed and true momenta is fitted by a Gaussian at each true
momentum bin, and taking the sigma of the fitted Gaussian as the resolution.
Figure 6.28 shows the momentum resolution as a function of the true momentum.
For the particle identification, the misidentification rates as a function of thetrue momentum are shown in Figure 6.29. Compared to APfit, the best single
ring fit of fiTQun has smaller misidentification rates and they are below 1%.
As discussed in Chapter 7, the most dominant background for p → e+ π 0
is νe CC events accompanied with a single π 0 in final state (νe CC1π 0 ). As
an example of the ring counting and PID performance, Figure 6.30 shows the78

Number of Events

Number of Events

CHAPTER 6. EVENT RECONSTRUCTION

4000
3500

Resolution (1σ):

3000

fiTQun (Dwall>1.5m): 19.6 cm
fiTQun (Dwall>2.0m): 19.5 cm

2500

APfit (Dwall>2.0m): 29.6 cm

2000

4000
fiTQun (Dwall>1.5m): 19.6 cm

3000

fiTQun (Dwall>2.0m): 19.6 cm

2500

APfit (Dwall>2.0m): 22.8 cm

2000

1500

1500

1000

1000

500

500

0
0

Resolution (1σ):

3500

10

20

0
0

30 40 50 60 70 80 90 100
Distance btw. true/recon. vertex (cm)1020

30 40 50 60 70 80 90 100
Distance btw. true/recon. vertex (cm)100

Vertex Resolution (cm)Vertex Resolution (cm)Figure 6.24: Distributions of the distance between the true and reconstructedvertex of FC true-fiducial CCQE single electron (left) and muon (right) eventsin the atmospheric neutrino MC. The green (red) histogram is for the best single
ring fit of fiTQun with the Dwall > 1.5 (2.0) m cut and black histogram is for
APfit with the Dwall > 2.0 m cut. The dotted lines show the vertex resolution
which is defined as 1σ deviation point.

90

fiTQun (Dwall>1.5m)
fiTQun (Dwall>2.0m)
APfit (Dwall>2.0m)

80
70
60
50
40

100
90
70
60
50
40

30

30

20

20

10
0
0

fiTQun (Dwall>1.5m)
fiTQun (Dwall>2.0m)
APfit (Dwall>2.0m)

80

10
200

400

600

0

800 1000 1200 1400
True Momentum (MeV/c)

200

400

600

800 1000 1200 1400
True Momentum (MeV/c)

Figure 6.25: Vertex resolution as a function of true momentum of FC truefiducial CCQE single electron (left) and muon (right) events in the atmospheric
neutrino MC. The green (red) points are for the best single ring fit of fiTQun
with the Dwall > 1.5 (2.0) m cut and black points are for APfit with the Dwall >2.0 m cut.

79

Number of Events

Number of Events

CHAPTER 6. EVENT RECONSTRUCTION

2500
Resolution (1σ):
2000

fiTQun (Dwall>1.5m): 3.04 deg.
fiTQun (Dwall>2.0m): 3.02 deg.

1500

APfit (Dwall>2.0m): 3.11 deg.

1000

4000
3500

Resolution (1σ):

3000

fiTQun (Dwall>1.5m): 1.86 deg.

2500

fiTQun (Dwall>2.0m): 1.84 deg.
APfit (Dwall>2.0m): 1.85 deg.

2000
1500
1000

500

500
0
0

2

0
0

4
6
8
10
12
14
Angle btw. true/recon. direction (deg.)

2

4
6
8
10
12
14
Angle btw. true/recon. direction (deg.)

10

Direction Resolution (deg.)

Direction Resolution (deg.)

Figure 6.26: Distributions of the angle between the true and reconstructed
direction of FC true-fiducial CCQE single electron (left) and muon (right) eventsin the atmospheric neutrino MC. The green (red) histogram is for the best single
ring fit of fiTQun with the Dwall > 1.5 (2.0) m cut and the black histogram is for
APfit with the Dwall > 2.0 m cut. The dotted lines show the direction resolution
which is defined as 1σ deviation point.

9fiTQun (Dwall>1.5m)
fiTQun (Dwall>2.0m)
APfit (Dwall>2.0m)

8
7
6
5
4
3
2

9

fiTQun (Dwall>1.5m)
fiTQun (Dwall>2.0m)
APfit (Dwall>2.0m)

8
7
6
5
4
3
2

1
0
0

10

1
200

400

600

0

800 1000 1200 1400
True Momentum (MeV/c)

200

400

600

800 1000 1200 1400
True Momentum (MeV/c)

Figure 6.27: Direction resolution as a function of true momentum of FC truefiducial CCQE single electron (left) and muon (right) events in the atmospheric
neutrino MC. The green (red) points are for the best single ring fit of fiTQun
with the Dwall > 1.5 (2.0) m cut and the black points are for APfit with the
Dwall > 2.0 m cut.

80

10

Momentum Resolution (%)

Momentum Resolution (%)

CHAPTER 6. EVENT RECONSTRUCTION

9fiTQun (Dwall>1.5m)
fiTQun (Dwall>2.0m)
APfit (Dwall>2.0m)

8
7
6
5
4
3
2

9

fiTQun (Dwall>1.5m)
fiTQun (Dwall>2.0m)
APfit (Dwall>2.0m)

8
7
6
5
4
3
2

1
0
0

10

1
200

400

600

0

800 1000 1200 1400
True Momentum (MeV/c)

200

400

600

800 1000 1200 1400
True Momentum (MeV/c)

5

Mis-PID Rate (%)

Mis-PID Rate (%)

Figure 6.28: Momentum resolution as a function of true momentum for FC
true-fiducial CCQE single electron (left) and muon (right) events in the atmospheric neutrino MC. The green (red) points are for the best single ring fit offiTQun with the Dwall > 1.5 (2.0) m cut and the black points are for APfit with
the Dwall > 2.0 m cut.

4.5

fiTQun (Dwall>1.5m)
fiTQun (Dwall>2.0m)
APfit (Dwall>2.0m)

4
3.5
3
2.5

5
4.5

3
2.5

2

2

1.5

1.5

1

1

0.5

0.5

0
0

fiTQun (Dwall>1.5m)
fiTQun (Dwall>2.0m)
APfit (Dwall>2.0m)

4
3.5

200

400

600

0

800 1000 1200 1400
True Momentum (MeV/c)

200

400

600

800 1000 1200 1400
True Momentum (MeV/c)

Figure 6.29: Mis-PID rate as a function of true momentum for FC true-fiducial
CCQE single electron (left) and muon (right) events in the atmospheric neutrino
MC. The green (red) points are for the best single ring fit of fiTQun with the
Dwall > 1.5 (2.0) m cut and the black points are for APfit with the Dwall > 2.0 m
cut.

81

700

3Reee

600

2Ree

Number of Events

Number of Events

CHAPTER 6. EVENT RECONSTRUCTION

1Re

500

others
400

800

3Reee

700

2Ree

600

1Re
others

500
400

300

300

200

200

100

100
-0.8 -0.6 -0.4 -0.20

0.2 0.4 0.6
Log (E

Number of Events

10

vis

0
-1

0.8 1
(GeV))

800
3Reee

700

2Ree

600

1Re

500

others40000.2 0.4 0.6
Log (E

900

vis

0.8 1
(GeV))

3Reee

800

2Ree

700

1Re

600

others

500
400

300

300

200

200

100
0
-1

-0.8 -0.6 -0.4 -0.2

10

Number of Events

0
-1

100
-0.8 -0.6 -0.4 -0.20

0.2 0.4 0.6
Log (E
10

vis

0
-1

0.8 1
(GeV))-0.8 -0.6 -0.4 -0.20

0.2 0.4 0.6
Log (E
10

vis

0.8 1
(GeV))

Figure 6.30: Event categories of νe CC single-π events with a single π 0 andno visible π ± in final state in the FC atmospheric neutrino MC for the best
multi-ring fit of fiTQun (left) and APfit (right). The top figures are for the
events with the Dwall > 2 m cut. The bottom figures are for the events with the
Dwall > 1.5 m cut. The horizontal axes are the log of the visible energy in GeV.
The colors correspond to three e-like ring events (red), two e-like ring events(blue), single e-like ring events (green) and the other events (black), respectively.

event categories of true νe CC1π 0 events in the FC atmospheric neutrino MC.
It can be seen that the fraction of three e-lke ring events of the best multiring fit of fiTQun is higher than APfit. For example, in the 1 GeV bin, thefraction of three (two) e-like ring events of fiTQun (Dwall > 1.5 m) is 52.8%(38.5%) while that of APfit (Dwall > 2.0 m) is 28.7% (53.4%). Since there
should be three shower rings in this neutrino mode, the fiTQun ring counting
gives more correct information. Figure 6.31 shows the distance between the true
vertex and reconstructed vertex and the π 0 invariant mass distributions in theνe CC1π 0 sample. In the document [97], there is more discussion about otherneutrino modes. As a result, it can generally be said that fiTQun has a higher
multi-ring fit performance than APfit.

82

250

Number of Events

Number of Events

CHAPTER 6. EVENT RECONSTRUCTION

Resolution (1σ):
fiTQun (Dwall>1.5m): 26.3 cm

200

fiTQun (Dwall>2.0m): 26.0 cm
APfit (Dwall>2.0m): 31.7 cm

150

300

fiTQun (Dwall>1.5m): 12.5 %

200

100

50

50
10

20

0
0

30 40 50 60 70 80 90 100
Distance btw. true/recon. vertex (cm)fiTQun (Dwall>2.0m): 12.4 %APfit (Dwall>2.0m): 12.9 %

150

100

0
0

Resolution (1σ):

250

50

100

150
200
250
300
π0 Invariant Mass (MeV/c2)

Figure 6.31: Distributions of the distance between the true and reconstructedvertex (left) and the π 0 invariant mass (right) for νe CC single-π events witha single π 0 and no visible π ± in final state in the FC atmospheric neutrino
MC. The green (red) histogram is for the best multi-ring fit of fiTQun with
the Dwall > 1.5 (2.0) m cut and the black histogram is for APfit with the
Dwall > 2.0 m cut. The dotted lines show the vertex resolution which is definedas 1σ deviation point.

6.5

Fiducial Volume

Since fiTQun is a new fitter and has better performance compared to APfit, it
is worth revisiting the definition of the fiducial volume (FV). In all the previous
proton decay analyses using APfit, FV was defined as the region where the vertex
distance from the nearest wall Dwall is greater than 2 m, which corresponds tothe mass of 22.5 kton. As described in Section 7.1, fiTQun has a same levelsignal efficiency for p → e+ π 0 even if the Dwall > 1.5 m cut is used, comparedto that of APfit with the Dwall > 2 m cut. The fiducial volume is expanded to24.7 kton by using the Dwall > 1.5 m cut. In order to confirm that the Dwall >1.5 m cut is safe for fiTQun, the reconstructed vertices of several samples forboth the data and MC are checked.
As for the first check, the reconstructed entering point of the stopping muon
events is shown in Figure 6.32. Two peaks in the figure corresponds to the
ID wall positions from the center of the tank (top wall: 1810 cm, barrel wall:
1690 cm). The peak position differences between the data and the MC for both
the z and R directions are about 1 cm and their errors are 0.07 ± 0.01%.
Next, Dwall of Michel electrons are shown in Figure 6.33. The agreement
between the data and the MC at the Dwall = 1.5 m bin is consistent. From theseresults, we can say that the vertex of single ring events is properly reconstructedby fiTQun.
In addition to the above, the atmospheric neutrino events which have multirings are checked, since we use multi-ring events in this analysis. Figure 6.34
shows the reconstructed vertices of the multi-ring atmospheric neutrino events
for both the z and R directions. For fiTQun, the reconstructed vertex of themost energetic ring is adopted for the determination of the vertex of the event.
83

CHAPTER 6. EVENT RECONSTRUCTION

Number of Events

Number of Events

The data distribution shapes are similar to the MC, and there are no strange
peaks in the whole region. Thus, we can use the events which are in the Dwall >1.5 m region.
9000
8000
7000

Data: 1806.191± 0.185
MC: 1807.625± 0.105
Error: 0.08± 0.01%

6000
5000
4000

3500

Data: 1683.943± 0.166
MC: 1685.210± 0.097
Error: 0.08± 0.01%

3000
2500
2000
1500

3000

1000

2000

500

1000
0
1400 1500 1600 1700 1800 1900 2000 2100 2200 2300 2400
Z (cm)

0
1400

1500

1600

1700

1800

1900
2000
R (cm)

Figure 6.32: Vertex distributions of the stopping muon events of fiTQun in thez and R directions for the data (black points) and the MC (red). Both of theerrors of the mean value are 0.08 ± 0.01%.

6.6

Energy Scale

The precise momentum determination of the data is necessary for the proton
decay analysis described Chapter 7, since the total momentum and invariant
mass are conclusive parameters. It is ideal to have a control data sample of the
proton decay with various final states, and the energy scale is then evaluated bycomparing it with the simulation. However, taking such a data is not realistic.
Four well-known control samples are instead used for the energy scale: high/low
energy stopping muon, Michel electron and neutral pion sample.
In this section, the absolute energy scale error, the time variation and thedetector uniformity are evaluated. We discuss only about fiTQun for the SK-IV
detector period. The data period used for the calibration of SKDETSIM, which
is April 2009, is selected for the absolute error estimation. Since the water
attenuation length Latt is changing in time as described in Section 3.13, timedependent Latt corrections are applied to both APfit and fiTQun.

6.6.1

High Energy Stopping Muon Sample

The high energy stopping muon events are used to check the energy scale ina range of 1 GeV to 10 GeV by taking ratios of the reconstructed momentumof muons to their track lengths (p/range), namely the energy loss dE/dx since
this is approximately constant. The range is obtained from the muon entering
position and the Michel electron vertex. It is expected to be independent on the
momentum reconstruction, since the muon entering position is defined as the ID
wall position estimated by the direction and the Michel electron vertex can be84

Number of Events

CHAPTER 6. EVENT RECONSTRUCTION

4500
4000
3500
3000
2500
2000
1500
1000
500
0
0

200 400

600

800 1000 1200 1400 1600 1800
Dwall of Michel Electron (cm)Number of Events

Number of Events

Figure 6.33: Dwall of Michel electron vertices of fiTQun for the data (black
points) and the MC (red).

140
120
100
80

350
300
250
200

60

150

40

100

20
0

400

50
-1500

-1000

-500

0

0
0

500
1000
1500
Vertex in Z, 2R or 3R (cm)

200

400

600

800

1000

1200 1400 1600 1800
Vertex in R, 2R or 3R (cm)Figure 6.34: Vertex distributions of the multi-ring (2 or 3-ring) atmospheric
neutrino events of fiTQun in the z and R directions for the data (black points)
and the MC (red).

85

CHAPTER 6. EVENT RECONSTRUCTION

well-reconstructed mainly by the hit time. Therefore, p/range depends almost
only on the momentum. The high energy stopping muon events are selected by
the following criteria:• Total observed charge < 70,000 p.e.s• Muon entering position is on the top (zent > 15 m and Rent < 15 m) or
barrel (Rent < 15 m, −15 m < zent < 15 m) wall
• Number of Michel electrons = 1• Decay time is within a range between 1.5 µs and 10 µs• Michel electron vertex is inside the fiducial volume
where zent and Rent are the muon entering position in the vertical and radius
axes. Figure 6.35 shows p/range distributions of the data and the MC between
the range of 5 m and 35 m with 5 m bin for fiTQun with Dwall = 1.5 m. The
errors are evaluated at each bin by comparing the peak position of a fitted
Gaussian of both the data and MC. All the errors are less than 1%.

6.6.2

Low Energy Stopping Muon Sample

The momentum of the stopping muons can be estimated by the Cherenkov
angle θ as the following formula, since the angle largely depends on the momentum especially the low energy region below around 500 MeV/c as shown in
Figure 6.36.
p(θ) = pmµn2 cos2 (θ) − 1(6.34)

where mµ is the mass of muons. However, no variable related to the Cherenkov
angle is available in fiTQun. Thus, the event selection is done mainly based on
APfit parameters as follow.
• Total observed charge is within a range between 500 p.e.s and 5000 p.e.s• zent > 17.2 m• Direction of the entering muon in the z-direction < −0.87• Goodness of the stopping muon fit > 0.6• Number of Michel electrons = 1• Decay time ≥ 1.2 µs
The ratios of the reconstructed momenta p(p.e.) to p(θ) at several p(p.e.) bins
are used for the error estimation. Figure 6.37 shows p(p.e.)/p(θ) distributions
for fiTQun with Dwall = 1.5 m. The maximum error is (−1.905 ± 0.402)%.

86

Number of Events

2200
2000
1800
1600
1400
1200
1000
800
600
400
200
0
0

Data: 2.219± 0.001
MC: 2.232± 0.001
Error: 0.58± 0.08%

2000
Data: 2.213± 0.002
MC: 2.217± 0.001
Error: 0.18± 0.08%

1800
1600
1400
1200
1000
800
600
400
200

0.5

1

0
0

1.5 2 2.5 3 3.5 4 4.5 5
Momentum/range of Muon (MeV/c/cm)

Number of Events

Number of Events

Number of Events

CHAPTER 6. EVENT RECONSTRUCTION

1600
Data: 2.238± 0.002
MC: 2.240± 0.001
Error: 0.06± 0.09%

1400
1200
1000
800

0.5

1

1.5 2 2.5 3 3.5 4 4.5 5
Momentum/range of Muon (MeV/c/cm)

1200
Data: 2.262± 0.002
MC: 2.254± 0.001
Error:-0.34± 0.10%

1000
800
600

600
400
400
200

200
0.5

1

0
0

1.5 2 2.5 3 3.5 4 4.5 5
Momentum/range of Muon (MeV/c/cm)

Number of Events

Number of Events

0
0

700
Data: 2.253± 0.003
MC: 2.245± 0.001
Error:-0.33± 0.13%

600
500
400

0.5

1

1.5 2 2.5 3 3.5 4 4.5 5
Momentum/range of Muon (MeV/c/cm)

120
Data: 2.226± 0.008
MC: 2.238± 0.003
Error: 0.54± 0.40%

100
80
60

300
40

200

20

100
0
0

0.5

1

0
0

1.5 2 2.5 3 3.5 4 4.5 5
Momentum/range of Muon (MeV/c/cm)

0.5

1

1.5 2 2.5 3 3.5 4 4.5 5
Momentum/range of Muon (MeV/c/cm)

Figure 6.35: Momentum over the range of the stopping muon events with various
ranges for the data (black points) and the MC (red). From the top left plot
to the bottom right plot, the ranges are defined as 5-10 m, 10-15 m, 15-20 m,
20-25 m, 25-30 m and 30-35 m, respectively.

87

Cherenkov Angle (degree)

CHAPTER 6. EVENT RECONSTRUCTION

40

30

20

10

0
0

200

400

600
800
1000
Muon Momentum (MeV/c)

Figure6.36: Cherenkov angle of muons as a function of their momentum (cos θ =q
1 + m2µ /p2 /n).

6.6.3

Michel Electron Energy Spectrum

The energy spectrum of Michel electrons emitted from the stopping muons iswell-known and useful for the energy scale check. Figure 6.38 shows the reconstructed momentum of Michel electrons selected by the same criteria as the highenergy stopping muons. The error is calculated by comparing the mean value ofthe spectrum for both the data and MC, and then is obtained as (−1.80±0.16)%
for fiTQun with Dwall = 1.5 m.

6.6.4

Neutral Pion Mass

Neutral current neutrino events often generate neutral pions, and their invariant
mass Mπ0 can be reconstructed by using two showering rings as follows;pMπ0 = 2p1 p2 (1 − cos Θ),(6.35)
where p1 and p2 are each reconstructed momentum and Θ is the opening angle
between the two rings. In order to evaluate the invariant mass of π 0 , FC eventsof the data and the atmospheric neutrino MC are used. The selection criteria
are as follows:• Vertex in the fiducial volume• Number of rings = 2• All rings are identified as showering ring
88

140

Data: 1.27±0.20
MC: 1.26±0.20
Error:-0.343±0.077%

120
100
80

Number of Events

Number of Events

CHAPTER 6. EVENT RECONSTRUCTION

100

Data: 1.19±0.18
MC: 1.17±0.18
Error:-1.904±0.402%

90
80
70
60
50

60

40

40

30
20

20

10
0.5

0
0

1
1.5
2
2.5
3
fiTQun-P(pe)/P(θ) in 200<P(θ)<280 MeV/c
Number of Events

0
0

80

0.5

1
1.5
2
2.5
3
fiTQun-P(pe)/P(θ) in 280<P(θ)<360 MeV/c

Data: 1.08±0.16
MC: 1.10±0.17
Error:1.779±0.384%

70
60
50
40
30
20
10
0
0

0.5

1
1.5
2
2.5
3
fiTQun-P(pe)/P(θ) in 360<P(θ)<440 MeV/c

Number of Events

Figure 6.37: Distributions of p(p.e.)/p(θ) at each p(θ) bin for the data (black
points) and the MC (red). From the top left plot to the bottom plot, the p(θ)
corresponds to 200-280 MeV/c, 280-360 MeV/c and 360-440 MeV/c, respectively.

1600
Data: 37.769± 0.052
MC: 37.088± 0.031
Error:-1.80± 0.16%

1400
1200
1000
800
600
400
200
0
0

10

20

30

40 50 60 70 80 90 100
Momentum of Michel Electron

Figure 6.38: Momentum of Michel electrons for the data (black points) and theMC (red).

89

Number of Events

CHAPTER 6. EVENT RECONSTRUCTION

Data: 139.8± 0.5
MC: 138.7± 0.1
Error:0.75±0.25%

180
160
140
120
100
80
60
40
20
0
0

50

100

150

200
250
300
π0 Invariant Mass (MeV/c2)

Figure 6.39: Invariant mass of neutral pions for the data (black points) and theMC (red). The solid lines show the Gaussian fitting lines.

• Number of Michel electrons = 0
For fiTQun with Dwall = 1.5 m, the distributions of the invariant π 0 mass
for the data and the MC are shown in Figure 6.39. The Gaussian fitted peak
positions are compared to both the data and MC, and the error is obtained as(0.75 ± 0.25)%.

6.6.5

Energy Scale Error

The final absolute energy scale error of fiTQun is defined as the most deviatederror in all the errors described as above. Figure 6.40 shows the summarized
errors, and then the absolute energy scale error is to be −1.90%. In addition
to this, the time variation of the energy scale should be considered. For that
purpose, the p/range and the momentum of Michel electrons are monitored intime as shown in Figure 6.41. The p/range has been almost stable ((−0.02 ±0.01)% per year), whereas the Michel electron momentum has been changing by(0.21 ± 0.03)% per year. This change is considered as due to the increasing the
PMT hit rate as shown in Appendix B. The time variation errors are defined as
the maximum value of the ratios of the RMS to the mean, which is 0.88% (0.28%
for the p/range, 0.88% for the momentum of Michel electrons). Then, the final
energy scale error is estimated as the square root of the final absolute and time
variation energy scale error, which is 2.1% for fiTQun with Dwall = 1.5 m (2.1%
for APfit with Dwall = 2.0 m).

90(MC-Data)/Data (%)

CHAPTER 6. EVENT RECONSTRUCTION

8
Max: 1.78, Min: -1.90
6
4
2
0
-2
-4
-6
-8

103
104
Momentum (MeV/c)

102

Figure 6.40: Absolute energy scale errors for fiTQun with Dwall = 1.5 m. The
horizontal axis is the estimated momentum. The black and red points show
the errors related to Michel electrons and neutral pions, respectively. The blue
and purple points show the errors related to the low and high energy stoppingmuons, respectively. The final absolute energy scale error is -1.90%.

6.6.6

Uniformity of Energy Scale Error

The detector uniformity of the energy scale is important especially for reconstructing the total momentum in proton decay analysis. In order to check theuniformity, Michel electrons are the good control sample since their vertices
and momenta are uniformly distributed in the detector and in all the directions,
respectively. Michel electrons are selected by the following criteria.
• Total observed charge > 1,000 p.e.s• N50 > 60
• Goodness of the stopping muon fit > 0.5• Number of Michel electrons = 1• Decay time > 1.2 µs• Michel electron vertex is inside the fiducial volume• | cos Θµe | < 0.25
Here, Θµe is the angle between parent muons and Michel electrons. The last
criterion is to select the Michel electron whose direction is perpendicular to the91Momentum/range of Muon (MeV/c/cm)

2.45

Momentum of Michel Electron (MeV/c)

CHAPTER 6. EVENT RECONSTRUCTION

42

2.4

-0.02±0.01% (First year=2008.83)
RMS/Mean=0.006/2.236=0.0028

2.35
2.3
2.25
2.2
2.15
2.1

41

2009 2010 2011 2012 2013 2014 2015 2016
Year

0.21±0.03% (First year=2008.83)
RMS/Mean=0.326/37.210=0.0088

40
39
38
37
36
35
34
33
32

2009 2010 2011 2012 2013 2014 2015 2016
Year

Figure 6.41: Time variation plots of the p/range (top) and the momentum ofMichel electron (bottom) for fiTQun with Dwall = 1.5 m. The red dotted lines
are linear fit functions. The green lines are the value of the MC.

92

Michel Electron Momentum (MC/Data)

CHAPTER 6. EVENT RECONSTRUCTION

1.02

1.01

1

0.99

0.98
-1

-0.5

0

0.5

1
cosθ

Figure 6.42: Averaged momentum of Michel electrons for the MC normalized
by the data as a function of cos θ. The most deviated error in all the cos θ bin
is 0.57%.

parent muon direction in order to avoid the muon polarization effect. Figure 6.42
shows the averaged momentum of Michel electrons for the MC normalized by
the data as a function of cos θ. In this thesis, the uniformity error is defined as
3σ in the vertical axis, which is 1.0% for fiTQun with Dwall = 1.5 m (0.9% for
APfit with Dwall = 2.0 m).

93Chapter 7

Data Analysis
Searches for p → e+ π 0 in SK-IV are performed by both APfit and fiTQun.
First, the selection criteria are defined by using some reconstruction parameters.
Then, the systematic errors of the signal efficiency and the background rate are
estimated. For the background rate estimation, each atmospheric neutrino MC
events are reweighed by considering the solar activity and the two flavor neutrino
oscillation effect (νµ ↔ ντ with ∆m2 = 2.5×10−3 eV2 and sin2 2θ = 1.0). Before
opening the data, we confirm the MC agreement with the data in the outside thesignal region sample. Finally, we discuss the result and set the partial lifetime
limit of protons by combining the result of SK-IV and the APfit results of SK-I
to SK-III.

7.1

Selection Criteria

In this decay mode, a positron and a neutral pion are emitted back-to-back in
the rest frame of a proton. The neutral pion decays promptly (∼ 10−16 sec) into
two γ-rays and these two create two e-like rings by the electromagnetic shower
process after traveling about 40 cm in water. A typical signal MC event which
has three e-like rings is shown in Figure 7.1. The p → e+ π 0 signal is selectedfrom the fully contained (FC) sample by following criteria for both APfit andfiTQun:1. Distance from the wall to the reconstructed vertex > 150 (200) cm for
fiTQun (APfit)
2. 2 or 3 rings
3. All rings should be identified as e-like
4. No Michel electron
5. 85 MeV/c2 < π 0 mass < 185 MeV/c2 for 3-ring events94

CHAPTER 7. DATA ANALYSIS

Figure 7.1: Typical p → e+ π 0 MC event display in SK-IV. A free proton
decayed at around the tank center. It emitted a positron on the right and aneutral pion on the left. As a result of the pion decay, three e-like rings areobserved. Each dot corresponds to a PMT and the color shows the amountof the observed charge. Solid lines represent the fiTQun reconstruction rings.
FiTQun found all the 3 e-like ring correctly.

6. 830 (800) MeV/c2 < total invariant mass < 1050 MeV/c2 andtotal momentum < 250 MeV/c for fiTQun (APfit)
7. No tagged neutrons for SK-IV
As described in Section 6.5, thanks to the better reconstruction of fiTQun,
the definition of the fiducial volume (FV) for fiTQun is expanded from thatof APfit. Figure 7.2 shows the signal efficiencies in each bin of the distancefrom the wall to the true vertex after applying the selection criteria 2-7. For
fiTQun, the reconstructed vertex of the most energetic ring is adopted for thedetermination of the vertex of the event. The efficiencies for fiTQun around near
the wall region are higher than that of APfit. The FV for fiTQun is selected as
1.5 m inside the wall where its efficiency (= (Number of final events)/(Number

95

Efficiency in each bin

CHAPTER 7. DATA ANALYSIS

1
0.8
0.6
0.4
0.2
0
0

200 400 600 800 1000 1200 1400 1600 1800

Distance from the wall to the true vertex (cm)Figure 7.2: Signal efficiencies in each bin of the distance from the wall to
the true vertex after applying the selection criteria 2-7. Black (blue) and red(purple) points show the efficiencies in all (free) protons for APfit and fiTQun,
respectively.

of FC true-FV events)) for all protons becomes similar to that of APfit while
the FV for APfit is defined as 2 m inside the wall.
Since the decaying π 0 has a momentum of 459.43 MeV/c, the two γ-rays are
Lorentz-boosted. Then, in the laboratory frame, one of them can take a very low
momentum which can not be reconstructed. Taking into account such an effect,
we choose the 2- or 3-ring events by the number of rings cut. These selectedrings should be e-like ring. The number of rings distributions after FCFV cut
are shown in Figure 7.3. The fraction of 2-ring event and 3-ring event for free
protons are 18.0% (37.9%) and 77.5% (60.5%) for fiTQun (APfit), respectively.
Figure 7.4 shows the fraction of the number of all e-like ring events and not-all
e-like ring events.
No muons appear in p → e+ π 0 . So, the selection criterion 3 requires thatthe number of Michel electrons is zero (Fig. 7.5).
The mass of π 0 must be reconstructed for the case of 3-rings event. Figure 7.6
shows the distributions of the reconstructed π 0 mass by using two of the threerings and the selection region between 85 to 185 MeV/c. The Gaussian fitted
π 0 mass peak and sigma are 138.0 (135.8) MeV/c2 and 18.6 (20.1) MeV/c2 for96

Number of Events

Number of Events

CHAPTER 7. DATA ANALYSIS

12000
10000
8000

8000
7000
6000
5000

6000

4000
3000

4000

2000
2000

1000

0
1

2

3

4

5

0
1

6
7
Number of Rings

2

3

4

5

6
7
Number of Rings

Figure 7.3: Number of rings distributions of fiTQun (left) and APfit (right) forthe p → e+ π 0 MC in SK-IV (black) and the atmospheric neutrino MC (red)after applying the selection criterion 1. The shaded histogram is for free proton
decay events. The green lines and arrows show the selected region.

Number of Events

Number of Events

16000

14000

14000

12000

12000

10000

10000

8000
8000
6000

6000

4000

4000

2000

2000
0
0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

0
0

2
PID

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2
PID

Figure 7.4: Distributions of a particle identification flag of fiTQun (left) and
APfit (right) for the p → e+ π 0 MC in SK-IV (black) and the atmospheric neutrino MC (red) after applying the selection criteria 1-2. The shaded histogram
is for free proton decay events. The green lines and arrows show the selected
region. The first and second bin correspond to events which is not all e-like andall e-like, respectively.

97

Number of Events

Number of Events

CHAPTER 7. DATA ANALYSIS

105

104

103

105

104

103

102

102

10

10

1
0

0.5

1

1.5

2

2.5

1
0

3
3.5
4
4.5
5
Number of Michel Electrons

0.5

1

1.5

2

2.5

3
3.5
4
4.5
5
Number of Michel Electrons

Figure 7.5: Number of Michel electrons distributions of fiTQun (left) and APfit(right) for the p → e+ π 0 MC in SK-IV (black) and the atmospheric neutrino
MC (red) after applying the selection criteria 1-3. The shaded histogram is forfree proton decay events. The green lines and arrows show the selected region.

Number of Events

Number of Events

fiTQun (APfit), respectively.

2500

2000

1500

1600
1400
1200
1000
800

1000

600
400

500
200
0
0

50

100

150

0
0

200
250
300
π0 Invariant Mass (MeV/c2)

50

100

150

200
250
300
π0 Invariant Mass (MeV/c2)

Figure 7.6: Distributions of π 0 mass of fiTQun (left) and APfit (right) for thep → e+ π 0 MC in SK-IV (black) and the atmospheric neutrino MC (red) after
applying the selection criteria 1-4. The shaded histogram is for free proton
decay events. The green lines and arrows show the selected region.
In order to check the reconstructed proton momentum and invariant mass,
the total momentum Ptot , energy Etot and invariant mass Mtot are calculated
as follows;
Ptot =

all
X

p~i ,

(7.1)

|~
pi | ,(7.2)

2 − P2 ,
Etot
tot

(7.3)

i=0

Etot =

allX
i=0

Mtot =q98

CHAPTER 7. DATA ANALYSIS

Table 7.1: Performance for the determination of the total momentum and invariant mass for all and free protons in SK-IV. The total momentum resolutions
are defined as 1σ deviation from the origin. The total invariant mass peak andwidth are obtained by Gaussian fitting.

Ptot (MeV/c)
Mtot (MeV/c2 )σµσ

fiTQun
all
free

all

161
929.5
32.6

157
909.2
34.7

54
941.7
27.6

APfit
free
66
921.4
30.2

where p~i is the momentum of each ring. Figure 7.7 and 7.8 show the total momentum and invariant mass distributions, respectively. The total momentum
resolutions defined as 1 σ deviation from the origin and the Gaussian fitted total invariant mass peak and sigma are summarized in Table 7.1. Thanks to thebetter momentum fit of fiTQun, the Ptot and Mtot distributions of fiTQun are
sharper than those of APfit for free protons. The difference of Mtot peaks between the distributions for all and free protons is due to the nuclear binding energy as discussed in Section 4.1. Since the observed charge inside the Cherenkov
cone with a half angle of 70◦ is only used in the momentum determination of
APfit (see Section 6.2.4) while fiTQun uses all observed charge, the Mtot peakof fiTQun is closer to the true proton invariant mass than that of APfit. In
the Mtot distribution for the background, there are two peaks; the largest peakaround 130 MeV/c2 is the π 0 mass produced by the delta resonance neutrino
interactions, the smaller peak around the origin is due to a mis-ring counting (CCQE interactions, mostly). The selection criterion 6 highly suppresses
the background from the previous criterion by 0.14 ± 0.01% (0.18 ± 0.02%) forfiTQun (APfit). Because most of the neutrino induced particles are emitted
to the forward region, neutrino events are not isotropic and have higher total
momenta and lower invariant mass.
[Two box analysis] In the previous study, the signal region is just a single box.
However, in order to improve the sensitivity, we divide the signal region definedby the selection criterion 6 into two regions by total momentum: Ptot < 100
MeV/c (low-Ptot box) and 100 ≤ Ptot < 250 MeV/c (high-Ptot box). In the lowPtot box, free proton decay events is concentrated while almost no atmospheric
neutrino background events exist as shown in Figure 7.7. On the other hand,in the high-Ptot box, bound proton decay events are dominant while larger
contamination from neutrino events can be seen. The low-Ptot box has muchbetter sensitivity compared to the higher one. Therefore, we can achieve higher
sensitivity by analyzing these two signal boxes separately if background becomes
non-negligible in the high-Ptot box.
Neutrons often are produced by atmospheric neutrino interactions and captured by hydrogen after traveling a few micro seconds in water. Then 2.2 MeV
99

Number of Events

Number of Events

CHAPTER 7. DATA ANALYSIS

800
700
600
500

700
600
500
400

400

300

300
200
200
100

100
100

200

300

400

500

0
0

600 700 800 900 1000
Total Momentum (MeV/c)

Number of Events

Number of Events

0
0

103

102

200

300

400

500

600 700 800 900 1000
Total Momentum (MeV/c)

100

200

300

400

500

600 700 800 900 1000
Total Momentum (MeV/c)

103

102

10

10

1

1

-1

10-1

10

10-2

100

0

100

200

300

400

500

10-2

600 700 800 900 1000
Total Momentum (MeV/c)

0

Figure 7.7: Distributions of the total momentum of fiTQun (left) and APfit(right) in the linear (top) and log (bottom) scale for the p → e+ π 0 MC in SKIV (black) and the atmospheric neutrino MC (red) after applying the selection
criteria 1-5 and the Mtot cut. The shaded histogram is for free proton decay
events. The green lines and arrows show the selected region.

100

Number of Events

Number of Events

CHAPTER 7. DATA ANALYSIS

1000

1200
1000
800

800

600
600
400
400
200

200

200

400

600

0
0

800
1000
1200
Total Invariant Mass (MeV/c2)

104

Number of Events

Number of Events

0
0

103
102

200

400

600

800
1000
1200
Total Invariant Mass (MeV/c2)

200

400

600

800
1000
1200
Total Invariant Mass (MeV/c2)

103
102

10

10

1

1

-1

-1

10

10

10-2

10-2

0

200

400

600

800
1000
1200
Total Invariant Mass (MeV/c2)

0

Figure 7.8: Total invariant mass distributions of fiTQun (left) and APfit (right)
in the linear (top) and log (bottom) scale for the p → e+ π 0 MC in SK-IV (black)and the atmospheric neutrino MC (red) after applying the selection criteria 1-5
and the Ptot cut. The shaded histogram is for free proton decay events. The
green lines and arrows show the selected region.

101

CHAPTER 7. DATA ANALYSIS

Number of Events

Number of Events

γ-ray is finally emitted from the de-excitation of a nucleus. But, the probabilityof the neutron emission from the de-excitation of nucleus after proton decay inoxygen is small and no neutron emission from hydrogen. Therefore, the neutron
tagging information is useful for proton decay analyses for further backgroundrejection. The neutron tagging algorithm is available only in SK-IV thanks to
the new electronics (QBEE) as described in [77]. It first searches for hit clusters
as 2.2 MeV γ-ray candidates with ≥ 7 hits within a 10 ns sliding window after
the prompt neutrino interaction. After that, 16 variables for a candidate eventare calculated and inputted to a neural network to separate the signal frombackground. The neutron tagging efficiency is 20.5 ± 2.1 % while a mis-tagging
rate is 1.8%. Thanks to this technique, the background rate for p → e+ π 0 is
reduced by almost 50% while only a few percent is lost in the signal efficiency.
Figure 7.9 shows the number of tagged neutrons distributions.

104

103

104

103

102

102

10

10

1

1

-1

-1

10

0

1

2

3

4

5

6

10

7
8
9
10
Neutron Multiplicity

0

1

2

3

4

5

6

7
8
9
10
Neutron Multiplicity

Figure 7.9: Number of tagged neutrons distributions of fiTQun (left) and APfit(right) for the p → e+ π 0 MC in SK-IV (black) and the atmospheric neutrino
MC (red) after applying the selection criteria 1-5, excluding the signal region(criterion 6). The shaded histogram is for free proton decay events. The green
lines and arrows show the selected region.
The signal efficiencies and background rates for p → e+ π 0 in SK-IV are summarized in Table 7.2. Figure 7.10 shows the signal efficiencies as a function of theselection criteria. The ratio between 2-ring and 3-ring events in the final signal
MC events in the total Ptot box is roughly 2:8 (4:6) for fiTQun (APfit). In orderto evaluate the background rate with a large statistics, we use the 3000 years atmospheric neutrino MC which is pre-selected for signal-like events as describedin Appendix D. The effect on the signal efficiencies caused by this pre-selection
is negligible small. Due to the tighter Mtot cut for fiTQun, the background rateof fiTQun is reduced by (33.2 ± 11.4)% in the total Ptot box compared to that ofAPfit while keeping the same level of the signal efficiency. The expected background rate in the 800 < Mtot < 1050 MeV/c2 and Ptot < 250 MeV region is+0.45
consistent with the experimental result of 1.63+0.42−0.33 (stat)−0.51 (syst) events per
Megaton·years obtained by the K2K’s one kiloton water Cherenkov detector
which measured beam neutrinos [98].

102CHAPTER 7. DATA ANALYSIS

Table 7.2: Signal efficiencies and background rates in SK-IV (statistical errors
only). all , free and bound denotes the signal efficiencies for all, free and boundprotons. As for the background, the 3000 years atmospheric neutrino MC which
is pre-selected for signal-like events is used.

Low-Ptotall (%)
free (%)
bound (%)
BG (/Mt/yr)

fiTQunHigh-Ptot

20.0 ± 0.3
82.1 ± 1.4
4.7 ± 0.2
0.028 ± 0.019

18.1 ± 0.3
3.3 ± 0.3
21.7 ± 0.3
0.778 ± 0.102

APfit
Low-Ptot

High-Ptot

19.0 ± 0.3
78.4 ± 1.4
4.2 ± 0.2
0.030 ± 0.021

19.0 ± 0.3
7.7 ± 0.4
21.8 ± 0.4
1.116 ± 0.132

The primary neutrino energy of the final background MC events is mostly
distributed in a range of 1 GeV to 2 GeV as shown in Figure 7.11. The ratio
between 2-ring and 3-ring events in the final background sample in the total Ptot
box is roughly 3:7 for both APfit and fiTQun. The breakdown of the neutrino
interaction mode in the background is summarized in Table 7.3. There are noνµ CC events in both APfit and fiTQun thanks to the good PID performance.
The dominant mode of neutrino interactions is νe charged current with a single π 0 (νe CC1π 0 ) interaction. A typical background νe CC1π 0 event display
is shown in Figure 7.12, which looks like the typical signal event (Figure 7.1).
There are two background events in the low-Ptot box for both APfit and fiTQun
including one common event. The common event is induced by νe CC1π 0 interaction leaving three visible e-like rings (three (two) rings are reconstructedby fiTQun (APfit)) and one proton ring, but the momentum of the proton isnear the Cherenkov threshold and thus this event mimics like a signal event.
The fiTQun-only event is a νe CC1π 0 event with three visible e-like rings which
are properly reconstructed by fiTQun (APfit found only two rings) and seems
an intrinsic background. The APfit-only event is induced by a NC single pion
interaction with multi-πs in final state, in which APfit found three e-like rings
and their reconstructed total mass and momentum are accidentally within thebox. All background event displays overlaying with the true and reconstructedrings are checked and then we find out that there is no serious reconstruction
failure in both APfit and fiTQun.

103

Signal Efficiency

CHAPTER 7. DATA ANALYSIS

1.1
1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2

tot

+ N Low
o n -P
eu tottron+ N Higo n h-P
eu tottron+ N Low
o n -P
eu tottron+ N Higo n h-P
eu tot
tro
n

-P
Hi

gh

tot

wP
Lo

ma

No

π0

ch
Mi

ss

e
el-

e
-lik
le

tot

Signal Efficiency

2o

Al

r3

rin

FV

0

gs

0.1

1.1
1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2

-P

tot

wP

Hi
gh

π0

Lo

ma

ss

e

No
M

ich

el-

e
-lik
Al
le

2o

r3

rin

FV

0

gs

0.1

Figure 7.10: Signal efficiencies after each selection criterion for the p → e+ π 0
MC in SK-IV. Top (bottom) figure is for fiTQun (APfit). Black (blue) histogram
is for all (free) proton decay. The signal efficiencies are normalized to the numberof events after the FC true-FV cut.

104

Number of Events

CHAPTER 7. DATA ANALYSIS

14
12
10
8
6
4
2
0
0

1000

2000

3000 4000 5000 6000 7000
Primary Neutrino Energy (MeV)

Figure 7.11: Primary neutrino energy distributions of the final events in the
total Ptot box in the 3000 years atmospheric neutrino MC. Black (red) histogram
is for APfit (fiTQun).

Table 7.3: Breakdown of neutrino interaction modes for the final events in the
total Ptot box in SK-IV (statistical errors only). The 3000 years atmospheric
neutrino MC which is pre-selected for signal-like events is used. Nev is thenumber of events without the reweighing factor. Numbers in parentheses are
Nev in the low-Ptot box.
Interactionνe CC single-πνe CCQE
νe CC coherent-πνe CC single-ηνe CC multi-πνe CC DIS
NC single-π
NC single-γ
NC multi-π
NC DIS

fiTQunFraction (%)
Nev

APfit
Fraction (%)

Nev

CommonNev

58.5 ± 9.9
18.4 ± 5.5
5.0 ± 2.9
5.0 ± 2.9
3.3 ± 2.4
3.3 ± 2.4
1.7 ± 1.7
1.7 ± 1.7
0060.3 ± 8.6
14.8 ± 4.3
4.9 ± 2.5
4.9 ± 2.5
3.7 ± 2.10
2.5 ± 1.70
3.7 ± 2.1
2.5 ± 1.7

49 (1)
12
4
4
3
0
2 (1)
0
3
2

21 (1)
6
2
1
0
0
0
0
0
0

35 (2)
11
3
3
2
2
1
1
0
0

105

CHAPTER 7. DATA ANALYSIS

Figure 7.12: Typical background MC event display in SK-IV. The type ofneutrino interaction of this event is νe CC1π 0 . The three solid blue lines are thefiTQun reconstructed e-like rings. FiTQun catches each edge of the Cherenkov
ring properly.

106CHAPTER 7. DATA ANALYSIS

7.2

Systematic Errors

In this section, the systematic error estimation of the fiTQun based analysis
in SK-IV is described, which is basically the same as that of the APfit based
analysis. The systematic errors of the APfit based analysis for SK-I to SK-IV
in the latest published paper [1] are used for the combined analysis discussed
in section 7.5. Since the observation time is well defined and the detector isalways full of water, the systematic error of the detector exposure is negligible.
However, we assign the systematic error of 1% for the exposure for just in case.

Pion Final State Interaction
A neutral pion generated by proton decay or a neutrino interaction in oxygen
can interact within that oxygen nucleus. By this nuclear effect, a pion canbe scattered, exchange its charge or even be absorbed (final state interaction,
FSI). It is also possible to interact with other nuclei in water after escaping theoriginal nucleus (secondary interaction, SI). Those interactions are simulated by
the NEUT pion cascade model as discussed in Section 4.1.
In order to evaluate this model uncertainty, the original model is weighted
by the 6 parameters based on some pion scattering data. The detailed weighting
method can be found in [61]. The parameters represent the interaction probabilities for quasi-elastic scattering and charge exchange in two energy regions (πmomentum is less or greater than 500 MeV/c), inelastic scattering and absorption. The 24 weighted parameter sets are prepared and the number of events inthe final sample is calculated for each set. For the pion final state interaction
error, RMS/mean of those 24 values is used.

Correlated Decay
For proton decays in oxygen, it is expected that 10% of nucleons decay underthe effect of other nucleons as described in Section 4.1. We conservatively assign a 100% source error for this effect because it is not well understood. The
systematic error of this decay is calculated by taking the fractional difference
between the following two numbers:
ndef = nH + n16 Onorm + n16 Ocorr
N16 Onorm ∓ αN16 Ocorr
n16 Onorm + (1 ± α)n16 Ocorr
nwei = nH +N16 Onorm

(7.4)
(7.5)

Here, n (N ) is the number of events in the final (FC true-FV) sample, H(16 Onorm, 16 Ocorr) corresponds to proton decays in hydrogen (normal decay in
16
O, correlated decay in 16 O) respectively. α is the source error and it is set to
be 1.

107

CHAPTER 7. DATA ANALYSIS

Fermi Motion
In the signal MC, the nucleon momentum is simulated by using the experimental
spectrum function as described in Section 4.1. On the other hand, the Fermi
gas model is used in the atmospheric neutrino MC. This model uncertainty is
considered by taking the fractional difference between the number of final events
under the spectral function (ndef ) and that of the Fermi gas model (nFermi ).
As shown in Figure 7.13, the proton momentum distribution for the Fermi
gas model is evaluated by weighting the spectrum function and its maximum
momentum is 225 MeV/c which is the Fermi surface momentum simulated in
NEUT. nFermi is defined as,nFermi = nH +11X
n16 O,ii=1

N16 O,iN16 OFermi,i(7.6)

where N16 OFermi,i is the number of events for the i-th momentum bin (the bin
width = 25 MeV) in the FC true-FV sample under the weighted momentum
distribution. N16 OFermi,i is normalized as follows.
11
X
i=1

N16 O,i =11X

N16 OFermi,i

(7.7)

i=1

The Fermi motion error of the APfit high-Ptot box is lower than that of fiTQun
since the fraction of free proton events in the high-Ptot box is larger.

Fiducial Volume
The systematic error of the fiducial volume affects the number of events, directly.
However, there is no control sample for the vertex reconstruction accuracy of
two or three ring events. For single non-shower-ring events, we can estimate thevertex reconstruction accuracy by using the cosmic ray muon events since theirentering points should be the ID wall position as shown in Figure 6.32. The
MC-data differences of the entering muon peak positions are smaller than 2 cm.
We can calculate the fiducial volume error based on this single non-shower-ringsample as 0.4% by shifting the fiducial volume boundary inward/outward with
2 cm.
Another method to evaluate the fiducial volume error is to use the FC subGeV multi-ring atmospheric neutrino MC and data whose event topologies are
more similar to p → e+ π 0 than that of single non-shower-ring events. Here,we check the overall vertex reconstruction accuracy in several volume regions;
Dwall > 50 cm and 100 cm (and also 150 cm for APfit). Figure 7.14 shows the
Dwall distributions and the MCs are normalized by area to the data. In each
distribution, the number of events whose vertices are within the fiducial volume
are calculated for both the data and MC, and an error is then estimated bycomparing that value for both the data and MC. Finally, the fiducial volume
error by this multi-ring sample is obtained as 0.6% (0.7% for APfit) by taking
108

CHAPTER 7. DATA ANALYSIS

1400
3500
1200
3000
1000
2500
800
2000
600

1500
1000

400

500

200

0
0

50

100

0
0

150
200
250
True Proton Momentum (MeV/c)

50

100

150
200
250
True Proton Momentum (MeV/c)

50

100

150
200
250
True Proton Momentum (MeV/c)

1200
3000
1000
2500
800
2000
600
1500
400

1000

200

500
0
0

50

100

0
0

150
200
250
True Proton Momentum (MeV/c)

Figure 7.13: True proton momentum distributions of the final events for the
low-Ptot (left) and the high-Ptot (right) box in the p → e+ π 0 MC. The top(bottom) figures are for fiTQun (APfit). The black (red) histograms show boundproton events under the experimental spectrum function (Fermi gas model) andthe blue histograms show free proton events.

109

CHAPTER 7. DATA ANALYSIS

the average and RMS of those errors. We use this 0.6% error as the systematic
error of the fiducial volume for both the low-Ptot and high-Ptot box.
DATA=4900.4,MC=4925.2,Error=-0.504
Number of Events

Number of Events

DATA=4900.4,MC=4930.2,Error=-0.606
500

400

300

500

400

300

200

200

100

100

0
0

200

400

600

800

0
0

1000 1200 1400 1600 1800
Distance from the wall (cm)

200

400

600

800

1000 1200 1400 1600 1800
Distance from the wall (cm)Figure 7.14: Dwall distributions of the FC sub-GeV multi-ring events for boththe atmospheric neutrino MC (lines) and data (dots) for fiTQun. The red(green) line shows the reconstructed (true) Dwall . The MCs are normalized byarea to the data.

Detector Non-Uniformity
The detector non-uniformity causes the imbalance of the reconstructed momentum which is important to reconstruct protons. As described in Section 6.6, thedetector non-uniformity error is assigned as 1.0% for fiTQun (0.9% for APfit).
The systematic error of the detector non-uniformity is evaluated by shifting thetotal momentum region upward/downward with a twice of the non-uniformity
error and see how the final number varies. Here, we multiply the original nonuniformity error by two since most of p → e+ π 0 events are observed by two
opposite sides each other.

Energy Scale
The energy scale error directly affects the total momentum and the invariant
mass selection, which is set to be 2.1% for fiTQun (2.1% for APfit) as describedin Section 6.6. The number of final events is compared by shifting the total
momentum and invariant mass regions upward/downward with this energy scaleerror, and its fluctuation gives the systematic error of the energy scale.

PID and Ring Counting
The agreement of likelihood distributions between the data and the MC is important for both the PID and ring counting uncertainty. As an example of thegood agreement of the PID likelihood distribution for single-ring events, the110CHAPTER 7. DATA ANALYSIS

mis-PID rates of both the data and MC for the stopping muon (its Michel electron) events are 0.10 ± 0.03% (0.61 ± 0.08%) and 0.15 ± 0.03% (0.48 ± 0.04%),
respectively (see Figure 6.19).
For more realistic estimation, we use the FCFV sub-GeV multi-ring sample
and FCFV sub-GeV e-like ring sample for the PID and ring counting, respectively. The systematic errors of the PID and ring counting likelihoods are evaluated by comparing each likelihood distribution of the atmospheric neutrino MC
and the data. If there are some systematic differences between the data andthe MC, the likelihood distribution would be wider/narrower and/or shifted upward/downward direction. Therefore, the likelihood distribution of the MC is
shifted and scaled so that the disagreement with the data becomes minimal (aχ2 test) as shown in Figure 7.15, and the fractional difference of events passed
through the threshold between the tuned MC and the data is used as the systematic error. The systematic errors for the PID and ring counting likelihoods
are obtained as 0.3% and 0.3%, respectively. We assign these errors for boththe low-Ptot and high-Ptot box.

Neutrino Flux and Cross Section
The systematic errors of the neutrino flux and the cross section are evaluatedby an event-by-event weighting method for each error based on atmospheric
neutrino analyses in SK, which can be found in [99, 100].
For the neutrino flux error, the following error sources are considered in thisanalysis.
• Energy dependent normalization• Neutrino flavor ratio• ν/ν ratio• Up/down asymmetry• Horizontal/vertical ratio• K/π production ratio• Neutrino flight length
The dominant source error in the neutrino flux is the energy dependent normalization which includes the uncertainty of the production of π and K in a
hadronic interaction model, hadronic interaction cross sections and atmosphericdensity profile. Those model uncertainties are evaluated in [65] by comparing
the observed data with the calculation of the atmospheric muon flux as shownin Figure 7.16. Table 7.4 shows each neutrino flux error.
For the neutrino cross section, the following items are considered.
• MA (axial vector mass) in quasi-elastic scattering and single-meson production111

CHAPTER 7. DATA ANALYSIS

400
350
300
250
200
150
100
50
0
-2000 -1500 -1000 -500

0

500

30

40

1000 1500 2000
ln(L e/L µ)

2000
1800
1600
1400
1200
1000
800
600
400
200
0
-10

0

10

20

60

70

sgn × ln(L50/ L 1Re)

2Rex

Figure 7.15: Distributions of thepPID likelihood ratio ln(Le /Lµ ) (top) andthe ring counting likelihood ratio | ln(L2Rex /L1Re )| (bottom) in the FCFV
sub-GeV sample for the data (black points), the nominal MC (black histogram)and the tuned MC (red histogram). In the top figure, the likelihood ratio of themost energetic ring is displayed.
112

Relative Uncertainty

CHAPTER 7. DATA ANALYSIS

0.30

δ tot
δπ
δK
δσ
δ air

0.25
0.20
0.15
0.10
0.05
0

0

10

10

1

10

2

10
E ν (GeV)

3

Figure 7.16: Uncertainties of the atmospheric neutrino flux as a function ofneutrino energy. δπ , δK , δσ , δair and δtot are the uncertainty of π and K
production, hadronic interaction cross sections, atmospheric density profile andtheir sum. Taken from [65].

• Quasi-elastic scattering for bound nucleons (total cross section)• Quasi-elastic scattering for bound nucleons (ν̄/ν)• Quasi-elastic scattering for bound nucleons (flavor ratio)• Single(multi)-meson production (total cross section)• Single(multi)-meson production (model dependence)• Coherent pion production• NC/CC ratio
The largest source error in the neutrino cross section is the single-meson production which is determined by comparing the form factor of the model [84]with the experimental data from bubble chamber experiments [85] as discussed
in 4.2. Table 7.5 shows each neutrino cross section error.

Michel Electron
The systematic error of the Michel electron tagging is evaluated by comparing
the number of tagged Michel electrons between the FCFV sub-GeV one µ-likering sample of the atmospheric neutrino MC and the data. Figure 7.17 shows113

CHAPTER 7. DATA ANALYSIS

Table 7.4: Systematic errors from uncertainties of the neutrino flux in SK-IV.
fiTQun

APfit

7.9%
4.5%

7.0%
4.0%

0.7%
1.1%
<0.1%

0.6%
0.8%
<0.1%

(Energy dependent)
Normalization(Eν <1GeV)(Eν >1GeV)νµ /νe ratio(Eν <1GeV)(1 < Eν <10GeV)(Eν >10GeV)ν¯e /νe ratio(Eν <1GeV)(1 < Eν <10GeV)(Eν >10GeV)ν¯µ /νµ ratio(Eν <1GeV)(1 < Eν <10GeV)(Eν >10GeV)
Up/down asymmetry
Horizontal/vertical ratio
K/π production ratio
Neutrino flight length

0.9%
2.2%
<0.1%

0.8%
1.6%
<0.1%

<0.1%
<0.1%
<0.1%
<0.1%
0.2%
<0.1%
<0.1%

<0.1%
0.2%
<0.1%
0.3%
<0.1%
<0.1%
<0.1%

Total

9.5%

8.3%

Table 7.5: Systematic errors from uncertainties of the neutrino cross sections in
SK-IV.
fiTQunAPfit

MA in QE and single-pion production
NC elastic and CCQE ratio
CCQE ν̄/ν ratio
CCQE µ/e ratio
Single-meson production
Single-pion production
DIS (model comparison)
DIS
Coherent pion production
NC/CC ratio

9.9%
7.5%
0.9%
0.1%
9.6%
14.5%
2.5%
0.4%
1.2%
0.7%

7.5%
6.2%
0.4%
0.1%
9.4%
22.4%
1.6%
<0.1%
1.5%
<0.1%

Total

21.6%

26.3%

114

CHAPTER 7. DATA ANALYSIS

the distributions of the number of Michel electrons. This systematic error forthe signal efficiency is negligible small.
DATA=82.7%,MC=82.2%,Error=0.602%
Number of Events

Number of Events

DATA=85.4%,MC=85.5%,Error=-0.053%
4000

4500
4000

3500

3500

3000

3000

2500

2500
2000
2000
1500

1500

1000

1000

500

500
0
0

0.5

1

1.5

2

2.5

3

0
0

3.5
4
4.5
5
Number of decay-e

0.5

1

1.5

2

2.5

3

3.5
4
4.5
5
Number of decay-e

Figure 7.17: Distributions of the number of Michel electrons for both fiTQun(left) and APfit (right). The dots (red lines) show the data (atmospheric neutrino MC) of the FCFV sub-GeV one µ-like ring events.

Neutron Tagging
In the SK-IV analysis, the neutron tagging information is used. The systematic
error of the neutron tagging is estimated by comparing the number of tagged
neutrons between the atmospheric neutrino MC and the data. The events outside the signal region are used in this check after applying the other selection
criteria. Figure 7.18 shows the distributions of the number of neutrons. Thissystematic error for the signal efficiency is negligible small.
DATA=451.0 evts,MC=465.1 evts,Error=-3.030 ± 6.408 (%)
Number of Events

Number of Events

DATA=426.0 evts,MC=414.7 evts,Error=2.724 ± 7.086 (%)
400
350
300
250

400
350
300
250

200

200

150

150

100

100
50

50
0
0

1

2

3

4

0
0

5
6
7
8
Number of Tagged Neutrons

1

2

3

4

5
6
7
8
Number of Tagged Neutrons

Figure 7.18: Distributions of the number of tagged neutrons for both fiTQun(left) and APfit (right). The dots (red lines) show the data (atmospheric neutrino MC) of the events outside the signal region after applying the other selection criteria.

115CHAPTER 7. DATA ANALYSIS

Table 7.6: Systematic errors for the signal efficiencies in SK-IV.

fiTQunLow-Ptot High-Ptot(Interaction)
Pion FSI/SI
Correlated decay
Fermi motion(Detector & Reconstruction)
Fiducial volume
Detector non-uniformity
Energy scale
PID
Ring counting
Total

APfit (reference)
Low-Ptot High-Ptot

2.9%
2.0%
9.0%

11.9%
10.0%
10.3%

2.6%
1.9%
8.3%

11.7%
9.4%
6.2%

0.6%
1.3%
1.5%
0.3%
0.3%

0.6%
0.5%
2.1%
0.3%
0.3%

0.7%
1.7%
1.8%
2.7%
0.1%

0.7%
0.3%
1.9%
2.7%
0.1%

9.9%

18.8%

9.7%

16.6%

Summary of the Systematic Errors
Table 7.6 and 7.7 show the systematic errors for the signal efficiencies and thebackground rate. The systematic errors of the signal efficiency for fiTQun (APfit) are 9.9% (9.7%) and 18.8% (16.6%) for the low-Ptot box and the high-Ptot
box, respectively. The systematic errors of the background rate for fiTQun(APfit) is 41.4% (40.6%) for the total Ptot box. Both fiTQun and APfit have
similar systematic errors in the SK-IV analysis.

116

CHAPTER 7. DATA ANALYSIS

Table 7.7: Systematic errors for the background rates for the total Ptot box in
SK-IV.

(Interaction)
Neutrino flux
Neutrino cross section
Pion FSI/SI(Detector & Reconstruction)
Fiducial volume
Detector non-uniformity
Energy scale
PID
Ring counting
Michel electron
Neutron tagging
Total

117

fiTQunAPfit (reference)

9.5%
21.6%
9.3%

8.3%
26.3%
6.6%

0.6%
10.3%
29.4%
0.3%
0.3%
0.1%
9.8%

0.7%
10.0%
25.5%
2.7%
0.1%
0.6%
9.4%

41.4%

40.6%

CHAPTER 7. DATA ANALYSIS

7.3

MC-Data Agreement in the Outside the Signal Region

Number of Events

Number of Events

Before opening the data in the signal region, the agreement between the SKIV data and the atmospheric neutrino MC is checked in the outside the signal
region data with a margin of 100 MeV. Figure 7.19 - 7.23 show the distributionsof the outside sample for the number of Michel electrons, π 0 mass, number oftagged neutrons, total momentum and invariant mass, respectively. In these
distributions, there are no significant differences between the data and the MC
can be found. The data was well reproduced by the MC.

2500

2000

2500

2000

1500

1500
1000

1000

500

500

0
0

0.5

1

1.5

2

2.5

0
0

3
3.5
4
4.5
5
Number of Michel Electrons

0.5

1

1.5

2

2.5

3
3.5
4
4.5
5
Number of Michel Electrons

Number of Events

Number of Events

Figure 7.19: Number of Michel electrons distributions of fiTQun (left) and APfit(right) for the outside the signal region sample in the SK-IV data (black) and the
500 years atmospheric neutrino MC (red) after applying the selection criteria
1-3. The MC is normalized by the live time of the data.

100

80

80
70
60
50

60
40
40

30
20

20
10
0
0

50

100

150

200

250

0
0

300 350 400 450 500
π0 Invariant Mass (MeV/c2)

50

100

150

200

250

300 350 400 450 500
π0 Invariant Mass (MeV/c2)

Figure 7.20: Distributions of π 0 mass of fiTQun (left) and APfit (right) for theoutside the signal region sample in the SK-IV data (black) and the 500 years
atmospheric neutrino MC (red) after applying the selection criteria 1-4. The
MC is normalized by the live time of the data.

118

Number of Events

Number of Events

CHAPTER 7. DATA ANALYSIS

103

102

102

10

10

1
0

103

1

2

3

4

5

6

1
0

7
8
9
10
Neutron Multiplicity

1

2

3

4

5

6

7
8
9
10
Neutron Multiplicity

Number of Events

Number of Events

Figure 7.21: Number of tagged neutrons distributions of fiTQun (left) and APfit(right) for the outside the signal region sample in the SK-IV data (black) and the
500 years atmospheric neutrino MC (red) after applying the selection criteria
1-5. The MC is normalized by the live time of the data.

120
100
80

100

80

60
60
40
40
20

20
0
0

100

200

300

400

500

0
0

600 700 800 900 1000
Total Momentum (MeV/c)

100

200

300

400

500

600 700 800 900 1000
Total Momentum (MeV/c)

Figure 7.22: Distributions of the total momentum of fiTQun (left) and APfit(right) for the outside the signal region sample in the SK-IV data (black) and the
500 years atmospheric neutrino MC (red) after applying the selection criteria
1-5. The MC is normalized by the live time of the data.

119

Number of Events

Number of Events

CHAPTER 7. DATA ANALYSIS

102

10

10

1
0

102

200

400

600

1
0

800
1000
1200
Total Invariant Mass (MeV/c2)

200

400

600

800
1000
1200
Total Invariant Mass (MeV/c2)

Figure 7.23: Total invariant mass distributions of fiTQun (left) and APfit (right)
for the outside the signal region sample in the SK-IV data (black) and the 500
years atmospheric neutrino MC (red) after applying the selection criteria 1-5.
The MC is normalized by the live time of the data.

120

CHAPTER 7. DATA ANALYSIS

7.4

ResultSearch for p → e+ π 0 in the 179.5 (163.0) kton·years SK-IV data has been
conducted by using fiTQun (APfit) with the selection criteria described in Section 7.1. The number of expected background events in the total Ptot box for
fiTQun (APfit) is about 0.14 (0.19) events. No candidate event was found in
the data for both fiTQun and APfit.
The scatter plots of the 2-dimensional space of the total invariant mass andmomentum for the signal MC, the atmospheric neutrino MC and the data in
SK-IV are shown in Fig. 7.24 - 7.26. Figure 7.26 demonstrates that fiTQun
can perform a very clean search for proton decay as well as APfit. In order tocompare the shape of the 2D distributions between the data and the MC, we usean estimator L which is defined as a minimum distance between an event point
and a vertical line (Mtot = 938.3 MeV/c2 , 0 ≤ Ptot ≤ 250 MeV/c) in the 2D
space. The good agreement in the 2D distribution for both fiTQun and APfit
can be seen in Figure 7.27. The event reduction plots are shown in Figure 7.28,
our MC reproduces the characteristics of the data well.

121

CHAPTER 7. DATA ANALYSIS

Figure 7.24: Total invariant mass and momentum scatter plots for the p →e+ π 0 MC in SK-IV after applying the selection criteria 1-5 and 7. The top(bottom) figure is for fiTQun (APfit). The light (dark) blue shows free (bound)proton decay events. The two boxes correspond to the Ptot < 100 MeV/c region(lower box) and 100 ≤ Ptot < 250 MeV/c region (upper box).

122

CHAPTER 7. DATA ANALYSIS

Figure 7.25: Total invariant mass and momentum scatter plots for the 500
years-live-time atmospheric neutrino MC in SK-IV after applying the selection
criteria 1-5 and 7. The top (bottom) figure is for fiTQun (APfit). The two boxes
correspond to the Ptot < 100 MeV/c region (lower box) and 100 ≤ Ptot < 250
MeV/c region (upper box).

123

Total Momentum (MeV/c)

CHAPTER 7. DATA ANALYSIS

1000
900
800
700
600
500
400
300
200
100

Total Momentum (MeV/c)

0
0

200

400

600

800
1000
1200
Total Invariant Mass (MeV/c2)

200

400

600

800
1000
1200
Total Invariant Mass (MeV/c2)

1000
900
800
700
600
500
400
300
200
100
0
0

Figure 7.26: Total invariant mass and momentum scatter plots for the FC data
set in SK-IV after applying the selection criteria 1-5 and 7. The top (bottom)
figure is the result of the 179.5 kton·years data (163.0 kton·years data) with
fiTQun (APfit). The two boxes correspond to the Ptot < 100 MeV/c region(lower box) and 100 ≤ Ptot < 250 MeV/c region (upper box).

124

Number of Events

CHAPTER 7. DATA ANALYSIS

103

102

10

Number of Events

1
0

200

400

600
800
1000
1200
Distance L in (M tot, P tot) Space

200

400

600
800
1000
1200
Distance L in (M tot, P tot) Space

103

102

10

1
0

Figure 7.27: Distributions of the distance L in SK-IV after applying the selection criteria 1-5 and 7. The top (bottom) figure is for fiTQun (APfit). The
black points, red and blue histograms correspond to the FC data, the 500 yearslive-time atmospheric neutrino MC and the p → e+ π 0 MC, respectively. The
atmospheric neutrino MC is normalized by the live time of the data. The signal
MC is scaled to the atmospheric neutrino MC.

125Number of Events

CHAPTER 7. DATA ANALYSIS

105
104
103
102
10
1

tot

+ N Low
o n -P
eu tottron+ N Higo n h-P
eu tottron+ N Low
o n -P
eu tottron+ N Higo n h- P
eu tot
tro
n

-P
Hi

gh

tot

Lo

wP

ss
ma

ich
No
M

π0

ele

e
-lik
le

tot

Number of Events

2o

Al

r3

FV

10-2

rin
gs

10-1

105
104
103
102
10
1

-P
Hi
gh

tot

wP
Lo

ma
π0

ch
No

Mi

ss

e
el-

e
-lik
le
Al

2o

r3

rin

FV

10-2

gs

10-1

Figure 7.28: Number of events for the FC data set (black points) and the 500
years-live-time atmospheric neutrino MC (red) in SK-IV after each selection
criterion. The top (bottom) figure is for fiTQun (APfit). The atmospheric
neutrino MC is normalized by the live time of the data.

126

CHAPTER 7. DATA ANALYSIS

7.5

Lifetime Limit

Although we could not find any candidate events for p → e+ π 0 , we can set thelower limit of the proton partial lifetime by using a Bayesian method. Since
we have 8 independent searches for p → e+ π 0 (i = 1 ∼ 4 correspond to the
low-Ptot box analyses for SK-I to SK-IV, i = 5 ∼ 8 correspond to the high-Ptot
box analyses for SK-I to SK-IV), we calculate the lower limit by combining theresults. As for the APfit based analysis, the signal efficiencies and the expectedbackground rates for SK-I to SK-III and the systematic errors for SK-I to SK-IV
are taken from [1]. The signal efficiencies and the expected background rates
for SK-IV and the systematic errors for the SK-IV fiTQun based analysis are
obtained by using the new MC described in Chapter 4.
The probability of detecting n events is given by a Poisson statistics if the
decay rate Γ (> 0), the exposure λi (> 0), the signal efficiency i (0 < i < 1)
and the number of background bi (> 0) are known.
P (ni |Γλi i bi ) =e−(Γλi i +bi ) (Γλi i + bi )ni
ni !

(7.8)

Applying Bayes’ theorem to this probability, we get the following formulas.
P (Γλi i bi |ni )P (ni ) = P (ni |Γλi i bi )P (Γλi i bi )= P (ni |Γλi i bi )P (Γ)P (λi )P (i )P (bi )(7.9)(7.10)

Here, we use that the decay rate, the exposure, the signal efficiency and thenumber of background are independent. The probability density function of the
decay rate is given as follows;
ZZZ
P (Γ|ni ) =P (Γλi i bi |ni )di dλi dbi(7.11)Z Z Z −(Γλi i +bi )e(Γλi i + bi )ni1
P (Γ)P (i )P (λi )P (bi )di dλi dbi ,
=
Ai
ni !
(7.12)
where the normalization constant Ai is,Z ∞Ai =P (Γ|ni )dΓ.

(7.13)

0

We assume the probability of the decay rate is uniform and the probabilities
P (λi ) and P (i ) are defined as a Gaussian;P (Γ) = 1,

(7.14)22

2

2

P (λi ) ∝ e−(λi −λ0i ) /2σλi(7.15)andP (i ) ∝ e−(i −0i ) /2σi ,127(7.16)

CHAPTER 7. DATA ANALYSIS

where λ0i (σλi ) and 0i (σi ) are the estimations (systematic errors) of the signal
efficiency and the exposure, respectively. In order to take into account thestatistical error of the number of background, we use a convolution of Poisson
and Gaussian distributions for the probability of the number of background as
follows;
Z ∞ −B nbi
22
e B(7.17)
P (bi ) ∝e−(Ci bi −B) /2σbi dB,
nbi !
0
where nbi is the number of background events in the 3000 (500) years MC for
SK-IV (SK-I to SK-III), Ci is a constant to normalize the MC live time to thedata live time and σbi is the systematic error of the number of background. The
lower limit of the decay rate, Γlimit , at 90% confidence level is given by,
0.9 =Z Γlimit NY=8
Γ=0

P (Γ|ni )dΓ.

(7.18)

i=1

Then, the lower limit of the partial lifetime of p → e+ π 0 is calculated as follows.
τ /Bp→e+ π0 =1Γlimit(7.19)

By using only the SK-IV data, we obtain the lower partial lifetime limit ofp → e+ π 0 for the APfit based analysis as,τ /Bp→e+ π0 = 8.58 × 1033 years at 90% C.L.,(7.20)

and for the fiTQun based analysis as,τ /Bp→e+ π0 = 9.47 × 1033 years at 90% C.L.

(7.21)The lower partial lifetime limit of p → e+ π 0 can be calculated by using theresults of the APfit based analyses for the SK-I to SK-IV data as,τ /Bp→e+ π0 = 1.79 × 1034 years at 90% C.L.

(7.22)

Finally, we conclude the lower partial lifetime limit of p → e+ π 0 by combining
the APfit based analysis for the SK-I to SK-III data with the fiTQun based
analysis for the SK-IV data as,τ /Bp→e+ π0 = 1.88 × 1034 years at 90% C.L.

(7.23)

This lifetime limit is 18% longer than that of the latest published paper (1.6 ×1034 years) [1].

128Chapter 8

Discussion
8.1

Comparison with GUT models

We discussed several attractive GUT models in Chapter 2, the predictions of
the proton partial lifetime for p → e+ π 0 are rewritten as follows [9, 14, 22, 23].
Minimal SU (5) : 1031±1 years
Minimal SUSY SU (5) : > 4.1 × 1033 years
Minimal SO(10)

:

32

10

− 10

39

years

(8.1)
(8.2)
(8.3)

In contrast to these predictions, we obtain the world best experimental limit byreplacing the reconstruction algorithm applied to the half of the data set from
APfit to fiTQun as,
τ /Bp→e+ π0 = 1.88 × 1034 years at 90% CL.

(8.4)

By this result, we can reconfirm that the minimal SU (5) model is completely
ruled out. Note that this conclusion can be derived only from the fiTQun result
and it is not trivial in case of a new reconstruction algorithm. This experimental
result excludes part of the predicted regions of both the minimal SUSY SU (5)and minimal SO(10).
In order to estimate how much our energy frontier is expanded, we canroughly write down the lifetime assuming the mass of a mediated boson MX isenough large as,
τp ∼

4
1 MX
.
2
5
αG mp

(8.5)

where αG is the coupling constant at the GUT scale. Assuming αG = 40−1 bylooking at Figure 2.2, the lower limit of the GUT energy scale is obtained asfollow.
MX & 4.5 × 1015 GeV
129

(8.6)

CHAPTER 8. DISCUSSION

This corresponds to the length of 4.4 × 10−32 m (the Planck scale is 1.62 ×
10−35 m). This experimental limit is the highest energy frontier for humanity
measured on the earth at this moment.

8.2

Further Improvements for fiTQun

In this section, we discuss about possibilities of further improvements for fiTQun.
The dominant background for p → e+ π 0 are νe CC events leaving a primary e±
and 2γs from a π 0 in final state. On the other hand, NC events can mimic suchan event topology by generating 2π 0 s (NC2π 0 ) when one of four gammas fromthem does not have enough energy to produce a visible ring for the fitter. The
ring counting of fiTQun tends to find more rings compared to that of APfit
and thus fiTQun can reject more NC2π 0 events by correctly finding all rings.
An example of a NC2π 0 event which is a background event for the APfit based
analysis in SK-IV but not for fiTQun is displayed in Figure 8.1. In this figure,it can be seen that fiTQun found all four e-like rings (APfit did not reconstruct
the larger ring in the top). However, as shown in Table 7.3, there are remaining
background NC events. Figure 8.2 shows an example of such events, the ringcounting of fiTQun did not work well for this event (a similar situation can befound in APfit). Therefore, it may be possible to achieve an analysis with no
NC background event by improving the ring counting algorithm in fiTQun.
Since fiTQun is the maximum likelihood fitter, it may be natural to construct a p → e+ π 0 hypothesis likelihood and conduct a likelihood analysis. For
example, if we use a three-e-like-ring hypothesis while fixing a π 0 mass for signal
events instead of the multi-ring hypothesis, it is naively expected to have better
momentum resolution and thus a higher signal efficiency. In order to establish alikelihood analysis by using such a dedicated likelihood fitter, we need a larger
statistics of the atmospheric neutrino MC to estimate the likelihood distribution
for background events because the background rate should be very small for a
p → e+ π 0 search. However, the process of the fitter will take unrealistic time
if we want to generate an over-thousand-years full MC. For that purpose, thepreselected method described in Appendix D will be useful to extract the eventsof interest.

8.3

Discovery Potential in Future

Although we may observe proton decay in SK in future, it is necessary to havea bigger detector in order to determine a correct GUT model by measuring
the branching ratios for many proton decay modes. The Hyper-Kamiokande(HK) project is proposed as a next generation water Cherenkov detector [101](Figure 8.3). The detector consists of a cylindrical stainless tank, 258 kton
ultra-pure light water and 40,000 20-inch PMTs (40% photocoverage). The
fiducial volume is 187 kton (178 kton) when we use a region of Dwall > 1.5(2.0) m. In HK, it is expected to use a new type of photosensor which has

130CHAPTER 8. DISCUSSION

Figure 8.1: Typical NC2π 0 event display in the SK-IV atmospheric neutrino
MC. The four solid blue (green) lines are the fiTQun reconstructed e-like rings(MC true shower rings from gammas). FiTQun found all four rings whereas
APfit did not reconstruct the larger ring in the top.

131

CHAPTER 8. DISCUSSION

Figure 8.2: A NC2π 0 event display in the SK-IV atmospheric neutrino MC.
The four solid blue (green) lines are the fiTQun reconstructed e-like rings (MCtrue shower rings from gammas). FiTQun failed to reconstruct two of four true
rings (APfit also failed).

132

CHAPTER 8. DISCUSSION

Figure 8.3:
from [101].

Schematic view of the Hyper-Kamiokande detector.

Taken

better performance than that of the SK PMT.
In principle, fiTQun is expected to work on HK since HK is a natural extension of SK. However, it should be taken into account of the difference of PMTs
and so on. In the following sentence, we discuss advantages to use fiTQun inboth SK and HK by comparing a discovery potential of p → e+ π 0 for both APfit(Dwall > 2.0 m) and fiTQun (Dwall > 1.5 m) assuming the HK performance isthe same as SK-IV whose photocoverage is the same as that of HK. Here, we
adopt the two-box analysis for both SK and HK.
The discovery potential for p → e+ π 0 is calculated by using the profile
likelihood ratio method [102]. A test statistic q0 is defined as the log likelihood
ratio between the null hypothesis (background-only) and the signal hypothesis(signal+bakcground);(ˆL(0,θ̂)−2 ln L(µ̂,µ̂ ≥ 0
θ̂)
q0 =
(8.7)
0
µ̂ < 0where µ̂ is the signal strength which means the decay rate Γ in this case. θ isa nuisance parameter while µ̂ is a parameter of interest. The null hypothesis
ˆlikelihood in the numerator is maximized by an estimator θ̂, whereas the signal

133

CHAPTER 8. DISCUSSION

hypothesis one in the denominator is maximized by an estimator θ̂. That is,
!
N
2
2
Y
θeff,iθbkg,i(Si + Bi )ni
1
L(Γ, θ̂) =×
exp − 2 − 2, (8.8)2 σ2ni !
2πσeff,i
2σeff,i2σbkg,ibkg,ii=0where N is the number of analysis bins for each total momentum box and eachSK, σeff,i (σbkg,i ) and θeff,i (θbkg,i ) are the systematic error and its nuisanceparameter for the signal efficiency (number of background) in the i-th bin. The
expected number ni is given as Si + Bi , where Si = Γλi i × (1 + θeff,i ) andBi = bi × (1 + θbkg,i ). The probability function of q0 for the µ = 0 case is given
as,
 q 
1
1
1
0
f (q0 |0) = δ(q0 ) + √ √ exp −
(8.9)
2
2
2 2π q0
Using the cumulative distribution of Equation 8.9, Φ, the p-value of the null
hypothesis is then calculated as,√
p0 = 1 − Φ( q0 )(8.10)
By using Equation 8.10, the median significance for the signal hypothesis canbe calculated.
The 3σ discovery corresponding proton partial lifetime as a function of live
time is obtained for both SK and HK as shown in Figure 8.4-8.5. As for the APfit
based analysis, the signal efficiencies and expected background rates for SK-I
to SK-III and the systematic errors for SK-I to SK-IV are taken from [1]. The
signal efficiencies and expected background rates for SK-IV and the systematic
errors for the SK-IV fiTQun based analysis are obtained by using the new MC
described in Chapter 4. At the point of the 30th year of exposure of SK (after 15
years of exposure from today), the APfit based analysis can approach at 3.14 ×1034 years while the combined analysis (APfit for SK-I to SK-III and fiTQun
for SK-IV) can approach at 3.48 × 1034 years which is 11% higher sensitivity.
After 30 years of exposure of HK (solid lines), the APfit based analysis canapproach at 1.46 × 1035 years while the fiTQun based analysis can approachat 1.56 × 1035 years which is 7% higher sensitivity. There may be a possibility
to construct a second HK tank which has the same configuration as the firstone after six years of the one tank operation. If we consider such a two tankoperation, the fiTQun (APfit) based analysis can approach at 2.30 (2.16) ×
1035 years after 30 years of exposure of HK (dotted lines).
Once fiTQun is able to work for the whole SK period properly, the discovery
potential of the full fiTQun based analysis is expected to much better. As an example of the full fiTQun based analysis, we calculate the 3σ discovery potential
assuming the performance (signal efficiency, background rate and systematic
errors) in SK-I to SK-III to be the same as that of SK-IV as shown in Figure 8.4(green line) get the value of 3.91 × 1034 years after 30 years of exposure.

134

CHAPTER 8. DISCUSSION

45 ×10τ/ β (years)

33

40
35
30
25
20
15
10
5
0
0

5

10

15

20
25
30
Live Time (years)

Figure 8.4: 3σ discovery potential for p → e+ π 0 in SK as a function of live
time. Each line shows the different analysis; Black: SK-I to SK-IV based on
APfit, Red: SK-I to SK-III based on APfit and SK-IV based on fiTQun, Green:
SK-IV based on fiTQun. In this thesis, we stand on the point of 14.9 years in
SK.

135

CHAPTER 8. DISCUSSION

250 ×10τ/ β (years)

33

200
150
100
50
0
0

5

10

15

20
25
30
Live Time (years)

Figure 8.5: 3σ discovery potential for p → e+ π 0 in HK as a function of live
time. Black (Red) line shows the HK based on APfit (fiTQun) assuming its
performance is the same as that of SK-IV. Dotted lines show the HK two tankoperation assumed to start after six years of the one tank operation.

136

CHAPTER 8. DISCUSSION

It may become a long journey for mankind to discover proton decay, but
it is just a moment for protons. If we really want to know the deep inside ofnature, we must keep observation and improve our detectors and techniques ina long-term strategy.

137

Chapter 9

Conclusion
A proton decay search for p → e+ π 0 using the improved event reconstruction
algorithm in Super-Kamiokande was discussed. We have newly developed the
new event reconstruction algorithm, fiTQun, based on the maximum likelihood
method. This is the first time to change the concept of the event reconstruction
for SK proton decay analyses. We choose p → e+ π 0 as a first application offiTQun based proton decay analyses since this decay mode is the most dominant
in many non-SUSY GUTs and the most sensitive in SK. In this thesis, we applied
fiTQun for the SK-IV data whose live time is almost half of the entire SK livetime.
In the signal MC study, fiTQun shows better performance compared to theconventional algorithm, called APfit. While keeping the same level of the signal efficiency as APfit, the fiducial volume of fiTQun is defined as 24.7 kton
which is 10% higher than that of APfit. One may expect signal candidates will
be found in this unexplored volume region by using fiTQun instead of APfit.
The three ring fraction of fiTQun (APfit) for free proton decay events is 77.5%(60.5%). The total invariant mass peak and its width of fiTQun (APfit) for
free proton decay events are 941.7 MeV/c2 (921.4 MeV/c2 ) and 27.6 MeV/c2(30.2 MeV/c2 ), respectively. Thanks to this better mass reconstruction performance of fiTQun, we can reduce the atmospheric neutrino background rate
about −30% compared to APfit by applying the tighter total invariant mass cut(800 ∼ 1050 MeV/c2 → 830 ∼ 1050 MeV/c2 ). After the proper calibrations forboth the data and fiTQun was done, fiTQun’s energy scale and detector nonuniformity error became similar to those of APfit (energy scale: 2.1% (2.1%),
non-uniformity: 1.0% (0.9%) for fiTQun (APfit)). The new electronics, QBEE,
employed in SK-IV enabled us to reduce the background by approximately halfusing the neutron tagging technique.
We have searched for p → e+ π 0 in the 179.5 kton·years FC data taken in
SK-IV from Sep. 2008 to Aug. 2016 by using fiTQun. Since the latest publishedpaper [1], 483 live days data were added from 4973 live days to 5456 live days
(including the SK-I to SK-III data). For more accurate background estimation,
we newly generated and used the 2500 years atmospheric neutrino MC with the138

CHAPTER 9. CONCLUSION

pre-selection for p → e+ π 0 , in addition to the conventional 500 years MC. In
the total signal region (=low-Ptot box + high-Ptot box), the signal efficiency
and number of expected background events are about 38% and 0.14 events,
respectively. Unfortunately, no signal candidate was found in our data. Then,the lower partial lifetime limit of p → e+ π 0 for the fiTQun based analysis canbe calculated by using Bayes’ theorem as 9.47 × 1033 years at 90 % C.L. This
limit is 10% longer than that of the APfit based analysis (8.58 × 1033 years),
mainly due to the larger fiducial volume of fiTQun. Finally, we combined theAPfit results for SK-I to SK-III (172.5 kton·years) [1] and the fiTQun result for
SK-IV, the lower lifetime limit is given as,
τ /Bp→e+ π0 = 1.88 × 1034 years at 90% C.L.
This lower lifetime limit is 18% longer than that of the latest published SK paper(1.6×1034 years [1]) and thus the most stringent constraint on non-SUSY GUTs
in the world.

139

Appendix A

Special Handling for the
Blasting in the Mine
During May 7th to Aug. 4th in 2015, the blasting work was conducted in the
Kamioka mine for constructing a room of a new water system for SK-Gd project.
The location of the blasting point is about 100 m away from the SK site. The
total number of blasting is 132 times and the blasting timing is recorded by aperson in charge of the Kamioka Mining and Smelting Company.
By eye-scanning the final FC data sample, we noticed that there are two
flasher events coincided with the blasting moments, which survived several
flasher rejection algorithms in the FC reduction (Fig. A.1). Figure A.2 showssome FC reduction output rates as a function of the time difference between
the SK DAQ time and the blasting time (top left: FC 3rd reduction output,
top right: Rejected events in FC 4th reduction, bottom left: FC 5th reduction output, bottom right: FCFV output). From these figures, it is obviousthat the blasting work caused flasher events although the mechanism is unknown (shaking the tank, fluctuating the grounding etc). In the top left figure,the peak position and width are fitted to be −39 s and 61 s by a function off (t) = p0 ×exp[−(t−p1 )2 /(2p2 )2 ]+p3 , respectively. The non-zero peak position
is considered to be due to the systematic time gap between SK and the mine
company.
Since the mechanism of flasher is uncertain, we decided not to use the data
which recorded at a time between −500 s and 400 s from the blasting time inany high energy analysis for the sake of cleanness of the data set. The deadtime due to this handling is about 1.4 days (∼ 11 FC neutrino events).

140

APPENDIX A. SPECIAL HANDLING FOR THE BLASTING IN THE
MINE

Figure A.1: Flasher events coincided with the blasting moments. The color
indicates the hit timing (red: early hit, blue: late hit). The top (bottom) event
happened at May 27th, 2015 at 22:29 JST (June 2nd, 2015 at 01:50 JST), it is
the same time as the nearby blasting moment.

141

APPENDIX A. SPECIAL HANDLING FOR THE BLASTING IN THE
MINE

hfc3tot
EntriesMean
RMS
χ2 / ndf
p0
p1
p2
p3

30
25
20

543
197.3
1972
162.9 / 72
26.08 ± 6.02
-38.79 ± 11.32
60.88 ± 12.69
3.781 ± 0.230

10

8

6
15
4

10

2

5
0
-4000

-3000

-2000

-1000

0

0
-4000

1000
2000
3000
4000
Super-K time - Blasting time (s)

20

-3000

-2000

-1000

0

1000
2000
3000
4000
Super-K time - Blasting time (s)

-3000

-2000

-1000

0

1000
2000
3000
4000
Super-K time - Blasting time (s)

5

18
16

4

14
12

3

10
8

2

6
4

1

2
0
-4000

-3000

-2000

-1000

0

0
-4000

1000
2000
3000
4000
Super-K time - Blasting time (s)Figure A.2: Some FC reduction output rates as a function of the time difference
between the SK DAQ time and the blasting time. Top left: FC 3rd reduction
output, top right: Rejected event in FC 4th reduction, bottom left: FC 5th
reduction output, bottom right: FCFV output.

142

Appendix B

Time Variation of PMT Hit
Rate
It is empirically known that the PMT hit rate has a time dependence andits variation varies depending on the detector period. Since the PMT gain is
increasing but the hit threshold does not change, there might include such athreshold effect in the change of the PMT hit rate. However, if we compare
the PMT hit distributions with different time periods, it appears that the distribution shape does not change but the number of hits is increased (Fig. B.1).
Thus, it seems that the dark rate of PMTs has been also changing. In any
case, those PMT hit rate change will affect to physics analysis in particular for
low energy events. As described Appendix C, the reconstructed momentum ofMichel electrons changed in time due to the PMT hit rate change. In order tocorrect such effect, the PMT hit rate plot as a function of time is necessary.
Also, it is worth to make it since there was no time variation plot of the hit ratewith a consistent estimation method for whole SK period.
For the time variation plot of the PMT hit rate, the stopping muon data
sample was used because it is well understand and available in whole SK period. First of all, the stopping muon fitter was applied to the data. Using thereconstructed entering point of muons, the time-of-flight subtracted hit timing
distribution can be made (Fig. B.2). Then, the early off timing hits in a 100 ns(250 ns) time window for SK-I to SK-III (for SK-IV) are counted up. The time
window is set to be a sufficiently flat area. Finally, the run-by-run averaged hit
rate is estimated by the following formula.
PMT hit rate (kHz) =

n
X
i

nHITi
× 106
T × nEVi × (nPMT − nBADi )(B.1)

where nHITi is the number of hits in the i-th run in the time window T , nEVi
is the number of events in the i-th run, nPMT is the number of ID PMTs andnBADi is the number of bad ID PMTs on average in the i-th run. Figure B.3
shows the time variation of the averaged PMT hit rate for whole SK period. The
143

APPENDIX B. TIME VARIATION OF PMT HIT RATE

Figure B.1: Hit charge distributions for a PMT set of a same production year(the hit charge (p.e.s) vs. the number of events). Black histogram is for theperiod from Dec. 2008 to Jan. 2009. Red histogram is for the period on Mar.
2015. Taken from a SK internal document.

PMT hit rate is increasing as a whole although there is variation. A significant
decrease in SK-II can be seen, it is thought that it is due to the decay of a flash
light rate caused by a chemical reaction in the FRP case. That flasher fromthe FRP case seems to be defeated by the other hit rate source at around 2008.
The gap between the end point of SK-II and the start point of SK-III might bedue to the tank opening, it activated the dark rate of PMTs. The replacementof electronic boards is the reason why the gap between the end point of SK-III
and the start point of SK-IV happened.
An application of the time dependent PMT hit rate for fiTQun is shown inAppendix C.

144

APPENDIX B. TIME VARIATION OF PMT HIT RATE

Averaged PMT hit rate (kHz)

Figure B.2: Time-of-flight subtracted hit timimig distribution for SK-II (the hit
timing (ns) vs. the number of events). Magenta line shows the peak position
and two blue lines indicate the time window. The peak around 1700 ns is dueto the signal reflection between a PMT and an ATM board.

6

5

4

3

2000

2005

2010

2015
Year

Figure B.3: Time variation of the averaged PMT hit rate from the beginningof SK to mid. 2015. The colors shows the detector periods; black, red, green
and blue show SK-I, -II, -III and -IV, respectively.

145

Appendix C

Development of FiTQun for
SK-I to SK-III
Almost half of the observed data in whole SK period at this moment consistsof SK-I, II and III data, and there might be proton decay events which were
not discovered by APfit. Towards full fiTQun-based analysis for SK-I to SKIV, we have developed fiTQun not only for SK-IV but also SK-I to SK-III. By
replacing the Cherenkov profiles, the timing distribution, etc., fiTQun can work
on SK-I to SK-III properly. We have slightly modified the code regarding thedata taking by the ATM during the development. After the basic development,
MC based performance for SK-I to SK-III became similar to that of SK-IV.
The detector quality had been largely changed over time as described in
Section 3.6 and Appendix B, such as the water attenuation length, the PMT
gain and dark rate. It affects the reconstructed momentum in data, and thenmade large data-MC discrepancies. Therefore, the time dependent corrections of
the detector variables were implemented to fiTQun. Figure C.1 shows the time
variation plots of the muon momentum over range (p/range) for the stopping
muon data in SK-I with a no correction case, with a PMT gain correction caseand with a PMT gain and water attenuation length correction case. Afterthe time dependent corrections, the time dependence of p/range was reduced
significantly. In addition to the above corrections, the dark rate correction
in time was useful for the low energy reconstruction such as Michel electron
momentum. Figure C.2 shows the time variation plots of Michel electrons inthe stopping muon sample with the no correction case and with the dark rate
correction case. Since the number of hits are low for Michel electrons, thereconstructed momentum is largely affected by the effect of the PMT dark
rate. After the dark rate correction, the time dependence of the Michel electron
momentum becomes flat.
However, there are still non-negligible data-MC discrepancies for SK-I to
SK-III. Figure C.3 shows the PID likelihood ratio (ln Le /Lµ ) distributions ofboth the data and MC for the FC sub-GeV single-ring atmospheric neutrino

146

APPENDIX C. DEVELOPMENT OF FITQUN FOR SK-I TO SK-III

Figure C.1: Time variation plots of the muon momentum over range for the
stopping muon data in SK-I. The top left, right and bottom plots are for theno correction case, with the PMT gain correction case and with the PMT gain
and water attenuation length correction case. The dotted red lines show liner
fitted slopes.

events. The data distributions for SK-I to SK-III shifted to the left (more µlike) while good agreement of data-MC can be seen in SK-IV. It is consideredthat these discrepancies come from the time distribution difference between thereal PMT response and the MC modeling in the ATM. Although the PMT time
distribution is an asymmetric Gaussian, SKDETSIM for SK-I to SK-III (ATMperiod) assume a symmetric Gaussian. In near future, fiTQun will be able toreconstruct all the SK data by improving the detector simulation.

147

APPENDIX C. DEVELOPMENT OF FITQUN FOR SK-I TO SK-III

Momentum of decay-e (MeV/c)

SK2-0.84±0.13% (First year=2003.08)
RMS/Mean=0.448/38.359=0.0117

40.5
40
39.5
39

38.5
38
37.5

2003

2003.5

2004

2004.5

2005

2005.52006
Year

Figure C.2: Time variation of the Michel electron momentum for the stopping
muon data in SK-II after the PMT gain and water attenuation length correction. Left (right) figure shows the momentum before (after) the PMT dark rate
correction. The dotted red lines show liner fitted slopes.

800
600
400
200
0

200
100
0

SK1 Data
νeCC 1e
νeCC oth.
νµ CC 1µνµ CC oth.
NC

-2000

400

200

0

SK3 Data
νeCC 1e
νeCC oth.
νµ CC 1µνµ CC oth.
NC

-2000

0

2000
ln(Le/Lµ )

1000

500

0

0

2000
ln(Le/Lµ )SK2 Data
νeCC 1e
νeCC oth.
νµ CC 1µνµ CC oth.
NC

-1000

0

1000
ln(Le/Lµ )02000
ln(Le/Lµ )SK4 Data
νeCC 1e
νeCC oth.
νµ CC 1µνµ CC oth.
NC

-2000

Figure C.3: Distributions of the PID log likelihood ratio (ln Le /Lµ ) of both
the data (dots) and MC (stacked histograms) for the FC sub-GeV single-ringatmospheric neutrino events. Taken from [97].

148Appendix D

Mass Production of
Atmospheric Neutrino MC
In order to obtain the background rate for the p → e+ π 0 analysis, we conventionally used the MC for 500 years of livetime for each SK which is 25 times
longer than the current SK observation time. Since this analysis is very clean,
the number of background events in the 500 years MC is low (around 10 events)
and thus the statistical error is large. In addition to that, it is worth to check the
difference of the background tendency between APfit and fiTQun because it is
the first time to apply fiTQun to this analysis. Therefore, we need a statisticallylarger atmospheric neutrino MC event set for SK-IV. Taking into account thestatical error size and our computer resource (Process speed of fiTQun multiring fit in our computer system is about 8.8 min/events for the preselected
events described below), the target livetime is set to be 3000 years. We produce
an additional 2500 years MC for SK-IV.
The following criteria are applied to the MC events before processing both
APfit and fiTQun since most of events are out of the signal region and it takes
amount of time to process the event reconstruction;
1. NHITAC < 16,
2. 6, 000p.e.s < PEtot < 13, 000p.e.s,|ri |i qi ~< 0.5,
PEtot

P

3.

where ~ri is an unit vector pointing to the i-th PMT which observes charge ofqi from a timing based fitted vertex. The first criterion is the one of the fullycontained cut. Since the energy region of proton decay is around 1 GeV, thesecond criterion is chosen so that both the signal and background events arecontained with a margin as shown in Figure D.1. The third criterion represents
the cut for an anisotropy in the observed charge of an event. Since a positron anda neutral pion are emitted back-to-back, the anisotropy of the signal events is149

APPENDIX D. MASS PRODUCTION OF ATMOSPHERIC NEUTRINO
MC

105

103

104
2

10

103

102
10
10
1
0

1
2000

4000

6000

8000

10000
12000
14000
Total Charge (PE)

0

2000

4000

6000

8000

10000
12000
14000
Total Charge (PE)Figure D.1: Total observed charge distributions for the p → e+ π 0 MC (left)
and the 500 years atmospheric neutrino MC (right). The histograms are theoriginal events (black), the events after the NHITAC cut (green), the final events
selected by APfit without both the fiducial volume and neutron cut (blue) andthe final events selected by fiTQun without both the fiducial volume and neutroncut (red), respectively.

expected to be small while neutrino induced particles are boosted to the incident
neutrino direction. Figure D.2 shows the anisotropy distributions. The fractionof survived background events is 3.5% after applying these criteria. After thepreselection for p → e+ π 0 , both APfit and fiTQun are applied and then the
resultant 3000 years MC is then used in this analysis.
The validation of the preselection is done by comparing the nominal 500
years MC and the preselected 3000 years MC for both APfit and fiTQun. Figure D.3 and D.4 show the reconstructed total momentum and invariant mass
distributions. In these figures, the 3000 years MC is normalized by the livetimeof the other MC. It can be seen the good agreement between both the MC
inside and around the signal region. Thus, the preselection does not affect the
distribution shape and the number of background events. Figure D.5 shows thetotal invariant mass and momentum scatter plots in the preselected 3000 years
MC for APfit and fiTQun after applying the selection criteria 1-5 for p → e+ π 0(not applied the neutron tagging cut).

150

APPENDIX D. MASS PRODUCTION OF ATMOSPHERIC NEUTRINO
MC

105

103

104
102

103
102

10

10
1
0

1
0.1

0.2

0.3

0.4

0.5

0.6

0.7
0.8
0.9
1
Charge Anisotropy

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7
0.8
0.9
1
Charge Anisotropy

Number of Events

Number of Events

Figure D.2: Charge anisotropy distributions for the p → e+ π 0 MC (left) andthe 500 years atmospheric neutrino MC (right). The histograms are the original
events (black), the events after the NHITAC cut (green), the events after the
NHITAC and the total observed charge cut (purple), the final events selectedby APfit without both the fiducial volume and neutron cut (blue) and the final
events selected by fiTQun without both the fiducial volume and neutron cut(red), respectively.

70
60
50

60
50
40

40
30
30
20

20

10

10
0
0

200

0
0

400
600
800
1000
1200
Total Momentum in 700<Mtot<1150MeV (MeV/c)

200

400
600
800
1000
1200
Total Momentum in 700<Mtot<1150MeV (MeV/c)

Figure D.3: Total momentum distributions in both the nominal 500 years(black) and preselected 3000 years (red) MC for fiTQun (left) and APfit (right).
The selection criteria 1-5 for p → e+ π 0 and the 700 MeV < Mtot < 1150 MeV
cut are applied.

151

Number of Events

Number of Events

APPENDIX D. MASS PRODUCTION OF ATMOSPHERIC NEUTRINO
MC

60
50
40

60
50
40

30

30

20

20

10

10

0
0

200

400

0
0

600
800
1000
1200
Total Mass in Ptot<350MeV/c (MeV/c2)

200

400

600
800
1000
1200
Total Mass in Ptot<350MeV/c (MeV/c2)

Figure D.4: Total invariant mass distributions in both the nominal 500 years(black) and preselected 3000 years (red) MC for fiTQun (left) and APfit (right).
The selection criteria 1-5 for p → e+ π 0 and the Ptot < 350 MeV cut are applied.

Figure D.5: Total invariant mass and momentum scatter plots in the preselected 3000 years MC for fiTQun (left) and APfit (right). The selection criteria
1-5 for p → e+ π 0 are applied.

152

Bibliography[1] K. Abe et al. “Search for proton decay via p → e+ π 0 and p → µ+ π 0
in 0.31megatonyears exposure of the Super-Kamiokande water Cherenkov
detector”. Phys. Rev., D95(1):012004, 2017.
[2] B. P. Abbott et al. “Observation of Gravitational Waves from a Binary
Black Hole Merger”. Phys. Rev. Lett., 116(6):061102, 2016.
[3] G. Aad et al. “Observation of a new particle in the search for the Standard
Model Higgs boson with the ATLAS detector at the LHC”. Phys. Lett.,
B716:1–29, 2012.
[4] S. Chatrchyan et al. “Observation of a new boson at a mass of 125 GeVwith the CMS experiment at the LHC”. Phys. Lett., B716:30–61, 2012.
[5] Y. Fukuda et al. “Evidence for oscillation of atmospheric neutrinos”. Phys.
Rev. Lett., 81:1562–1567, 1998.
[6] H. F. Dylla and J. G. King. “Neutrality of Molecules by a New Method”.
Phys. Rev., A7:1224–1229, 1973.
[7] C. Patrignani et al.
C40(10):100001, 2016.

“Review of Particle Physics”.

Chin. Phys.,[8] H. Georgi and S. L. Glashow. “Unity of All Elementary Particle Forces”.
Phys. Rev. Lett., 32:438–441, 1974.
[9] P. Langacker. “Grand unification and the standard model”. arXiv, hepph/9411247.
[10] P. Langacker and N. Polonsky. “Uncertainties in coupling constant unification”. Phys. Rev., D47:4028–4045, 1993.
[11] U. Amaldi et al. “Comparison of grand unified theories with electroweak
and strong coupling constants measured at LEP”. Phys. Lett., B260:447–
455, 1991.
[12] J. R. Ellis, S. Kelley, and D. V. Nanopoulos. “Probing the desert using
gauge coupling unification”. Phys. Lett., B260:131–137, 1991.

153

BIBLIOGRAPHY[13] P. Langacker and M. Luo. “Implications of precision electroweak experiments for Mt , ρ0 , sin2 θW and grand unification”. Phys. Rev., D44:817–
822, 1991.
[14] J. Hisano, H. Murayama, and T. Yanagida. “Nucleon decay in the minimal
supersymmetric SU(5) grand unification”. Nucl. Phys., B402:46–84, 1993.
[15] H. Murayama and A. Pierce. “Not even decoupling can save minimal
supersymmetric SU(5)”. Phys. Rev., D65:055009, 2002.
[16] S. P. Martin. “A Supersymmetry primer”. 1997. [Adv. Ser. Direct. HighEnergy Phys.18,1(1998)].
[17] P. Minkowski. “µ → eγ at a Rate of One Out of 109 Muon Decays?”.
Phys. Lett., B67:421–428, 1977.
[18] T. Yanagida. “Horizontal Symmetry and Masses of Neutrinos”. Prog.
Theor. Phys., 64:1103, 1980.
[19] M. Gell-Mann, P. Ramond, and R. Slansky. “Complex Spinors and Unified
Theories”. Conf. Proc., C790927:315–321, 1979.
[20] H. Georgi. “The State of the ArtGauge Theories”. AIP Conf. Proc.,
23:575–582, 1975.
[21] H. Fritzsch and P. Minkowski. “Unified Interactions of Leptons andHadrons”. Annals Phys., 93:193–266, 1975.
[22] N. T. Shaban and W. James Stirling. “Minimal left-right symmetry andSO(10) grand unification using LEP coupling constant measurements”.
Phys. Lett., B291:281–287, 1992.
[23] D. Lee, R. N. Mohapatra, M. K. Parida, and M. Rani. “Predictions
for proton lifetime in minimal nonsupersymmetric SO(10) models: An
update”. Phys. Rev., D51:229–235, 1995.
[24] H. Weyl. “Electron and Gravitation. I.”. Z. Phys., 56:330–352, 1929.
[25] E. C. G. Stueckelberg. “Interaction forces in electrodynamics and in thefield theory of nuclear forces”. Helv. Phys. Acta, 11:299–328, 1938.
[26] E. P. Wigner. “Invariance in Physical Theory”. Am. Phil. Soc., 93:521,
1949.
[27] F. Reines, C. L. Cowan, and M. Goldhaber. “Conservation of the numberof nucleons”. Phys. Rev., 96:1157–1158, 1954.
[28] S. P. Rosen. “Another Method of Searching for Proton Decay”. Phys.
Rev. Lett., 34:774, 1975.

154

BIBLIOGRAPHY[29] G. N. Flerov, D. S. Klochkov, V. S. Skobkin, and V. V. Terent’ev. “Spontaneous fission of Th*232 and the stability of nucleons”. Sov. Phys. Dokl.,
3:79, 1958.
[30] J. C. Evans, Jr. and R. I. Steinberg. “Nucleon Stability: A Geochemical
Test Independent of Decay Mode”. Science, 197:989–991, 1977.
[31] E. L. Fireman. “Depth Dependence Of Ar-37 From K-39: Consequences
For Solar Neutrino Experiments And The Baryon Conservation Principle.”. pages 389–393, 1979.
[32] C. L. Bennett. “Limitations on Proton Decay Modes from a Passive Detection Scheme”. pages 120–124, 1981.
[33] F. Reines, C. L. Cowan, and H. W. Kruse. “Conservation of the numberof nucleons”. Phys. Rev., 109:609–610, 1958.
[34] G. K. Backenstoss et al. “An investigation of the stability of nucleons”.
Nuovo Cimento, 16:749–755, 1960.
[35] C. C. Giamati and F. Reines. “Experimental test of the conservation of
nucleons”. Phys. Rev., 126:2178–2188, 1962.
[36] W. R. Kropp and F. Reines. “Experimental test of baryon conservationand neutrino flux limits”. Phys. Rev., 137:B740–B751, 1965.
[37] H. S. Gurr, W. R. Kropp, F. Reines, and B. Meyer. “Experimental testof baryon conservation”. Phys. Rev., 158:1321–1330, 1967.
[38] L. Bergamasco and P. Picchi. “An Experimental Lower Limit on the
Proton Lifetime”. Nuovo Cim. Lett., 11:636, 1974.
[39] F. Reines and M. F. Crouch. “Baryon conservation limit”. Phys. Rev.
Lett., 32:493–494, 1974.
[40] J. Learned, F. Reines, and A. Soni. “Limits on Nonconservation of
Baryon Number”. Phys. Rev. Lett., 43:907, 1979. [Erratum: Phys. Rev.
Lett.43,1626(1979)].
[41] G. Battistoni et al. “Fully Contained Events in the Mont Blanc Nucleon
Decay Detector”. Phys. Lett., B118:461, 1982.
[42] M. R. Krishnaswamy, M. G. K. Menon, N. K. Mondal, V. S. Narasimham,
B. V. Sreekantan, Y. Hayashi, N. Ito, S. Kawakami, and S. Miyake. “Results From the Kgf Proton Decay Experiment”. Nuovo Cim., C9:167–181,
1986. [Conf. Proc.C850418,97(1985)].
[43] John E. Bartelt et al. “Monopole Flux and Proton Decay Limits From
the Soudan-i Detector”. Phys. Rev., D36:1990, 1987. [Erratum: Phys.
Rev.D40,1701(1989)].

155

BIBLIOGRAPHY[44] Christoph Berger et al. “Results from the Frejus experiment on nucleon
decay modes with charged leptons”. Z. Phys., C50:385–394, 1991.
[45] M. L. Cherry, M. Deakyne, K. Lande, C. K. Lee, R. I. Steinberg, and B. T.
Cleveland. “Experimental Test of Baryon Conservation: A New Limit on
the Nucleon Lifetime”. Phys. Rev. Lett., 47:1507–1510, 1981.
[46] C. McGrew et al. “Search for nucleon decay using the IMB-3 detector”.
Phys. Rev., D59:052004, 1999.
[47] K. S. Hirata et al. “Experimental Limits on Nucleon Lifetime for Lepton+ Meson Decay Modes”. Phys. Lett., B220:308–316, 1989.
[48] M. Shiozawa et al. “Search for proton decay via p → e+ π 0 in a large water
Cherenkov detector”. Phys. Rev. Lett., 81:3319–3323, 1998.
[49] H. Nishino et al. “Search for Nucleon Decay into Charged Anti-lepton plus
Meson in Super-Kamiokande I and II”. Phys. Rev., D85:112001, 2012.
[50] K. Abe et al. “Search for proton decay via p → νK + using 260kilotonyear
data of Super-Kamiokande”. Phys. Rev., D90(7):072005, 2014.
[51] Y. Fukuda et al. “The Super-Kamiokande detector”. Nucl. Instrum.
Meth., A501:418–462, 2003.
[52] R. Becker-Szendy et al. “IMB-3: A Large water Cherenkov detector for nucleon decay and neutrino interactions”. Nucl. Instrum. Meth., A324:363–
382, 1993.
[53] H. Ikeda et al. “Front end hybrid circuit for Super-Kamiokande”.
Nucl.Instrum.Meth., A320:310–316, 1992.
[54] H. Nishino et al. “High-speed charge-to-time converter ASIC for theSuper-Kamiokande detector”. Nucl. Instrum. Meth., A610:710–717, 2009.
[55] K. Abe et al. “Calibration of the Super-Kamiokande Detector”. Nucl.
Instrum. Meth., A737:253–272, 2014.
[56] R. Pope and E. Fry. “Absorption spectrum (380-700 nm) of pure water. II.
Integrating cavity measurements”. Applied Optics, 33:8710–8723, 1997.
[57] K. Nakamura et al. “The Reaction C-12 (e, e’ p) at 700-MeV and DWIA
Analysis”. Nucl. Phys., A268:381–407, 1976.
[58] T. Yamazaki and Y. Akaishi. “Nuclear medium effects on invariant mass
spectra of hadrons decaying in nuclei”. Phys. Lett., B453:1–6, 2000.
[59] Y. Hayato. “NEUT”. Nucl. Phys. Proc. Suppl., 112:171–176, 2002.
[60] L. L. Salcedo et al. “Computer Simulation of Inclusive Pion Nuclear
Reactions”. Nucl. Phys., A484:557–592, 1988.

156

BIBLIOGRAPHY[61] P. de Perio. “NEUT pion FSI”. AIP Conf. Proc., 1405:223–228, 2011.
[62] G. Rowe, M. Salomon, and R. H. Landau. “An Energy Dependent Phase
Shift Analysis of Pion - Nucleon Scattering Below 400-MeV”. Phys. Rev.,
C18:584–589, 1978.
[63] R. D. Woods and D. S. Saxon. “Diffuse Surface Optical Model for NucleonNuclei Scattering”. Phys. Rev., 95:577–578, 1954.
[64] H. Ejiri. “Nuclear deexcitations of nucleon holes associated with nucleon
decays in nuclei”. Phys. Rev., C48:1442–1444, 1993.
[65] M. Honda et al. “Calculation of atmospheric neutrino flux using theinteraction model calibrated with atmospheric muon data”. Phys. Rev.,
D75:043006, 2007.
[66] M. Honda et al. “Improvement of low energy atmospheric neutrino
flux calculation using the JAM nuclear interaction model”. Phys. Rev.,
D83:123001, 2011.
[67] J. Alcaraz et al. “Cosmic protons”. Phys. Lett., B490:27–35, 2000.
[68] J. Alcaraz et al. “Helium in near Earth orbit”. Phys. Lett., B494:193–202,
2000.
[69] T. Sanuki et al. “Precise measurement of cosmic ray proton and helium
spectra with the BESS spectrometer”. Astrophys. J., 545:1135, 2000.
[70] S. Haino et al. “Measurements of primary and atmospheric cosmic - rayspectra with the BESS-TeV spectrometer”. Phys. Lett., B594:35–46, 2004.
[71] K. Niita et al. “PHITS: A particle and heavy ion transport code system”.
Radiat. Meas., 41:1080–1090, 2006.
[72] S. Roesler et al. “The Monte Carlo event generator DPMJET-III”. arXiv,
hep-ph/0012252:1033–1038, 2000.
[73] G. Battistoni et al. “The FLUKA atmospheric neutrino flux calculation”. Astropart. Phys., 19:269–290, 2003. [Erratum: Astropart.
Phys.19,291(2003)].
[74] G. D. Barr et al. “A Three - dimensional calculation of atmospheric
neutrinos”. Phys. Rev., D70:023006, 2004.
[75] J. A. Formaggio and G. P. Zeller. “From eV to EeV: Neutrino Cross
Sections Across Energy Scales”. Rev. Mod. Phys., 84:1307–1341, 2012.
[76] S. Jadach et al. “The tau decay library TAUOLA: Version 2.4”. Comput.
Phys. Commun., 76:361–380, 1993.

157

BIBLIOGRAPHY[77] T. J. Irvine. “Development of Neutron-Tagging Techniques and Application to Atmospheric Neutrino Oscillation Analysis in Super-Kamiokande”.
Ph. D. Thesis, University of Tokyo, 2014.
[78] C. H. Llewellyn Smith. “Neutrino Reactions at Accelerator Energies”.
Phys. Rept., 3:261–379, 1972.
[79] R. A. Smith and E. J. Moniz. “NEUTRINO REACTIONS ON NUCLEAR TARGETS”. Nucl. Phys., B43:605, 1972. [Erratum: Nucl.
Phys.B101,547(1975)].
[80] R. Gran et al. “Measurement of the quasi-elastic axial vector mass inneutrino-oxygen interactions”. Phys. Rev., D74:052002, 2006.
[81] A. A. Aguilar-Arevalo et al. “First Measurement of the Muon Neutrino
Charged Current Quasielastic Double Differential Cross Section”. Phys.
Rev., D81:092005, 2010.
[82] J. Nieves, I. Ruiz Simo, and M. J. Vicente Vacas. “Inclusive Charged–
Current Neutrino–Nucleus Reactions”. Phys. Rev., C83:045501, 2011.
[83] D. Rein and L. M. Sehgal. “Neutrino Excitation of Baryon Resonancesand Single Pion Production”. Annals Phys., 133:79–153, 1981.
[84] K. M. Graczyk and J. T. Sobczyk. “Form Factors in the Quark Resonance Model”. Phys. Rev., D77:053001, 2008. [Erratum: Phys.
Rev.D79,079903(2009)].
[85] C. Wilkinson et al. “Reanalysis of bubble chamber measurements of muonneutrino induced single pion production”. Phys. Rev., D90(11):112017,2014.
[86] D. Rein and L. M. Sehgal. “Coherent pi0 Production in Neutrino Reactions”. Nucl. Phys., B223:29–44, 1983.
[87] M. Glck, E. Reya, and A. Vogt. “Dynamical parton distributions revisited”. Eur. Phys. J., C5:461–470, 1998.
[88] A. Bodek and U. K. Yang. “Modeling neutrino and electron scatteringcross-sections in the few GeV region with effective LO PDFs”. AIP Conf.
Proc., 670:110–117, 2003. [,110(2003)].
[89] M. Derrick et al. “Properties of the Hadronic System Resulting fromanti-Muon-neutrino p Interactions”. Phys. Rev., D17:1, 1978.
[90] S. Barlag et al. “Charged Hadron Multiplicities in High-energy Anti-muon
Neutrino n and Anti-muon Neutrino p Interactions”. Z. Phys., C11:283,
1982. [Erratum: Z. Phys.C14,281(1982)].
[91] T. Sjostrand. “PYTHIA 5.7 and JETSET 7.4 physics and manual”. 1994.

158

BIBLIOGRAPHY[92] R. Brun et al. “GEANT3”. 1987.
[93] C. Zeitnitz and T. A. Gabriel. “The GEANT - CALOR interface andbenchmark calculations of ZEUS test calorimeters”. Nucl. Instrum. Meth.,
A349:106–111, 1994.
[94] R. B. Patterson and others. “The Extended-track reconstruction for MiniBooNE”. Nucl. Instrum. Meth., A608:206–224, 2009.
[95] M. Shiozawa. “Reconstruction algorithms in the Super-Kamiokande large
water Cherenkov detector”. Nucl. Instrum. Meth., A433:240–246, 1999.
[96] H. Nishino. “Search for Nucleon Decay into Charged Antilepton plus
Meson in Super-Kamiokande”. Ph. D. Thesis, University of Tokyo, 2009.
[97] S. Tobayama. “An Analysis of the Oscillation of Atmospheric Neutrinos”.
Ph. D. Thesis, University of British Columbia, 2016.
[98] S. Mine et al. “Experimental study of the atmospheric neutrino backgrounds for p → e+ π 0 searches in water Cherenkov detectors”. Phys.
Rev., D77:032003, 2008.
[99] Y. Ashie et al. “A Measurement of atmospheric neutrino oscillation parameters by SUPER-KAMIOKANDE I”. Phys. Rev., D71:112005, 2005.
[100] G. Mitsuka. “Study of Non-Standard Neutrino Interactions with Atmospheric Neutrino Data in Super-Kamiokande”. Ph. D. Thesis, Universityof Tokyo, 2009.
[101] K.
Abe
et al.
“Hyper-Kamiokande
DesignReport”.
KEK-PREPRINT-2016-21, ICRR-REPORT-701-2016-1, 2016.
[102] G. Cowan et al. “Asymptotic formulae for likelihood-based tests of
new physics”. Eur. Phys. J., C71:1554, 2011. [Erratum: Eur. Phys.
J.C73,2501(2013)].

159