Doctor thesis
Development of Neutrino Energy Reconstruction
Algorithm for Neutrino Mass Hierarchy Determination

Tokyo Institute of Technology, Department of Physics

Yuji Okajima

May 26, 2020

Abstract
From the discovery of neutrino oscillation, its mixing angles and squared diﬀerences between
three neutrino mass eigenvalues have been measured. However, the ordering of the neutrinomass eigenstates, which is called neutrino mass hierarchy, has not been decided yet. To decide
the mass hierarchy, the search for the oscillation probability enhancement owing to MSW eﬀect
is important. To detect this enhancement eﬃciently, the neutrino energy resolution is important.
In Super-Kamiokande, the neutrino interaction event is reconstructed based on the Cherenkov
light distribution and intensity. The Cherenkov light intensity depends on the particle charge
and track length. Because of the existence of strong interaction in the water, hadron’s track
length tends to be shorter than leptons. Because the particle momentum is reconstructed assuming it is lepton, the reconstructed hadron momentum tends to be smaller than the true value.
Especially, the momentum of neutral particles, for example neutron, cannot be reconstructedbecause they do not emit any Cherenkov lights. Therefore, the reconstructed neutrino energy
based on the Cherenkov light tends to be smaller than the true value. This deficit is called
the missing energy. Because the uncertainty of the missing energy makes the resolution of thereconstructed neutrino energy worse, its correction is required.
This thesis describes the correction method of this missing energy based on the numbers oftagged Michel electrons and neutrons. This correction improved the neutrino energy resolution
in the range from 4% to 26% for SK-I, II, III, and from 12% to 39% for SK-IV. This improved
energy resolution enables us to set the energy binning width for oscillation analysis narrower.
This new energy reconstruction method and binning improved the sensitivity of the neutrino oscillation. The sensitivity was evaluated by χ2 test between MC outputs based onsome candidates of oscillation parameters and two true patterns whose parameters are set
to (∆m212 , sin2 θ23 , sin2 θ12 , sin2 θ13 , δCP ) = (7.65 × 10−5 [eV2 ], 0.575, 0.309, 4.19), and ∆m223 =2.5 × 10−3 (−2.5 × 10−3 )[eV2 ] as examples of normal (inverted) hierarchy parameters. As aresult, the diﬀerence of minimum χ2 values between normal hierarchy and inverted hierarchy
candidates (|χ2min (N H) − χ2min (IH)|) was improved from 3.04 (1.68) to 4.05 (2.05) when thetrue hierarchy is normal (inverted) hierarchy.

iContents
1 Introduction
1.1 Author Contribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 Thesis Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1
1
2

2 Neutrinos
2.1 Neutrino Oscillation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 The Oscillation of Solar Neutrinos . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3 The Oscillation of Atmospheric Neutrino . . . . . . . . . . . . . . . . . . . . . . .
2.4 Mixing Parameter Measurements . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.4.1 θ12 , ∆m221 Measurement . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.4.2 θ23 , ∆m232 Measurement . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.4.3 θ13 , ∆m231 Measurement . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.5 Atmospheric Neutrino Oscillation in Matter . . . . . . . . . . . . . . . . . . . . .
2.6 Mass Hierarchy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.7 MSW Eﬀect in Atmospheric Neutrinos . . . . . . . . . . . . . . . . . . . . . . . .

3
3
4
5
7
7
7
8
8
10
11

3 Super-Kamiokande
3.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 SK periods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3 Cherenkov Radiation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4 Components of Super-Kamiokande Detector . . . . . . . . . . . . . . . . . . . . .
3.5 Water Tank . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.6 Photomultiplier Tubes for Inner Detector . . . . . . . . . . . . . . . . . . . . . .
3.7 Outer Detector Systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.8 Purification Systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.8.1 Air Purification System . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.8.2 Water Purification System . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.9 Data Acquisition Systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.9.1 Systems for SK-I, II, and III . . . . . . . . . . . . . . . . . . . . . . . . .
3.9.2 Electronics for SK-IV . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

15
15
16
16
18
18
19
19
21
22
22
23
23
26

4 Calibration
4.1 PMT Calibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.1.1 Gain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.1.2 Relative Diﬀerence of QE . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.1.3 Timing Calibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

28
28
28
31
33

ii

CONTENTS
4.2

iii

Calibration for Photon Propagation . . . . . . . . . . . . . . . . . . . . . . . . .
4.2.1 Water Transparency . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2.2 Light Reflection at PMT . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2.3 Light Reflection at Black Sheet . . . . . . . . . . . . . . . . . . . . . . . .
Energy Scale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3.1 Calibration Sources . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3.2 Diﬀerence Between MC and Data . . . . . . . . . . . . . . . . . . . . . . .
4.3.3 Time Variation of Energy Scale . . . . . . . . . . . . . . . . . . . . . . . .
4.3.4 Evaluation Result of Energy Scale Uncertainty . . . . . . . . . . . . . . .

33
33
40
40
40
41
45
47
47

5 Monte Carlo Simulation of Atmospheric Neutrino Events
5.1 Atmospheric Neutrino Flux . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2 Neutrino Interaction Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2.1 NEUT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2.2 GENIE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.3 Detector Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.3.1 Background Data Combination for Neutron Tagging Analysis . . . . . . .
5.3.2 π − Capture . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

49
49
52
52
58
58
60
60

6 Event Reconstruction Based on The Primary Signal
6.1 Vertex and Direction Fitting . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.1.1 Point Fit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.1.2 Ring Edge Search . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.1.3 TDC-Fit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.2 Ring Counting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.3 Particle Identification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.4 Precise Single Ring Fitting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.5 Momentum Reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

61
61
61
61
62
63
66
72
72

7 Event Selection
7.1 Fully Contained (FC) Events . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.1.1 First Selection (FC1, Cosmic Ray and Low Energy Events Cut) . . . . . .
7.1.2 Second Selection (FC2, Cosmic-ray Muon and Flasher Events Cut) . . . .
7.1.3 Third Selection (FC3, Cosmic Muon Cut and Low Energy Event Cut) . .
7.1.4 Fourth Selection (FC4, Flasher Cut) . . . . . . . . . . . . . . . . . . . . .
7.1.5 Fifth Selection (FC5, Remaining Background Cut) . . . . . . . . . . . . .
7.1.6 Fiducial Volume Selection . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.1.7 FC Data Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.2 Partially Contained Events . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.2.1 First Selection (PC1, Through Going Muon and Low Energy Events Cut)
7.2.2 Second Selection (PC2, Cosmic Ray Muon Cut) . . . . . . . . . . . . . .
7.2.3 Third Selection (PC3, PMT Flasher Event Cut) . . . . . . . . . . . . . .
7.2.4 Fourth Selection (PC4, Remaining Cosmic Ray Muon Event Cut) . . . .
7.2.5 Fifth Selection (PC5, Remaining Background Event Cut) . . . . . . . . .
7.2.6 Fiducial Volume Selection . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.2.7 PC Selection Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . .

75
75
75
76
77
80
80
84
84
85
85
86
87
87
90
94
94

4.3

CONTENTS
7.3

7.4

Upward-Going Muon Events . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.3.1 First Selection (UPMU1, Too High or Low Energy Event Cut) . . . . . .
7.3.2 Second Selection (UPMU2, Downward-Going Muon Event Cut) . . . . . .
7.3.3 Third Selection (UPMU3, Event Type Identification) . . . . . . . . . . . .
7.3.4 Fourth Selection (UPMU4, Corner Clipping Muon Event Cut) . . . . . .
7.3.5 Upward-Going Muon Event Type Criteria . . . . . . . . . . . . . . . . . .
7.3.6 Background Subtraction . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Selection Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

iv
95
95
95
96
96
97
97
99

8 Energy Correction Based on Secondary Interactions
100
8.1 Hadron Tagging Algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
8.1.1 Neutron Tagging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
8.1.2 Decay Electron Search . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
8.2 Missing Energy Correction Method . . . . . . . . . . . . . . . . . . . . . . . . . . 108
8.2.1 Function Type Separation . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
8.2.2 Energy Correction as a Function of the Primary Ring Momentum . . . . 109
8.2.3 Energy Bias Rejection from the Resonance Energy Region . . . . . . . . . 110
8.2.4 Energy Resolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
8.2.5 Systematic Error Estimation . . . . . . . . . . . . . . . . . . . . . . . . . 115
9 Oscillation Analysis
127
9.1 Analysis Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
9.1.1 χ2 Calculation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
9.1.2 FC Event Type Criteria . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128
9.1.3 Separation Method of Multi-GeV Multi-ring Electron Like Events . . . . 128
9.1.4 PC Event Types . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129
9.1.5 Binning Setting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129
9.2 Systematic Errors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129
9.2.1 Systematic Errors Related to Neutrino Flux . . . . . . . . . . . . . . . . . 129
9.2.2 Systematic Errors Related to Neutrino Interaction . . . . . . . . . . . . . 134
9.2.3 Systematic Errors Related to Event Selection and Reconstruction . . . . . 135
9.3 Sensitivity Study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136
10 Conclusion and Outlook

142

Acknowledgement

149

List of Figures
2.1

2.2

2.3

2.4

2.5

2.6

The comparison of B8 solar neutrino flux with µ and τ flavor and the electron
neutrino flux. The solid band shows the fluxes measured with the neutrino reactions in SNO, and the area between the broken straight lines are assigned to theflux predicted by the standard solar model[9]. . . . . . . . . . . . . . . . . . . . .
The comparison of zenith angle distribution between measured data, the expected
flux without oscillation and that with best-fit oscillation parameters for data sets
separated by the produced particle flavor and its energy. Except for the figure forthe partially contained events, all figures show the distribution for the fully contained events. The definition of fully-contained and partially contained events aredescribed in Sec. 7. Sub-GeV and multi-GeV events indicate the events with lessand more than 1.33 GeV visible energy, which is the energy reconstructed assuming the detected Cherenkov light to be produced by an electron. The Sub-GeV
e-like and µ-like fully-contained events are shown separately for the events withless and more than 400 MeV/c momentum. The multi-GeV e-like fully-contained
events are shown separately for the events with less and more than 400 MeV/cmomentum. The hatched areas are assigned to the expected distribution with nooscillation, the dots are assigned to the data, and the solid line are assigned tothe expectation based on the best-fit νµ → ντ oscillation parameters.[21] . . . . .
The comparison of neutrino mass ordering between normal and inverted hierarchy.
The flavor compositions of mass eigenstates are indicated by color bars. ∆m232 has
positive and negative values for the normal and inverted hierarchies, respectively.
The relationship between the lightest neutrino mass and the eﬀective Majorana
neutrino mass (mββ ) (left), and the upper limit of the eﬀective Majorana neutrinomass based on the measurement result of each nucleus decay as a function of themass number (right) [34]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Schema of the atmospheric neutrinos passing through the earth. The path lengths
in mantle and core of atmospheric neutrino can be estimated from zenith angleof the atmospheric neutrino detected in Super-Kamiokande. The core has larger
electron density than mantle, and it changes the oscillation probability owing to
the MSW eﬀect. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Oscillation probabilities for neutrinos and anti neutrinos when the mass hierarchy
is normal hierarchy. The left two figures show νµ and ν̄µ survival probabilities,
and the right two figures show νe and ν̄e appearance probabilities. These figures
assume (∆m232 , sin2 θ23 , sin2 θ13 ) = (2.5 × 10−3 , 0.5, 0.0219). These figures aretaken from [35]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

v

5

6

11

12

13

14

LIST OF FIGURES

vi

3.1
3.2

Super-Kamiokande detector . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
Schema of Cherenkov radiation. The angle of cone depends on the particle velocity
and refractive index. n is refractive index of medium, c is light speed in vacuum.
In the time interval t, Cherenkov light moves t(c/n) and the charged particle
moves βct. Because the opening angle of Cherenkov radiation from the particle
with relativistic momentum is almost constant, the PMT signal distribution shows
the ring shape. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
3.3 Schema of Super-Kamiokande and its components. . . . . . . . . . . . . . . . . . 18
3.4 Design of the photomultiplier tube in the inner detector [40] . . . . . . . . . . . . 19
3.5 The quantum eﬃciency of the inner detector PMT photocathode as a function oflight wavelength[40] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
3.6 The shape of wavelength shifter plate on the surface of the outer detector PMT[43]. 20
3.7 The left figure shows the outer detector PMT quantum eﬃciency as a function oflight wavelength, and the right figure shows the absorption (shorter wavelength
curve) and emission (longer wavelength curve) spectra of the wave shifting plate.
The absorption spectrum uses the molar extinction coeﬃcient µ as Y-axis on theleft. The absorption probability follows Pabs = 1 − 10−µρd , where ρ is the flour’sconcentration in mol/l, and d is the thickness of the material. The flour emission
spectrum uses the Y-axis on the right [43]. . . . . . . . . . . . . . . . . . . . . . . 21
3.8 Schematic views of the detector with the OD segmentation and its corner. As
shown in the right figure, the Cherenkov light which is produced in the barrel andendcap region of OD are separated to reject the corner clipping muon. . . . . . . 22
3.9 Water purification system[42] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
3.10 ATM structure[46] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.11 DAQ system of SK-I,II,III . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.12 DAQ system of SK-IV . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
4.1

4.2
4.3

4.4

4.5

Location of the standard PMTs in the inner detector (Left figure). The red point
indicates the standard PMT. The standard and normal PMTs are separated tosome groups based on the geometrical relationship to the light source (Right
figure). These figures are taken from [44]. . . . . . . . . . . . . . . . . . . . . . .
Pulse shape of the signal from each PMT. the signal is detected at the timing
when the voltage exceeds the threshold. . . . . . . . . . . . . . . . . . . . . . . .
The single photoelectron (p.e.) charge distributions in SK-III with diﬀerent charge
range (left: 5 pC, right: 50 pC). The green line in the left figure is the distribution
estimated based on the result of the measurement with half threshold and double
gain. The red line is the region approximated by linear function. This figure istaken from [44]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Normalized hit probability as a function of PMT position. the upper figure shows
the dependence on z (cm). Below left and right figures show the position dependence of top and bottom PMTs, where the x-axis is the square of the distancefrom the center to the PMT (cm2 ). The red plot is the data and the blue plot is
MC. This figure is taken from [44]. . . . . . . . . . . . . . . . . . . . . . . . . . .
Normalized hit probability comparison between PMTs used in SK-II and newly
installed from SK-III, similar to Fig. 4.4. The green line shows the newer PMTs,
and the red line shows the older one. This figure is taken from [44]. . . . . . . . .

29
30

31

32

34

LIST OF FIGURES
Schematic view of the timing calibration system. the light intensity can bechanged by a variable optical filter. This figure is taken from [44]. . . . . . . . .
4.7 The (time-of-flight corrected hit timing)-(charge) distribution. The sign of Yaxis value is reversed from the hit timing. As Y-axis value becomes larger, thehit timing becomes earlier. This figure is taken from [44]. . . . . . . . . . . . . .
4.8 Schema of water transparency calibration. The barrel PMTs are separated to
5 groups, which is labeled B1 to B5. By comparing the timing distributions of
PMTs in the top and these barrel regions, detector parameters are tuned. This
figure is taken from [44]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.9 TOF subtracted timing distributions for the laser beam with 405 nm wave length.
The red and black lines are assigned to the data and MC with best tune. The top
figure is for the PMTs at the top of the tank, and the second to bottom figures
is for the PMTs in B1 to B5 regions, respectively. The definition of B1 to B5
regions are written in Fig. 4.8. The time region between the earlier two blue lines
are regarded as hits from the scattered light, and the time region between later
two lines are regarded as the reflected light. This figure is taken from [44]. . . . .
4.10 Inverse water transparency based on April 2009 data. Parameters are fitted with
(4.9), (4.10), and (4.11). This figure is taken from [44]. . . . . . . . . . . . . . . .
4.11 The schema of the module for the black sheet reflection measurement[44] . . . .
4.12 Comparison of the black sheet reflection between MC and data. The top, middle, bottom figure show the measurement result with the 337, 400, and 420 nm
wavelength light sources, respectively. The left figure shows the charge ratio (redcircle: MC, black cross: data). The ratio between data and MC for the left plots
value is shown in the right plots.[44] . . . . . . . . . . . . . . . . . . . . . . . . .
4.13 The reconstructed invariant mass distributions for SK-IV period. The black dots
and red boxes are assigned to the data and MC. . . . . . . . . . . . . . . . . . .
4.14 The schema of the definition of the time regions for the cluster detection algorithm [49]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.15 The reconstructed momentum of decay electrons for SK-IV period. The blackand red lines are assigned to the data and MC, respectively. . . . . . . . . . . .
4.16 Diﬀerence of the energy scale between MC and data. The colors of data points
indicate the calibration sources, and each source’s comparison target is defined inSec. 4.3.2. The error bars to x-axis are assigned to the deviation of the momentum
distributions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.17 Time variation of reconstructed decay electron momentum. Dotted lines indicate±1% from the mean value in each period. The shaded regions in SK-IV indicate
T2K run period, in which neutrino beam comes from J-PARC. . . . . . . . . . .
4.18 Time variation of ratio of reconstructed stopping muon momentum to its track
length. Dotted lines indicate ±1% from the mean value in each period. The
shaded regions in SK-IV indicate T2K run period, in which neutrino beam comesfrom J-PARC. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

vii

4.6

5.1

The azimuthal angle dependence of atmospheric neutrino flux at 1 GeV. This
figure is taken from [50]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

35

36

37

38
39
40

41
42
43
44

46

47

48

50

LIST OF FIGURES

viii

Comparisons of µ± fluxes between the data taken at Tsukuba and the calculation
based on the modified JAM [50]. The left figure shows their momentum distribution, and the right one shows the ratio of data to MC. The red dots solid lines
in the left figure show the observed and calculated µ+ fluxes, and the blue ones
indicate their µ− fluxes, respectively. The dots in the right figure show the ratioof data to the introduced model, which is based on JAM, and broken lines show
the ratio of this model and DPMJET-III. . . . . . . . . . . . . . . . . . . . . . . 50
5.3 Comparisons of µ± fluxes between the data taken at balloon altitude at FortSumner and the calculation based on the modified JAM [50]. The meaning of dots
and lines is the same as that of Fig. 5.2. . . . . . . . . . . . . . . . . . . . . . . . 51
5.4 Atmospheric neutrino fluxes based on the calculation result of Honda used for SK
Monte Carlo, older calculation result of Honda (HKKMS06), the model of Bartol
group[67], and that of Fluka[68]. This plot is taken from[50]. . . . . . . . . . . . 51
5.5 Comparison of the charged current quasi-elastic scattering cross sections between
data and NEUT[74]. The solid and broken lines are assigned to the calculated
cross sections for the interactions with free nucleons and bound nucleons in oxygen, respectively. These figures are taken from [74] . . . . . . . . . . . . . . . . . 53
5.6 Comparison of the cross sections of single pion productions between data andNEUT[74]. The solid red and black lines indicate the cross sections based on
NEUT5.3.6 and older version of NEUT, respectively. These figures are written
based on [74] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
5.7 Comparison between charged current inclusive cross sections summarized by PDG[84]and the cross section estimated based on NEUT 5.3.6. The red line is assigned
to the cross section based on NEUT5.3.6. . . . . . . . . . . . . . . . . . . . . . . 55
5.8 Comparison of the coherent π production cross sections between results of theexperiments and models. The left and right figures show the cross sections of thecharged and neutral current interactions with π production, respectively. The
red and blue points in the left figure are assigned to the measurement results ofK2K[86] and SciBooNE[87]. The circle, box, and cross points in the right figure
are assigned to the measurement results of MiniBooNE[88], Aachen-Padova[89],and Gargamelle[90]. The solid and broken lines in both figures are assigned tothe cross sections based on Rein-Sehgal model with and without lepton mass
correction, respectively. The dotted and dash-dotted lines in the left figure areassigned to the cross section based on other models, which are the models of
Kartavtsev et al.[91] and Alvarez-Ruso et al.[92], respectively. These figures aretaken from [74]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
5.9 π + and π − interaction cross sections on 12 C[96]. . . . . . . . . . . . . . . . . . . 57
5.10 The cross sections of coherent pion production as a function of the transmitted
momentum [113]. The black and blue lines are assigned to the output of ReinSehgal and Berger-Sehgal models, respectively. . . . . . . . . . . . . . . . . . . . 58
5.2

LIST OF FIGURES
5.11 The cross section comparison of π + ’s absorption and charge exchange with carbon particles [114]. The red and black colors are assigned to the cross sectionsof absorption and charge exchange, respectively. The filled and hollow dots are
assigned to DUET and the other experimental data[116][117][118][119], respectively. The solid, broken, and dotted lines are assigned to the model prediction
based on GEANT4, FLUKA and NEUT before the tuning based on the DUET
data in this figure, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.12 The comparison of the number of produced neutrons between models. Blue, redand green lines are assigned to the calculation results of NEUT5.3.6, GENIE2.12.2(with hA model), and NEUT5.4.0, respectively. The average numbers of the
produced neutrons are largely diﬀerent between the models. . . . . . . . . . . . .
5.13 The structure of the background for MC as a function of the elapsed time. The
dummy trigger data is introduced from 18−ttrg [µs] to 535−ttrg [µs] after the SHE
trigger for the neutrino interaction. ttrg is the triggered time from the start of thesimulation, which is about 1 µs. Until 18 − ttrg [µs], the background is simulated
assuming it comes from the dark noises of the PMTs, which are independent ofeach other. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.1
6.2

6.3

6.4

6.5
6.6

6.7

6.8

7.1

ix

59

59

60

Schema of Cherenkov ring opening angle distribution with θpeak and θedge positions. 62
The schema of the ring search idea. By setting each PMT position to conical vertex, the cones with 42◦ opening angle with antiparallel direction with Cherenkov
ring direction are generated. By searching for the position where the edges ofthese cones are concentrated. This figure is taken from [121]. . . . . . . . . . . . 64
A map of the charge weighted based on the opening angle dependence of theCherenkov light. Two peaks in this figure are regarded as the candidates of theCherenkov ring. This figure is taken from [121]. . . . . . . . . . . . . . . . . . . . 65
The distribution of ring counting likelihood diﬀerence between 1-ring and 2-ringfor multi-GeV events. The black dots and red line are assigned to the distributionsof the data and MC with the same live time. . . . . . . . . . . . . . . . . . . . . 67
Schema of Cherenkov light emission area. . . . . . . . . . . . . . . . . . . . . . . 69
Particle identification likelihood distribution of single ring multi-GeV even of SKIV. The black dots and red line are assigned to the distributions of the data andMC with the same live time. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
Particle identification likelihood distribution of all the rings in multi ring multiGeV events. The black dots and red line are assigned to the distributions of thedata and MC with the same live time. . . . . . . . . . . . . . . . . . . . . . . . . 71
The distributions for SK-IV multi-GeV single ring electron like FC events basedon MC. (Left:The distance between true and reconstructed vertex. Right: The
opening angle between true and reconstructed direction.) The green regions contain 68% of all events in the distribution, and its upper limit is defined as theirresolutions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
The schema of the target events of each selection types. The solid arrows andbroken lines show the behaviors of charged particles and neutrinos, respectively. .

76

LIST OF FIGURES
7.2

7.3

7.4

7.5

8.1

8.2

8.3

8.4

The schema of the cluster detection algorithm [121]. The size of circles represents
the total charge observed in each patch. The arrow shows the vector of thecharge gradient which head for the patches with the maximum charge in theneighborhood patches. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
The schema of the candidate selection. The yellow circles are assigned to thecandidate PMTs selected by the first process of the PMT selection, and the pink
circles are assigned to the PMTs with hits. The number of hits in the candidate
itself and eight of the nearest PMTs, which the red and blue boxes are assigned
to, are required to be the largest in the detector. . . . . . . . . . . . . . . . . . .
The scatter plots of azimuthal and zenith angles for stopping muon (left plot)
and non-showering (middle plot) and showering (right plot) through-going muon
events based on the data. The particles with θ = 0 head for vertically upperdirection. The distribution of downward and horizontal cosmic-ray muons seemto have some clusters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
The zenith angle distributions of stopping muon and showering and non-showering
through-going muon events in region 1 and 2. The particles with θ = 0 head forvertically upper direction. The black and white points show the distributions inregion 1 and region 2. The solid curves show the region 2’s fitted functions. . . .

x

81

88

97

98

Schema of the energy reconstruction for ν + p → (lepton) + n based on thecurrent and new method and their energy spectrum with and without MSW
eﬀect. New method corrects the missing energy taken by hadrons, and it isexpected to improve the energy resolution. It is expected that the width of theregion where the flux enhancement from the MSW eﬀect becomes narrower, as theenergy resolution becomes better due to the improvements in energy correction
method. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
Missing energies and their fitting function of events with one decay electron andno tagged neutrons based on MC. In this analysis, the events are separated tosome event types. The data points show the median value of the missing energy
distribution in their primary ring momentum bins. . . . . . . . . . . . . . . . . . 101
Output of the neural-network for the data and MC of atmospheric neutrino events.
The dashed line shows the threshold to separate the neutron capture event andthe background events. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
The schema of the event gates. The event gates are produced when any triggers
described in Table 3.3 and Table 3.4 are satisfied. The width of the event gates
until SK-III and from SK-IV are set to 30 µs and 40 µs from the start of the
primary event gate, respectively. This event gates contain one primary event gate
and some sub-event gate. The primary event gate is defined as the gate which
passes the event selection in Sec. 7.1, and the sub-events gate are the gate which
passed the threshold of the global trigger. The width of the primary-event andsub-event gates are 1.3 µs. The data in shaded region are not used. Until SK-III,
only the region where the primary and sub-event gates contain are registered. . 105

LIST OF FIGURES
8.5

xi

Energy bias of SK-IV charged current events in multi ring electron like event
group as a function of the primary ring momentum (right figure), and the bias asa function of true neutrino energy (left figure). Black dots and lines are assigned
to the median values of each energy region. The regions between two broken lines
are the resonance region. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
8.6 Energy biases of SK-IV charged current events with the target neutrino in eachfunction group as a function of the primary ring momentum. Red, black, andgreen line shows the bias of Evis,ID , Eν,corr1 , Eν,corr2 , respectively. These figurescontain not only the events with the energy in the resonance region, but also outof the region. Because the correction from Eν,corr1 to Eν,corr2 takes only the eventswith the energy in the resonance region into account, Eν,corr2 shows larger bias
than Eν,corr1 in this figure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
8.7 Energy biases of SK-IV charged current events with the target neutrino in eachfunction group as a function of the true neutrino energy. The target neutrino andthe line colors definitions is the same as Fig. 8.6. The region between two brokenlines is the resonance region. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112
8.8 Energy biases of SK-IV charged current events with the target neutrino in eachfunction group as a function of the number of tagged neutrons. The target neutrino and the line colors definitions is the same as Fig. 8.6. . . . . . . . . . . . . . 113
8.9 Energy biases of SK-IV charged current events with the target neutrino in eachfunction group as a function of the number of tagged decay electrons. The target
neutrino and the line colors definitions is the same as Fig. 8.6. . . . . . . . . . . . 114
8.10 Captured neutron multiplicity of the events in the resonance energy region basedon the toy MC output. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
8.11 Generated muon multiplicity of the events in the resonance energy region basedon the toy MC output. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118
8.12 Reconstructed energy ratios of the SK-IV outputs based on NEUT 5.4.0 and
GENIE 2.12.6 to the SK-IV output based on NEUT 5.3.6 output. The red lines
show Evis,ID = 1.33 GeV. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
8.13 Reconstructed energy ratios of the SK-IV events with high (23.9%) and low(19.5%) tagging eﬃciencies to the SK-IV events with normal tagging eﬃciency(21.7%). The red line definition is the same as Fig. 8.12. . . . . . . . . . . . . . . 121
8.14 The cross section of π + − O interactions based on data and theoretical calculation
which is used in MC. The experimental data is taken from [131]. This figure istaken from [121] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122
8.15 The π + decay probability as a function of its momentum[121]. The black solidline shows the π + decay probability in our MC simulation output, and the red
dotted and blue dashed lines show the probabilities with the longer and smaller
mean free path of π + . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123
8.16 The π + decay probability as a function of π + momentum. The probability falls
from 0.1 GeV to 1 GeV because the cross section of π + absorption on oxygen
have a peak around 200 MeV/c. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123
8.17 Reconstructed energy ratios of the events with high and low π + decay probabilities
to the events with normal probability for SK-IV. The red line definition is thesame as Fig. 8.12. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124

LIST OF FIGURES

xii

8.18 The energy spectrum of the γ-ray produced by the π − capture by 12 C nucleus[135]. The histogram colored in magenta shows the data, and the other color lines
show the expected energy spectra based on the diﬀerent models. The systematic
error from π − capture is estimated by using the energy spectrum of the data itself
and the diﬀerence between the energy spectra of the data and SF+RPA model. . 125
8.19 Reconstructed energy ratios of the events based on the normal MC to the events
based on the MC which has the neutron emission process from the π − capture.
The red line definition is the same as Fig. 8.12. . . . . . . . . . . . . . . . . . . . 126
9.1

9.2

9.3

9.4

9.5

9.6

9.7

Schema of the energy binning for oscillation analysis. The energy region in thewhite boxes are separated to 10 bins based on the zenith angle. On the other
hand, the gray boxes have only one bin. . . . . . . . . . . . . . . . . . . . . . . . 130
Binning for 1-ring Multi-GeV events and Multi-ring e-like events and the oscillation probability from νµ to νe based on normal hierarchy. The oscillation
probability table is taken from [35]. . . . . . . . . . . . . . . . . . . . . . . . . . . 131
The upper figure shows relative uncertainties of the error sources for the atmospheric neutrino flux[137]. The total uncertainty (δtot ) includes the eﬀects of all
the error sources in this figure, which are the pion and kaon production (δπ andδK ), hadronic interaction cross section (δσ ) and the atmospheric density profile(δair ). The lower figure shows the uncertainty of the absolute normalization whichis calculated based on the sum of δπ and δσ [136]. . . . . . . . . . . . . . . . . . . 132∆χ2 (≡ χ2 − χ2min ) distributions as a function of ∆m223 , sin θ23 and δCP whose
true hierarchy is the normal hierarchy. The red and blue lines are assigned to thesensitivities of the current and new analysis method, respectively. The solid andbroken lines are assigned to the sensitivities of the normal and inverted hierarchy
candidates, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138∆χ2 (≡ χ2 − χ2min ) distributions as a function of ∆m223 , sin θ23 and δCP whose truehierarchy is the inverted hierarchy. The red and blue lines are assigned to thesensitivities of the current and new analysis method, respectively. The solid andbroken lines are assigned to the sensitivities of the normal and inverted hierarchy
candidates, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139
Reconstructed Energy spectrum of upward-going (−1 < cos Θ < −0.9) 1-ringMulti-GeV νe like samples. The left figure shows the spectrum based on thecurrent binning setting and reconstruction method, and the right figure shows
the spectrum based on the new binning and method. . . . . . . . . . . . . . . . . 140
χ2 distribution of upward-going(−1 < cos Θ < −0.9) 1-ring Multi-GeV νe likesamples. χ2 in this figure is calculated by setting En and On following to theparameters in Table 9.2 with normal and inverted hierarchies, respectively. The
blue and red lines are assigned to the distributions based on current and newreconstruction method and binning, respectively . . . . . . . . . . . . . . . . . . 141

List of Tables
3.1
3.2
3.3

3.4
4.1
4.2
4.3

6.1

7.1

7.2

7.3
7.4
7.5
7.6
7.7

8.1

SK phases and the features . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Feature of the inner detector PMTs . . . . . . . . . . . . . . . . . . . . . . . . .
Global trigger thresholds. Because the width of a rectangular pulse is 200 ns, thisthreshold corresponds to the number of hits in 200 ns. The pulse height is 11 mV
per hit. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Triggers threshold and event width . . . . . . . . . . . . . . . . . . . . . . . . . .

16
19Absolute gain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Fitted parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Systematic error table about energy scale. Total energy scale errors are calculated
combining the absolute energy scale and time variation errors in quadrature. . .

31
39Vertex and angular resolutions of fully contained single ring events after the precise single ring fitting (MS-Fit). . . . . . . . . . . . . . . . . . . . . . . . . . . .

24
27

47

73

The eﬃciencies of neutrino events estimated with MC atmospheric neutrino events.
To evaluate the event rate, the distance between the detector wall and true andreconstructed vertices are defined as dwallv and dwall, respectively. From FC1
to FC5, the estimation is done by using only the neutrino events that satisfy
dwallv>200 cm, Evis > 30 MeV and NHITAC<16 hits. FCFV means the events
which pass all FC selections and fiducial volume selection. . . . . . . . . . . . . . 84
The background contamination rate table of the single ring event types. Thisrate is based on the result of the eye-scan. Sub-GeV and Multi-GeV events
satisfy Evis <1.33 GeV and Evis ≥1.33 GeV, respectively. . . . . . . . . . . . . . 85
The table of the selection criteria applied to SK-I and SK-II period and SK-III
and SK-IV period. ⃝ and △ indicate hard and soft cut criteria, respectively. . . 90
The remaining event rate at each selection step. . . . . . . . . . . . . . . . . . . . 94
Signal eﬃciencies, background rates and event rates after the PC selection. The
background contamination rates are estimated based on the eye-scanning result[124]. 95
The remaining event rate at each of the upward-going muon selection steps. . . . 99
The number of remaining events after the upward-going muon selections, and thenumber of background events to be subtracted. . . . . . . . . . . . . . . . . . . . 99
variables for neutron tagging based on neural-network analysis. The detailed
description is written in [1]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103

xiii

LIST OF TABLES
8.2

8.3
8.4
9.1
9.2
9.3
9.4
9.5
9.6

xiv

Eﬃciency and mis-tagging rate before and after the neural-network selection. The
contamination per events means the number of mis-tagging neutron candidates
in a neutron searching window, which have 517 µs width. . . . . . . . . . . . . . 104
Correspondence table of event type and function group. The definition of theevent types is described in Sec. 9. . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
Energy resolution of each event type for each period. . . . . . . . . . . . . . . . . 115
The event types and their separation criteria. . . . . . . . . . . . . . . . . . . . . 128
The table of ”true” parameters and Fixed/Fitted types for the oscillation sensitivity study. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137
Oscillation sensitivities (in χ2 diﬀerence) when the true oscillation parameters
follows Table 9.2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137
Allowed region of ∆m232 at 68% and 90% C.L., where 68% (90%) C.L. is definedas χ2 − χ2min = 1.0(2.7). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140
Allowed region of sin2 θ23 at 68% and 90% C.L., where 68% (90%) C.L. is definedas χ2 − χ2min = 1.0(2.7). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140
Allowed region of δCP at 68% and 90% C.L., where 68% (90%) C.L. is defined as
χ2 − χ2min = 1.0(2.7). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141

Chapter 1

Introduction
The theories that describe the elemental particles and their interactions are summarized tostandard model. If particles are produced following to standard model, almost the same numberof anti-particles is expected to be produced. However, the universe contains the large numberof baryons and few anti-baryons in reality. To solve this problem, physics beyond the standard
model have to be considered. Neutrino oscillation is also occurred by the physics beyond thestandard model. It is sensitive to the lepton CP violation, neutrino mixing angles, and thesquared diﬀerences of neutrino invariant masses. The CP violation is considered to contributeto the baryon asymmetry[4]. Based on the various measurement results, these parameters havebeen measured. However, some unsolved problems still remain. The neutrino mass hierarchy
problem is an example. The determination of the hierarchy is important for Majorana neutrino
search, δCP measurement, and so on. The atmospheric neutrino oscillation, which is sensitive
to the mass hierarchy, is measured by Super-Kamiokande. To identify the mass hierarchy, themethod to improve the sensitivity has been considered from many aspect[1][2][3]. This study
contributed as a part of this improvement.

1.1

Author Contribution

To make the energy reconstruction more precise and improve the sensitivity for the mass hierarchy, I have developed a new neutrino energy reconstruction algorithm. In this thesis, the resultof sensitivity study with this new energy reconstruction algorithm is described.
In addition to the works written in this thesis, author also contributed for the calibrationin Super-Kamiokande. To check the time variation of the background which related to themis-tagging rate of neutrons, new calibration process was introduced. This process measures
PMT signals without any neutrino events or cosmic-ray muons periodically. In addition, author
have calibrated the neutron tagging eﬃciency in 2016. The neutron tagging eﬃciency had been
already measured in 2009 with Am/Be neutron source[1]. I have measured the eﬃciency tocheck if the time variation of neutron tagging eﬃciency is negligible or not. As a result, thediﬀerence between tagging eﬃciencies measured in 2009 and 2016 was within error margin.

1

CHAPTER 1. INTRODUCTION

1.2

2

Thesis Overview

Chapter 2 will describe basic physics about neutrino oscillation. The method how the neutrino
oscillation parameters are measured and the importance of the mass hierarchy determination
are also included. Chapter 3 contains the description of Super-Kamiokande, which consists fromthe method how to detect the signal from neutrino interaction, the photo-multiplier tube, watertank, water and air purification systems, and data acquisition systems. Chapter 4 will describe
the calibration methods of these systems. Chapter 5 will describe the Monte Carlo simulationof atmospheric neutrino events, which is used for the comparison with the data. Chapter 6
will describe the method how to classify the candidate events and reject the background events.
Chapter 7 will describe the event reconstruction method based on the signal from the primaryparticles. Chapter 8 will describe the new method of the missing energy correction based on
the number of hadrons. This chapter will also describe how to detect the hadrons. Chapter 9
describes the method how to measure the oscillation parameters and mass hierarchy, and theresult of sensitivity study.

Chapter 2

Neutrinos
2.1

Neutrino Oscillation

Neutrino oscillation is a phenomenon that the neutrino flavors, which are the eigenstates of
weak interactions, are changed. It happens because neutrino flavors are not the same as mass
eigenstates. The state of each neutrino flavor να can be written as follows:∑|να ⟩ =
Uαi |νi ⟩(2.1)
iwhere U is the mixing matrix and νi is neutrino mass eigenstate. From the measurement resultof Z boson resonance[5], the number of lepton flavors is assumed to be three, and other flavor
leptons with weak interaction do not exist. Based on this, U is defined as the following 3 × 3
matrix:




1
0
0
cos θ130 sin θ13 e−iδCPcos θ12 sin θ12 0 − sin θ12 cos θ12 0
010
U = 0 cos θ23 sin θ23  iδ
0 − sin θ23 cos θ23001− sin θ13 e CP 0
cos θ13(2.2)
where θ12 , θ23 and θ13 are called mixing angles, and δCP is the charge-parity(CP) symmetry
violating phase. When sin δCP ̸= 0, the CP symmetry is violated. This matrix is known as
Pontecorvo-Maki-Nakagawa-Sakata mixing matrix [6]. Based on this matrix, the time propagation of the neutrino based on Schrödinger equation can be written as follows:|να (t)⟩ =∑Uαi |νi (t)⟩ =

i

3
∑
i=1

=

∑

Uαi e−iEi t |νi (0)⟩√2 2Uαi e−i p +mi t |νi (0)⟩

(2.3)

i

where p and mi are the neutrino momentum and i-th mass eigenvalue, and t is the elapsed time.
Assuming the neutrino propagates nearly at light speed, να → νβ oscillation probability can be3

CHAPTER 2. NEUTRINOS

4

written as follows:
P (να → νβ ) = |⟨νβ |να (t)⟩|2 =

∑∑
i

=

∑∑
i

≃
=

√2

−i(

p +m2i −√p2 +m2j )t∗ −i((p+mi /2p)−(p+mj /2p))t∗Uαi Uβj Uαj
e
Uβi
2

2

(2.4)

j

∑∑
i

∗
∗
Uβi
Uαi Uβj Uαj
e

j

∑∑
i

j

∗
∗ −i(Ei −Ej )t
Uβi
Uαi Uβj Uαje∗∗ −i(mi −mj )t/2p
Uβi
Uαi Uβj Uαj
e
2

2

j

In standard model, neutrinos are considered to be massless particles. In this case, the mass
eigenstates are degenerated to only one state and neutrino oscillation does not happen. However,
the existence of neutrino oscillation had been proved by the studies in Super-Kamiokande[7],SNO [8] [9] [10] and the other experiments. These results suggest physics beyond the standard
model.

2.2

The Oscillation of Solar Neutrinos

The nuclear fusion interaction in the Sun produces electron neutrinos with O(10) MeV or lessenergy. These neutrinos are called solar neutrinos. In 1960s, the Homestake experiment measured the solar neutrino flux by using a detector with 615 tons of tetrachlorothylene (C2 Cl4 ).
This detector measured the neutrino flux by detecting the quantity of 37 Ar gas produced by theinteraction below:
37
Cl + νe → 37 Ar + e−
(2.5)
However, the measured neutrino flux was only one third of the predicted flux based on the Standard Solar Model (SSM) [13]. This result agreed with that of Kamiokande[14], GALLEX[15],and SAGE[16]. This diﬀerence from the prediction had been called solar neutrino problem. This
problem is caused by the neutrino oscillation. Super-Kamiokande searched for the signature of
the oscillation based on the measurement of elastic scattering between the neutrinos and electrons in the water[17]. Elastic scattering is mainly sensitive to electron neutrinos. Therefore, if
the electron neutrino is oscillated to the other neutrinos, the number of elastic scattering events
will decrease. As a result, its measured solar neutrino flux was much smaller than the predictionof SSM without neutrino oscillation. On the other hand, Sudbury Neutrino Observatory (SNO)
also measured the neutrino oscillation based on the interaction between heavy water[8][9]. SNO
searched for the solar neutrino produced by 8 B decay in pp-chain(8 B →8 Be∗ + e+ + νe ). Owing
to D2 O liquid in its detector, this experiment could detect not only the elastic scattering, butalso the neutral and charged current interactions as follows:
Charged current: νe + d → p + p + e−

(2.6)

Neutral current: νx + d → p + n + νx

(2.7)
(2.8)

The charged current interaction occurs only for the electron neutrino. On the other hand,the neutral current interaction and elastic scattering occur for all neutrino flavors. Fig. 2.1

CHAPTER 2. NEUTRINOS

5

Figure 2.1: The comparison of B8 solar neutrino flux with µ and τ flavor and the electron
neutrino flux. The solid band shows the fluxes measured with the neutrino reactions in SNO,
and the area between the broken straight lines are assigned to the flux predicted by the standard
solar model[9].
shows the µ and τ neutrinos flux as a function of the electron neutrino flux deduced from those
three interactions. Combining this result with the flux measured with the elastic scatteringin Super-Kamiokande, the total flux of µ and τ flavor components was evaluated to ϕµτ =62 −1
3.45+0.65−0.62 [10 cm s ], which is 5.3σ above zero. In this evaluation, the solar neutrino energy
spectrum is assumed to follow the standard spectrum shape of 8 B decay[23]. In addition, the total flux of 8 B solar neutrinos measured with NC interaction in SNO was evaluated allowing the so+0.55
+1.57
6
2 −1
lar neutrino energy spectrum distortion. The result is ϕSNO
NC = 6.42−1.57 (stat)−0.58 (syst) [10 cm s ],
+1.01
6
2
−1
which is consistent with the prediction based on SSM model (ϕSSM = 5.05−0.81 [10 cm s ])[18].
Based on these results, the existence of solar neutrino oscillation was proved. Soon after, theKamLAND experiment[24] observed reactor neutrino oscillation caused by the same oscillation
parameter (See Sec. 2.4.1).

2.3

The Oscillation of Atmospheric Neutrino

When a cosmic ray particle approaches the earth, it interacts with the atmosphere, and producescharged mesons (π ± , K ± ). They decay into leptons. For example, π + decays as follows:
π + → µ + + νµ
→ e+ + νµ + ν̄µ + νe

(2.9)

Based on the interaction of Eq. (2.9), the neutrino flavor ratio R = (νµ + ν̄µ )/(νe + ν̄e ) is expected
to be about two. However, the detected neutrino flavor ratio in Kamiokande was much smaller
than the prediction. When the detected flavor ratio (Rdata ) is compared with the ratio estimatedCHAPTER 2. NEUTRINOS

6

based on the MC (RMC ), their ratio was estimated to be R = Rdata /RMC = 0.60+0.07−0.06 ± 0.05%
on average[19], depending on the incident direction and energy[20]. This diﬀerence is occurredbecause atmospheric neutrino is also oscillated as well as solar neutrino. This oscillation existence
was proved by Super-Kamiokande study[21]. Fig. 2.2 shows the zenith angle distribution of the

Figure 2.2: The comparison of zenith angle distribution between measured data, the expected
flux without oscillation and that with best-fit oscillation parameters for data sets separated by
the produced particle flavor and its energy. Except for the figure for the partially contained
events, all figures show the distribution for the fully contained events. The definition of fullycontained and partially contained events are described in Sec. 7. Sub-GeV and multi-GeV
events indicate the events with less and more than 1.33 GeV visible energy, which is the energy
reconstructed assuming the detected Cherenkov light to be produced by an electron. The SubGeV e-like and µ-like fully-contained events are shown separately for the events with less andmore than 400 MeV/c momentum. The multi-GeV e-like fully-contained events are shown
separately for the events with less and more than 400 MeV/c momentum. The hatched areas
are assigned to the expected distribution with no oscillation, the dots are assigned to the data,
and the solid line are assigned to the expectation based on the best-fit νµ → ντ oscillation
parameters.[21]atmospheric neutrino events. As shown in this figure, the flux of µ-like events largely depends
on the neutrino zenith angle asymmetrically, and the number of upward-going µ-like events ismuch smaller than the expectation with no-oscillation. This discrepancy between the zenith
angle distributions based on MC and data is the evidence of the existence of the νµ oscillationto ντ . Because the atmospheric neutrinos are produced at the surface of the earth, the neutrino
travel length changes from O(10) km to O(104 ) km depending on the zenith angle of the neutrino
incident direction to the detector. The neutrino oscillation probability depends on the travel
length as written in Eq. (2.28). Therefore, the magnitude of the discrepancy between the fluxes
based on MC and data depends on the zenith angle. Currently, Not only the νµ disappearance
measurement, but also ντ appearance is measured for the neutrino oscillation search. The

CHAPTER 2. NEUTRINOS

7

hypothesis of no-tau-appearance was excluded with a significance of 4.6σ[22].

2.4

Mixing Parameter Measurements

Based on the results of the various experiments, the mixing parameters and complex phaseof PMNS matrix (θ12 , θ13 , θ23 , δCP ) and the squared diﬀerences between three neutrino mass
eigenvalues (∆m221 , ∆m231 , ∆m232 , where ∆m2ij ≡ ∆m2i − ∆m2j ) are evaluated. In this section,
the measurement methods of these parameters are described.

2.4.1

θ12 , ∆m221 Measurement

Because the oscillation probability of electron neutrinos to the other flavor neutrinos dependson θ12 , ∆m221 , these parameters are evaluated based on the results of the solar neutrino and long
baseline reactor neutrino experiments, which produce electron neutrinos. The solar neutrino
experiment has been already described in Sec. 2.2. The reactor neutrinos are produced by the βdecays in the nuclear fission decay chain in reactors. In this decay chain, some ν̄e s with O(1) MeV
are produced. The reactor neutrino measurements search for the oscillation of these neutrinos.
KamLAND is a reactor antineutrino measurement experiment with a long baseline (typically
180 km)[24]. This experiment uses a 1 ton liquid scintillator detector to detect the reactor ν̄evia ν¯e + p → e+ + n, and measures the deficit of the reactor ν̄e . The produced neutrons are also
detected by using the γ-ray from the neutron capture on a hydrogen (n + p → d + γ). Owing
to this neutron tagging, the solar νe can be rejected from the detected reactor antineutrino
candidates. In the global analysis result summarized by PDG in 2016[25], the θ12 and ∆m221
values with 3σ errors are estimated as follows:2−5∆m221 = 7.37+0.59−0.44 × 10 [eV ]2sin θ122.4.2= 0.297+0.057−0.047(2.10)(2.11)

θ23 , ∆m232 Measurement

The oscillation between νµ and ντ largely depends on θ23 , ∆m232 . These parameters are estimated based on the number of νµ (ν¯µ ) of atmospheric or accelerator neutrinos. The oscillationof atmospheric neutrinos has been already described in Sec. 2.3. The accelerator neutrinos areproduced by the collision between a target and accelerated protons. K2K and MINOS experiments are its examples. K2K is the first accelerator based experiment for neutrino oscillation. It
uses the neutrinos with about 1.3 GeV on average which are produced by the collision between
12 GeV proton beam and a target. The neutrinos travelled over a distance of 250 km, and weredetected by the Super-Kamiokande water Cherenkov detector. This experiment could detect
the oscillation if ∆m2 is O(10−3 ) eV2 , which corresponds to ∆m232 . K2K ran for a period of 5
years from 1999. MINOS is the experiment which uses the neutrinos produced by the 120 GeV
proton beam at Fermilab. These neutrinos traveled over a distance of 735 km, and they weredetected by the 5.4 kton iron-scintillator in the Soudan mine[26]. Both K2K and MINOS detected the deficit in the observed νµ events, which comes from the νµ oscillation to the otherflavor neutrinos[27][28]. Based on their deficits, the nonexistence of the neutrino oscillation was
rejected in a diﬀerent way from the atmospheric neutrino observation. In the global analysis
result summarized by PDG in 2016[25], the θ23 and ∆m232 values with 3σ errors are evaluatedCHAPTER 2. NEUTRINOS

8

as follows:|∆m232 | = (2.54 ± 0.12) × 10−3 [eV2 ](2.12)

2= 0.425+0.190−0.044 (∆m32 > 0)2
= 0.589+0.047−0.205 (∆m32 < 0)(2.13)

2sin θ23

2.4.3(2.14)

θ13 , ∆m231 Measurement

As well as θ12 , ∆m221 , the oscillation probability from electron neutrino to the other flavor neutrinos depends on θ13 , ∆m231 . Because θ13 is smaller than θ12 , and ∆m231 is much larger than ∆m221 ,the electron neutrino source with smaller L/E than ∆m12 measurement is required. The shortbaseline reactor antineutrino and long baseline accelerator neutrino is suitable for this requirement. Comparing to the long baseline reactor antineutrino experiment in KamLAND, which is
described in Sec. 2.4.1, these experiments have smaller L/E value, owing to their shorter baseline and higher energy. Double Chooz[29], Daya Bay[30], and RENO[31] experiments measure
the short baseline reactor antineutrinos. Their measurement methods are basically the sameas the long baseline reactor antineutrino measurement, which is described in Sec. 2.4.1. Theyuse gadolinium-doped liquid scintillator. Gadolinium captures neutrons which are produced byν̄e + p → e+ + n, and reject the background νe events like solar neutrinos. These experiments
measure the ν̄e disappearance. Because reactor produces large amount of neutrinos, θ13 , ∆m231
can be estimated with small statistical error. In 2012-2013, These experiments reported sin2 2θ13
is about 0.1[29][30][31]. Especially, Daya Bay excluded θ13 = 0 at 5.2σ. On the other hand, theνe appearance measurement at T2K, which is the successor experiment of K2K experiment, alsoestablished the nonzero θ13 at 7.3σ in 2013[32]. This experiment uses the neutrinos with theenergy spectrum peaked at 0.6 GeV produced by the collision between accelerated proton beam
and a target. They travel over a distance of 295 km and are detected by Super-Kamiokande
detector. In the global analysis result summarized by PDG in 2016[25], the θ13 and ∆m231 values
with 3σ errors are evaluated as follows:|∆m231 |2−3= (2.56+0.13−0.11 ) × 10 [eV ]2.5(2.15)sin θ13 = 0.0215 ± 0.0025 (∆m231 > 0)

(2.16)= 0.0216 ± 0.0026 (∆m231 < 0)(2.17)

2

Atmospheric Neutrino Oscillation in Matter

When a neutrino passes through matter, the neutrino oscillation eigenstate is changed by the
weak potential that causes coherent forward scattering with electrons. This eﬀect is called
Mikeyev-Smirnov-Wolfenstein (MSW) eﬀect[11][12]. The charged current interaction with electrons occurs only for the electron neutrinos among all the atmospheric neutrinos. Therefore,
the potentials in matter for electron neutrino and the other neutrinos can be written as follows:
Ve = VCC + VN C
Vµ,τ= VN C(2.18)(2.19)

where VN C and VCC are the potentials of neutral and charged current interactions, respectively.
VN C value is not changed by the neutrino flavor. These potentials changes Hamiltonian in

CHAPTER 2. NEUTRINOS

9vacuum as follows:
H0 → H0 + V(2.20)

where H0 is the Hamiltonian in vacuum, and V is the potential which satisfies the following
equation:
  
 
νe
VCC + VN C
0
0
νe





0
VN C
0
νµ 
V νµ =(2.21)
ντ
0
0
VN C
ντ
Based on this Hamiltonian, Schrödinger equation is calculated as follows:
 
  
 
 
Ve 0
0
νe
ν1
νe
νe
d
0  νµ i νµ  = (H0 + V ) νµ  = H0 U ν2  +  0 Vµτ
dt
ντ
ντ
ν3
0
0 Vµτ
ντ






(p + m21 /2E)ν1
VCC + VN C
0
0
νe
0
VN C
0  νµ = U (p + m22 /2E)ν2  + 2ντ(p + m3 /2E)ν3
0
0
VN C

    2
m1 /(2E)
00
p + VNC + VCC
00νe U† + 0
m22 /(2E)00p + VNC0  νµ = U 
00
m23 /(2E)00p + VNCντ  

2−∆m21 /(2E) 00
VCC 0 0νe
 U† +  0000
0 0 νµ = CI3 + U 0
0 ∆m232 /(2E)00 0ντ(2.22)
where C = p+VN C +m22 /(2E) and U follows Eq. (2.2). For the analysis of atmospheric neutrino,
which have > 1 GeV energy, ∆m221 ∼ 10−5 eV2 is neglected. Then, Eq. (2.22) can be written as
follows:
  



  
νe
0 0
0
VCC 0 0νed   
†






νµ = CI3 + U 0 0
0
0
0 0
νµ 
i
U +
dt
2
ντ
0 0 ∆m32 /2E
0
0 0
ντ



 
2
2
∆m32
∆m32
2
−iδCP
sin
θ
+
V
0
sin
θ
cos
θ
e
νe
13
13
13
CC
2E
2E


 †  000= CI3 + R23  R23  νµ
∆m232
∆m232
+iδ
2
ντ
CP
0
2E sin θ13 cos θ13 e
2E cos θ13



 
sin2 θM,13
0 sin θM,13 cos θM,13 e−iδCPνe∆m2M,32†  



0
0
0
νµ
= CI3 +
R23
R23
2E
ντ
sin θM,13 cos θM,13 e+iδCP 0
cos2 θM,13(2.23)
where R23 , ∆mM,32 and θM,13 are defined as follows:


1
0
0
R23
≡ 0 cos θ23 sin θ23 
0 − sin θ23 cos θ23√
2
2∆m32,M ≡ (∆m32 sin 2θ13 )2 + (2EVCC − ∆m232 cos 2θ13 )2
cos 2θM,13∆m232 cos 2θ13 −2EVCC
2(∆m32 sin 2θ13 )2 +(∆m232 cos 2θ13 −2EVCC )2≡√(2.24)(2.25)(2.26)

CHAPTER 2. NEUTRINOS
VCC is defined as follows:

10√
VCC = ± 2GF Ne(2.27)

where GF is the Fermi coupling constant and ne is the number density of electrons. The VCC
becomes positive and negative values for neutrino and anti-neutrino, respectively. When Ne = 0,
which corresponds to the vacuum environment, ∆m2M,32 and θM,13 become ∆m232 and θ13 .
By using these parameters, the oscillation formula is translated to the same style as that in
vacuum. In vacuum, the neutrino oscillation probability between νµ → νe can be expressed asfollows:1 :∆m232
L(2.28)
P (νµ → νe ) = sin2 2θ23 sin2 2θ13 sin2
4E
In this equation, the elapsed time t is translated to its path length (t → L) assuming the velocity
of light. The oscillation probability in the matter can be calculated by changing ∆m32 →∆mM,32 , θ13 → θM,13 as follows:
P (νµ → νe ) = sin2 2θ23 sin2 2θM,13 sin2

∆m2M,32
4E

L(2.29)As shown in Eq. (2.24), sin 2θM,13 depends on the neutrino energy. It has a peak aroundE ∼ cos 2θ · ∆m232 /(2VCC ), and the oscillation probability is proportional to sin2 2θM as shown in
Eq. (2.29). Therefore, the neutrino oscillation probability will increase when the neutrino energy
satisfies E ∼ cos 2θ · ∆m232 /(2VCC ). Because neutrino energy is a positive value, the oscillation
probability has a peak depending on the mass hierarchy. If ∆m232 is positive (negative), only theoscillation probability of neutrino (anti-neutrino) shows this peak. This MSW eﬀect is manifest
in the solar and atmospheric neutrinos.

2.6

Mass Hierarchy

Based on the measured mixing parameters described in the previous sections, the ordering of
the neutrino mass eigenstates are narrowed to two candidates, which are called normal hierarchy(NH) and inverted hierarchy (IH). As shown in Fig. 2.3, the diﬀerence between NH and IH exists
in the sign of ∆m232 . This hierarchy is important to check if neutrino is a Majorana particle ornot, i.e. the anti-neutrino and neutrino are the same particles or not. The existence of Majorana
particle is possible to describe the baryon asymmetry[4]. Currently, the Majorana neutrino issearched by neutrino-less double beta decay experiments, for example GERDA[33], KamLANDZen[34], and CUORE[36]. KamLAND-Zen measures the decay of 136 Xe. In standard model,
the double beta decay produces 2 neutrinos (136 Xe →136 Ba + 2e− + 2ν̄e ). However, if neutrino
is a Majorana particle, the produced ν̄e is possible to interact with the other hadron as νe .
As a result, electrons are produced without neutrino emission. The neutrino-less double beta
decay has not been discovered yet, and the lower limit of the decay half time has been extended
to longer. Based on this life time limit, the upper limit of the Majorana neutrino mass canbe calculated. The 90% C.L upper limit was evaluated to 61-165 meV from the measurement
result of KamLAND-Zen, [34]. As shown in Fig. 2.4, the eﬀective Majorana neutrino mass as afunction of the lightest neutrino mass depends on the mass hierarchy. If the mass hierarchy isproved to be inverted hierarchy, the nonexistence of Majorana neutrino can be proved when the
1

2

2

[eV ]
For convenience, this equation is often expressed as P (νµ → νe ) = sin2 2θ23 sin2 2θ13 sin2 1.27∆m
L[m]
E [MeV]CHAPTER 2. NEUTRINOS

m23
Flavor''Composi9ons'
νe'
νμ'
ντ

11

m22
m21

⇠ 2.4 ⇥ 10 3 [eV2 ]

m22
m21

⇠ 7.6 ⇥ 10 5 [eV2 ]⇠ 2.4 ⇥ 10 3 [eV2 ]⇠ 7.6 ⇥ 10 5 [eV2 ]m2
3Normal'Hierarchy'(NH) Inverted'Hierarchy'(IH)
Figure 2.3: The comparison of neutrino mass ordering between normal and inverted hierarchy.
The flavor compositions of mass eigenstates are indicated by color bars. ∆m232 has positive andnegative values for the normal and inverted hierarchies, respectively.
upper limit of the mass is analyzed to be less than about 0.01 eV. Therefore, the mass hierarchy
determination is important.

2.7

MSW Eﬀect in Atmospheric Neutrinos

Fig. 2.5 is the schema of the path of atmospheric neutrinos. As shown in this figure, the zenith
angle of the neutrino incident direction to Super-Kamiokande indicates the path lengths in coreand mantle of the neutrino. Therefore, the oscillation probability can be written as a functionof the neutrino energy and the zenith angle. The MSW eﬀect in atmospheric neutrinos areshown in the oscillation probability of neutrinos which pass through the earth. Fig. 2.6 shows
the oscillation probabilities as a function of the zenith angle of neutrino incident direction to
Super-Kamiokande and neutrino energy. The earth core and mantle have the diﬀerent densities.
Therefore, the oscillation probability distribution seems to have two layers. As shown in (b) ofFig. 2.6, the νe appearance probability shows significant enhancement in 2-10 GeV. If the mass
hierarchy is normal (inverted) hierarchy, this enhancement is only shown in (anti-)neutrino. The
mass hierarchy is evaluated by searching for this enhancement.

CHAPTER 2. NEUTRINOS

12

Figure 2.4: The relationship between the lightest neutrino mass and the eﬀective Majorana
neutrino mass (mββ ) (left), and the upper limit of the eﬀective Majorana neutrino mass basedon the measurement result of each nucleus decay as a function of the mass number (right) [34].

CHAPTER 2. NEUTRINOS

13

Figure 2.5: Schema of the atmospheric neutrinos passing through the earth. The path lengths inmantle and core of atmospheric neutrino can be estimated from zenith angle of the atmospheric
neutrino detected in Super-Kamiokande. The core has larger electron density than mantle, andit changes the oscillation probability owing to the MSW eﬀect.

CHAPTER 2. NEUTRINOS

14

Figure 2.6: Oscillation probabilities for neutrinos and anti neutrinos when the mass hierarchy isnormal hierarchy. The left two figures show νµ and ν̄µ survival probabilities, and the right two
figures show νe and ν̄e appearance probabilities. These figures assume (∆m232 , sin2 θ23 , sin2 θ13 ) =(2.5 × 10−3 , 0.5, 0.0219). These figures are taken from [35].

Chapter 3

Super-Kamiokande

Figure 3.1: Super-Kamiokande detector

3.1

Overview

Super-Kamiokande(SK) is Kamiokande’s successor experiment with a larger water Cherenkov
detector. Originally, Kamiokande experiment searched for the nucleon decay, but it could also
detect neutrino interaction. This experiment did not find nucleon decay event and rejected
SU(5) model, which is the simplest nucleon decay model. On the other hand, this experimentobserved neutrinos from supernova SN1987A [37] and the solar neutrino disappearance which
was known as ”Solar neutrino problem” [17]. In addition, this experiment observed unexpected
atmospheric neutrino disappearance [38]. This disappearance was considered to come fromneutrino oscillation. To prove this, more statistics was required.
To search for new physics, SK experiment started from April.1996. This experiment uses
50 kton detector, which is 17 times larger detector than that of Kamiokande. In 1998, theevidence of neutrino oscillation was detected as described in Sec. 2.3. This detector is located15

CHAPTER 3. SUPER-KAMIOKANDE

16

in the 1000 m underground (2700 m water equivalent) of Mt. Ikenoyama, in Gifu Prefecture,
Japan, which is the same position as Kamiokande. Owing to the deep underground location,
most part of cosmic ray muons, which is the background of neutrino interaction, is stopped in
the rock. This detector is filled with pure water, and Cherenkov light from charged particle
is measured. This light is detected by using photomultiplier tubes (PMT). Their signals aresent to data acquisition (DAQ) system and their information of charge and hit time are stored.
Because of the directivity of Cherenkov light, this detector can reconstruct the direction of
produced particles and identify the flavor by using the distribution of PMTs with signal.

3.2

SK periods

The measurement term is separated into four periods; SK-I, SK-II, SK-III, SK-IV. SK-I continued from April 1996 to July 2001. In November 2001, PMT chain reaction implosion happened
and about half of PMTs were destroyed. By using the remaining PMTs and reducing the photocathode coverage from 40% to 19%, the detector was rebuilt. Measurement was resumed in
October 2002 as SK-II period. From SK-II, all PMTs were put into fiber reinforced plastic(FRP) case to reduce the risk of chain reaction implosion. This period continued until October
2005. After that, new PMTs were introduced and the coverage returned to 40%. The next period, SK-III, continued from July 2006 to September 2008. To satisfy requirements to the DAQ
system such as the detection of neutron capture signal, new front-end electronics was developedand introduced after the SK-III period. After the introduction, SK-IV started, and it continues
until May 2018. Features of four periods are summarized in Table 3.1.

PeriodBegin
EndLive time (day)Number of ID PMTs
Photocathode coverage
Front-end electronics

SK-I
Apr. 1996
Jul. 2001
1489.2
11146
40%
ATM

SK-II
Oct. 2002
Oct. 2005
798.6
5182
19%
ATM

SK-III
Jul. 2006
Sep. 2008
518.1
11129
40%
ATM

SK-IV
Sep. 2008
May. 20183244.4
11129
40%
QBEE

Table 3.1: SK phases and the features

3.3

Cherenkov Radiation

When charged particle passes though dielectric medium, the molecules in the medium are polarized. When the charged particle is faster than the light speed in the medium, which satisfy
Eq. (3.1), it polarizes the medium molecules asymmetrically.
v≥

c
n

(3.1)

When these medium atoms are depolarized, it radiates electromagnetic wave, which is called
Cherenkov light[39]. This light is emitted as a cone, and the opening angle of Cherenkov emission

CHAPTER 3. SUPER-KAMIOKANDE

17

follows Eq. (3.2).
cos θC =

c/n
1
=
cβ
nβ

(3.2)

v is the absolute velocity of charged particle, n is the refractive index of the medium and c isthe light speed in vacuum. The refraction rate of pure water is about 1.34. So the opening angle
is 42◦ when the particle satisfy β ∼ 1. Because the opening angle does not directly depend on
the particle’s momentum, but its speed, the opening angle of the Cherenkov radiation from therelativistic particles is constant. Therefore, the distribution of Cherenkov light on the detector
wall shows ring shape as shown in Fig. 3.2.

Figure 3.2: Schema of Cherenkov radiation. The angle of cone depends on the particle velocity
and refractive index. n is refractive index of medium, c is light speed in vacuum. In the timeinterval t, Cherenkov light moves t(c/n) and the charged particle moves βct. Because the opening
angle of Cherenkov radiation from the particle with relativistic momentum is almost constant,
the PMT signal distribution shows the ring shape.
The radiation intensity follows the following equation:()
d2 N
2παq 212παq 2=1−=
sin2 θCdxdλ
λ2
β 2 n(λ)2
λ2

(3.3)

where α is fine structure constant, λ is the wavelength of the photon, x is the length along thecharged particle path, and q is the charge of the particle. Because the sign of charge does not
aﬀect to the flux, the Cherenkov light from lepton and anti-lepton can not be distinguished.

CHAPTER 3. SUPER-KAMIOKANDE18In Super-Kamiokande, the neutrino and anti-neutrino identification uses the number of Michelelectrons and so on.

3.4

Components of Super-Kamiokande Detector

After this section, each component of Super-Kamiokande is described. Fig. 3.3 shows the schema

Figure 3.3: Schema of Super-Kamiokande and its components.
of Super-Kamiokande component positions, and the chapter number in which they are described.

3.5

Water Tank

The Super-Kamiokande detector has a cylindrical shape with 39.3 m diameter, 41.4 m heightand 50 kton volume. This detector is separated into two regions, inner detector (ID) and outerdetector (OD). ID has a cylindrical shape with 33.8 m diameter, 36.2 m height and 32 ktonvolume. The OD wraps ID with approximately 2 m thickness (2.05 m for top and bottom, and2.20 m for barrel). Both detectors are filled with pure-water. By using Cherenkov light, physics
measurements are done in ID. OD works as veto for incoming cosmic ray. To detect the light,
20 inch and 8 inch diameter photomultiplier tubes are introduced in the ID and OD, respectively.

CHAPTER 3. SUPER-KAMIOKANDE

3.6

19

Photomultiplier Tubes for Inner Detector

Photomultiplier tubes in the inner detector (R3600 model, developed by Hamamatsu Photonics
K.K) have 20 inch diameter. Fig. 3.4 shows their detailed structure. The inner surface of the

Figure 3.4: Design of the photomultiplier tube in the inner detector [40]front glass is coated with bi-alkali (Sb-K-Cs) photocathode. Fig. 3.5 shows its quantum eﬃciency
as a function of light wavelength. This photocathode has sensitivity to 300-600 nm wavelength
light, and its quantum eﬃciency is 22% at maximum. This PMT has 11 stages Venetian blind
type dynodes[41]. By applying high voltage to this PMT, and providing electric potential
diﬀerence between each dynode, generated photoelectrons are multiplied by the collision withthese dynodes. To set their gain to ∼ 107 , high voltages of about 2 kV are applied to PMTs.
These performances are summarized in Table 3.2
Shape
PMT Diameter
Gain
Dynodes
Voltage
Quantum Eﬃciency

Hemispherical
20 inch∼ 107
11 stages venetian blind type
∼ 2000 V
22% at peak

Table 3.2: Feature of the inner detector PMTs

3.7

Outer Detector Systems

In the outer detector, 1885 8 inch PMTs (R1408) are installed. To improve the sensitivity toparticles passing through the OD, the wavelength shifting plate with 1.3 cm × 60 cm × 60 cm
is equipped around each PMT. Fig. 3.6 shows the shape of the plate. At the start of SK-I,
IMB-3 PMTs and wavelength shifter were reused from the IMB experiment. However, because

CHAPTER 3. SUPER-KAMIOKANDE

20

Figure 3.5: The quantum eﬃciency of the inner detector PMT photocathode as a function oflight wavelength[40]Figure 3.6: The shape of wavelength shifter plate on the surface of the outer detector PMT[43].

CHAPTER 3. SUPER-KAMIOKANDE

21of the accident described in Sec. 3.2, half of the PMTs were replaced with newly produced
ones for SK-II. Fig. 3.7 shows the OD PMT’s quantum eﬃciency and the wavelength shifter’sFigure 3.7: The left figure shows the outer detector PMT quantum eﬃciency as a functionof light wavelength, and the right figure shows the absorption (shorter wavelength curve) andemission (longer wavelength curve) spectra of the wave shifting plate. The absorption spectrum
uses the molar extinction coeﬃcient µ as Y-axis on the left. The absorption probability follows
Pabs = 1 − 10−µρd , where ρ is the flour’s concentration in mol/l, and d is the thickness of thematerial. The flour emission spectrum uses the Y-axis on the right [43].
absorption and emission spectrum. The wavelength is shifted to the region where the PMT has
higher quantum eﬃciency. In addition, the OD wall is lined with reflective sheets (Type 1073B
Tyvek) to suppress the light absorption. The PMTs in the tank are connected to the front-end
electronics with cables, which pass through twelve cable holes on the top of the tank. OD PMTs
are not installed around the four holes out of the twelve holes. If a cosmic-ray muon entersfrom this no OD PMT region, the OD veto system possibly fails. To compensate this insensible
OD area, plastic scintillator veto counters are installed on the four holes. The OD hit clusters
indicates the points of entrance to and exit from the inner detector. However, when the cosmic
ray muon passes through the corner of the inner detector, it is diﬃcult to identify the photons
produced at the muon’s entrance and exit points. This muon is called corner clipping muon. To
separate these photons, OD segmentation is added as shown in Fig. 3.8 at the period between
SK-II and SK-III run. It improves the eﬃciency of the event selection, which is described inSec. 7.2.

3.8

Purification Systems

The water quality is important for Cherenkov light detection. If the water is contaminated,
Cherenkov light will be frequently absorbed or scattered before their incidence on the PMTs,
and the performance of the event reconstruction will be worsened. Especially, the removal ofradioactive contaminants is important. It is the source of not only the background for low energy
events, but also radiation exposure of the mineworkers. The surrounding rock contains uranium,
which generates radon gas in the detector. To remove these contamination sources and maintainthe water quality, air and water purification systems were constructed. In this section, theseCHAPTER 3. SUPER-KAMIOKANDE

22

Figure 3.8: Schematic views of the detector with the OD segmentation and its corner. As shownin the right figure, the Cherenkov light which is produced in the barrel and endcap region of
OD are separated to reject the corner clipping muon.
systems are described briefly. The detail is described in [42].

3.8.1

Air Purification System

In Super-Kamiokande experimental area, Rn concentration becomes 1500 Bq/m3 in summer and30 Bq/m3 in winter. To keep the Rn concentration at low level, air purification is done with
two steps. As a first step, air circulation system supplies flesh air from outside with a rate of
70 m3 /min. This system keeps the Rn contamination in whole experimental area at ∼50 Bq/m3 .
As a next step, the air in the detector tank is purified by using Rn free air system. This system
uses activated carbon to absorb Rn gas. It reduces the Rn concentration level at ∼30 mBq/m3 .

3.8.2

Water Purification System

Water purification system consists of a number of components. Fig. 3.9 shows schematic viewof the water purification system. As a first step, mesh filters remove the dust larger than 1 µm
in size. After that, several components remove diﬀerent kinds of contaminants. High molecularand heavy ion particles are rejected by ion exchanger and inverse osmosis. Radon and oxygengases are removed by vacuum degasifier. The growth of bacteria is suppressed by sterilizing
with UV light. After these purifications, this water passes through ultra filter, which removes
particles larger than 10 nm. From SK-II, membrane degasifier had been used to remove Rn
gas from the water after ultra filter. However, this component was bypassed on 2014 in SK-IV
because it was found that the Rn concentration increased when the water passed through thisdegasifier. This system keeps the water temperature with 0.01 ◦ C accuracy to reduce convection
inside the tank.

CHAPTER 3. SUPER-KAMIOKANDE

23

Figure 3.9: Water purification system[42]This system circulates the water in the tank with 60 ton/hour flow rate and reduces the Rnconcentration to 1.83±0.31 mBq/m3 .

3.9

Data Acquisition Systems

The PMT signals are sent to the front-end electronics passing through cables with the uniform
length of 70 m. For SK-I, II, and III period, Analog Timing Modules (ATM) were used for thefront-end electronics. In 2008, the data acquisition system was upgraded for SK-IV. The new
system uses QTC-Based Electronics with Ethernet (QBEE). In this section, the feature of these
old and new data acquisition systems are described.

3.9.1

Systems for SK-I, II, and III

Fig. 3.10 shows the structure of ATM module. ATM has pairs of QACs (Charge to Analog
Converter) and TACs (Time to Analog Converter). QAC integrates the charge of input signal
in a 400 ns gate. On the other hand, TAC measures the time interval between two triggers by
integrating a constant current. When the first trigger inputs in TAC, the charge integration isstarted, and it is stopped when the second trigger inputs. The integrated charge indicates thetime interval. The outputs of QAC and TAC are sent to ADC (Analog to Digital Converter),but it takes 5.5 µs to digitize an event. The pairs of QACs and TACs enable to store the input
during the digitization process. However, when a strong signal is input to ATM, a part of thesignal is reflected to PMT. The reflected signal is input to ATM again 700 ns after the first
input. Therefore, the signals in 400 ns - 900 ns from the first input are neglected to remove thisnoise.
When the PMT signal is input to ATM, it splits into several signals. The discriminatorreceives one of these signals. When the input pulse height exceeds the threshold which corresponds to a quarter of single photoelectron, QAC and TAC start to integrate the input charge
and the constant current, respectively. At the same time, a rectangular signal with 200 ns widthCHAPTER 3. SUPER-KAMIOKANDE

24

Figure 3.10: ATM structure[46]is generated. The rectangular signals of all the channels in the ATM are summed and sent toglobal trigger system. This signal is called HITSUM signal. The global trigger system sums up
the HITSUM signals of all the PMTs, and a global trigger is generated when the pulse heightexceeds the threshold. Because Super-Kamiokande measures many kinds of particles, including
solar and atmospheric neutrinos, three kinds of triggers with diﬀerent HITSUM threshold, which
are called the high energy trigger (HE) , the low energy trigger (LE), and the super low energy
trigger (SLE) are prepared. The SLE trigger is introduced in 1997 to analyze the solar neutrino
with lower energy. OD trigger is also generated based on the similar procedure. When thenumber of OD PMT hits in 200 ns time window exceeds 19 hits, the trigger is generated.
Their thresholds are shown in Table 3.3. When the charge integrations of TAC and QAC

SLE
LE
HE

SK-I
186 mV
320 mV
340 mVSK-II
110 mV152 mV
180 mV

SK-III
186 mV
302 mV
320 mV

Table 3.3: Global trigger thresholds. Because the width of a rectangular pulse is 200 ns, thisthreshold corresponds to the number of hits in 200 ns. The pulse height is 11 mV per hit.

continue for 1300 ns, their charges are cleared. It means the timing dynamic range is uniformly
1300 ns. The threshold of the SLE trigger corresponds to the signal from ∼ 5 MeV electron. As
shown in Fig. 3.11, the digitized data are sent to the online host computer passing through the
VME memory module which is called SMP (Super Memory Partner) and the slave computer.

CHAPTER 3. SUPER-KAMIOKANDE

Figure 3.11: DAQ system of SK-I,II,III

25

CHAPTER 3. SUPER-KAMIOKANDE

3.9.2

26

Electronics for SK-IV

QBEE measures the charge and hit time of input signal with QTC (Charge to Time Converter)and multi-hit TDC (Time to Digital Converter). When a signal is input to QTC, it outputs
the rectangular pulse whose width is proportional to the charge. The pulse width is measured
by getting its rise and fall time with TDC. QTC has three channels with diﬀerent gains. Thesechannels enable to have a wide dynamic range from 0.2 pC to 2500 pC. The QBEE dynamic
range covers higher energy than the upper limit of the ATM dynamic range, which is 450 pC [40].
It improves the performance of the high energy event reconstruction. QBEE process speed is so
fast that all PMT signals can be stored. The data are registered in the disk passing through the

Figure 3.12: DAQ system of SK-IV
system as shown in Fig. 3.12. The QBEE data are read out by front-end PCs, and their data aremerged by merger PCs in parallel. At this process, the merger PC applies software triggers to
reduce the data size. Unlike the hardware trigger of ATM, the gate width of the software triggers
can be set one by one. Table 3.4 shows the threshold and gate width of software triggers in SKIV. This new electronics has many advantages for atmospheric neutrino measurement. First,
owing to the software trigger, the detection of neutron capture, which has very long capturetime(205 µs), was enabled. Second, owing to the nonexistence of dead time between 400-900 ns
from the previous signal, the decay electron tagging eﬃciency was improved. Third, owing to
the wider dynamic range of charge, the energy resolution for high energy neutrino events was
improved.

CHAPTER 3. SUPER-KAMIOKANDE

SK-IV triggers
OD
SLE
LE
HE
SHE
AFT
(After)Threshold (Number of hits)
22
34→31 (Changed in 2015)47
50
70→58 (Changed in 2011)
SHE and no OD
is generated

27

Time width [µs]
(VETO)-0.5 → 1.0-5 → 35-5 → 35-5 → 35
35 → 535

Table 3.4: Triggers threshold and event widthChapter 4

Calibration
4.1

PMT Calibration

The Cherenkov light information is reconstructed from the charge and detection time of PMTs.
To reconstruct this precisely, the PMTs performances have to be known, and their individual
diﬀerences have to be calibrated. PMT calibration consists of three factors, gain, quantum
eﬃciency and detection timing. These factors are described in the following sections.

4.1.1

Gain

PMT gain can be adjusted by changing its high voltage. The high voltages are set for all PMTs
to produce the same charge from the same intensity light. This setting uses Xe lamp in a
scintillator ball. This light source is set in the center of the inner detector. The each PMT’s
incident light intensity from the scintillator ball depends on the distance and direction from thescintillator ball to the PMT. Therefore, the inner detector has 420 ”standard” PMTs, which are
the PMTs calibrated before the installation. As shown in Fig. 4.1, the PMTs in the detectors
are separated to some groups based on the their geometrical relationship to the light source.
The high voltage of each PMT is set to produce the same charge as the standard PMTs in thesame group. Because their gains are set to correct the individual diﬀerence of their quantum
eﬃciency, adjusted gains also have the individual diﬀerences. To calibrate the gain one by one,their averaged absolute gain of all PMTs and relative diﬀerence are checked, respectively.
4.1.1.1

Relative Diﬀerence of Gain

This calibration procedure uses the stable light source (nitrogen-laser-driven dye laser light)that is deployed in the center of the detector. Relative diﬀerences of gain between PMTs areevaluated by comparing two kinds of the PMT outputs that are produced from diﬀerent intensity
lights. The average charge Qobs (i) is measured with higher intensity light for each PMT i. On
the other hand, the number of each PMT’s single photoelectron detection (Nobs (i)) is measured
with lower intensity light. The light intensity in the Nobs (i) measurement is weak enough for
all PMTs to detect only single photoelectron at maximum, then Qobs (i) and Nobs (i) depends on
the similar factors:
Qobs (i) ∝ Is × a(i) × ϵqe (i) × G(i)(4.1)

Nobs (i) ∝ Iw × a(i) × ϵqe (i)

(4.2)

28

CHAPTER 4. CALIBRATION

29

Figure 4.1: Location of the standard PMTs in the inner detector (Left figure). The red point
indicates the standard PMT. The standard and normal PMTs are separated to some groups
based on the geometrical relationship to the light source (Right figure). These figures are takenfrom [44].

CHAPTER 4. CALIBRATION

30

where Is and Iw are the average intensities of high and low intensity flashes, a(i), ϵqe (i), and G(i)are acceptance, quantum eﬃciency, and one-photoelectron gain of each ID-PMT, respectively.
The relative gain is derived by taking the ratio of (4.1) to (4.2):
G(i) ∝
4.1.1.2

Qobs (i)Nobs (i)(4.3)

Absolute Gain Conversion Factor

Absolute gain calibration is required to convert the detected charge to the number of photoelectrons. This calibration uses the cumulative charge distribution of single photoelectron after the
relative gain correction based on Sec. 4.1.1.1. To measure the charge distribution, Ni-Cf source
is deployed at the center of the detector. This source consists of the nickel source with hollow
structure and 252 Cf radiation source in the center of the nickel. When 252 Cf decays, it produces
neutrons. These neutrons are captured by the surrounding Ni, and produce γ-ray with 9 MeV
at maximum. γ-ray can be detected because it interacts with the water and produce electrons,
which emits Cherenkov light.

Figure 4.2: Pulse shape of the signal from each PMT. the signal is detected at the timing whenthe voltage exceeds the threshold.
In usual, the electronics does not record the PMT signals when its pulse height is less than
25% of the single photoelectron’s height as shown in Fig. 4.2. To convert the charge to thenumber of photoelectrons, the signal with small pulse height also has to be considered. To takethe signal with low pulse height into account, this measurement is also done with the halved
threshold and doubled PMT gain. The region with less than 0.3 pC, in which the pulse heightis less than the halved threshold, is approximated by linear function. To reject the backgroundeﬀect, the measurement without Ni-Cf source has been also done. Combining these measurement
and approximation results, the single photoelectron charge distribution is evaluated as shownin Fig. 4.3. Its mean value is defined as the absolute gain. The gain has been continuously
increasing, but this reason has not been identified yet. To correct this change, the measured
gains in each period are used as a conversion factor from charge to the number of photoelectrons.
They are determined as shown in Table 4.1, respectively.

CHAPTER 4. CALIBRATION

31

Figure 4.3: The single photoelectron (p.e.) charge distributions in SK-III with diﬀerent charge
range (left: 5 pC, right: 50 pC). The green line in the left figure is the distribution estimatedbased on the result of the measurement with half threshold and double gain. The red line is theregion approximated by linear function. This figure is taken from [44].

SK-I
SK-II
SK-III
SK-IV

Absolute Gain (pC/1pe)
2.055
2.297
2.243
2.645

Table 4.1: Absolute gain

4.1.2

Relative Diﬀerence of QE

As mentioned in Section 4.1.1, PMTs have individual diﬀerences in their quantum eﬃciencies(QE). The relative diﬀerence of quantum eﬃciency is analyzed by comparing the number ofdetected photoelectrons and the expected incident photons of each PMT. To count the numberof photoelectrons precisely, the Ni-Cf source is used again. The number of incident photons isestimated by MC, which takes the eﬀect of the light attenuation and reflection into account.
Owing of the light attenuation in the water and the reflection on the surface of the wall, thephoton detection probability has position dependence even if the diﬀerence of the acceptance iscorrected. Fig. 4.4 shows the position dependence of the hit probability whose counts of each
PMT hits are corrected based on:
Nobs,corr (i) = Nobs (i) × R(i)2 /a(θ(i))(4.4)

where i is the index of ID-PMT, R(i) is the distance from the light source to i-th PMT, anda(θ(i)) is the acceptance as a function of the angle between the direction which i-th PMT
head for and the direction from the light source to i-th PMT. As shown in Fig. 4.4, it shows

CHAPTER 4. CALIBRATION

32

Figure 4.4: Normalized hit probability as a function of PMT position. the upper figure shows
the dependence on z (cm). Below left and right figures show the position dependence of top andbottom PMTs, where the x-axis is the square of the distance from the center to the PMT (cm2 ).
The red plot is the data and the blue plot is MC. This figure is taken from [44].

CHAPTER 4. CALIBRATION

33

the consistent result with the data. By comparing the number of each PMT hits based on
MC and data individually, the relative QE is evaluated. As a result, it was proved that QE
largely depends on the PMT production time. As shown in Fig. 4.5, the PMTs which are newly
installed in SK-III have higher quantum eﬃciencies than PMTs used in SK-II. In addition,
the diﬀerence is also shown between the PMTs installed in SK-III, which are produced in thediﬀerent production run. The PMTs in the bottom region were produced in 1996-1997. On
the other hand, the top and barrel PMTs are produced earlier, in 1992-1995. Because of thisdiﬀerence, the QE of bottom PMT is higher than the other region by 4.5%. This diﬀerence comesfrom the improvement of glass transparency by Hamamatsu Photonics K.K., the manufacturerof the ID PMTs. Therefore, this measurement shows consistent result with the expectation.

4.1.3

Timing Calibration

To translate the detection time to the photon incident time, the individual diﬀerence of thetransition time from PMTs to their electronics and the dependence on their charge have to
be considered. As mentioned in Section. 4.1.1.2, the PMT signal is detected when the output
voltage exceeds the threshold. The detection time is also registered at this timing. When thelight intensity is strong, many photoelectrons are produced and the rise time of its pulse becomesshorter. Therefore, the detection time depends on the signal strength. To correct this eﬀect, thecharge and detection time relationship is analyzed. As shown in Fig. 4.6, nitrogen-laser-driven
dye laser with 398 nm wavelength is used for this measurement, and its light is guided to thediﬀuser ball in the center of the detector passing through the optical fiber. A monitor PMT isused to define the light injection time. Based on this measurement, the 2-dimensional correlationtable of timing and charge is produced as shown in Fig. 4.7. In this analysis, the time-of-flightfrom the diﬀuser ball to each PMT is subtracted from its detection time. In this distribution,
the axis of the charge has 180 bins, which is called QBins. The QBin is defined as:
QBin = 0.2 × Q (pC) (Q < 10 pC)(4.5)

QBin = 50 × log10 (Q (pC)) (10 pC < Q < 3981 pC)(4.6)

After the event separation with QBins, the transit time of photoelectrons is analyzed by getting
the peak of the timing distribution of each QBin. The variation of the timing peak is fitted bypolynomial function, and their detection times are corrected by subtracting the timing oﬀset
based on this function.

4.2

Calibration for Photon Propagation

To select the suitable model of photon propagation in the tank for MC, the properties of SK
tank components have to be calibrated. In this section, the calibration methods of the propertiesare described.

4.2.1

Water Transparency

Basically, photon with wavelength λ attenuates following to:I(l, λ) = I0 (λ)e−l/L(λ)(4.7)

CHAPTER 4. CALIBRATION

34

Figure 4.5: Normalized hit probability comparison between PMTs used in SK-II and newly
installed from SK-III, similar to Fig. 4.4. The green line shows the newer PMTs, and the redline shows the older one. This figure is taken from [44].

CHAPTER 4. CALIBRATION

35

Figure 4.6: Schematic view of the timing calibration system. the light intensity can be changed
by a variable optical filter. This figure is taken from [44].

CHAPTER 4. CALIBRATION

36

Figure 4.7: The (time-of-flight corrected hit timing)-(charge) distribution. The sign of Y-axis
value is reversed from the hit timing. As Y-axis value becomes larger, the hit timing becomes
earlier. This figure is taken from [44].
where l, I0 (λ), and I(l, λ) are the travel length, the initial light intensity, and the intensity atl, respectively. L(λ) is the total attenuation length caused by scattering and absorption. Thisvariable is defined as:
1
L(λ) =(4.8)
αabs (λ) + αsym (λ) + αasym (λ)
where αabs (λ), αsym (λ), and αasym (λ) are the amplitude of absorption, symmetric scattering,
and asymmetric scattering, respectively. Scattering eﬀect is separated to the components withforward-backward symmetry and that with asymmetry. The angular distribution of the symmetric component has 1 + cos2 θ factor, where θ is the scattering angle of a photon. The
symmetric scattering comes from Rayleigh and symmetric component of Mie scattering, and theasymmetric scattering comes from the forward component of Mie scattering. To calibrate theseamplitudes, collimated laser beam is injected as shown in Fig. 4.8, changing the wavelength.
Fig. 4.9 shows the distribution of the hit time after the time of flight (TOF) subtraction. The
TOF is calculated based on the length between the target of the laser at the bottom and each
PMT. Owing to this TOF subtraction, the signal from the light reflected on the bottom of thedetector produces a peak around 1850 ns, regardless of the PMT position. The earlier signals
are considered to come from the light scattered in the water. The absorption and scatteringparameters are tuned for MC by using the hit timing data which is taken earlier than the line
around 1830 µs in Fig. 4.9. To do this calibration, the PMTs on the top and barrel wall of theinner detector are separated to six groups based on their position. To evaluate amplitude of
the scattering and absorption for the light with each wavelength, the timing distributions for
the PMTs in each group are produced with MC based on various patterns of αabs , αsym , andαasym . The amplitudes which minimize χ2 value for the diﬀerence between the MC output anddata are selected as the most suitable amplitudes. Fig. 4.10 shows the relationship between the
wavelength and the evaluated amplitudes of absorption and scattering. Each amplitude is fitted

CHAPTER 4. CALIBRATION

37

Figure 4.8: Schema of water transparency calibration. The barrel PMTs are separated to 5groups, which is labeled B1 to B5. By comparing the timing distributions of PMTs in the top
and these barrel regions, detector parameters are tuned. This figure is taken from [44].

CHAPTER 4. CALIBRATION

38

Figure 4.9: TOF subtracted timing distributions for the laser beam with 405 nm wave length.
The red and black lines are assigned to the data and MC with best tune. The top figure is for
the PMTs at the top of the tank, and the second to bottom figures is for the PMTs in B1 to
B5 regions, respectively. The definition of B1 to B5 regions are written in Fig. 4.8. The time
region between the earlier two blue lines are regarded as hits from the scattered light, and thetime region between later two lines are regarded as the reflected light. This figure is taken from[44].

CHAPTER 4. CALIBRATION

39

Figure 4.10: Inverse water transparency based on April 2009 data. Parameters are fitted with
(4.9), (4.10), and (4.11). This figure is taken from [44].
with the empirically determined functions as follows:
αabs= P0 ×

αsym

=

αasym

P1
+C
4
λ
(

(4.9)

)

P4
P5
× 1.0 + 2
λ4
λ
(
)
P7
2
= P6 × 1.0 + 4 × (λ − P8 )λ(4.10)(4.11)The value C uses the amplitude based on the experimental result for λ ≥ 464 nm [47] [48]. For
λ ≤ 464 nm, it follows the formula:
C = P0 × P2 × (λ/500)P3(4.12)As a result of the fitting, the parameters are determined to the value shown in Table 4.2
P0
0.624

P1
2.96 × 107

P2
3.24 × 10−2

P3
10.9

P4
8.51 × 107

P5
1.14 × 105

Table 4.2: Fitted parameters

P6
1.00 × 10−4P7
4.62 × ×106

P8392

CHAPTER 4. CALIBRATION

4.2.2

40

Light Reflection at PMT

The light reflection at the PMT surface is evaluated by using the same measurement result as
Sec. 4.2.1. The peak around 1850 ns in Fig. 4.9 shows the light reflected on the bottom. Basedon the intensity of this peak, the complex refractive index (Real part: nreal , Imaginary part;nimg ) of the photocathode is tuned. The relationship between the light absorption coeﬃcient(α) and nimg can be written as follows:
α = 4πnimg /λ(4.13)

where λ is the wavelength of the incident light. The best fit value were nimg = 1.667 andnreal = 2.31, 2.69, 3.06, 3.24 at λ = 337, 365, 400, 420 µm

4.2.3Light Reflection at Black Sheet

The reflection at the black sheet is evaluated by using the module shown in Fig. 4.11. The

Figure 4.11: The schema of the module for the black sheet reflection measurement[44]reflected charges (Qreflected ) are measured by changing the incident angle (30, 45, 60 degree)
and wavelength (337, 400, 420 nm). As a reference, the direct charges without black sheet
reflection (Qdirect ) were also measured. As shown in the Fig. 4.12, the reflectivity was compared
between data and MC. The agreement between data and MC is better than the ±1% level at
all wavelength and incident angle range.

4.3

Energy Scale

The energy, flavor and direction of the charged particles produced in the neutrino interaction arereconstructed based on the event reconstruction tools described in Sec 6. The particle energy
largely depends on the signal intensity of ID PMTs. To convert the charges of ID PMTs to theCHAPTER 4. CALIBRATION

41

Figure 4.12: Comparison of the black sheet reflection between MC and data. The top, middle,
bottom figure show the measurement result with the 337, 400, and 420 nm wavelength lightsources, respectively. The left figure shows the charge ratio (red circle: MC, black cross: data).
The ratio between data and MC for the left plots value is shown in the right plots.[44]particle momentum, energy scale is calibrated based on the control samples, which are decay
electron momentum, π 0 invariant mass, and Cherenkov angle and track length of cosmic raymuons. These measurements cover diﬀerent momentum ranges. In this section, these calibration
sources and the systematic error of energy scale are described briefly. The detailed description
is given in [45].

4.3.1

Calibration Sources

4.3.1.1

π 0 Invariant Mass

The atmospheric neutrino interaction with nuclei often produces π 0 . It decays to two γ-rays
just after its production, which produce two of the electron like Cherenkov rings. In two body
decay events, the invariant mass of the parent particles can be reconstructed based on produced
particles’ momenta as follows:√
Mπ0 = 2Pγ1 Pγ2 (1 − cos θ)(4.14)
where Pγ1 and Pγ2 are the momenta of produced γ-rays, and cos θ is the opening angle between
two γ-rays. Fig. 4.13 shows the distribution of reconstructed invariant mass based on the data
and MC. These distributions are fitted with a Gaussian function, and the energy calibration isdone by using their peak position. In this figure, the diﬀerence of peaks between data and MC
is less than 2%. This calibration covers ∼100 MeV/c region, which is the same level as π 0 mass.

CHAPTER 4. CALIBRATION

42

Figure 4.13: The reconstructed invariant mass distributions for SK-IV period. The black dots
and red boxes are assigned to the data and MC.

CHAPTER 4. CALIBRATION
4.3.1.2

43

Decay Electron Momentum

Muon decay events are used for energy scale calibration because the energy spectrum of decay
electron is well known. The decay electrons are searched based on the signal distributions in
the triggered event gate with 1.3 µs width. The basic definition of event gate and the method
how to search for the decay electron is described in Sec 8.1.2. For this calibration, the decay
electron candidates are required to satisfy the additional criteria below:• 2.0 µs<(Muon decay time following the stopping muon detection)< 8.0 µs• (The number of ID PMTs with hits in 50 ns time window)>60 hits (30 hits for SK-II)• Glow >0.5
• The reconstructed muon decay vertex is in the fiducial volume, which is the region inside
ID and more than 2 m away from the wall.
Muon decay vertex is reconstructed based on the algorithm for low energy events[49]. Glow is
the goodness of the reconstruction. To search for the decay electrons, and calculate Glow andreconstruct vertex, this algorithm separate the triggered event gate to a signal time window andbackground regions as shown in Fig. 4.14. SK-IV events have wider gate than 1.3 µs, but thisFigure 4.14: The schema of the definition of the time regions for the cluster detection algorithm [49].
analysis uses the same width gate as the ATM trigger gate. The signal region is set to have themaximum significance S, which is defined as follows:S≡

Nsig − Nbg√
Nbg(4.15)

CHAPTER 4. CALIBRATION

44

where Nsig is the number of hits in the selected signal time window, and Nbg is the expected
number of the background hits in the signal time window, which is calculated as follows:
Nbg =

∆Tsig
(NATM − Nsig )∆TATM − ∆Tsig(4.16)

where ∆TATM and ∆Tsig are the width of the ATM trigger gate and the selected signal time
window, and NATM is the number of hits in the ATM trigger gate. By using this time window,
the goodness of reconstruction is calculated as follows:
)
(
Nsig
∑ 1
(tres,i − t0 )2 )1
(4.17)
Glow = ∑N
exp −
sig 1
σi2
2σi2
2 i=1
i=1 σiwhere tres,i is the time-of-flight subtracted time of i-th hit, and σi is i-th PMT time resolution.
Fig. 4.15 shows the comparison of the reconstructed momentum of decay electrons for SK-IV

Figure 4.15: The reconstructed momentum of decay electrons for SK-IV period. The black andred lines are assigned to the data and MC, respectively.
between data and MC. The discrepancy of the mean of their momenta was 1.5%, and that forthe other periods were also less than 2% as shown in Fig. 4.16.
4.3.1.3

High Energy Stopping Muon

High-energy muon’s path length in the detector largely depends on its momentum at entrance
point. By comparing the muon momentum that is estimated based on its track length and theCHAPTER 4. CALIBRATION

45

detected charge in the inner detector, the energy scale can be calibrated. This calibration uses
the events that satisfy the criteria below:• The muon entrance point is on the top of the detector.
• cos Θ > 0.94 (Downward-going muon)• One decay electron is detected.
• Reconstructed muon track length L satisfies 7 < L < 30 m.
where Θ is the zenith angle of muon direction. The track length is calculated by using the
distance between the muon’s entrance point and the reconstructed decay electron vertex.
4.3.1.4

Low Energy Stopping Muon

As described in Sec. 3.3, the cone angle of Cherenkov light depends on the particle speed.
When muon’s energy is non-relativistic and its speed is suﬃciently slower than light speed, its
Cherenkov light cone angle depends on its momentum as follows:√11m2
cos θ ==
1+ 2
(4.18)nβnpwhere n is the refraction index of water and m is muon mass. This calibration uses the events
that satisfy the criteria below:• The muon entrance point is on the top of the detector.
• cos Θ > 0.9 (Downward-going muon)• One decay electron is detected.
• P Etot <1500 p.e (750 p.e.s for SK-II)
where Θ is the zenith angle of muon direction, and P Etot is the total charge in the inner detector.

4.3.2

Diﬀerence Between MC and Data

To estimate the systematic error of energy scale after the calibration, this calibration data withthe scaled energy is compared with the expected performance based on MC. The comparedvalues are shown as follows:• Reconstructed π 0 invariant mass (Sec. 4.3.1.1)• Reconstructed decay electron momentum (Sec. 4.3.1.2)• (Reconstructed momentum)/(Muon track length) of High energy stopping muon (Sec. 4.3.1.3)• (Reconstructed momentum based on normal method)/(Reconstructed momentum basedon Cherenkov cone angle) of Low energy stopping muon (Sec. 4.3.1.4)
These evaluation covers the events with less than 10 GeV momentum. Fig. 4.16 shows thediﬀerence of these values between data and MC. The largest diﬀerence of all calibration sources
between MC and data is defined as the uncertainty of absolute energy scale.

CHAPTER 4. CALIBRATION

46

Figure 4.16: Diﬀerence of the energy scale between MC and data. The colors of data points
indicate the calibration sources, and each source’s comparison target is defined in Sec. 4.3.2.
The error bars to x-axis are assigned to the deviation of the momentum distributions.

CHAPTER 4. CALIBRATION

4.3.3

47

Time Variation of Energy Scale

SK>I%

43
4.3

42
4.2

41
4.1

Jul.%96

Decay Electron Momentum [ MeV/c ]Decay Electron Momentum [ MeV/c ]

44
4.4

Jul.%97

Jul.%98

Jul.%99

Jul.%00

41
4.1

SK>II%

40
4.0

39
3.9

3.8
38

Jul.%01

Apr.%03

Oct.%03

Apr.%04

Oct.%04

Apr.%05

Oct.%05

07/01/1996
01/01/1997
07/01/1997
01/01/1998
07/01/1998
01/01/1999
07/01/1999
01/01/2000
07/01/2000
01/01/2001
07/01/200

01/01/0304/01/0307/01/0310/01/0301/01/0404/01/0407/01/0410/01/0401/01/0504/01/0507/01/0510/01/05

date

date

43
4.3

SK>III%

Decay Electron Momentum [ MeV/c ]Decay Electron Momentum [ MeV/c ]Because of the change of detector environment, for example the PMT gain, dark noise rate,and so on, the reconstructed energy has time dependence. This eﬀect is evaluated by using the
decay electron and stopping muon measurement results. Fig. 4.17 and Fig. 4.18 show the time
variation of energy scale based on the values for comparison which are mentioned in Sec. 4.3.2.
The largest shift from the mean value in all term of each SK period is assigned to the systematic
error for the time variation of energy scale.

42
4.2

41
4.1

40
4.0

39
3.9

SK>IV%

38
3.8

37
3.7

36
3.6

Nov.%06

Mar.%07

Jul.%07

Nov.%07

Mar.%08

Jul.%08

09/01/06
11/01/06
01/01/07
03/01/07
05/01/07
07/01/0709/01/07
11/01/07
01/01/08
03/01/08
05/01/08
07/01/08

date

01/01/09

Jan.%10 01/01/11

01/01/10

Jan.%12 01/01/13

01/01/12

date

Jan.%14 01/01/15

01/01/14

Jan.16

01/01/16

Figure 4.17: Time variation of reconstructed decay electron momentum. Dotted lines indicate±1% from the mean value in each period. The shaded regions in SK-IV indicate T2K run period,in which neutrino beam comes from J-PARC.

4.3.4

Evaluation Result of Energy Scale Uncertainty

Based on Sec. 4.3.2 and Sec.4.3.3, the total energy uncertainties are evaluated as shown in
Table 4.3. The total energy scale errors in this table are defined as the systematic errors of

Absolute Energy Scale
Time Variation
Total Energy Scale Error

SK-I
3.13%
0.90%
3.26%

SK-II
2.70%
0.56%
2.75%

SK-III
1.62%
1.81%
2.43%

SK-IV
2.08%
0.47%
2.13%

Table 4.3: Systematic error table about energy scale. Total energy scale errors are calculated
combining the absolute energy scale and time variation errors in quadrature.
energy scale. This error does not contains the eﬀect from the position dependence As written
in [136], the uniformity of energy scale in the detector is checked based on the decay electron
measurement.

2.4
2.4

SK:I%

2.35
2.35

2.3
2.3
2.25
2.25
2.2
2.2
07/01/1996
Jul.%9601/01/199707/01/1997
Jul.%9701/01/199807/01/1998
Jul.%9801/01/199907/01/1999
Jul.%9901/01/200007/01/2000
Jul.%0001/01/200107/01/2001
Jul.%01

48

Stop Muon Momentum / Range [ MeV/c/cm ]Stop Muon Momentum / Range [ MeV/c/cm ]CHAPTER 4. CALIBRATION

2.45
2.45

2.4
2.4

2.35
2.35

2.3
2.3

Apr.%03

Stop Muon Momentum / Range [ MeV/c/cm ]Stop Muon Momentum / Range [ MeV/c/cm ]SK:III%

2.45
2.45

2.4
2.4
2.35
2.35

2.3
2.3
10/01/06
01/01/07 04/01/07 07/01/07
10/01/07 01/01/08 04/01/08 07/01/08
Oct.%06
Apr.%07
Oct.%07
Apr.%08

date

Oct.%03

Apr.%04

Oct.%04

Apr.%05

Oct.%05

01/01/0304/01/0307/01/03 10/01/0301/01/0404/01/0407/01/04 10/01/0401/01/0504/01/0507/01/0510/01/05

date

2.5
2.5

SK:II%

2.4
2.4

date

SK:IV%

2.35
2.35

2.3
2.3

2.25
2.25
2.2
2.2
01/01/09

Jan.%10 01/01/11

01/01/10

Jan.%12 01/01/13

01/01/12

date

Jan.%14 01/01/15

01/01/14

Jan.16

01/01/16

Figure 4.18: Time variation of ratio of reconstructed stopping muon momentum to its track
length. Dotted lines indicate ±1% from the mean value in each period. The shaded regions in
SK-IV indicate T2K run period, in which neutrino beam comes from J-PARC.

Chapter 5Monte Carlo Simulation of
Atmospheric Neutrino Events
For neutrino oscillation analysis, the measurement result is compared with the output of Monte
Carlo simulation. This simulation consists of an atmospheric neutrino flux model, a neutrino
interaction model, and a detector simulator. This section describes these models and the simulator.

5.1

Atmospheric Neutrino Flux

The atmospheric neutrino flux simulation follows the calculation result of M. Honda et al. at
2011[50]. This calculation uses the several kinds of models to estimate the flux of the incoming
primary cosmic ray, and the neutrino yield from the interaction between atmosphere and theprimary cosmic ray particles. As a first step, the primary cosmic ray flux is calculated basedon AMS (Alpha Magnetic Spectrometer [51], [52]) and BESS (Balloon-borne Experiment witha Superconducting Spectrometer [53], [54]) data[55], [56]. The charged particles incoming to
the earth are bent by the geo-magnetic field. Because the magnetic field around the earth isdirected north, the cosmic-ray particles with positive and negative charges are bent to eastward
and westward, respectively. As a result, the azimuthal angle dependence of atmospheric neutrino
flux is not symmetric. Fig. 5.1 shows the azimuthal angle dependence of incident neutrino flux.
To take this eﬀect into account, IGRF2005 geo-magnetic model[58] was introduced into thecalculation. Atmospheric neutrinos are produced by the interaction between the cosmic ray andthe atmosphere. The atmosphere density follows US-standard’76 model[59]. Their hadronicinteractions follow JAM and DPMJET-III models. DPMJET-III[60] is the hadron interaction
event generator based on Dual Parton Model[61], and JAM is the interaction model used in
PHITS[62]. DPMJET-III is used for events whose primary cosmic-ray has more than 32 GeVenergy, and JAM is used for events with less than 32 GeV. These models were modified toagree with the experimental data based on [63] [50]. The atmospheric muon fluxes based onthese interaction models are compared with the data, which are taken at Tsukuba[64], Mt.
Norikura[65], and balloon altitude at Fort Sumner [66]. As shown in Fig. 5.2 and Fig. 5.3, thismodel agrees with the data taken at balloon altitude and Tsukuba at low energy region where
JAM is used. Based on this model, the neutrino flux at Kamioka is calculated. The calculated
neutrino energy spectra are shown in Fig. 5.4.

49

CHAPTER 5. MONTE CARLO SIMULATION OF ATMOSPHERIC NEUTRINO EVENTS50

Figure 5.1: The azimuthal angle dependence of atmospheric neutrino flux at 1 GeV. This figure
is taken from [50].

Figure 5.2: Comparisons of µ± fluxes between the data taken at Tsukuba and the calculation
based on the modified JAM [50]. The left figure shows their momentum distribution, and theright one shows the ratio of data to MC. The red dots solid lines in the left figure show theobserved and calculated µ+ fluxes, and the blue ones indicate their µ− fluxes, respectively. The
dots in the right figure show the ratio of data to the introduced model, which is based on JAM,
and broken lines show the ratio of this model and DPMJET-III.

CHAPTER 5. MONTE CARLO SIMULATION OF ATMOSPHERIC NEUTRINO EVENTS51

Figure 5.3: Comparisons of µ± fluxes between the data taken at balloon altitude at FortSumner
and the calculation based on the modified JAM [50]. The meaning of dots and lines is the sameas that of Fig. 5.2.

Figure 5.4: Atmospheric neutrino fluxes based on the calculation result of Honda used for SK
Monte Carlo, older calculation result of Honda (HKKMS06), the model of Bartol group[67], andthat of Fluka[68]. This plot is taken from[50].

CHAPTER 5. MONTE CARLO SIMULATION OF ATMOSPHERIC NEUTRINO EVENTS52

5.2

Neutrino Interaction ModelSK observes signals from the interactions between neutrinos and water. In this thesis, thisinteraction is simulated based on NEUT version 5.3.6. Because the simulation result largely
depends on the interaction model, the uncertainty has to be considered. Especially, there are nodata about the number of neutrons produced by neutrino interactions. Therefore, the corrected
neutrino energy based on NEUT 5.3.6 is compared with the corrected energies based on theother models, which are NEUT version 5.4.0 and GENIE version 2.12.6. The comparison result
is described in Sec. 8. In this section, these interaction models are described.

5.2.1

NEUT

NEUT is a simulation program library of neutrino interactions, which is originally developed
for Kamiokande study[69]. NEUT does not take only the neutrino interactions, but also nuclear
eﬀect to the produced particles into account. In this section, these interactions are describedbased on NEUT version 5.3.6 as the first half. After that, the diﬀerence between version 5.3.6
and 5.4.0 is described.
5.2.1.1

Neutrino Interactions Based on NEUT 5.3.6

1. (Quasi-)Elastic Scattering (νl + N → l + N ′ )The cross section of charged current quasi-elastic scattering with a free proton follows
Llewellyn-Smith’s model[70]. On the other hand, the interaction with oxygen cannot be
described by this model because of the nuclear in-medium eﬀect. NEUT 5.3.6 calculates
this interaction following the global relativistic Fermi gas model of Smith and Moniz[71].
This cross section depends on the parameter for the axial vector form factor, which ismodeled as a dipole form factor parameterized by the axial vector mass MAQE . MAQE of
NEUT5.3.6 is set to 1.21 GeV/c2 based on the results of K2K[72] and MiniBooNE[73].
Fig. 5.5 shows the cross sections of the charged current quasi-elastic scattering estimatedwith NEUT and data. The simulated cross section seems to be diﬀerent from data. As
mentioned in the following section, this diﬀerence is suppressed in NEUT version 5.4.0 bychanging the interaction model.
2. Interaction Producing Single Meson or γ (νl + N → l + N ′ + meson or γ)
Single meson and γ production in the region of hadronic invariant mass (W ) lower than
2 GeV is simulated using a resonant model based on Rein and Sehgal’s method with
Garczyk-Sobczyk form factor[75][76]. This method considers the interaction with single
pion production, and assumes this interaction occurs via nucleon resonance and its decay(νN → lN ∗ , N ∗ → N π). Fig. 5.6 shows the relationship between neutrino energy and thecross sections of the interactions producing single pion. The interactions producing single
photon, K or η production are also modeled based on this method, changing the decay
probability of the resonance nucleons.
3. Deep Inelastic Scattering (νl + N → l + N ′ + hadrons)
NEUT models the deep inelastic scattering considering the parton distribution in the
nucleon. The cross section of the charged current deep inelastic scattering is calculated as

CHAPTER 5. MONTE CARLO SIMULATION OF ATMOSPHERIC NEUTRINO EVENTS53

Figure 5.5: Comparison of the charged current quasi-elastic scattering cross sections between
data and NEUT[74]. The solid and broken lines are assigned to the calculated cross sections forthe interactions with free nucleons and bound nucleons in oxygen, respectively. These figures
are taken from [74]Figure 5.6: Comparison of the cross sections of single pion productions between data andNEUT[74]. The solid red and black lines indicate the cross sections based on NEUT5.3.6 andolder version of NEUT, respectively. These figures are written based on [74]

CHAPTER 5. MONTE CARLO SIMULATION OF ATMOSPHERIC NEUTRINO EVENTS54
follows[77]:
)
((
)
(
)
( 2)
( 2)
G2F MN Eν
d2 σ
y2
y
=
1−y+
+ C1 F2 x, q ± y 1 − + C2 xF3 x, q
dxdy
π
2
2
2
2
2
yMl
m
ml
xyMN
(5.1)
C1 =
−
− l2 −
4MN Eν x
2Eν
4Eν
2MN Eν x
m2l
C2 = −
4MN Eν x
where q 2 is the 4 momentum transfer from the interacted neutrino to the hadronic system,
MN and Ml are the masses of the interacted nucleon and lepton, Eν and El are their energy,
and x and y satisfy x = −q 2 /(2M (Eν − El )) and y = (Eν − El )/Eν . F2 and xF3 are called
structure functions, and they are taken from GRV98 parton distribution function[78]. To
improve the agreement with the real data in the low q 2 region, this function is correctedbased on the method of Bodek and Yang[79]. On the other hand, the cross sections of
the neutral current deep inelastic scatterings are estimated based on the relation withthat of this charged current deep inelastic scattering, which is evaluated based on theexperimental results[80] [81]. Two kinds of the interaction models are prepared for the
events whose interacted nucleon’s invariant masse (W ) are heavy and light. As a firststep, the interaction with light invariant mass (1.3 GeV < W < 2.0 GeV) is considered.
The interactions with single meson production are a part of the deep inelastic scattering,but the interaction with W < 2 GeV has been introduced independently as describedin the previous item. To remove the single meson production with W < 2 GeV fromthe deep inelastic scattering process, NEUT uses the probability function of the pion
multiplicitiy calculated based on the result of FermiLab 15-foot hydrogen bubble chamber
experiment[82]. In this region, all of the outgoing mesons are considered to be π. As a next
step, W > 2 GeV interactions are described. In this region, this process is modeled basedon PYTHIA and JetSet[83]. Fig. 5.7 shows the comparison of the inclusive cross sectionsof the charged current interactions between data and modeling which includes quasi-elastic
scattering, single meson production, and deep inelastic scattering. The interaction with
W < 1.3 GeV are not considered in this process.
4. Neutral (Charged) Current Coherent π Production (ν +16 O → ν(l± ) +16 O + π 0 (π ∓ ))This interaction produces π without breaking the interacted oxygen nucleus. This is modeled based on the Rein and Sehgal method[85]. However, the measurement results of K2K
and SciBooNE indicate the charged current coherent π production cross section is muchsmaller than this model as shown in Fig. 5.8.
5. Meson Exchange Current (νl + N N ′ → l + N ′′ N ′′′ )In the neutrino interactions, the correlation between correlated pair of nucleons is notconsidered as an approximation. To take the interaction between neutrino and the multiple
nucleons into account, Meson Exchange Current (MEC) model [93] is introduced for the
charged current interaction. This process considers the interaction with nucleons as follows:νl + N N ′ → l + N ′′ N ′′′where N , N ′ , N ′′ , and N ′′′ are nucleons.

(5.2)

CHAPTER 5. MONTE CARLO SIMULATION OF ATMOSPHERIC NEUTRINO EVENTS55

Figure 5.7: Comparison between charged current inclusive cross sections summarized by
PDG[84] and the cross section estimated based on NEUT 5.3.6. The red line is assigned tothe cross section based on NEUT5.3.6.

CHAPTER 5. MONTE CARLO SIMULATION OF ATMOSPHERIC NEUTRINO EVENTS56

Figure 5.8: Comparison of the coherent π production cross sections between results of theexperiments and models. The left and right figures show the cross sections of the charged andneutral current interactions with π production, respectively. The red and blue points in the left
figure are assigned to the measurement results of K2K[86] and SciBooNE[87]. The circle, box,
and cross points in the right figure are assigned to the measurement results of MiniBooNE[88],Aachen-Padova[89], and Gargamelle[90]. The solid and broken lines in both figures are assigned
to the cross sections based on Rein-Sehgal model with and without lepton mass correction,
respectively. The dotted and dash-dotted lines in the left figure are assigned to the cross section
based on other models, which are the models of Kartavtsev et al.[91] and Alvarez-Ruso et al.[92],respectively. These figures are taken from [74].
5.2.1.2

Nuclear Eﬀects

The hadrons produced by the neutrino interaction in oxygen nucleus may interact with its
nucleons before leaving it. NEUT calculates this process based on the cascade model, in which
the produced hadrons propagate step by step in the nucleus. In each step, the interaction with
nucleus is possible to happen following their interaction probabilities. In this process, the nucleon
density is based on Woods-Saxon type distribution[94]. This model is slightly diﬀerent between
the produced particles. The nuclear eﬀect for pion is modeled based on the four types of pion
interactions, which are inelastic scattering, charge exchange, absorption, and particle production.
When pion momentum is less than 500 MeV, its interaction probability are calculated based on
the method of [95]. The calculated probabilities are tuned by comparing the calculated crosssections with various experimental data, which are summarized in [96]. The interacted pion’s
momentum and direction are calculated based on the phase shift analysis result of πN scatteringexperiments[97] with correction[98]. When the pion momentum is higher than 500 MeV/c, its
interaction probability is calculated as follows:( )ρPi = ρ0 σi
dx
(5.3)
ρ0
where σi is each cross section of elastic scattering, charge exchange, and hadron production. dx
is the step size of calculation. ρ0 and ρ are the nucleon densities at the center and the positionof the pion. The σi for π ± scattering is modeled based on the interaction with free proton anddeuteron, which are given by the fitting result of SAID group [99] [100] [101]. The cross sections

CHAPTER 5. MONTE CARLO SIMULATION OF ATMOSPHERIC NEUTRINO EVENTS57

Figure 5.9: π + and π − interaction cross sections on 12 C[96].
calculated based on these models show good agreement with measured data as shown in Fig. 5.9.
The nuclear eﬀects for K, η, and nucleons also use the cascade model based on the experiment
data for each particle [102] [103] [104] [105] [108].
5.2.1.3

NEUT Improvement

In this thesis, the oscillation analysis is done by using NEUT5.3.6 as an interaction model.
However, NEUT version 5.4 has been already produced, adding some improvements. This new
version NEUT is used for the systematic error estimation, as mentioned in Sec. 8.2.5. In thissection, the updates of NEUT5.3.6 to NEUT5.4.0 are described.
Quasi-Elastic Scattering Model
The interaction model is changed from global relativistic Fermi-Gas model to local FermiGas model. This new model takes the position dependence of the nucleon density into
account. In addition, the axial mass is changed from 1.21 GeV/c2 to 1.05 GeV/c2 . As aresult of these changes, the calculated cross section decreased.
Deep Inelastic Scattering
For the low W model, Bodek-Yang’s correction method was updated based on [110]. In
addition, CKM matrix elements are newly introduced in the calculation of the structure
functions to rotate the mass eigenstate of quarks to the flavors for weak interaction. The
neutral current interaction also uses the same structure function without CKM matrix
elements. The detailed description is given in [111].
Coherent π Production
As mentioned in Sec. 5.2.1.1, the current model has much larger cross section than themeasurement result. To correct this, the cross section model is changed to Berger andSehgal’s model[112]. Fig. 5.10 shows the coherent pion cross sections as a function of thetransmitted momentum from the interacted neutrino before and after the modification.
The calculated cross section is reduced after the modification.
Nuclear Eﬀect
As mentioned in Sec. 5.2.1.2, the cross section related to the nuclear eﬀects are calculated

CHAPTER 5. MONTE CARLO SIMULATION OF ATMOSPHERIC NEUTRINO EVENTS58

Figure 5.10: The cross sections of coherent pion production as a function of the transmitted
momentum [113]. The black and blue lines are assigned to the output of Rein-Sehgal andBerger-Sehgal models, respectively.
based on the experiment data. As an additional data of π ± A scattering, the result ofDual Use Experiment at TRIUMF (DUET)[114] [115] is introduced to the cross section
evaluation. DUET measured the interaction between carbon target and the π + beams
with five diﬀerent momenta from 201.6 MeV/c to 295.5 MeV/c, and their cross sectionsof absorption and charge exchange are analyzed.

5.2.2

GENIE

The systematic error evaluation uses not only NEUT 5.3.6 and NEUT 5.4.0, but also GENIE 2.12.2 for comparison. GENIE is a neutrino Monte Carlo generator developed by theexperimental neutrino physics community [120]. GENIE has some final state interaction (FSI)models[134]. In this analysis, ”hA” FSI model is used for the simulation with GENIE. Its output
is diﬀerent from the output of NEUT. For example, its neutron multiplicity tends to be larger
than NEUT models as shown in Fig. 5.12. This figure shows the comparison of the neutron
multiplicities between NEUT5.3.6, NEUT5.4.0 and GENIE. Because the neutrino energy reconstruction uses the number of tagged neutrons, this diﬀerence is possible to change the corrected
neutrino energy. The systematic error of the corrected neutrino energy is estimated in Sec.8.2.5.

5.3

Detector Simulation

The particles produced by the neutrino interaction are detected based on the signals of PMTs
in the tank. The distribution of signals is estimated by simulating the behavior of the produced
particles and the Cherenkov light production. This simulation is done by using SKDetSim,
which is the detector simulation tool based on GEANT3[106]. SKDetSim simulates hadron in-CHAPTER 5. MONTE CARLO SIMULATION OF ATMOSPHERIC NEUTRINO EVENTS59

Event Rate

Figure 5.11: The cross section comparison of π + ’s absorption and charge exchange with carbon
particles [114]. The red and black colors are assigned to the cross sections of absorption andcharge exchange, respectively. The filled and hollow dots are assigned to DUET and the otherexperimental data[116][117][118][119], respectively. The solid, broken, and dotted lines are assigned to the model prediction based on GEANT4, FLUKA and NEUT before the tuning basedon the DUET data in this figure, respectively.

0.5
NEUT5.3.6
GENIE2.12.2

0.4

NEUT5.4.0
0.3
0.2
0.1
0
0

1

2

3

4

5

6
7
8
9 10
Number of Neutrons

Figure 5.12: The comparison of the number of produced neutrons between models. Blue, redand green lines are assigned to the calculation results of NEUT5.3.6, GENIE2.12.2 (with hA
model), and NEUT5.4.0, respectively. The average numbers of the produced neutrons are largely
diﬀerent between the models.

CHAPTER 5. MONTE CARLO SIMULATION OF ATMOSPHERIC NEUTRINO EVENTS60
teractions based on GCALOR[107] and NEUT. The pions with less than 500 MeV/c momentum
are simulated based on the cascade model of NEUT, and the other particles are simulated basedon GCALOR.

5.3.1

Background Data Combination for Neutron Tagging Analysis

As described in Sec. 8.1.1, in SK-IV period, neutron is searched from not only SHE trigger
region but also AFT trigger region by searching for the signal from neutron capture, which iscalled neutron tagging. However, the source of the background for this process has not been
known quantitatively yet. Therefore SKDetSim cannot simulate this background. Instead of
the simulation, the background data is taken by using the dummy trigger data. This trigger
is produced periodically even when it does not contains any neutrino signal. By combining
the signals taken by this trigger and the neutron capture signals which is simulated without
any background, the PMT signals in the neutron search region can be reproduced. Fig. 5.13

Figure 5.13: The structure of the background for MC as a function of the elapsed time. The
dummy trigger data is introduced from 18 − ttrg [µs] to 535 − ttrg [µs] after the SHE trigger forthe neutrino interaction. ttrg is the triggered time from the start of the simulation, which isabout 1 µs. Until 18 − ttrg [µs], the background is simulated assuming it comes from the darknoises of the PMTs, which are independent of each other.
shows the schema of the MC and dummy trigger data combination. The dummy trigger data is
introduced after 17 µs from the neutrino interaction.

5.3.2

π − Capture

When the produced π − loses its energy, it is captured by water. In SKDetSim, this process
is simulated without neutron emission. In practice, this process is possible to produce someneutrons. This diﬀerence changes the hadron multiplicities and is possible to change the reconstructed neutrino energy. This eﬀect is estimated in Sec. 8.2.5.4.

Chapter 6

Event Reconstruction Based on The
Primary Signal
The neutrino oscillation probability depends on its flavor, energy and travel length of the neutrino. Ideally, the oscillation analysis should be done by using the number of events in thebins separated based on these physical quantities, but they cannot be reconstructed directly.
Alternatively, the current oscillation analysis uses the reconstructed momenta, directions andflavors of the produced charged particles, which emit the Cherenkov light in the water. Theseobservables reflect the physical quantities of the interacted atmospheric neutrino. They are reconstructed based on the charge and timing information of the ID PMTs. The reconstruction
tool is called APFit. In this section, the method of the event reconstruction is described.

6.1

Vertex and Direction Fitting

6.1.1

Point Fit

First of all, the vertex is fitted assuming all Cherenkov lights are emitted from the same point.
The distributions of the time-of-fright subtracted hit time, which is called residual time, arecomputed for various hypothetical vertices. When the assumed vertex is close to the neutrino
interaction vertex, the sharp peak appears. The reconstructed vertex is defined as the position
with the maximum goodness that is defined as following:()
1 ∑(ti − t0 )2
Gpfit =
(6.1)
exp −
Nhit
2(1.5σPMT )2iwhere Nhit is the number of hit PMTs, ti is the residual time of the i-th PMT, σPMT is the
typical PMT timing resolution of 2.5 ns, and t0 is a free parameter, which is shifted for Gp tohave the maximum value. The hypothetical vertex is searched in the grid whose intervals areabout 30 cm.

6.1.2

Ring Edge SearchThis step searches for the edge of Cherenkov rings. This information is used for more precisevertex fitting. In this step, the directions of produced particles are not known yet. To estimate
the direction and edge of Cherenkov rings, this step follows the procedure below.
61

CHAPTER 6. EVENT RECONSTRUCTION BASED ON THE PRIMARY SIGNAL

62• Make the opening angle distribution of PMT signals. q(θ) is assigned to the angular
distribution of the observed charge as a function of the opening angle.
• Evaluate the opening angle of Cherenkov ring edge θedge . θedge is defined as follows:θedge > θpeak2
d q(θ)
dθ2

θedge

(6.2)

=0

(6.3)

where θpeak is the peak of q(θ). Fig. 6.1 shows the schema of the opening angle distribution
and the position of θpeak and θedge . When multiple angles satisfy these criteria, the nearest
angle to θpeak is set to θedge

Figure 6.1: Schema of Cherenkov ring opening angle distribution with θpeak and θedge positions.
• Calculate the goodness below:∫ θedgeGθ (θedge ) =0q(θ)dθ
×sin θedge(dq(θ)
dθ θ=θedge)2)((θedge − θexp )2
× exp −
2σθ2

(6.4)

where θexp is the expected Cherenkov angle based on the total charge inside the Cherenkov
cone, σθ is its resolution, and q(θ) is the angular distribution of the detected charge.
These steps are done repeatedly until maximizing Gedge , the direction and edge angle with the
maximum Gedge are chosen as the parameters of the primary Cherenkov ring.

6.1.3

TDC-FitAfter the ring edge search, the vertex is reconstructed again with the method called TDC-fit.
This method takes the particle’s track length and scattered photons into account. The PMTs
inside and outside the cone are separated by comparing the opening angle with θedge , which is
calculated in the previous section. The residual hit time of PMTs inside and outside the cone

CHAPTER 6. EVENT RECONSTRUCTION BASED ON THE PRIMARY SIGNAL

63

is calculated assuming the light is emitted along the track and at the vertex, respectively. It
follows the equation below:

1
n
0

⃗
⃗
⃗
⃗

 ti − c × |Xi − O| − c × |Pi − Xi | (inside the Cherenkov ring)
ti =
(6.5)


 t0 − 1 × |P⃗i − O|⃗(outside the Cherenkov ring)ic⃗ is the vertex position, X⃗ i is the emission position of the Cherenkov light toward thewhere Oi-th PMT, n is the refractive index of the water, P⃗i is the position of the i-th PMT, and t0i is
the hit time of the i-th PMT. Based on this residual time, the vertex is fitted by searching forthe position with the maximum goodness. In this analysis, the contributions of PMTs insideand outside the Cherenkov cone to the goodness are calculated with diﬀerent methods. The
contribution of the PMTs inside the Cherenkov cone is defined as:
)
(

∑
1
(ti − t0 )2(for ti ≤ t0 )
exp −

i 2

2 × (1.35 × σ)2σi (qi )
GI =
(6.6)
(
)

∑
1
(ti − t0 )2



exp −
(for ti > t0 )i 2
2 × (1.5 × σ)2σi (qi )where σi (qi ) is the timing resolution of the i-th PMT depending on its charge qi , t0 is the meanof the ti , and σ is the timing resolution for the charge averaged over all PMTs with any signals.
When the PMT is outside the Cherenkov ring, the definition is changed. When the hit time islater than the mean (=t0 ), the contribution from the scattered photons is taken into account.
(
(
)
)

∑
1
(ti − t0 )2


exp −
−1 ×2−1(for ti ≤ t0 )
i 2
2(1.5 × σ)2σi (qi )
GO =
(6.7)
(
[
(
)
]
)

∑
1
(ti − t0 )2max exp −, Gscatt (ti , t0 ) × 2 − 1(for ti > t0 )i 2
2(1.5 × σ)2σi (qi )Gscatt (ti , t0 ) follows:
)
(
) (
)
(
Rq
Rq
(ti − t0 )2
ti − t0
Gscatt (ti , t0 ) =× exp −
+ 1−
× exp −
1.52
2(1.5 × σ)2
1.52
60 nsec(6.8)

where Rq is the fraction of the charge inside the Cherenkov ring. The numerical factors in theseequations are optimized based on the Monte-Carlo study. Based on these contributions, thegoodness GT is calculated as follows:
GI + GOGT = ∑2i (1/σi (qi ))(6.9)

After that, single ring events are reconstructed with a more precise method, which is called
MS-fit. It is described in Sec. 6.4.

6.2

Ring Counting

To count the number of Cherenkov rings, Hough transform is executed. In this analysis, theCherenkov opening angle is assumed to be 42◦ , which is the Cherenkov angle of particles travelingCHAPTER 6. EVENT RECONSTRUCTION BASED ON THE PRIMARY SIGNAL

64

through water with the speed of light in vacuum. If the PMT signal detects the photons fromthe Cherenkov radiation, the light source should be on the surface of the 42◦ cone whose cone
axis is antiparallel to the Cherenkov cone, and conical vertex is the PMT position. Fig. 6.2
shows the schema of this idea. If the Cherenkov cone axis is set to the true direction, the cones

42 deg. ringhit PMT(possible center)

centerCherenkov ring

(most probable)

Figure 6.2: The schema of the ring search idea. By setting each PMT position to conical vertex,
the cones with 42◦ opening angle with antiparallel direction with Cherenkov ring direction are
generated. By searching for the position where the edges of these cones are concentrated. This
figure is taken from [121].
of the PMTs will intersect at the vertex. In practice, the search of the intersect points of thecones are done by summing up the PMT charges weighted by the expected charge distributionfunction f (θi ). where θi is the opening angle between the particle direction and the direction
from the vertex to i-th PMT. The total weighted charge is calculated for each zenith angle Θ andazimuthal angle Φ of the particle, and it is mapped on a (Θ, Φ) plane. Fig. 6.3 is an example.
This map contains two peaks. They show the ring candidates and their directions. When N
rings has been already selected, the new candidates are tested by comparing the likelihoods with
N and N + 1 rings, assuming the candidate ring as (N + 1)-th ring. The likelihood with N rings
is defined as:
N
∑
∑
exp
obs
LN =
log(P (qi ,αn qi,n))(6.10)in=1where the likelihood of hits inside any Cherenkov rings are summed up, qiobs is the observed
exp
charge of i-th PMT, αn qi,nis the expected charge of i-th PMT from n-th ring without particle
obs
exp
type assumption, P (q , q ) is a probability function. The value of scale factor αn (n =1, 2, · · · , N ) are set to have maximum LN with a constraint of lower momentum limit. If theq exp is less than 20 p.e.s, P (q obs , q exp ) follows the probability density function based on theCHAPTER 6. EVENT RECONSTRUCTION BASED ON THE PRIMARY SIGNAL

65

Figure 6.3: A map of the charge weighted based on the opening angle dependence of theCherenkov light. Two peaks in this figure are regarded as the candidates of the Cherenkov
ring. This figure is taken from [121].

CHAPTER 6. EVENT RECONSTRUCTION BASED ON THE PRIMARY SIGNAL

66

convolution of a Poisson distribution and the single photoelectron charge distribution. Whenq exp is more than 20 p.e.s, P (q obs , q exp ) follows:()
1(q obs − q exp )obs exp
P (q , q ) = √exp −(6.11)2σ 2
2πσ
When there are not any candidates which satisfy LN +1 ≥ LN , the number of rings is determined
to be N . When the requirement is satisfied, another likelihood evaluation is executed. In thisevaluation, the values below are compared with the expected value from MC.
1. the diﬀerence between L(N + 1) and L(N )2. the average of the expected charge near the (N + 1)-th Cherenkov ring edge (Qedge ),excluding the expected charge from the other rings.
3. the average of the expected charge outside the (N + 1)-th ring and Qedge , excluding the
expected charge from the other rings.
4. the diﬀerence between each PMT’s observed and expected charge assuming only N rings
exist.
5. the diﬀerence between the average of the expected charge around the peak of (N + 1)-thring and the average of the expected charge inside and outside the detector, excluding the
expected charge from the other rings.
Based on each parameter, the likelihood of (N + 1)-ring pattern is output respectively. Whenthe total log-likelihood is larger than N ring pattern, the new ring candidate is regarded asa Cherenkov ring. Fig. 6.4 shows the ring counting likelihood distribution of the FC eventswith multi-GeV momentum as an example. The Data/MC ratio of the multi-GeV single andmulti-ring events in this figure are 92 ± 2% and 98.5% ± 3%. The diﬀerence between the single
ring events based on data and MC comes from the unsuitable interaction model for MC, which
is described in Sec. 5.2.1.1.

6.3

Particle Identification

The flavor determination of the lepton is crucial in neutrino oscillation analysis. For this analysis,
the shapes of Cherenkov rings are classified into 2 types, electron like and muon like.
Because of the electromagnetic shower and the multiple scattering eﬀect, the Cherenkov ringof electron tends to be smeared. On the other hand, the ring of muon tends not to be smearedby these eﬀects, owing to its heavy mass. As a result, the Cherenkov ring of muon shows sharperedge than that of electron. To classify the rings, each PMT’s observed charge is compared withthe expected charge from the Cherenkov radiation of electron and muon, and the candidate withhigher likelihood is selected. Assuming electrons to point like source of Cherenkov radiation,
the expected charge of i-th PMT from an electron is calculated as follows:
qiexp (e) = αe × Qexp (pe , θi ) ×
where the variables are defined as follows:( )1.5R
1
×× f (Θi ) + qiscattri
exp(ri /L)(6.12)

CHAPTER 6. EVENT RECONSTRUCTION BASED ON THE PRIMARY SIGNAL

67

Figure 6.4: The distribution of ring counting likelihood diﬀerence between 1-ring and 2-ring for
multi-GeV events. The black dots and red line are assigned to the distributions of the data andMC with the same live time.

CHAPTER 6. EVENT RECONSTRUCTION BASED ON THE PRIMARY SIGNAL

68• αe : a normalization factor• ri : the distance from the fitted vertex to the i-th PMT• θi : the opening angle between the i-th PMT and the ring direction• pe : the momentum of electron, which is reconstructed• L: light attenuation length in water• f (Θi ): PMT acceptance correction factor as a function of the photon incident angle Θi• Qexp (pe , θi ): the expected charge as a function of the opening angle and the electron
momentum.
• R: radius of the virtual sphere (16.9m)
The light intensity at the surface of the PMT is calculated by a factor of (R/ri )1.5 exp(−ri /L).
the power index of (R/ri ) is determined based on a Monte Carlo simulation.
Taking the long and straight track into account, the expected charge of i-th PMT from amuon is calculated as follows:()
sin2 θxi
1exp
qi (µ) = αµ ×× f (Θi ) + qiscatt(6.13)+ qiknock ×dθx
exp(ri /L)ri (sin θxi + ri dx |x=xi )where the variables are defined as follows:• αµ : a normalization factor• xi : muon path length from the vertex to the points of Cherenkov light emission toward
the i-th PMT• θxi : opening angle between muon direction and Cherenkov light emitted at track position
with xi path length from the vertex.
• qiscatt , qiknock : expected charge for the i-th PMT from scattered photons and knock-on
electrons, which depend on the charge of each ring and its distribution, and also depends
on the vertex.
As mentioned in Eq. (3.3), the intensity of Cherenkov light emitted during a unit track length
depends on sin2 θxi . As shown in Fig. 6.5, the Cherenkov light is emitted from the track withdx length to the area of 2πr sin θ(dx sin θ + rdθ). To take this into account, the terms ofri (sin θxi + ri (dθx /dx)|x=xi ) is introduced into Eq. (6.13). By using these expected charges, theelectron and muon likelihoods are evaluated as follows:


∏
∑
exp
exp 
L(e or µ) =P qiobs , qi,n(e or µ) +qi,n(6.14)′
θi <(1.5×θC )n′ ̸=nexp
where P (q obs , q exp ) is the same probability function shown in Sec. 6.2, qi,n(e or µ) is the expectedexp
charge from n-th ring, which is the identification target, assuming the particle is e or µ, and qi,n′′is the expected charge from the n -th ring, which is not the identification target. This likelihood

CHAPTER 6. EVENT RECONSTRUCTION BASED ON THE PRIMARY SIGNAL

69

Figure 6.5: Schema of Cherenkov light emission area.
is translated to the probabilities of electron and muon which is evaluated from the ring pattern(Pnpattern (e or µ)). It follows:()(χ2n (e or µ) − min[χ2n (e), χ2n (µ)])2
patternPn(e or µ) = exp −(6.15)2σχ2 2nwhere
χ2n (e or µ) = −2 log Ln (e or µ) + constant(6.16)√and σχ2n is the resolution of the χ2n distribution, which is approximated by σχ2n = 2N , where N
is the number of PMTs used in this calculation. On the other hand, the probabilities from theCherenkov opening angle are also estimated for the single ring event reconstruction as follows:()(θnobs − θnexp (e or µ))2
angle
Pn (e or µ) = exp −(6.17)2(δθn )2
where θnobs is the reconstructed opening angle of n-th ring, and δθ is its fitting error. By using
these probabilities, the particle identification (PID) likelihood is defined as:√√
PP ID = − log P (µ) − − log P (e)(6.18)
where P (e or µ) is defined as:{ patternP(e, µ) + P angle (e, µ) (single-ring event)P (e, µ) =P pattern (e, µ)(multi-ring event)(6.19)

CHAPTER 6. EVENT RECONSTRUCTION BASED ON THE PRIMARY SIGNAL

70

Figure 6.6: Particle identification likelihood distribution of single ring multi-GeV even of SK-IV.
The black dots and red line are assigned to the distributions of the data and MC with the same
live time.

CHAPTER 6. EVENT RECONSTRUCTION BASED ON THE PRIMARY SIGNAL

71

Figure 6.7: Particle identification likelihood distribution of all the rings in multi ring multi-GeV
events. The black dots and red line are assigned to the distributions of the data and MC withthe same live time.

CHAPTER 6. EVENT RECONSTRUCTION BASED ON THE PRIMARY SIGNAL

72

Fig. 6.6 and Fig. 6.7 show the particle identification likelihood distributions of all detectedparticles in the multi-GeV events with single ring and multiple rings, respectively. This e/µclassification is applied even when the particle is not e or µ. For example, charged pions areclassified to µ-like because it has heavy mass and produces a Cherenkov ring with sharp edge.

6.4

Precise Single Ring Fitting

The vertex fitting in Section 6.1 uses only the timing information, but it has worse resolution in
the particle direction compared to that in the perpendicular direction. For more precise fitting,the ring pattern is used for the fitting of single ring event. This fitting software is called MSFit. This uses the same likelihood as Eq. (6.14), which is the likelihood for PID. The vertex isadjusted in parallel to the particle direction, and the position with the maximum likelihood isdefined as the vertex. To check the performance, the distribution of the distance between true
and reconstructed vertex is calculated for the vertex resolution analysis. The angular resolution
analysis is also done with the similar method. As shown in Fig. 6.8, the upper limit of the

Figure 6.8: The distributions for SK-IV multi-GeV single ring electron like FC events basedon MC. (Left:The distance between true and reconstructed vertex. Right: The opening angle
between true and reconstructed direction.) The green regions contain 68% of all events in the
distribution, and its upper limit is defined as their resolutions.
distance and angle regions with 68% events of all are defined as their resolutions. They aresummarized in Table 6.1.

6.5

Momentum Reconstruction

The charged particle momentum is reconstructed based on the total charge in each ring. Multiring events distribute the charge of each PMT to the rings based on the rate of the expectedCHAPTER 6. EVENT RECONSTRUCTION BASED ON THE PRIMARY SIGNAL
Vertex resolutions [cm]SK-I SK-II
Single ring subGeV
e-like 31.3
35.4
(Plepton <1.33 GeV) µ-like 23.8
30.4
Single ring multiGeV e-like 35.5
36.0
(Plepton >1.33 GeV) µ-like 24.9
26.9

SK-III
30.9
23.6
35.6
24.9

SK-IV
31.1
22.8
35.1
23.6

Angular resolutions [deg]SK-I SK-II
Single ring subGeV
e-like3.1
3.2(Plepton <1.33 GeV) µ-like
1.9
2.2
Single ring multiGeV e-like
1.3
1.4(Plepton >1.33 GeV) µ-like
0.9
1.0

SK-III
3.1
1.9
1.3
0.8

SK-IV
3.1
1.9
1.3
0.8

73

Table 6.1: Vertex and angular resolutions of fully contained single ring events after the precise
single ring fitting (MS-Fit).

charge:

exp
qi,n
obs
obs
qi,n = qi × ∑ exp
n′ qi,n′(6.20)

exp
where qi,nis the charge of i-th ring expected to come from n-th ring. To suppress the eﬀect of
the background such as decay electron, only the PMTs which passed the cuts below is used forthe momentum reconstruction.

θi,n < 70◦ , −50 ns < t < 250 ns(6.21)

where θi,n is the opening angle between the direction of n-th ring and the i-th PMT, and t is theresidual time interval between the peak and the i-th PMT hit. The following equation shows
the n-th ring charge with the correction of the water attenuation, the PMT acceptance, and thescattered light contribution:
(
)
) N∑
N∑
PMT (
PMT
cos Θi
GMC
obs
ri /L
qi,n × e
×RT OTn =α×−Si(6.22)Gdata
f (Θi )ii

where:• α: a normalization factor• Gdata , GMC : relative PMT gains for data and MC• L: light attenuation length in the water• NPMT : The number of PMTs with signal which passed the cuts of Eq. (6.21)• ri : distance from the vertex to the i-th PMT• Θi : Photon incident angle on the i-th PMT

CHAPTER 6. EVENT RECONSTRUCTION BASED ON THE PRIMARY SIGNAL

74• f (Θi ): Correction factor of the i-th PMT acceptance as a function of Θi• Si : expected charge of the i-th PMT from scattered photons
RT OTn is converted to the charged particle momentum based on the result of the energy scalecalibration in Sec. 4.3.

Chapter 7

Event Selection
The DAQ system in Super-Kamiokande acquires the data from not only the neutrino interactions, but also a large amount of the backgrounds. It registers O(106 ) events with the HE
trigger per day. The definition of the HE trigger is described in Sec. 3.9. On the other hand,the number of the atmospheric neutrino events in this data is only about eight events per day.
To study the atmospheric neutrino properties, the background event rate has to be suppressed
suﬃciently smaller than the neutrino event rate. The acquired events are subject to three kindsof event selection as follows.
Fully Contained (FC) Event Selection
This selection searches for events whose primary interaction occurred in the inner detector
and the generated charged particles stopped in the inner detector.
Partially Contained (PC) Event Selection
This selection searches for events whose primary interaction occurred in the inner detector,
but the generated particles went out of the inner detector and the outer detector detected
their Cherenkov light.
Upward-Going Muon (UPMU) Event Selection
This selection searches for events with muons entering from the bottom side of the tank.
These UPMU events are considered to be generated from the neutrino interaction with therock surrounding the detector. The cosmic ray muons generated on the other side of theearth interact with the earth before their arrival to the detector. Therefore, UPMU events
do not contain cosmic-ray muons. The events are separated again to upward stopping
muon and upward through-going muon based on whether or not the muon stopped in the
inner detector.
Fig. 7.1 shows the schema of the selection target behavior. The events that passed any of theselections are used for the oscillation analysis.

7.1

Fully Contained (FC) Events

7.1.1

First Selection (FC1, Cosmic Ray and Low Energy Events Cut)This step rejects the cosmic ray muon events and the low energy events.

75

CHAPTER 7. EVENT SELECTION

76

Figure 7.1: The schema of the target events of each selection types. The solid arrows and brokenlines show the behaviors of charged particles and neutrinos, respectively.
The low energy events cut removes the solar neutrino events, radiation background and soon. The solar neutrino events are studied with the diﬀerent selection criteria, which is notdiscussed in this thesis. Because the Cherenkov light intensity is an indicator of the neutrinoenergy, the low energy events are rejected based on the maximum total charge of ID PMT in a
300 ns timing window. The candidate events are required to satisfy the criteria below:• PE300 ≥ 200 p.e.s (for SK-I, III, and IV),
• PE300 > 100 p.e.s (for SK-II)
where PE300 is the maximum total charge collected by ID PMTs in a 300 ns timing window
sliding between -500 ns and 500 ns from the mean hit time of the signals in ATM trigger gate.
The cosmic ray muon is produced outside the detector, for example in the atmosphere, andit produces Cherenkov light in the outer detector. Therefore, the cosmic ray muon and FC
neutrino events can be identified by checking if the PMTs in the outer detector have many hits
or not. The candidates which do not satisfy the criteria below are rejected as the cosmic-ray
muons:• NHITA800 ≤ 50 hits (for SK-I, II, III)
• NHITA800 ≤ 55 hits (for SK-IV)
where NHITA800 is the number of OD PMT hits in a timing window from 500 ns to 1300 ns.

7.1.2

Second Selection (FC2, Cosmic-ray Muon and Flasher Events Cut)This step rejects the cosmic-ray muons with lower energy than the muons rejected by the firststep selection. In addition, the backgrounds from flashing PMTs and from PMTs with electricalnoise are also rejected.

CHAPTER 7. EVENT SELECTION

77

As an indicator of the muon energy, the total charge in ID PMT is used for the selection.
When the total charge exceed the criteria, the number of OD PMT hits is checked again withthe stricter criteria than the first selection. The selection criteria are set as follows:• PEtot < 100000 p.e.s(50000 p.e.s for SK-II)• OD trigger is generated• NHITA800 > 25 hits (30 hits for SK-IV)
where PEtot is the total charge collected from ID PMTs in the ATM trigger gate. In thisselection, the threshold of NHITA800 is set lower than that of first selection, and OD trigger,
which is described in Sec. 3.9, is also utilized. The candidates which satisfy all of these criteria
are rejected as the cosmic-ray muons.
Super-Kamiokande PMTs sometimes emit light due to electrical discharges. When a ID
PMT flashed, the ID PMTs close to the flashed PMT can detect its photons, and OD PMT
does not detect the light. Therefore, this event cannot be rejected based on OD PMT signal.
This event is called flasher event. To reject this noise, this step selects the events without the
charge concentration to a single PMT. This selection is also possible to reject the events that
contains a single PMT with large electronics noise. The selection criteria are defined as follows:• Flasher cut– PEmax /PE300 > 0.5
where PEmax is the maximum charge observed by a single ID PMT in the ATM trigger gate.
The candidate events that satisfy the criteria are rejected as the flasher events.

7.1.3

Third Selection (FC3, Cosmic Muon Cut and Low Energy Event Cut)This step rejects the cosmic ray muon events with small number of hits in OD. The selection
criteria are set for each background source as follows. If any of the following criteria are satisfied,the event is rejected.
7.1.3.1

Hard Muon Cut

The high energy cosmic ray muons are identified based on the number of OD PMT hits withthe shorter time window than FC1 and FC2 as follows:• NHITA500 ≥ 40 hits
where NHITA500 is the maximum number of OD PMT hits in the 500 ns time window sliding
until 1300 ns. The candidates which satisfy this criteria are rejected as the cosmic ray muon
event.
7.1.3.2

Through Going Muon Cut

Through going muons enter and go out passing through the OD. Therefore two signal clusters
are generated in the OD. These two clusters enable to reconstruct the muon track assuming
their positions to be the entrance and exit points of the muon. Based on this muon track, thecriteria of the through-going muon identification are set as follows:

CHAPTER 7. EVENT SELECTION

78• PEmax ≥ 231 p.e.s• Number of ID PMTs with hit in ATM trigger gate ≥ 1000 hits
• NHITAin ≥ 10 hits or NHITAout ≥ 10 hits
• Gµ > 0.75
where NHITAin and NHITAout are the number of OD PMTs hits within 8 m from the entrance
and exit points. The hits are collected from the time window whose range is from 500 ns to1300 ns. Gµ is the goodness defined as follows:
))
(
∑( 1
1
(ti − Ti )2
Gµ = ∑ 1 ×
exp −
(7.1)
2 × (1.5 × σi )2
σi2
i σ2
i

i

where ti and σi are the observed hit time of the i-th PMT and its resolution, and Ti is the hit
time expected from the muon track and its entrance time. The candidates which satisfy all of
the criteria are rejected as the cosmic ray muon events.
7.1.3.3

Stopping Muon Cut

Unlike the through-going muon, the stopping muon generates one cluster in the OD around theentrance point. To analyze the relationship between the ID PMT signal distribution and themuon track, the muon direction has to be reconstructed. It is reconstructed by maximizing thetotal charge of ID PMTs inside the Cherenkov cone with 42◦ opening angle, fixing the conicalvertex to the entrance point. Based on the reconstructed direction, the expected hit time ofeach PMT can be estimated. The goodness is defined in the same way as Eq. (7.1). By using
this goodness, the stopping muon identification criteria are defined as follows:• NHITAin ≥ 10 hit (for all periods)or• Gµ ≥ 0.5 and NHITAin ≥ 5 hit (only for SK-I)
The candidates which satisfy these criteria are rejected as the cosmic ray muon events.
7.1.3.4

Cable Hole Muon Cut

As mentioned in Sec. 3.7, The plastic scintillator veto counters are installed on the four cableholes to compensate the insensible area of OD. By using the signals of the scintillator vetocounters which are described in Sec.3.7, the reduction criteria are defined as follows.
• Any scintillator veto counter detects signal• The vertex is reconstructed within 4 m from the position of the scintillator veto counter
with signal
The vertex reconstruction method is the same as that for stopping muon cut in the previous
section. The candidates which satisfy these criteria are rejected as the cosmic ray muon.

CHAPTER 7. EVENT SELECTION
7.1.3.579Flasher Cut

The PMT flashing event tends to generate longer timing distribution than that of neutrino
interaction events. Based on this diﬀerence, the flasher identification criteria are defined asfollows:• NMIN100 ≥ 15 hits (for SK-I)
NMIN100 ≥ 20 hits (for SK-II, III, and IV)
or• NMIN100 ≥10 hits and the number of hit ID PMTs ≤ 800 hits (only for SK-I)
where NMIN100 is the minimum number of hits in 100 ns time window sliding from 1200 ns
to 1700 ns. The candidates which satisfy these criteria are rejected as the flasher event. The
remaining flasher events with low intensity are rejected based on the goodness of event reconstruction. The reduction criteria are defined as follows:• Nhit ≤ 250 and Glow ≤ 0.4
where Nhit is the number of hits in the ATM trigger gate. Glow is the goodness calculated
based on the vertex reconstruction algorithm described in Sec. 4.3.1.2. Flasher events haveworse goodness than neutrino interaction events because of the wider time distribution. The
candidates which satisfy these criteria are rejected as the flasher event.
7.1.3.6

Coincidence Muon Cut

When a cosmic ray muon enters just after low energy events triggered by SLE trigger, it will not
be rejected by OD signals of the delayed muon because the OD hit timing is too late comparing
to the normal event. On the other hand, owing to the muon, this event has a large amountof charge enough to pass the low energy cut. To reject such events, the coincidence event
identification criteria is defined as follows:• NHITAoﬀ ≥ 20
• PEoﬀ ≥ 5000 p.e.s (2500 p.e.s for SK-II)
where NHITAoﬀ is the number of OD PMTs hits in the oﬀ timing window from 1300 ns to 1800
ns, and PEoﬀ is the total charge of ID PMTs collected in the same oﬀ timing window. The
candidates which satisfy these criteria are rejected as the coincidence event.
7.1.3.7

Low Energy Events Cut

This rejection removes low energy events with diﬀerent criteria from the 1st selection. The
criteria are defined as follows:• The number of ID PMT hits in the ATM trigger gate < 500 hits
• NHIT50 < 50 hits (25 hits for SK-II)
where NHIT50 is the maximum number of hits in the 50 ns width time-of-flight subtracted hit
time window sliding in the ATM trigger gate. The candidates which satisfy these criteria arerejected as the low energy events.

CHAPTER 7. EVENT SELECTION

7.1.4

80

Fourth Selection (FC4, Flasher Cut)In this step, the PMT flasher events are rejected. Because the flashing events occur repeatedly
from the same PMTs, the same signal distribution tends to be produced. Therefore, the events
having similar signal patterns to the past events are rejected in this step. To do this rejection,
each event remaining after the third selection is compared the charge distribution with the otherremaining events. After the separation of the inner detector to 1450 patches, the correlation
value r is calculated as follows:
r=

A
B
B
1 ∑ (QAi − ⟨Q ⟩) × (Qi − ⟨Q ⟩)N
σA × σB(7.2)

iwhere the signs of A and B indicate the diﬀerent FC event candidates which passed the thirdB
selection, respectively. N is the number of patches, QA
i , Qi are the total observed charges in theAB
i-th patch, ⟨Q ⟩ and ⟨Q ⟩ are the average charge and σA and σB are their RMS of all patches.
If the distance between the ID PMTs with the maximum charge of the two events is less than
0.75 m, 0.15 is added to r value. The threshold rth is set by using the linear function of thetotal ID PMT charge P Etot . The oﬀset and slope of the function is set by fixing two specificpoints, which are (0.5, 200 p.e.s) and (1.0, 30000 p.e.s) as (rth , P Etot ) values. The function is
calculated as follows:{
rlow = a0 log10 Qlow + a1((rth , P Etot ) = (rlow , Qlow ))(7.3)
rhigh = a0 log10 Qhigh + a1 ((rth , P Etot ) = (rhigh , Qhigh ))⇒ rth =

rlow log10 Qhigh − rhigh log10 Qlow
(P Etot − Qlow ) + rlow
log10 Qhigh − log10 Qlow

(7.4)

where (rlow , Qlow ) and (rhigh , Qhigh ) are the fixed two points, and a0 and a1 is the slope andoﬀset of the linear function, which are calculated by solving Eq. (7.3). When r exceeds rth , Boththe target and comparison events are marked as matched events, and they are rejected.

7.1.5

Fifth Selection (FC5, Remaining Background Cut)This step rejects remaining backgrounds from the cosmic ray muon, flasher, and DAQ systemproblems. The selection criteria are defined for each background source.
7.1.5.1

Invisible Muon

When a low energy muon enters the detector, the delayed signal from its decay electrons is
possible to produce the trigger, instead of the parent muon. In this case, the muon’s Cherenkov
light is detected in OD before the trigger. To reject this event, additional criteria to the early
OD signals are applied. These rejection criteria are defined as follows:• PEtot <1000 p.e.s (500 p.e.s for SK-II)• NHITACearly ≥5 hit (if DISTclust <500 cm)
NHITACearly ≥10 hit (if DISTclust >500 cm)• NHITACearly + NHITAC500 ≥10 hit (if DISTclust <500 cm)

CHAPTER 7. EVENT SELECTION

81

where NHITACearly is the maximum number of hits in the OD hit cluster, which is detected in
the 200 ns time window sliding between -8000 ns and +800 ns, NHITAC500 is the maximum
number of hits in the OD hit cluster between +800 ns and +1300 ns, and DISTclust is the distance
between OD hit clusters with NHITACearly and NHITAC500 hits. The candidates which satisfy
these criteria are rejected as the invisible muon event.
The hit clusters are formed by using the charge gradient between neighborhood patches. To
search for the hit cluster, the ID and OD walls are separated to 21×21 and 11×11 patches withthe same solid angles whose center is the detector center, respectively. Fig. 7.2 shows the schema
more charge

Cluster 1

less charge

Cluster 2

Cluster 3

Figure 7.2: The schema of the cluster detection algorithm [121]. The size of circles represents
the total charge observed in each patch. The arrow shows the vector of the charge gradient
which head for the patches with the maximum charge in the neighborhood patches.
of the cluster setting.
7.1.5.2

Coincidence Muon

In addition to the selection in Section 7.1.3.6, the remaining coincidence muons are rejectedbased on the criteria as follows:• PE500 < 300 p.e.s (150 p.e.s for SK-II)• PElate ≥ 20 p.e.s
where PE500 is the total charge of ID PMTs in the time window from +800 ns to +1300 ns, andPElate is the maximum charge of OD PMTs in the 200 ns time window sliding between +1300
ns and +2500 ns. The candidates which satisfy these criteria are rejected as the coincidence
muon event.
7.1.5.3

Long Tail Flasher

The remaining flasher events are rejected based on the long tail of the hit time distribution andthe goodness of the reconstructed vertex. The rejection criteria are defined as follows:

CHAPTER 7. EVENT SELECTION

82• NMIN100 ≥ 6 and Gp <0.4 (for all period)or• NMIN100 < 6 and Gp <0.3 (for SK-II, III, and IV)
where NMIN100 is the minimum number of ID PMT hits in the 100 ns time window sliding
between +1200 ns and +1700 ns from the trigger, and Gp is the goodness of the point fit, whichis the simple vertex fitter described in Section 6. The candidates which satisfy these criteria arerejected as the flasher event.
7.1.5.4

Electronics Noise

Electric noise may be produced from the high voltage systems or electronics board. Because
this noise has diﬀerent pulse shape from the normal signal, and its charge tends to be smaller
than single photoelectron level, the electric noise identification criteria are defined as follows:• N0 ≥ 250 hits (for SK-I, III, IV)
N0 ≥ 125 hits (for SK-II)• N0 ≥ N1 + 100 hits (for SK-I, III, IV)
N0 ≥ N1 + 50 hits (for SK-II)
where N0 is the number of ID PMT hits with less than 1 p.e.s charge signal, and N1 is thenumber of ID PMT hits with more than 1 p.e.s charge. The candidates which satisfy thesecriteria are rejected as the background event caused by electronics noise.
7.1.5.5

Dead Electronics

The PMT signals are input to the electronics installed in four rooms on the detector. Theserooms are called electronics huts. To remove the events whose any electronics huts are notworking correctly, the selection criteria are defined as follows:• All electronics huts have one or more OD hits.
Because PMTs produce dark noise, all electronics huts should have one or more hits in the event
gate if the electronics works correctly. The candidates which do not satisfy these criteria arerejected as the events with dead electronics.
7.1.5.6

Calibration Related Events

In order to monitor the change of water transparency, PMT gain, and so on, the automatedcalibration with some light sources is executed periodically during the measurement. The calibration events are taken with the trigger for the calibration, which is diﬀerent from the trigger
for the normal measurement. But it is still possible that the normal measurement events arecontaminated with these calibration events. When the auto-calibration is executed, the intensityof the light source is measured by its monitor PMT. Therefore, by checking if the monitor PMT
have any signals or not, this contamination is removed.

CHAPTER 7. EVENT SELECTION
7.1.5.783

Stopping Muon

The additional criteria to identify the remained stopping muon events are defined based on the
OD PMT hits in the clusters as follows:• NHITAC≥ 10 hit (for SK-I)
or• NHITA′in ≥ 10 hit
or• NHITA′in ≥ 5 hit and Gµ >0.5
where the definition of Gµ is the same as Eq. (7.1), NHITA′in is the number of hits in the OD
hit cluster around the entrance point whose hit time is in the time window from +500 ns to1300 ns.
7.1.5.8

Cable Hole Muon

To reject the case where cosmic ray muon enters from the no OD area without any signals inthe scintillators on the cable holes, additional identification criteria of the cable hole muon are
defined as follows:• Gµ ≥ 0.4
and P Etot > 4000 p.e.s
and cos Θ < −0.6
and the distance between the muon entrance point and the nearest hole position<250 cm
or• NHITAin ≥ 4 hit
where Θ is the zenith angle of the produced particle’s direction, NHITAin,APFIT is the numberof hits in the OD PMTs within 8 m from the muon entrance point reconstructed with APFit,
which is described in Sec 6. The reconstruction process of APFit is described in Sec. 6.
7.1.5.9

Wrong OD Timing Hit

Because of the problem of the TDC modules used for OD PMTs until SK-III, irregular hit times
are sometimes registered as the OD PMT signals. In that case, the OD hit cluster is generated
at the early timing, and abnormal number of sub-events is also generated. The definition ofsub-event is described in Sec 8.1.2.1. To cut this eﬀect, the identification criteria are defined asfollows:• PEtot ≥ 1000 p.e.s (for SK-I and III)
PEtot ≥ 500 p.e.s (for SK-II)• NHITACearly ≥ 20 hit (for SK-I, II, and III)• Seven subevents exist(for SK-I, II, and III)

CHAPTER 7. EVENT SELECTION

84

where PEtot is the total charge collected from ID PMTs in the ATM trigger gate, and NHITACearlyis the number of hits in the OD cluster, which are detected in the 200 ns time window slidingfrom -8000 ns to +800 ns. The candidates which satisfy all of these criteria are rejected as theevents with wrong OD timing hit. Because this problem is caused by the problem of the OD
DAQ system used for SK-I, II, and III, this reduction is not applied to SK-IV events.

7.1.6

Fiducial Volume Selection

The final FC data set is obtained by selecting the events that satisfy the criteria below:• dwall>200 cm• NHITAC<16 hits (10 hits for SK-I)• Evis >30 MeV
where dwall is the distance between the detector wall and the reconstructed vertex, NHITAC
is the number of hits in the OD hit cluster with the largest charge, and Evis is the sum of thereconstructed energies of all the Cherenkov rings, assuming the rings come from electrons. The
vertex and energy are reconstructed as written in Sec. 6.

7.1.7

FC Data Summary
Event Rate [/day]
Before Selection
FC1
FC2
FC3
FC4
FC5
FCFV

∼ O(106 )∼ O(103 )∼ O(102 )∼ O(10)∼ O(10)∼ O(10)∼8

Remaining neutrino event rate [%]SK-I SK-II SK-III
SK-IV
100.0 100.0 100.0
100.0
100.0 99.96 100.0
100.0
99.97 99.92 99.98
99.99
99.67 99.78 99.81
99.83
99.02 99.36 98.52
98.87
98.85 99.29 98.46
98.83
98.18 99.31 98.07
98.07

Table 7.1: The eﬃciencies of neutrino events estimated with MC atmospheric neutrino events.
To evaluate the event rate, the distance between the detector wall and true and reconstructed vertices are defined as dwallv and dwall, respectively. From FC1 to FC5, the estimation is done byusing only the neutrino events that satisfy dwallv>200 cm, Evis > 30 MeV and NHITAC<16 hits.
FCFV means the events which pass all FC selections and fiducial volume selection.
Table 7.1 shows the event rate and the neutrino event selection eﬃciencies at each step of the
reduction, which is estimated based on MC output. As shown in this table, these selections reject
the most triggered events, but they keep more than 98% of all the neutrino events. After this selection, the event rate of the samples which pass all the selections becomes ∼8 events/day. From
SK-II, the remaining background events are searched by comparing with the signal distributionsof the typical background events, for example flasher events, based on eye-scan. Table 7.2 is thetable of the background rate of single ring events based on the eye-scan result. As shown in thistable, it is confirmed that the number of the background events becomes much fewer than theneutrino events.

CHAPTER 7. EVENT SELECTION
SK-I
1-ring Sub-GeV e-like
1-ring Sub-GeV µ-like
1-ring Multi-GeV e-like
1-ring Multi-GeV µ-like

Not
Scanned85
SK-II
0.11% ± 0.08%
0.30%±0.13%<0.30%<0.30%

SK-III
0.08% ± 0.08%
0.32%±0.16%
0.36%±0.36%<0.42%

SK-IV<0.02%<0.02%<0.08%<0.08%

Table 7.2: The background contamination rate table of the single ring event types. This rateis based on the result of the eye-scan. Sub-GeV and Multi-GeV events satisfy Evis <1.33 GeV
and Evis ≥1.33 GeV, respectively.

7.2

Partially Contained Events

As mentioned in Sec. 3.7, OD segmentations were installed in the top and bottom regions
between SK-II and SK-III periods. Therefore, the selection criteria of the partially contained
events are largely diﬀerent between SK-I and SK-II events and SK-III and SK-IV events. The
selection is separated to five steps and fiducial volume cut. In this section, the selection criteria
are described.

7.2.1

First Selection (PC1, Through Going Muon and Low Energy Events
Cut)

This step rejects the through going muon and low energy events. Basically, these backgrounds
are rejected by using the number of OD hit clusters and the intensity of the Cherenkov light.
• Selection Criteria for SK-I and II
Because SK-I and SK-II do not have the OD segmentations, the OD signals around theentrance and exit points of the through going muon are possible to be taken as a single OD
hit cluster when the points are near. To reject such events, the OD hit time distributionis also used as a criterion. If the hit time distribution is broad, the event is rejected as thethrough going muon event. The selection criteria are defined as follows:– PEtot ≥ 1000 p.e.s (500 for SK-II)– TWIDA≤ 260 ns (170 ns for SK-II)– NCLSTA ≤ 1 (SK-I only)
where PEtot is the total charge collected from the ID PMTs in the ATM trigger gate,
TWIDA is the width of the time region that satisfies that the number of OD hits in
a 100 ns sliding time window have more than 14 hits, and NCLSTA is the number of
OD hit clusters. The candidates which do not satisfy these criteria are rejected as thethrough-going muon events or low energy events.
• Selection Criteria for SK-III and SK-IV
OD segmentation enables to separate the OD hits produced at the entrance and exit of thecosmic-ray muon more easily. To reject the through-going muon, some criteria about theOD hits in each location are introduced. The candidates which do not satisfy the criteria
below are rejected as the through-going muon events or low energy events.

CHAPTER 7. EVENT SELECTION

86– SK-III∗ PEtot ≥ 1000 p.e.s∗ NHITAtop <10 hits or NHITAbottom <10 hits
∗ NHITAendcap <25 hits or NHITAside <70 hits
∗ ODRmean <2100 cm if NHITA′500 <20 hits
– SK-IV∗ PEtot ≥ 1000 p.e.s∗ NHITAtop <11 hits or NHITAbottom <11 hits
∗ NHITAendcap <29 hits or NHITAside <84 hits
∗ ODRmean <2100 cm if NHITA′500 <22 hits
where NHITAlocation is the number of the OD PMT hits in each location. ”Endcap”
corresponds to the top and bottom regions. NHITA′500 is the maximum number of OD
PMT hits in the 500 ns time window sliding in the whole region of the trigger gate.
ODRmean is the mean of the distance between the all patterns pairs of the OD PMTs withhits, which is defined by:
ODRmean =1
Npair

NHITA-1∑∑ NHITAi=1|⃗xi − ⃗xj |

(7.5)

j=i

where ⃗xi is the position of the i-th OD PMT.
The performances of the detector, for example the dark noise rate and gain of the PMTs, arechanged as time goes on, and the PC event rate has been increasing. In the data measured in2016, ∼30000 events/day are remaining at this step.

7.2.2

Second Selection (PC2, Cosmic Ray Muon Cut)This step rejects the cosmic-ray muons by detecting hit clusters at the entrance point to thedetector. This selection uses the number of clusters in OD. The method how to search for thecluster is described in Sec. 7.1.5. In this selection, the ID and OD walls are separated to 21×21
and 11×11 patches with the same solid angles whose center is the detector center, respectively.
The hit clusters are formed by using the charge gradient between neighborhood patches. The
candidates which satisfies the selection criteria below are selected.
• SK-I– There is at most one OD hit cluster with more than 6 hits.
– NHITACmin <7 hits
– PE200 >1000 p.e.s (if OD hit cluster with more than 6 hits exists)• SK-II– There is at most one OD hit cluster with more than 6 hits.
– NHITAendcap < min[20 hits,MAX(NHITAside )]– NHITAC<12 hits + PE200 ×0.085 (hits/p.e.s)CHAPTER 7. EVENT SELECTION

87• SK-III, and SK-IV– There is at most one OD hit cluster with more than 10 hits.
– NHITAendcap < min[20 hits,MAX(NHITAside )]where NHITAlocation is the number of OD PMTs with hit in the location, NHITACmin is thesmaller number of hits among the endcap and barrel regions in the largest OD hit cluster. PE200
is the charge detected within 200 cm from the PMT with the largest charge in the ID hit cluster
nearest to the OD hit cluster with the largest charge, NHITAC is the number of hits in the OD
hit cluster with the largest charge, and MAX(NHITAside ) follows:
MAX(NHITAside ) = exp(a − b × NHITAside ) hits

(7.6)

(a, b)(NHITAside < 75 hits)(7.7)(a, b) = (4.675, 0.008) (NHITAside ≥ 75 hits)

(7.8)

where
= (5.8, 0.023)

In the data measured in 2016, ∼10000 events/day remained at this step.

7.2.3

Third Selection (PC3, PMT Flasher Event Cut)In this step, we aim to reject the ID PMT flasher events for all periods. To do this work, the longtail of its hit time distribution, which is expected in the flasher events, is searched. For the SK-I
and SK-II period, the stopping muon events are also rejected by searching for their entrance
point. The candidates which satisfy the criteria below are rejected as background events:• NMIN100 > 14 hits
or NMIN100 > 9 hits (the number of hits in the ID PMTs<800 hits)• NHITAin >10 hits (Only for SK-I and SK-II)
where NMIN100 is the minimum number of ID PMTs with hits in a 100 ns timing windows
sliding from 1200 ns to 1700 ns, and NHITAin is the number of OD PMTs which are locatedwithin 8 m from the reconstructed vertex and have hits in the time window from 800 ns to1300 ns. The event rate does not changed significantly from the second step.

7.2.4

Fourth Selection (PC4, Remaining Cosmic Ray Muon Event Cut)In this step, we aim to reject the remaining cosmic ray muon events with some criteria. Whenthrough-going muon enters to the ID, the earliest signal tends to be detected by the PMT aroundthe muon entrance point. On the other hand, the earliest hit of the partially contained event willbe the exit point of the particle, because the muon has faster speed than the Cherenkov light in
the water. Therefore, by comparing the particle direction and the earliest hit position, we canreject the through-going muon events. As the other through-going muon rejection, the track
length is also checked. The event whose track length is very long is rejected as the throughgoing muon events with high energy. In addition, this step also rejects the corner clipping muon
events. These are the events whose muon skims the ID. To reject the corner clipping muon,
the events whose vertices are near any corners of the ID are rejected. As a final criterion, thegoodness value evaluated by the event reconstruction tool is used. The detailed criteria before
and after the OD segmentation introduction are described, respectively.

CHAPTER 7. EVENT SELECTION
7.2.4.1

88

SK-I and SK-II

The vertex and particle direction are reconstructed by Point-fit, which is described in Sec. 6.1.1.
The candidates which satisfy the criteria below are selected:• d⃗point-fit · d⃗PMT > −0.8
where d⃗point-fit is the particle direction reconstructed with Point-fit, d⃗PMT is the direction
from the reconstructed vertex to the ID PMT around the particle entrance point. This
PMT is selected as follows:
1. Search for the PMTs which have maximum number of hits in the eight of the nearest
PMTs.
2. If plural PMTs satisfy the criteria 1, select the PMT with the earliest hit.
Fig. 7.3 shows the schema of the first process of the PMT selection.

Figure 7.3: The schema of the candidate selection. The yellow circles are assigned to thecandidate PMTs selected by the first process of the PMT selection, and the pink circles are
assigned to the PMTs with hits. The number of hits in the candidate itself and eight of thenearest PMTs, which the red and blue boxes are assigned to, are required to be the largest in
the detector.
• DCORN<150 cm
where DCORN is the distance between the reconstructed vertex and the nearest corner ofthe inner detector.

CHAPTER 7. EVENT SELECTION

89• TLMU<30 m if Gµ >0.85
where TLMU is the estimated track length assuming it is a through-going muon, and Gµ is
the goodness of the through-going muon fit, which follows the same definition as Eq. (7.1).
7.2.4.2

SK-III and SK-IV

After the OD segmentation introduction, the reduction is done based on the information of themuon fitter ”muboy”[122]. This fitter can classify the events as either through-going muon,
stopping muon, multiple muons or corner clipper muon. When the event is classified as stoppingmuon, additional criteria are applied. The neutrino event selection criteria are defined as follows:
1. angleµ < 90◦
where angleµ is the opening angle between the reconstructed direction of the particle andthe direction from the reconstructed vertex to the center of the OD hit cluster, which arereconstructed by muboy.
2. d⃗muboy · d⃗PMT > −0.8
where d⃗muboy is the muon direction reconstructed by the muboy fitting, and d⃗PMT is thesame definition as that of SK-I and SK-II.
3. lengthµ <17.5 m
where lengthµ is the muon track length reconstructed by the muboy fitting.
4. goodnessµ <0.52
where goodnessµ is the goodness returned by the muboy fitting. This goodness becomes
higher when the event is cosmic ray muon.
5. cornerµ >300 cm
where cornerµ is the length between the nearest inner detector corner and the entrance
point of the muon reconstructed by muboy.
To select the neutrino event, it is required for the events which muboy classifies as throughgoing or stopping muon events to satisfy four or more criteria above. In addition, for the events
identified as stopping muon, the 2nd criteria are required to exist in the satisfied criteria, andthe events are also required to satisfy the additional criteria below:• goodnessµ <0.5 or ehitµ <10 hits
where ehitµ is the number of hits in the time window from 800 ns to 1300 ns detected by
OD PMTs within 8 m from the entrance point. The entrance point reconstruction is doneby using muboy fitter.
On the other hand, the corner clipping and multiple muon events are required to satisfy only
two criteria to be passed. After this selection, the additional criteria are defined to reject low
energy background event:• PEtot <3000 p.e.s
where PEtot is the total charge of the ID PMTs detected within the ATM trigger gate. In thedata measured on 2016, ∼500 events/day remain at this step.

CHAPTER 7. EVENT SELECTION

7.2.5

90

Fifth Selection (PC5, Remaining Background Event Cut)This step rejects the remaining backgrounds. The criteria are diﬀerent between events in SK-I
and SK-II periods and SK-III and SK-IV periods. Table 7.3 shows the table of the cuts to be
applied. For SK-I and SK-II events, this step requires to satisfy all the applied criteria. On theNo.

Selections

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

Low Energy Event Cut
Through-going Muon Cut (1st)Through-going Muon Cut (2nd)Through-going Muon Cut (3rd)Through-going Muon Cut (4th)
Stopping Muon Cut (1st)
Stopping Muon Cut (2nd)
Stopping Muon Cut (3rd)
Stopping Muon Cut (4th)
Stopping Muon Cut (5th)Cable Hole Muon Cut
Corner Clipping Muon Cut (1st)Corner Clipping Muon Cut (2nd)Corner Clipping Muon Cut (3rd)
Decay Electron Cut

Periods
SK-I and SK-II SK-III and SK-IV
⃝
⃝
⃝
⃝
△
⃝
△
⃝
⃝
△
⃝
⃝
⃝
△
△
△
⃝
⃝
⃝
△
⃝
⃝
△

Table 7.3: The table of the selection criteria applied to SK-I and SK-II period and SK-III and
SK-IV period. ⃝ and △ indicate hard and soft cut criteria, respectively.
other hand, SK-III and SK-IV have two kinds of the cuts, called hard cut and soft cut. Whencandidates satisfy any hard cut criteria, it is rejected as background events. On the other hand,
this step allows the candidate to satisfy at most one soft cut criterion. The selection criteria aredescribed as follows:
1. Low Energy Event Cut
This cut rejects the low energy background events that satisfy the following criteria:• PEtot < 3000 p.e.s
where PEtot is the total charge within the ATM trigger gate detected by ID PMTs. Thiscut is applied only to SK-I and SK-II samples as a hard cut. For SK-III and SK-IV, thesame criterion is applied at the fourth reduction.
2. Through-going Muon Cut (1st)This selection rejects the through-going muons entering from top or bottom and exitingfrom the other part of the endcap. The candidates which satisfy the criteria below arerejected as through-going muon:• NHITAtop ≥ 7 hits and NHITAbottom ≥ 7 hits

CHAPTER 7. EVENT SELECTION

91• PEAtop ≥ 10 p.e.s and PEAbottom ≥ 10 p.e.s• 0.75 < (TDIFFA/ns)×c/40 m <1.5
where NHITAlocation and PEAlocation is the number of hits and total charge of OD PMTs
within 8 m from the entrance or exit point on the location, and TDIFFA is the diﬀerenceof the average hit time of OD PMTs in the two regions within 8 m from the entrance andexit point on the location. This cut is applied for all periods as a hard cut criterion.
3. Through-going Muon Cut (2nd)This cut rejects the candidates based on similar criteria to the 1st through-going muoncut, but the entrance and exit point, and the track length of the muon are reconstructedwith MS-fit, which is described in Sec. 6.4. The candidates which satisfy the criteria arerejected as through-going muons.
• NHITAin ≥ 5 hits and NHITAout ≥ 5 hits
• 0.75 < TRACK/(TDIFFA ×c) <1.5
where NHITAlocation is the number of hits of OD PMTs within 8 m from the location,
and TRACK is the distance between the entrance and exit point of the muon which isreconstructed by MS-fit. For SK-I and SK-II, these criteria are hard cut criteria, but they
are soft cut for SK-III and SK-IV.
4. Through-going Muon Cut (3rd)This cut rejects the events with two large clusters, which correspond to the entrance andexit points of the cosmic-ray muon. The criteria are defined as follows:• NCLSTA5(1)≥ 10 hits and NCLSTA5(2) ≥ 17 hits
where NCLSTA5(x) is the number of OD PMTs with hits in the cluster with the x-th
highest charge. For SK-I and SK-II, these criteria are hard cut criteria, but they are soft
cut for SK-III and SK-IV.
5. Through-going Muon Cut (4th) This step rejects the through-going muon events
by counting the number of OD hit clusters when the particle energy is high. This step
searches for the OD hit clusters by using a diﬀerent algorithm from Sec. 7.2.2, which iscalled ”nncluster” algorithm. This algorithm defines an OD hit cluster as the area where
includes one seed OD PMT with high charge and 80 of the nearest OD PMTs to theseed PMT. When the number of hits in the 80 PMTs exceeds a threshold, the group isregarded as the OD hit cluster. The candidates which satisfy the criteria below are rejected
as through-going muons.
• Ncluster > 1 if Evis > 20 GeV
where Ncluster is the number of the OD hit clusters detected by nncluster algorithm. These
criteria are applied only to SK-III and SK-IV as hard cut criteria.
6. Stopping Muon Cut (1st)This cut rejects the stopping muon by using the signal of the OD cluster produced atits entrance point. The candidates which satisfy the criteria below are rejected The cut
criteria are defined as follows:

CHAPTER 7. EVENT SELECTION

92• NHITAin ≥ 10 hits
where NHITAin is the number of hits in the OD PMTs within 8 m from the entrance point.
The entrance point is reconstructed based on MS-fit.
7. Stopping Muon Cut (2nd)This cut checks the opening angle between the particle direction and the direction from thevertex to OD hit cluster. If the detected particle is the stopping muon, the OD hit cluster
corresponds to the signals produced around the muon entrance point. On the other hand,
the OD hit cluster of the partially contained neutrino events corresponds to the signals
produced around the particle exit point. Considering this eﬀect, the cut criteria are defined
as follows:• θTDC > 90◦ or θMS > 90◦
where θTDC and θMS are the opening angles between the particle direction and the direction
from the vertex to OD hit cluster which are calculated by using TDC-fit and MS-fit,
respectively. These fitters are described in Sec. 6.1.3 and Sec. 6.4 This cut is applied forall periods as a hard cut criterion.
8. Stopping Muon Cut (3rd)This cut uses the stopping muon fit. This fit assumes the position of the earliest hit
cluster in the ID to the particle entrance point. Because of the diﬀerence of the vertex
between stopping muon and partially contained neutrino events, neutrino events tend to
be reconstructed with lower goodness than stopping muon events. The cut criteria aredefined as follows.
• Stopping muon fit goodness > 0• PEcone /PEtot ≥ 0.6• NHITAin > 6 hits
PEcone is the total charge of PMTs within a 42◦ Cherenkov opening angle, which is reconstructed by the stopping muon fitter. NHITAin is the number of hits in the OD PMTs
within 8 m from the entrance point. For SK-I and SK-II, these criteria are hard cut
criteria, but they are soft cut for SK-III and SK-IV.
9. Stopping Muon Cut (4th)As shown in the third step reduction for SK-I and SK-II events, the stopping muon events
can be rejected by using the number of OD hits around the muon entrance point. As asoft cut for SK-III and SK-IV, the cut criteria are defined as follows:• NHITAin > 6 hits (if |P⃗point-fit − P⃗MS | < 15 m)
where NHITAin is the number of hits in the OD PMTs within 8 m from the entrance
point, and P⃗fitter is the vertex reconstructed by the fitter. These criteria are applied only
to SK-III and SK-IV as soft cut criteria.
10. Stopping Muon Cut (5th)This cut is the same as one of the fourth step PC selection criteria for SK-III and SK-IV.
It is also applied in PC5 as a soft cut for SK-III and SK-IV. The cut criteria is as follows:

CHAPTER 7. EVENT SELECTION

93• angleµ < 90◦
where angleµ is the Cherenkov opening angle reconstructed by muboy.
11. Cable Hole Muon Cut
The cable hole muon is identified by using the opening angle between direction of theCherenkov ring and the direction from the vertex to the plastic scintillator veto counteron the cable holes. The criteria are defined as follows:• a veto counter is hit• d⃗TDC · d⃗veto-vertex > −0.8
where d⃗TDC is the Cherenkov ring direction reconstructed by TDC-fit, and d⃗veto-vertex isthe direction from the veto counter with hit to the reconstructed vertex.
12. Corner Clipping Muon Cut (1st)This step rejects corner clipping muon events by using the track length. When the particle
speed is near to the light speed, the intensity of the Cherenkov light is nearly constant. And
the relative muons lose ∼ 2 MeV/cm at least in the water [123]. Therefore, by comparing
the path length and the reconstructed energy, the consistency is checked as follows:• TRACK>15 m• Evis /(2 [MeV/cm]) < TRACK − 15 m
where TRACK is the track length from the reconstructed vertex to the exit point, whichis reconstructed with MS-fit.
13. Corner Clipping Muon Cut (2nd)This cut uses the same criteria as a part of the fourth reduction for SK-I and SK-II. At
this step, it is only applied to SK-III and SK-IV.
• DCORN> 1.5 m
where DCORN is the distance between the reconstructed vertex and the nearest corner ofthe inner detector.
14. Corner Clipping Muon Cut (3rd)
To cut the corner clipping muons with low energy, the OD hit clusters in the segmented
regions are searched separately. By searching for the positions of the hit clusters on the
OD endcap and barrel walls, and the ID wall with the largest charge in each the region,and by checking the distance and the number of hits in the clusters, corner clipping muonis searched.
• Evis < 1 GeV• DISTsegclust < 15 m
• NHITAsegclust > 30 hits
where DISTsegclust is the maximum distance in the all possible pairs selected from the threeclusters with the maximum charges in the ID wall and OD endcap and barrel walls, and
NHITAsegclust is the number of hits in the OD hit cluster with the smaller charge of theselected pair for DISTsegclust evaluation.

CHAPTER 7. EVENT SELECTION

94

15. Decay Electron Cut
When the high energy neutrino interacts, many pions tends to be produced around thevertex. Decay electrons are produced by their decay. Therefore, the cosmic-ray muons aresearched by checking if the decay electron exists or not. The criteria are defined as follows:• Evis >25 GeV• no decay electrons are detected.
where Evis is the total reconstructed energy which assumes the Cherenkov rings come fromelectrons.
After this step, about 1 events/day remains at this step.

7.2.6

Fiducial Volume Selection

As well as the FC selection, the fiducial volume cut is applied to PC selected events to rejectradiation background from the detector wall.
1. dwall>200 cm
2. NHITAC≥16 hits (10 hits for SK-I)
where dwall is the distance between the detector wall and the reconstructed vertex, and NHITAC
is the number of hits in the OD hit cluster with the largest charge. The second criterion is set
to reject the FC samples, which are required to satisfy NHITAC<16 hits as shown in Sec. 7.1.6.

7.2.7

PC Selection Summary
Selection Step
Before Selection
PC1
PC2
PC3
PC4
PC5
PC5 with FV cut

Event rate [day−1 ]O(106 )
O(104 )O(103 )O(103 ) (Almost the same as PC2)O(102 ) ∼ O(103 )∼ O(1)∼ O(0.1)

Table 7.4: The remaining event rate at each selection step.

The remaining event rate at each PC selection step is summarized in Fig. 7.4. The signal
eﬃciencies, background contamination rates, and PC5 event rates before and after the fiducial
volume (FV) cut are summarized in Table 7.5. This signal eﬃciencies are calculated based on
MC output, and the event rates after PC5 selection are corrected by the eﬃciency. The signal
eﬃciency is significantly improved by the installation of OD segmentation just before the startof SK-III, which is described in Sec. 3.7.

CHAPTER 7. EVENT SELECTION

Signal Eﬃciency in MC
Background rate
PC5 Event Rate (/day)(Corrected by the signal eﬃciency)
Event rate after the FV cut (/day)(Corrected by the signal eﬃciency)

95
SK-I
86%
0.2%
1.21±0.03

SK-II
77%
0.7%
1.16±0.04

SK-III95%
1.8%
1.20±0.05

SK-IV
96%
0.6%
1.22±0.02

0.69±0.02

0.68±0.03

0.70±0.04

0.65±0.02

Table 7.5: Signal eﬃciencies, background rates and event rates after the PC selection. The
background contamination rates are estimated based on the eye-scanning result[124].

7.3

Upward-Going Muon Events

This upward-going muon event selection is separated to four steps. After this selection, theexpected number of background events is subtracted from the remaining events. In this section,
the selection criteria are described briefly. The detailed description is given in [125].

7.3.1

First Selection (UPMU1, Too High or Low Energy Event Cut)This step rejects too high or low energy events based on the total charge in the ID. The
downward-going muons are also rejected by using muboy fitter outputs. The selection criteria are defined as follows:• PEtot > 8000 p.e.s (3000 p.e.s for SK-II)• PEtot < 1.75 × 106 p.e.s (8 × 105 p.e.s for SK-II)• goodnessµ > 0.3• cos Θµ < 0 (upward-going)• muboy succeeds to categorize the event type.
where Θµ is the zenith angle reconstructed by muboy fit, and goodnessµ is its goodness value.
As mentioned in Sec. 7.2.4.2, muboy fitter categorizes the applied events to some events types,
for example stopping muon and through-going muon. However, it is possible to fail to categorize
the event. The events that do not satisfy these criteria are removed.

7.3.2

Second Selection (UPMU2, Downward-Going Muon Event Cut)

To separate the upward-going and downward-going muon more precisely, the events are reconstructed and rejected based on muboy and some fitters which are specialized for stopping muonand through-going muon event reconstruction. The applied fitters are diﬀerent between stopping
and through-going muon events. The applied fitters are described in [126]. When any applied
fitters judge the event to downward-going muon or fail to fit, this event is rejected.

CHAPTER 7. EVENT SELECTION

7.3.3

96

Third Selection (UPMU3, Event Type Identification)

Before the third reduction, the events are reconstructed by using the fitter called precise-fitter.
This fitter consists of three fitters, which are called MS-fit and TDC-fit, and OD-fit. MS-fit andTDC-fit are described in Sec. 6.1.3 and Sec. 6.4. These fitters are adapted for upward-going
muon reconstruction, assuming the events have only one muon which comes from the outside ofthe detector. OD-fit is the reconstruction tool based on OD PMT hits. The entrance and exit
points of muon and its direction are reconstructed based on the OD cluster information. The
detailed description of these fitters are written in [127]. The most appropriate parameters forthese fitters are selected as output of precise-fitter. Based on its output, the events are separated
to the background, stopping muon, and through-going muon events. The selection criteria ofthrough going muons are defined as follows:• UMEHIT8M>10 hits (16 hits for SK-II)• UMOHIT8M≥10 hits (16 hits for SK-II)• Track length (Reconstructed by precise-fit) ≥ 700 cm• cos Θ (Reconstructed by precise-fit) < 0
On the other hand, the criteria of stopping muon events are defined as follows:• UMEHIT8M>10 hits (16 for SK-II)• UMOHIT8M<10 hits (16 for SK-II)• Fit momentum (Reconstructed by precise-fit) ≥ 1.6 GeV/c• cos Θ (Reconstructed by precise-fit) < 0
where UMEHIT8M and UMOHIT8M are the number of OD PMTs with hits within 8 m fromthe entrance and exit point of the muon, and Θ is the zenith angle of muon direction. Whenthe event does not satisfy either criteria of stopping muon or through-going muon, it is rejected
as background events.

7.3.4

Fourth Selection (UPMU4, Corner Clipping Muon Event Cut)

As a final step, corner clipping muon events are rejected by counting the number of OD hit
clusters. When an event has two or more OD hit cluster, the directions reconstructed withmuboy and with precise-fit are compared. If they do not match well, the goodness of muboyreconstruction is checked. The selection criteria of the corner-clipping muon is shown as follows:• Ncluster > 1• cos Θ > 0.035• The opening angle between the particle directions reconstructed by muboy and precise-fit
is more than 5◦or
The muboy’s goodness of fit is less than 0.35.
where Ncluster is the number of OD hit clusters reconstructed with nncluster, which is describedin 4th through-going muon cut of Sec. 7.2.5 and Θ is the zenith angle of the particle direction.
When the event satisfies all of these criteria, it is rejected as a corner-clipping muon event.

CHAPTER 7. EVENT SELECTION

7.3.5

97

Upward-Going Muon Event Type Criteria

The remaining Upward-going muon events are separated to stopping muon and through-going
muon events. The through-going muon is possible to lose its energy via an electromagnetic
shower when it has very large energy (≳ 1 TeV). To separate the upward-going muon events withand without the electromagnetic shower, the showering algorithm was developed. This algorithm
compares expected charges of non-showering muons and the observed charges, separating themuon track to the segments with 50 cm length. As a first step, the parameters for the separation
are defined as follows:∑ −2 ( (Qicorr −⟨Qcorr ⟩)2 (⟨Qcorr ⟩−⟨Qexp ⟩)2 )+χ2 = N
(7.9)
i=3
σ2
σ2
Qicorr⟨Qiexp ⟩∆ = ⟨Qcorr ⟩ − ⟨Qexp ⟩ − δQ

(7.10)

where Qicorr and ⟨Qcorr ⟩ are the sum of the observed charges in the Cherenkov cone from ith segment and the average of all segments, ⟨Qexp ⟩ is the expected average charge from eachsegment assuming the event to be a non-showering event. δQ is a tuning parameter to set thepeak value in the ∆ distribution produced by the MC to zero. The observed charge is correctedconsidering the eﬀect of the light attenuation in water and the PMTs’ acceptance for the incident
Cherenkov lights. If ∆ and χ2 satisfy the criteria below, the event is identified to showeringmuon.
(χ2 > 50.0)(7.11)∆ > 4.5 − 0.08χ2 (χ ≤ 50.0)(7.12)∆ > 0.52

7.3.6

Background Subtraction

Because of the non-uniformity of rock thickness, the quantity of cosmic-ray muon events depends
on the detected muon direction. Fig. 7.4 shows the scatter plots of zenith angle versus azimuthal

Figure 7.4: The scatter plots of azimuthal and zenith angles for stopping muon (left plot) andnon-showering (middle plot) and showering (right plot) through-going muon events based on
the data. The particles with θ = 0 head for vertically upper direction. The distribution of
downward and horizontal cosmic-ray muons seem to have some clusters.
angle, which contain not only upward and horizontal muon, but also downward muon events.
The through-going muon events are separated to showering and non-showering events to identify

CHAPTER 7. EVENT SELECTION

98

whether the detected muons mainly lose its energy via an electromagnetic shower or ionization.
As shown in these figures, the background distributions have azimuthal angle dependence. To
calculate the amount of background events, the muon events are separated into two groups basedon their azimuthal angles as follows:• Region 1: ϕ < 60◦ and ϕ > 310◦
• Region 2: 60◦ < ϕ < 310◦
where ϕ is the azimuthal angle of the detected muon, which is defined as the anti-clockwise angle
between a certain direction and the muon direction. ϕ = +220.583◦ is assigned to the true north.
Region 1 corresponds to the region with few cosmic-ray muons, and region 2 has much more
downward and horizontal going cosmic-ray muons than region 1. To remove these background
events from oscillation analysis, the zenith angle distribution of muon events in region 2 iscompared with the distribution in region 1. Fig. 7.5 shows the zenith angle distributions inregion 1 and 2, which are normalized by the run time and the coverage of azimuthal angle. The

Figure 7.5: The zenith angle distributions of stopping muon and showering and non-showering
through-going muon events in region 1 and 2. The particles with θ = 0 head for vertically upperdirection. The black and white points show the distributions in region 1 and region 2. The solid
curves show the region 2’s fitted functions.
zenith angle distributions in region 2 are fitted as follows:
F (cos θ) = P0 + eP1 +P2 cos θ(7.13)

where θ is the zenith angle, Pn is the n-th parameter of fitting. In this fitting, P0 value is fixed
to the averaged rate in region 1, assuming the exponential part of the fitted curve corresponds
to the diﬀerence of the distribution normalized by the width of azimuthal angle range between
region 2 and region 1. Therefore, the background rates are estimated as follows:∫ cos θ2NBG =eP1 +P2 cos θ d(cos θ)(7.14)cos θ1

where θ1 and θ2 are the lower and upper limits of this subtraction. The particles with θ1,2 = 0
head for vertically upper direction. This subtraction range is defined as follows:• Through-going muons: −0.1 < cos θ < 0.0CHAPTER 7. EVENT SELECTION

99• Stopping muons: −0.2 < cos θ < 0.0
For the oscillation analysis, the number of UPMU events is evaluated with the subtraction of
the expected number of background events.

7.4

Selection Summary

Table 7.6 shows the remaining event rate at each step, and Table 7.7 shows the number of
selected UPMU events and estimated background events to be subtracted.
Selection Step
Before Selection
UPMU1
UPMU2
UPMU3
UPMU4

Event rate [day−1 ]O(106 )
O(104 )
O(102 )O(1)O(1)

Table 7.6: The remaining event rate at each of the upward-going muon selection steps.

Stopping Muon
Through-going Muon
Total
(Before BG subtraction)
Background (Stopping)
Background (Through-going)

eventsevents/day
eventsevents/day
eventsevents/day
eventsevents

SK-I

SK-II

SK-III460
0.28±0.01
1863
1.13±0.032323
1.41±0.03
29.6
26.8

221
0.27±0.02848
1.02±0.041069
1.29±0.04
14.0
33.1

206
0.32±0.02728
1.14±0.04934
1.47±0.05
15.2
10.8

SK-IV
(Until 2017 Sep)676
0.27±0.013091
1.16±0.023809
1.44±0.02
56.8
28.0

Table 7.7: The number of remaining events after the upward-going muon selections, and thenumber of background events to be subtracted.

Chapter 8

Energy Correction Based on
Secondary Interactions

ParEcles/from/hadron/interacEons
Ehadron//(Missing/Energy)

Neutrino

Neutron
Correct/Missing/EnergyEvis,/IDcharged/lepton
protonEν,corr

Schema/of/the/reconstructed/energy/distribuEons/
w//MSW/Eﬀect/
w/o/MSW/Eﬀect

New/
MethodCurrent//
Method

90

90

80

80

70

70

60

60

50

50

40

40

30

30

20

20
10

10
0

1

Evis,/ID
MSW/eﬀect/is/shown
2

3

4

5

6

7

8

9

0

1

Eν,/corr
MSW/eﬀect/is/shown
2

3

4

5

6

7

8

9

Figure 8.1: Schema of the energy reconstruction for ν + p → (lepton) + n based on the current
and new method and their energy spectrum with and without MSW eﬀect. New method corrects
the missing energy taken by hadrons, and it is expected to improve the energy resolution. It is
expected that the width of the region where the flux enhancement from the MSW eﬀect becomesnarrower, as the energy resolution becomes better due to the improvements in energy correction
method.
100

CHAPTER 8. ENERGY CORRECTION BASED ON SECONDARY INTERACTIONS 101
Fig. 8.1 shows the schema of energy reconstruction based on the current and new reconstruction method. Super-Kamiokande detects the Cherenkov light of particles generated by theneutrino interaction. This interaction possibly generates some hadrons. After the generation,
they interact with the water via the strong interaction during their transmission in the tank.
Because their momenta are reconstructed based on the result of the energy scale calibrations
with leptons as shown in Sec. 4.3 and Sec. 6.5, where strong interactions are not taken into
account, hadrons’ reconstructed momenta tend to be smaller than their true energy. Also, thechargeless particles, like neutrons, do not emit Cherenkov light, and their energy cannot be reconstructed directly. This eﬀect is shown as the diﬀerence between the neutrino energy and thesum of the generated particles’ reconstructed momenta. This diﬀerence is called missing energy.
The neutrino interaction produces non-negligible missing energy from the produced hadrons.
For example, single ring electron like events with more than 1 GeV primary ring momentum
have more than 500 MeV missing energy, and its value depends on the primary Cherenkov ring
momentum as shown in Fig. 8.2. In oscillation analysis, the mass hierarchy is searched by com-

2500

3000

Eν ,true-Evis,ID [MeV]Eν ,true-Evis,ID [MeV]

3000

2500

Single Ring Electron Like

2000
1500

1500

1000

1000

500

500

0

0

-500

-500

-1000

-1000

104 P [MeV]pri

104 P [MeV]
pri

3000

Eν ,true-Evis,ID [MeV]Eν ,true-Evis,ID [MeV]

3000
2500

Single Ring Muon Like

2000

2500

Multi Ring Electron Like

2000

Multi Ring Muon Like

2000

1500

1500

1000

1000

500

500

0

0

-500

-500

-1000

-1000

104 P [MeV]pri

104 P [MeV]
pri

Figure 8.2: Missing energies and their fitting function of events with one decay electron and notagged neutrons based on MC. In this analysis, the events are separated to some event types.
The data points show the median value of the missing energy distribution in their primary ring
momentum bins.
paring the number of events in each bin separated by the reconstructed neutrino direction andenergy between data and MC output. As shown in Fig. 8.1, as the energy resolution becomesworse, the energy region with MSW eﬀect becomes wider but the diﬀerence between the total
fluxes with and without MSW eﬀect is not changed. Therefore, worse energy resolution makes

CHAPTER 8. ENERGY CORRECTION BASED ON SECONDARY INTERACTIONS 102
it diﬃcult to identify if the diﬀerence comes from neutrino oscillation or statistical flutuations.
Because FC events contains whole track, all of the produced Cherenkov light and secondary
particles are produced in the inner detector. Therefore, the eﬀect of the missing energy to
the reconstructed energy can be corrected by using the number of the detected hadrons. By
doing this, neutrino energy resolution can be improved, thus potentially contributing to better
performance in oscillation analysis. In this chapter, the correction method and the performance
evaluation result are described.

8.1

Hadron Tagging Algorithms

The energy correction is done based on the number of neutrons and decay electrons. In thissection, the detection methods of these particles are described.

8.1.1

Neutron Tagging

When a neutron is thermalized in the water, it slows down (thermalizes) and is captured byhydrogen predominantly. The neutron capture cross section of hydrogen is 1000 times larger thanthat of oxygen [129]. When the neutron is captured by a hydrogen, it emits a 2.2 MeV γ-ray via
the reaction n+p → d+γ, which can generate electrons by their electromagnetic interaction. By
using the Cherenkov light from this γ-ray, the captured neutron is tagged. The neutron lifetime
in water (thermalization followed by hydrogen capture) is measured to be 204.8±0.4 µs [130].
Because of the long lifetime, SK-I, II, and III events, whose gate width is only 1.3 µs, cannot
detect the neutron capture events. In SK-IV, the neutron capture events are searched fromthe event gates generated by both SHE trigger and the following AFT trigger. AFT trigger isgenerated just after the SHE trigger generation without any OD triggers. Combining the gates,
neutron capture events are searched from the gate with 535 µs width, which is slightly morethan twice the exponential lifetime. The intensity of the Cherenkov light from a 2.2 MeV γ-ray
is very weak. Therefore, there are many fake neutron events with comparable signal intensity to
the neutron capture events. For example, the accidental coincidence of PMT dark noise, particle
emission from nuclear spallation by cosmic ray and radiation isotopes in the detector are possible
to be the fake neutron events. To reject them, the neutron tagging algorithm is developed. Thisalgorithm consists of two processes. In this section, these processes are described briefly. The
detailed description is shown in [1].
8.1.1.1

Initial Candidate Selection

The 2.2 MeV γ-ray is searched based on the hit time distribution. As a first step, the hit time
distribution is corrected with their time-of-flight, assuming their light source is at the same
position as the primary interaction vertex. After that, the neutron candidates are searched bycounting the number of hits in a sliding 10 ns time window. This window slides from 18 µs
to 535 µs, where the origin of the time is defined as the primary interaction time. Neutron
candidates are not searched until 18 µs because there are large amount of fake neutrons fromafter pulses of the PMTs and decay electrons from muons. When seven or more hits are detected
in 10 ns window, the faster window edge is defined as t0 , and the following candidate selection

CHAPTER 8. ENERGY CORRECTION BASED ON SECONDARY INTERACTIONS 103
is applied.
7≤

≤ 50

(8.1)

N200 ≤ 200

(8.2)

N10

where N10 is the maximum number of PMT hits in the 10 ns time window sliding between t0 andt0 + 30 ns, and N200 is the number of hits in the 200 ns time window whose center is the sameas that of the 10 ns time window with the maximum number of hits. The upper limits of N10
and N200 are applied to reject high energy events like cosmic-ray muons. After this selection,
32.7% of neutron capture events are tagged. On the other hand, this selection gets 4.4 of fakeneutrons as the candidates per primary event.
8.1.1.2

Neural-Network Analysis

To remove the fake neutrons, a neural-network based signal-background classification is performed. Neural-network is a structure used for machine learning. By inputting the variablesof many signal and fake neutron events, their patterns are studied. In this analysis, sixteen
variables, which are listed on Table 8.1, are input to the network. These variables are related
Variables
N10
Nc
Nlow
N300
ϕrms
θmean
trms
min − trms (3)min − trms (6)BSenergy
BSwallN Fwall
∆N10

∆trms

(N F − BS)dis(N F − AP )dis

Description
Number of hits in 10 ns
Number of cluster hits
Number of hits on low probability PMTs
Number of hits in 300 ns
Root-mean-square of azimuthal angles of hits
Mean opening angle
Root-mean-square of hit timing
Minimum Root-mean-square of hit timing of 3 hits in 10 ns window
Minimum Root-mean-square of hit timing of 6 hits in 10 ns window
Energy reconstructed by Bonsai-fit
Distance between the detector wall and reconstructed vertex by Bonsai-fit
Distance between the detector wall and reconstructed vertex by Neut-Fit
The diﬀerence between N10 and the number of hits in 10 ns after
the subtraction of ToF which is based on secondary vertex
reconstructed by Neut-Fit.
The diﬀerence between the RMS values of the ToF subtracted hit time
distributions in 10 ns time windows calculated based on
the primary and secondary vertices.
Distance between the reconstructed vertices with NF-fit and Bonsai-fit
Distance between the reconstructed vertices with NF-fit and APFit

Table 8.1: variables for neutron tagging based on neural-network analysis. The detailed description is written in [1].

to the hit timing and position of PMTs with signal, the reconstructed vertex, and energy of theneutron capture events. The detailed description of the variables is written in [1]. This network

CHAPTER 8. ENERGY CORRECTION BASED ON SECONDARY INTERACTIONS 104

Neutron Candidates before neural net

was trained with the neutron capture events and the fake neutron events in the 500 years of theatmospheric neutrino MC events, which is taken by the dummy trigger, which is described inSec. 5.3.1. Based on the trained neural-network, the candidates in the data and MC events areclassified. Their outputs are shown in Fig. 8.3. As the value of the neural-network output inthe figure is larger, the candidate is more likely to be the signal. Currently, the threshold of the

105

Data
Total MC

104

Signal MC
103

Bckg MC

2

10

10

1

10−1
−0.2

0

0.2

0.4

0.6

0.8

1

1.2

Neural Net Output

Figure 8.3: Output of the neural-network for the data and MC of atmospheric neutrino events.
The dashed line shows the threshold to separate the neutron capture event and the background
events.
neural-network output is set to 0.832. When the output exceeds the threshold, the candidate isregarded as a neutron capture signal. As a result, the tagging eﬃciency and mis-tagging rateis evaluated as shown in Table 8.2. It is also seen that the distribution of the data agrees with
MC very well.

After Initial Candidate Selection
After Neural-NetworkEﬃciency
32.7%
21.7%

Contamination per event
4.4
0.018

Table 8.2: Eﬃciency and mis-tagging rate before and after the neural-network selection. The
contamination per events means the number of mis-tagging neutron candidates in a neutron
searching window, which have 517 µs width.

8.1.2

Decay Electron Search

Decay electrons are produced by muon decay. The muons are produced by not only the neutrino
interactions, but also by the decays of π + . When a π − is produced in the detector, it will be
captured by a nucleus before its decay. By using the number of decay electrons as an indicatorof the number of generated π + , the neutrino energy is corrected.

CHAPTER 8. ENERGY CORRECTION BASED ON SECONDARY INTERACTIONS 105
8.1.2.1

Gate Type Classification

Because of the trigger system diﬀerence between the old and new DAQ systems, the decay
electrons are searched with diﬀerent ways between the events until SK-III and the events in
SK-IV. Because the muon lifetime of about 2.2 µs [132] is longer than the old DAQ system’s
event gate width of 1.3 µs, the decay electrons in the events until SK-III are searched not only
from the primary event gates, but also in the diﬀerent gates. Fig. 8.4 shows the schema of therelationship between a primary event gate, and sub event gates. Based on the gate types, theFigure 8.4: The schema of the event gates. The event gates are produced when any triggers
described in Table 3.3 and Table 3.4 are satisfied. The width of the event gates until SK-III andfrom SK-IV are set to 30 µs and 40 µs from the start of the primary event gate, respectively.
This event gates contain one primary event gate and some sub-event gate. The primary event
gate is defined as the gate which passes the event selection in Sec. 7.1, and the sub-events gate
are the gate which passed the threshold of the global trigger. The width of the primary-eventand sub-event gates are 1.3 µs. The data in shaded region are not used. Until SK-III, only theregion where the primary and sub-event gates contain are registered.
decay electron candidates in the events until SK-III are classified to three categories as follows:
Sub-event Type (SK-I, II, III)The candidates detected in the diﬀerent gate from the primary event gate. This gate isrecorded as a sub-event.
Primary-event Type (SK-I, II, III)The candidates detected in the same gate as the primary event.

CHAPTER 8. ENERGY CORRECTION BASED ON SECONDARY INTERACTIONS 106
Split Type (SK-I, II, III)The candidates detected in the diﬀerent gate and recorded as a diﬀerent event from theprimary event.
On the other hand, SK-IV events have the larger gates with 40 µs width, which is enough todetect the decay electron without any sub-event and the other event gates. In practice, to takethe eﬀect of primary events into account, decay electron candidates in SK-IV period are alsoclassified to two types based on their detection time. The criteria are defined as follows:
Sub-Event Type (SK-IV)
The candidates that satisfy the criteria below:• ∆t > 600ns• No other decay electron candidates within ±500 ns from the event trigger
where ∆t is the time interval between the primary event signal and decay electron signal.
Primary-Event Type (SK-IV)
The candidates which do not satisfy the sub-event type criteria.
8.1.2.2

Decay Electron Search Algorithm

As a first step, the decay electron candidates are searched from the time-of-flight subtracted
hit time distribution. Because the primary-event type decay electrons exist near the primary
event signals, the time regions with the decay electron candidates contain the background hits
from the primary event. On the other hand, SK-IV does not have the trigger of the sub-event.
The candidate selection criteria take these problems into account. To perform this search, thedistribution is separated to bins with 10 ns width. The number of hits in the i-th bin is definedas N (i). Using this definition, the decay electron candidates are selected using the criteria as
follows:• SK-I, II, III– Sub-Event Type∗ To generate sub-event gate, the pulse height of the HITSUM signal is required to
exceed the global trigger threshold on Table 3.3. All of the produced subevents
are regarded as the candidates.
– The Primary and Split Types∗ The time region that satisfies all the criteria below:
N (i0 ) > N (i0 − 1) + 5 hits
N (i0 ) + N (i0 + 1) > N (i0 − 1) + 10 hits
N (i0 ) > NBG + 3σBG
NBG = (N (i0 − 1) + N (i0 − 2) + N (i0 − 3))/3√
σBG = NBG
• SK-IV(8.3)

CHAPTER 8. ENERGY CORRECTION BASED ON SECONDARY INTERACTIONS 107– ∆ti0 >450 ns∗ The time region which satisfies the criteria as follows:−j+5∑

N (i0 + i) > 9 (j = 1, 2, 3, 4, 5)

(8.4)

i=−j

– ∆ti0 < 450 ns∗ The same criteria as Eq. (8.3).
where i0 is a free parameter to search for the candidate, and ∆ti0 is the time interval between
the primary event and the faster edge of the i0 -th bin. To set SK-IV Sub-Event type event to thesame situation as the other periods data, 1.3 µs gates are generated around the Sub-Event typedecay electron candidates. After that, the decay electron candidates which pass the additional
selection below are selected.
• For All Periods– Sub-Event Type∗ Glow ≥ 0.5∗ PEtot,sub <2000 p.e.s (1000 p.e.s for SK-II)∗ NHIT50,sub ≥ 30 hits (16 hits for SK-II)– Primary-Event Type and Split Type
∗ EXCESS30 > 50 hits (25 hits for SK-II)′
or EXCESS30 /σBG
> 19.9/3. (15.5/3. for SK-II)• Only for SK-I, II, and III– Sub-Event Type∗ ∆t < 20 µs∗ NHITsub > 50 hits (25 hits for SK-II)– Primary-Event and Split Type∗ ∆t < 20 µs
where the goodness Glow follows the same definition as Section 7.1.3.5, PEtot, sub is the total
charge of ID PMTs in the sub-event gate, NHITsub is the number of hits in the trigger gate, and
NHIT50,sub is the maximum number of hits in a 50 ns time window sliding in the sub event gate.
EXCESS30 is the background subtracted value of the maximum number of hits in the 30 ns′
window sliding in the 60 ns width time region around the decay electron candidate, and σBG
isthe deviation of the number of background hits. The subtracted background hits of EXCESS30
is estimated by using the time region just before the 30 ns time window with maximum hits.
When the time window with maximum hits is set from t0 to t0 +30 ns, the number of background′ ) is estimated by using the hits in the time window of t − thits (NBG0win < ti < t0 , where twin is
the variable width changing from 20 ns to 100 ns based on the time interval between the primary

CHAPTER 8. ENERGY CORRECTION BASED ON SECONDARY INTERACTIONS 108
event and the time window with maximum hits . When the time window with maximum hits′
has N30 hits, EXCESS30 and σBG
is calculated as follows:
30 [ns]′EXCESS30 = N30 − NBG×twin√
30 [ns]
′
′
σBG
=
NBG
×
twin

(8.5)
(8.6)

For oscillation analysis, the additional selection criteria are applied as follows:• For All Periods– Sub-Event Type∗ NHIT50,sub ≥ 60 hits (30 hits for SK-II)– Primary-Event Type and Split Type∗ NHIT50,sub > 40 hits (20 hits for SK-II)• Only for SK-I, II, and III– 0.1 µs< ∆t < 0.8 µs, 1.2 µs< ∆t (if Evis < 1330 MeV)
– 1.2 µs< ∆t (if Evis ≥ 1330 MeV)
• Only for SK-IV– 0.1 µs< ∆t
The timing cut between 0.8 µs and 1.2 µs from the primary event is only applied to SK-I, II,
and III, because old DAQ system has the reflection noise from the primary events, as mentioned
in Section 3.9. The reflected signals are possible to generate the sub event trigger. The decay
electron tagging eﬃciency is evaluated with SubGeV muon MC events. In SK-III, the tagging
eﬃciency was 79% for µ+ and 64% for µ− . On the other hand, In SK-IV, the tagging eﬃciency
was 96% for µ+ and 83% for µ− . Because of the µ− capture by nuclei, the tagging eﬃciency ofµ− is smaller than that of µ+ .

8.2

Missing Energy Correction Method

In order to make the neutrino energy resolution better and improve the sensitivity of mass
hierarchy, The missing energy functions of the most energetic ring momentum, the number ofdecay electrons and tagged neutrons are prepared. To evaluate the correction function, the MC
output is used. This function will be applied not only to the MC output, but also to the data.
The correction function setting is separated to three processes below:• Event Separation Based on the Event Type• Evaluation of Missing Energy as a Function of the Primary Ring Momentum for Each
Event Types• Correction Function Tuning Based on the Events in the Resonance Energy Region
In this section, the method how to set the correction function is described.

CHAPTER 8. ENERGY CORRECTION BASED ON SECONDARY INTERACTIONS 109

8.2.1

Function Type Separation

The missing energy depends not only the number of the tagged neutrons and decay electrons,
but also the event type. For example, νµ and ν̄µ can produce muons in the neutrino interaction,and its decay electron can be tagged. On the other hand, νe and ν̄e do not produce any muonsin the neutrino interaction. Therefore, e-like and µ-like events have diﬀerent missing energy
dependence on the number of the tagged neutrons and decay electrons. In addition, whenhadrons are produced in the neutrino interaction, they or their daughter particles are possible
to produce Cherenkov rings. Based on the event types, the events are separated to the fourfunction types as shown in Table 8.3.
Event typeMulti-GeV single-ring νe -like
Multi-GeV single-ring ν̄e -like
Multi-GeV single-ring µ-like
Multi-ring νe -like
Multi-ring ν̄e -like
Multi-ring µ-like

Function group
Multi-GeV single-ring e-like
Multi-GeV single-ring µ-like
Multi-ring e-like
Multi-ring µ-like

Table 8.3: Correspondence table of event type and function group. The definition of the event
types is described in Sec. 9.

8.2.2

Energy Correction as a Function of the Primary Ring Momentum

The missing energy largely depends on the parent neutrino energy, but it is an unobservablevalue. Instead of the neutrino energy, the primary ring momentum is used for the missing energy
correction because the momentum of the lepton generated from a neutrino interaction, whichtends to be the origin of the primary ring, is correlated with the parent neutrino energy and it
is observable. Except for the multi ring muon like events, the correction function is defined asfollows:
Eν,corr1 = Evis, ID + f0 (Nn , NDcye ) + f1 (Nn , NDcye )Ppri
Evis, ID =

nring
∑

Ei

(8.7)
(8.8)

i=1

Because the multi ring muon like events did not follow linear function of the primary ring
momentum at the high momentum region, additional term is added to the correction function.
The function for the multi ring muon like events is defined as follows:Eν,corr1 = Evis, ID + f0 (Nn , NDcye ) + f1 (Nn , NDcye )Ppri+ fstep (Nn , NDcye , Ppri )(Ppri − f2 (Nn , NDcye ))2{f3 (Nn , NDcye ) (f2 (Nn , NDcye ) > Ppri )fstep (Nn , NDcye ) =0(f2 (Nn , NDcye ) < Ppri )where the variables are defined as follows:(8.9)(8.10)

CHAPTER 8. ENERGY CORRECTION BASED ON SECONDARY INTERACTIONS 110• Eν,corr1 : The corrected energy in this procedure• Nn : The number of tagged neutrons• NDcye : The number of tagged decay electrons• Ppri : Momentum of the most energetic ring
• Evis,ID : The sum of the reconstructed momentum of detected leptons with consideration
for their flavors.
• Ei : Energy of each Cherenkov ring
• nring: Number of the Cherenkov rings• fi (Nn , NDcye ): Fitting parameters
The red lines in Fig. 8.2 shows the fitted function of the missing energy as a function of the
primary ring momentum. These figures use only the events of the charged current interaction
with the target neutrinos generating no neutrons and one decay electron. The target neutrinos
are defined as νe and ν̄e for e-like samples, and νµ and ν̄µ for µ-like samples.

8.2.3

Energy Bias Rejection from the Resonance Energy Region

In the previous section, the missing energy is corrected to set the median values of the missing energy distribution of each primary momentum bins to zero. In practice, the neutrino
oscillation probability depends on the neutrino energy. The oscillation probability between
2.5 GeV< Eν,true <11.2 GeV (resonance region) shows large diﬀerence between normal and inverted hierarchy. Therefore, the energy bias in the resonance region should be small. As shown

Figure 8.5: Energy bias of SK-IV charged current events in multi ring electron like event groupas a function of the primary ring momentum (right figure), and the bias as a function of true
neutrino energy (left figure). Black dots and lines are assigned to the median values of eachenergy region. The regions between two broken lines are the resonance region.
in Fig. 8.5, Eν,corr1 of the events in the resonance region does not match to the true neutrino

CHAPTER 8. ENERGY CORRECTION BASED ON SECONDARY INTERACTIONS 111
energy. Therefore, Ecorr1 is scaled to remove the energy bias from the events in the resonanceregion. The correction function is defined as follows:Eν,corr2 = C(Nn , NDcye )Eν,corr1(8.11)

ν ,true

0.3
0.1

0.1
0
-0.1

-0.2

-0.2
-0.3
Single Ring Electron Like

0.5
0.3
0.2

(E

ν ,true

10P [MeV]
pri

Evis,IDEν ,corr1
Eν ,corr2

0.4

-0.5

ν ,true

10

0.1

10

104P [MeV]
pri

Evis,IDEν ,corr1
Eν ,corr2

0.4
0.3
0.2
0.1
0

-0.1

-0.1

-0.2

-0.2
-0.3

Multi Ring Electron Like

-0.4
-0.5

3

0.5

0

-0.3

Single Ring Muon Like

-0.4

4

-Erecon)/Eν ,true3

ν ,true-0.4-Erecon)/E

0.2

-0.1

-0.5

Evis,IDEν ,corr1
Eν ,corr2

0.4

0

-0.3

(E

0.5
0.3

(E

(E

0.2

ν ,true

Evis,IDEν ,corr1
Eν ,corr2

0.4-Erecon)/E

0.5

ν ,true-Erecon)/Eν ,true

where Eν,corr2 is the corrected energy with the additional energy scaling and C(Nn , NDcye ) is afitting parameter. The figures from Fig. 8.6 to Fig. 8.9 show the examples of the energy biases

3

10

4

10P [MeV]
pri

Multi Ring Muon Like

-0.4
-0.5

3

10

104P [MeV]
pri

Figure 8.6: Energy biases of SK-IV charged current events with the target neutrino in eachfunction group as a function of the primary ring momentum. Red, black, and green line shows
the bias of Evis,ID , Eν,corr1 , Eν,corr2 , respectively. These figures contain not only the events withthe energy in the resonance region, but also out of the region. Because the correction fromEν,corr1 to Eν,corr2 takes only the events with the energy in the resonance region into account,Eν,corr2 shows larger bias than Eν,corr1 in this figure.
as a function of neutrino true energy Eν,true , the primary ring momentum Ppri , as well as thenumber of tagged neutron and decay electrons, respectively. The energy bias ∆E is defined asfollows:∆E = (Eν,true − Erec )/Eν,true(8.12)
where Erec is the reconstructed energy. The red, black, and green lines show the biases when
Erec is set to Evis,ID , Eν,corr1 , and Eν,corr2 , respectively. As shown in Fig. 8.6, the energy
biases as a function of the visible energy are removed with the first step correction describedin the previous section. Fig. 8.7 shows the remained biases as a function of the true neutrino

0.2

true

-Erecon )/E

true

0.3

0.5
0.4

0.1

0.2

0

0

-0.1

-0.1

-0.2

-0.2

-0.3

-0.3

-0.4
-0.5

Evis,IDEν ,corr1
Eν ,corr2

0.3

true

0.1

Evis,IDEν ,corr1
Eν ,corr2(E-Erecon )/E

0.4

Single Ring Muon Like

Single Ring Electron Like

0.5

(E

true

CHAPTER 8. ENERGY CORRECTION BASED ON SECONDARY INTERACTIONS 112

-0.4
3

10

-0.5

4

10
Etrue [MeV]3

10Multi Ring Muon Like

0.2

true

-Erecon )/E

true

0.3

0.5
0.4

0.1

0.2

0

0

-0.1

-0.1

-0.2

-0.2

-0.3

-0.3
-0.4

-0.4
-0.5

Evis,IDEν ,corr1
Eν ,corr2

0.3

true

0.1

Evis,IDEν ,corr1
Eν ,corr2(E-Erecon )/E

0.4

(E

true

Multi Ring Electron Like
0.5

104
Etrue [MeV]

3

10

4

10
Etrue [MeV]

-0.5

3

10

104
Etrue [MeV]

Figure 8.7: Energy biases of SK-IV charged current events with the target neutrino in eachfunction group as a function of the true neutrino energy. The target neutrino and the line colors
definitions is the same as Fig. 8.6. The region between two broken lines is the resonance region.

ν ,true

0.3
0.1

0
-0.2

(2.5 GeV < E

ν ,true

3

4

5

6

-0.5
0

7
8
9
10
Tagged neutrons

0.5

Evis,IDEν ,corr1
Eν ,corr2

0.4
0.3

(E

0.2

ν ,true

2

0.1

3

4

5

6

< 11.2 GeV)

7
8
9
10
Tagged neutrons

Evis,IDEν ,corr1
Eν ,corr2

0.2
0.1
0

Multi Ring Muon Like

-0.2
-0.3

(2.5 GeV < E

-0.4
-0.5
0

2

0.3

-0.1
Multi Ring Electron Like

1

0.4

-0.1
-0.3

ν ,true

0.5

0
-0.2

(2.5 GeV < E-0.4-Erecon)/E1< 11.2 GeV)

ν ,trueν ,true-0.5
0Single Ring Muon Like

-0.3

-0.4

-Erecon)/E

0.1
-0.1

-0.3

ν ,true

0.2-0.1
Single Ring Electron Like

Evis,IDEν ,corr1
Eν ,corr2

0.4

0
-0.2

(E

0.5
0.3

(E

(E

0.2

ν ,true

Evis,IDEν ,corr1
Eν ,corr2

0.4-Erecon)/E

0.5

ν ,true-Erecon)/Eν ,true

CHAPTER 8. ENERGY CORRECTION BASED ON SECONDARY INTERACTIONS 113

ν ,true

1

2

3

4

5

6

< 11.2 GeV)

7
8
9
10
Tagged neutrons(2.5 GeV < E

-0.4
-0.5
0

ν ,true

1

2

3

4

5

6

< 11.2 GeV)

7
8
9
10
Tagged neutrons

Figure 8.8: Energy biases of SK-IV charged current events with the target neutrino in eachfunction group as a function of the number of tagged neutrons. The target neutrino and theline colors definitions is the same as Fig. 8.6.

ν ,true

0.3
0.1

0

(2.5 GeV < E

ν ,true

3

4

5

6

-0.5
0

7
8
9
10
Tagged Decay-e

0.5

Evis,IDEν ,corr1
Eν ,corr2

0.4
0.3

(E

0.2

ν ,true

2

0.1

3

4

5

6

< 11.2 GeV)

7
8
9
10
Tagged Decay-e

Evis,IDEν ,corr1
Eν ,corr2

0.2
0.1
0

Multi Ring Muon Like

-0.2
-0.3

(2.5 GeV < E

-0.4
-0.5
0

2

0.3

-0.1
Multi Ring Electron Like

1

0.4

-0.1
-0.3

ν ,true

0.5

0
-0.2

(2.5 GeV < E-0.4-Erecon)/E1< 11.2 GeV)

ν ,trueν ,true-0.5
0Single Ring Muon Like

-0.2
-0.3

-0.4

-Erecon)/E

0.1
-0.1

-0.3

ν ,true

0.2-0.1
Single Ring Electron Like

Evis,IDEν ,corr1
Eν ,corr2

0.4

0
-0.2

(E

0.5
0.3

(E

(E

0.2

ν ,true

Evis,IDEν ,corr1
Eν ,corr2

0.4-Erecon)/E

0.5

ν ,true-Erecon)/Eν ,true

CHAPTER 8. ENERGY CORRECTION BASED ON SECONDARY INTERACTIONS 114

ν ,true

1

2

3

4

5

6

< 11.2 GeV)

7
8
9
10
Tagged Decay-e(2.5 GeV < E

-0.4
-0.5
0

ν ,true

1

2

3

4

5

6

< 11.2 GeV)

7
8
9
10
Tagged Decay-e

Figure 8.9: Energy biases of SK-IV charged current events with the target neutrino in eachfunction group as a function of the number of tagged decay electrons. The target neutrino andthe line colors definitions is the same as Fig. 8.6.

CHAPTER 8. ENERGY CORRECTION BASED ON SECONDARY INTERACTIONS 115
energy. The low energy events have large biases which comes from the visible energy threshold at
1.33 GeV. Fig. 8.8 and Fig. 8.9 show the biases of the events in the resonance region. As shownin these figures, the biases remained at the first step correction are removed with the second
step correction. From next section, Eν,corr2 is defined as the reconstructed neutrino energy.

8.2.4

Energy Resolution

The resolutions of the reconstructed energies before and after the correction are evaluated as aperformance study. The oscillation analysis compares the two-dimensional distribution of thereconstructed zenith angle and energy of the interacted neutrino between MC and data. In thisanalysis, the neutrino direction is defined as the average of all particles’ directions weighted by
RT OT value, which is described in Sec. 6.5. This evaluation uses the 500-year MC output of
the atmospheric neutrino interaction events. The result is summarized in the Table 8.4, where
Evis,ID and Ecorr are the reconstructed energy before and after the missing energy correction.
The energy resolution σE in this table is defined as follows:()(Eν,rec − Eν,true ) − Mean(Eν,rec − Eν,true )σE = RMS(8.13)
Eν,recwhere Eν,rec is the neutrino energy reconstructed based on each method. As shown in thisEvent typeMulti-GeV single-ring νe -like
Multi-GeV single-ring ν̄e -like
Multi-GeV single-ring µ-like
Multi-ring νe -like
Multi-ring ν̄e -like
Multi-ring µ-like

SK-I
Evis,ID Eν,corr2
19.0
16.9
14.9
14.1
16.3
15.0
20.7
18.9
19.1
17.6
30.1
22.5

Energy Resolution σE [%]SK-II
SK-III
Evis,ID Eν,corr2 Evis,ID Eν,corr2
18.4
16.3
18.8
17.0
14.5
13.8
14.8
14.2
16.4
15.1
16.3
14.9
20.6
19.0
20.8
18.9
18.1
16.7
18.7
17.2
30.2
22.3
29.8
22.2

SK-IV
Evis,ID Eν,corr2
17.0
12.5
12.7
11.2
16.0
12.9
19.6
14.5
16.8
13.6
30.1
18.4

Table 8.4: Energy resolution of each event type for each period.
table, the energy resolution is largely improved. Especially, the improvement of SK-IV is moresignificant than the other period owing to the information of the tagged neutron multiplicity.

8.2.5

Systematic Error Estimation

The energy correction functions are set based on the MC If the interaction model is not correct
and the relationship between the neutrino energy and the number of decay electrons, tagged
neutrons and primary ring momentum is diﬀerent from the real data, the energy correction
does not work correctly. Therefore, the diﬀerence of the number of tagged decay electrons andneutrons between data and MC may smear the corrected energy distribution. To perform theoscillation analysis with the corrected energy, four of new systematic errors are introduced. In
this section, the estimation methods of these systematic errors are described.

CHAPTER 8. ENERGY CORRECTION BASED ON SECONDARY INTERACTIONS 116
8.2.5.1

Energy Shift from the Diﬀerence of Neutrino Interaction Models

The fraction of each neutrino interaction mode depends on the neutrino interaction model.
In addition, the multiplicities of the captured neutrons and produced decay electrons in eachinteraction mode also depend on the model. As mentioned in Section 5, the MC events use
NEUT5.3.6 and its interaction model. Because there are no data about the relationship between
the neutron multiplicity and the neutrino energy, the uncertainty is evaluated by comparing
between events based on the diﬀerent neutrino interaction models, NEUT 5.3.6, NEUT 5.4.0,
and GENIE 2.12.6[133] [134]. These interaction models are described in Sec. 5.2. Currently,
the oscillation analysis is done by using MC output of 500 year atmospheric neutrino events.
Because it takes long time to newly simulate it, the MC output based on NEUT 5.3.6 arereweighted to have the same neutron, muon, interaction mode, and invariant mass distributions
as the estimated distributions of the other interaction models. The distribution estimation isdone by using toy MC. The method is shown as follows:
1. Make the PDF Tables of Neutron and µ+ Multiplicities
The neutron and decay electrons can be produced via the decay of the hadrons producedin the neutrino interaction. For example, π + produces decay electrons as follows:
π + → µ+ νµ→ e+ νµ νe ν̄µ(8.14)

To estimate how many decay electrons and neutrons are produced by the hadrons, theprobability density function of the produced number of neutron and decay electrons haveto be evaluated for each hadron. In practice, PDF tables of neutron and µ+ multiplicities
as a function of the momentum are prepared. The number of µ+ indicates that of decay
electrons. Because of the π − capture, few µ− is produced and the eﬀect against the toy MC
is negligible. To make the PDF table of the neutron and µ+ multiplicities for each hadron
momentum region, the behaviors of the various hadrons in the detector are simulated bychanging their momenta. The simulated hadrons are proton, neutron, π ± , K ± , KL0 , KS0 ,and η.
2. Estimate Neutron and Decay Electron Multiplicities
By using the neutrino interaction generators with diﬀerent interaction models, the types
and momenta of the particles that are produced in the atmospheric neutrino interactions
are stored event by event. From this information and the PDF tables which is made in theprevious step, the number of neutrons and µ+ can be estimated without MC simulation.
Fig. 8.10 and Fig. 8.11 show the estimated captured neutron and produced µ+ multiplicitiesof each flavor neutrino events based on the three models in the resonance region.
3. Estimate the Neutrino Events Based on Each Interaction ModelIn the previous process, the number of decay electrons and neutrons can be estimated. To
estimate the remaining parameters to be used for the missing energy correction, for example the reconstructed momenta of Cherenkov rings, the output of Monte Carlo simulation
which is described in Sec. 5 is used. Toy MC have the information of the interaction mode,
invariant mass, neutrino energy, the number of decay electrons, and the number of captured neutrons. By weighting the MC events to make the five dimensional distributions of
the interaction mode, invariant mass, neutrino energy as well as number of decay electrons

CHAPTER 8. ENERGY CORRECTION BASED ON SECONDARY INTERACTIONS 117

ν e CC Events
0.25

Models and Its Mean Values

0.2

NEUT (Current version): 5.03
GENIE: 5.26

0.15

NEUT (New version): 5.37

Entries (Normalized)

Entries (Normalized)

0.25

Models and Its Mean Values

0.2

0.1

0.1

0.05

0.05

NEUT (Current version): 5.23
GENIE: 5.43

0.15

NEUT (New version): 5.57

Entries (Normalized)

Entries (Normalized)

0.25

Models and Its Mean Values

0.2ν µ CC Events
Models and Its Mean Values

0.2

NEUT (Current version): 4.22
GENIE: 4.67

0.15

0.1

0.1

0.05

0.05

0
0
1
2
3
4
5
6
7
8
9
10
The Number of Captured Neutrons (Based on Toy MC)

NEUT (New version): 4.48

0
0
1
2
3
4
5
6
7
8
9
10
The Number of Captured Neutrons (Based on Toy MC)ν µ CC Events
0.25

NEUT (Current version): 4.14
GENIE: 4.58

0.15

0
0
1
2
3
4
5
6
7
8
9
10
The Number of Captured Neutrons (Based on Toy MC)

ν e CC Events

NEUT (New version): 4.57

0
0
1
2
3
4
5
6
7
8
9
10
The Number of Captured Neutrons (Based on Toy MC)

Figure 8.10: Captured neutron multiplicity of the events in the resonance energy region basedon the toy MC output.

CHAPTER 8. ENERGY CORRECTION BASED ON SECONDARY INTERACTIONS 118

0.8

Models and Its Mean Values

0.7

NEUT (Current version): 0.67

0.6

GENIE: 0.66

0.5

NEUT (New version): 0.68

0.4

Entries (Normalized)

Entries (Normalized)

ν e CC Events
0.9

0.8

NEUT (Current version): 0.20

0.6

GENIE: 0.17

0.5

NEUT (New version): 0.20

0.4
0.3

0.2

0.2

0.1

0.1
1

2

0
0

3
4
5
6
7
8
9
10
The Number of Muons (Based on Toy MC)

Models and Its Mean Values

0.7

0.3

0
0

ν e CC Events

0.9

1

2

0.8

Models and Its Mean Values

0.7

NEUT (Current version): 1.68

0.6

GENIE: 1.67

0.5

NEUT (New version): 1.69

0.4

Entries (Normalized)

Entries (Normalized)

ν µ CC Events
0.9

0.8

NEUT (Current version): 1.20

0.6

GENIE: 1.18

0.5

NEUT (New version): 1.19

0.4
0.3
0.2

0.1

0.1
2

3
4
5
6
7
8
9
10
The Number of Muons (Based on Toy MC)

Models and Its Mean Values

0.7

0.2

1

ν µ CC Events

0.9

0.3

0
0

3
4
5
6
7
8
9
10
The Number of Muons (Based on Toy MC)

0
0

1

2

3
4
5
6
7
8
9
10
The Number of Muons (Based on Toy MC)

Figure 8.11: Generated muon multiplicity of the events in the resonance energy region based on
the toy MC output.

CHAPTER 8. ENERGY CORRECTION BASED ON SECONDARY INTERACTIONS 119
and captured neutrons to the same as the toy MC output based on each interaction model,
the relationship between the visible energy and the reconstructed energy based on eachinteraction model is estimated.

1.1

1.08
1.06

1 ring e-likeMaximum shift
GENIE/NEUT5.3.6
NEUT5.4.0/NEUT5.3.6

Ecorr2 rate

Ecorr2 rate

Based on this relationship between the visible energy and the reconstructed energy, the systematic error is evaluated by comparing the rate of the reconstructed energy based on NEUT 5.3.6
and the other interaction models as a function of the visible energy. Fig. 8.12 shows the recon1.1

1.08
1.06

1.04

1.04

1.02

1.02

1

1

0.98

0.98

0.96

0.96

0.94

0.94

0.92

0.92

0.9

10 E

1.06

104E

[MeV]

mring e-likeMaximum shift
GENIE/NEUT5.3.6
NEUT5.4.0/NEUT5.3.6

vis,ID

Ecorr2 rate

Ecorr2 rate

1.08

0.9

4

vis,ID

1.1

1 ring µ -likeMaximum shift
GENIE/NEUT5.3.6
NEUT5.4.0/NEUT5.3.6

1.1

1.08
1.06

1.04

1.04

1.02

1.02

1

1

0.98

0.98

0.96

0.96

0.94

0.94

[MeV]

mring µ -likeMaximum shift
GENIE/NEUT5.3.6
NEUT5.4.0/NEUT5.3.6

0.92

0.92
0.9

0.9

4

10 E

vis,ID

[MeV]

104E

vis,ID[MeV]

Figure 8.12: Reconstructed energy ratios of the SK-IV outputs based on NEUT 5.4.0 and
GENIE 2.12.6 to the SK-IV output based on NEUT 5.3.6 output. The red lines show Evis,ID =1.33 GeV.
structed energy distribution’s median value ratios of the SK-IV outputs based on NEUT 5.4.0
and GENIE 2.12.6 to that based on NEUT 5.3.6 as a function of the visible energy. The blackline shows the absolute value of the maximum energy shift of NEUT 5.4.0 and GENIE 2.12.6,
which is defined as follows:
Rmax = 1 +|Median(Ecorr2,max ) − Median(Ecorr2, NEUT 5.3.6 )|Ecorr2, NEUT 5.3.6(8.15)

where Ecorr2,max and Ecorr2, NEUT 5.3.6 are the corrected energy based on interaction model withthe largest energy shift and NEUT 5.3.6. It is defined as the systematic error value. The redline shows the position of Evis,ID = 1.33 GeV. Because single ring events have the threshold at
1.33 GeV visible energy assuming all particles are electron like, the energy shift of the single
ring events with Evis,ID tends to be large below the red line. The diﬀerence of the reconstructedCHAPTER 8. ENERGY CORRECTION BASED ON SECONDARY INTERACTIONS 120
energies between the interaction models in the resonance energy region is 4.8% at maximum.
This shift was introduced into the oscillation analysis as the uncertainty of the reconstructed
momentum. This uncertainty does not change the total number of neutrino events, but change
the reconstructed momentum distribution. As mentioned in 9.1.5, the oscillation analysis uses
the number of events in some bins separated by the reconstructed energy and the direction.
Owing to the reconstructed energy uncertainty, the number of events in each bin is possible tobe changed.
8.2.5.2

Energy Shift from the Uncertainty of Neutron Tagging Eﬃciency

Currently, the systematic error of the neutron tagging eﬃciency is set to 10% based on the resultof the calibration with Am/Be neutron source [1]. The eﬀect to the systematic error of thereconstructed energy is evaluated. This study also uses the MC outputs of 500 year atmospheric
neutrino events which is weighted based on NEUT 5.3.6, NEUT 5.4.0, and GENIE 2.12.6.
As mentioned in Table 8.2, the neutron tagging eﬃciency is evaluated to 21.7%. Instead of
applying the neutron tagging algorithm, the tagged neutrons are selected from the true captured
neutrons by applying the probability of the tagging eﬃciency. The dependence of eﬃciency on
the number of neutrons is not taken into account. This is because the capture time of thermalized
neutrons (∼ 2.0 × 105 ns) is much longer than the hit time distributions of signals from eachneutron capture event (∼ 10 ns), and its eﬀect is considered to be negligible. In addition, bychanging the eﬃciency, the eﬀect of the tagging eﬃciency uncertainty to the number of tagged
neutrons and the reconstructed energy can be estimated. Considering the 10% uncertainty ofthe tagging eﬃciency, the number of neutrons is estimated with three tagging eﬃciencies of
21.7% (Standard), 21.7% × (1 − 0.1) = 19.5% (Low), and 21.7% × (1 + 0.1) = 23.9% (High).
By using the estimated number of the tagged neutrons based on these eﬃciencies, the neutrinoenergy is reconstructed. Fig. 8.13 shows the reconstructed energy ratios of events with the”Low” and ”High” eﬃciencies to that with the ”Standard” eﬃciency. The systematic errorof the reconstructed energy from the neutron tagging eﬃciency uncertainty is defined as the
maximum shift value of the six conditions. The systematic error of the reconstructed energy inthe resonance energy region is 1.1% at maximum.
8.2.5.3

Energy Shift from the Uncertainty of Decay Electron Multiplicity

As mentioned in Section 8.1.2, π + decay produces µ+ , which is the origin of decay electron.
However, π + is possible to interact with 16 O. Based on the diﬀerence of its cross section between
data and MC, the uncertainty of the π + decay probability is evaluated [121]. This evaluation isbased on the cross section of π + −16 O interaction.
Fig. 8.14 is the comparison between the π + − O cross sections based on the theoretical
calculation and data. The uncertainty of the cross section as a function of π + momentum isevaluated based on this figure. By using this uncertainty as a diﬀerence of the mean free pathof π + in the water, π + decay probability is evaluated.
Fig. 8.15 shows the π + decay probability as a function of π + momentum. The probability
uncertainty is changed from 3% to 13% depending on the π + momentum. For this systematic
error estimation, the uncertainty is fixed to 10%. Taking this uncertainty into account, thenumber of decay electrons is estimated based on the π + decay probability and π + multiplicity.
As a first step, the π + decay probability is evaluated by using the MC output of the 100
year atmospheric neutrino events which generate only one π + . Fig. 8.16 shows the π + decay

CHAPTER 8. ENERGY CORRECTION BASED ON SECONDARY INTERACTIONS 121

1 ring µ -likeMaximum shift
NEUT5.3.6 (+10% Eff)
NEUT5.3.6 (-10% Eff)
GENIE2.12.6 (+10% Eff)
GENIE2.12.6 (-10% Eff)
NEUT5.4.0 (+10% Eff)
NEUT5.4.0 (-10% Eff)

1.08
1.06

1.04

Ecorr2 rate

Ecorr2 rate

1 ring e-like

1.1

1.08
1.06

1.04

1.02

1.02

1

1

0.98

0.98

0.96

0.96

0.94

0.94

0.92

0.92

0.9

0.9

4

10 E

vis,ID

vis,ID

1.06

1.04

Maximum shift
NEUT5.3.6 (+10% Eff)
NEUT5.3.6 (-10% Eff)
GENIE2.12.6 (+10% Eff)
GENIE2.12.6 (-10% Eff)
NEUT5.4.0 (+10% Eff)
NEUT5.4.0 (-10% Eff)

1.1

1.08
1.06

1.04

1.02

Maximum shift
NEUT5.3.6 (+10% Eff)
NEUT5.3.6 (-10% Eff)
GENIE2.12.6 (+10% Eff)
GENIE2.12.6 (-10% Eff)
NEUT5.4.0 (+10% Eff)
NEUT5.4.0 (-10% Eff)

1.02

1

1

0.98

0.98

0.96

0.96

0.94

0.94

0.92
0.9

[MeV]

mring µ -like

Ecorr2 rate

Ecorr2 rate

1.08

104E

[MeV]

mring e-like

1.1

Maximum shift
NEUT5.3.6 (+10% Eff)
NEUT5.3.6 (-10% Eff)
GENIE2.12.6 (+10% Eff)
GENIE2.12.6 (-10% Eff)
NEUT5.4.0 (+10% Eff)
NEUT5.4.0 (-10% Eff)

1.1

0.92
0.9

4

10 E

vis,ID

[MeV]

104E

vis,ID[MeV]

Figure 8.13: Reconstructed energy ratios of the SK-IV events with high (23.9%) and low (19.5%)
tagging eﬃciencies to the SK-IV events with normal tagging eﬃciency (21.7%). The red line
definition is the same as Fig. 8.12.

CHAPTER 8. ENERGY CORRECTION BASED ON SECONDARY INTERACTIONS 122

Figure 8.14: The cross section of π + − O interactions based on data and theoretical calculation
which is used in MC. The experimental data is taken from [131]. This figure is taken from [121]probability as a function of its momentum. The red line shows the fitting function defined asfollows:
Probπ+ →Dcye (Pπ+ ) = C0 × Erfc(C1 × (Pπ+ − C2 ))()(Pπ+ − µ)2+ C3 × exp −2σ 2(8.16)

where Ci , µ, and σ are the fitting parameters, and Pπ+ is the π + momentum. In this figure,there is a peak around 2 GeV/c. It is considered to be caused by the peak of the cross sectionof the scattering with proton [132]. The interaction is possible to make the π + energy smaller
than the peak of the π + capture cross section, which is around 300 MeV/c [121]. After that,the probabilities are multiplied by (1.0+0.1) (High) and (1.0-0.1) (Low). By using the numberof estimated muon and the neutrino flavor, the number of tagged decay electrons is estimated.
Because of the µ± existence in νµ and ν̄µ events and the decay electron tagging eﬃciency
dependence on the number of the generated decay electrons, the tagging eﬃciency is evaluatedfor each pattern of the number of decay electrons and the flavor of the generated lepton in theneutrino interaction. Because both e+ and e− do not produce any muons , they use the sameprobability table. Fig. 8.17 shows the reconstructed energy ratios of the SK-IV events with highand low π + decay probabilities to the SK-IV events with the normal probability as an exampleof the energy shift from the uncertainty of the number of decay electrons. The systematic errorof the reconstructed energy in the resonance energy region is 2.8% at maximum.

CHAPTER 8. ENERGY CORRECTION BASED ON SECONDARY INTERACTIONS 123(Generated µ +)/(Generated π+)

Figure 8.15: The π + decay probability as a function of its momentum[121]. The black solidline shows the π + decay probability in our MC simulation output, and the red dotted and blue
dashed lines show the probabilities with the longer and smaller mean free path of π +

1.4
1.2
1
0.8
0.6
0.4
0.2
0
10-2

10-1

1

102
10
π+ Momentum [GeV]

-1.5
-2
-0.5
-1
0.5
0
1.5
12

Figure 8.16: The π + decay probability as a function of π + momentum. The probability falls
from 0.1 GeV to 1 GeV because the cross section of π + absorption on oxygen have a peak around
200 MeV/c.

CHAPTER 8. ENERGY CORRECTION BASED ON SECONDARY INTERACTIONS 124

1 ring µ -likeMaximum shift
NEUT5.3.6 (+10% Eff)
NEUT5.3.6 (-10% Eff)
GENIE2.12.6 (+10% Eff)
GENIE2.12.6 (-10% Eff)
NEUT5.4.0 (+10% Eff)
NEUT5.4.0 (-10% Eff)

1.08
1.06

1.04

Ecorr2 rate

Ecorr2 rate

1 ring e-like

1.1

1.08
1.06

1.04

1.02

1.02

1

1

0.98

0.98

0.96

0.96

0.94

0.94

0.92

0.92

0.9

104E

vis,ID

0.9

vis,ID

1.06

1.04

1.1

1.08
1.06

1.04

1.02

Maximum shift
NEUT5.3.6 (+10% Eff)
NEUT5.3.6 (-10% Eff)
GENIE2.12.6 (+10% Eff)
GENIE2.12.6 (-10% Eff)
NEUT5.4.0 (+10% Eff)
NEUT5.4.0 (-10% Eff)

1.02

1

1

0.98

0.98

0.96

0.96

0.94

0.94

0.92

0.92

0.9

[MeV]

mring µ -likeMaximum shift
NEUT5.3.6 (+10% Eff)
NEUT5.3.6 (-10% Eff)
GENIE2.12.6 (+10% Eff)
GENIE2.12.6 (-10% Eff)
NEUT5.4.0 (+10% Eff)
NEUT5.4.0 (-10% Eff)

Ecorr2 rate

Ecorr2 rate

1.08

104E

[MeV]

mring e-like

1.1

Maximum shift
NEUT5.3.6 (+10% Eff)
NEUT5.3.6 (-10% Eff)
GENIE2.12.6 (+10% Eff)
GENIE2.12.6 (-10% Eff)
NEUT5.4.0 (+10% Eff)
NEUT5.4.0 (-10% Eff)

1.1

104E

vis,ID

0.9
[MeV]

104E

vis,ID[MeV]

Figure 8.17: Reconstructed energy ratios of the events with high and low π + decay probabilities
to the events with normal probability for SK-IV. The red line definition is the same as Fig. 8.12.

CHAPTER 8. ENERGY CORRECTION BASED ON SECONDARY INTERACTIONS 125
8.2.5.4

Energy Shift from the Neutrons Produced by π − Capture

Most π − s at rest are captured by nucleon before their decay. The captured π − s are possible toproduce neutrons as follows:
π− + n + p → n + nAA(8.17)

Z +π−→ γ+n+

Z +π−→ γ+n+n+

A−1(Z − 1)
A−2(Z − 1)(8.18)(8.19)

98% of π − is captured by Eq. (8.17) process and the residual 2% of π − s emit γ-rays. The γ-ray

Figure 8.18: The energy spectrum of the γ-ray produced by the π − capture by 12 C nucleus[135]. The histogram colored in magenta shows the data, and the other color lines show the
expected energy spectra based on the diﬀerent models. The systematic error from π − capture isestimated by using the energy spectrum of the data itself and the diﬀerence between the energy
spectra of the data and SF+RPA model.
produced by 12 C’s radiative absorption of π − follows the energy spectrum shown in Fig. 8.18.
In this figure, the red line shows the expected energy spectrum of the events with the single
neutron emission[135], and magenta line shows the measurement result. As shown in this figure,the γ-ray energy spectrum of the single neutron emission events covers only a part of the whole
distribution. The diﬀerence between energy spectra of data and SF+RPA model at low energy is
considered to be caused by multiple neutron emission events. In Super-Kamiokande, the π − does
not tend to be captured by 12 C but 16 O. However, the number of neutrons from the π − capture isestimated with this distribution as an approximation. From Fig. 8.18, the two neutrons emission
events seem to exist below 100 MeV region. To estimate the eﬀect of the neutron emission fromthe π − capture, the number of produced neutrons in the MC is set depending on the γ energy
as follows:• Single neutron emission: Eγ > 100 MeV

CHAPTER 8. ENERGY CORRECTION BASED ON SECONDARY INTERACTIONS 126• two neutron emission: Eγ < 100 MeV

1.1

1.08
1.06

1.04

1 ring e-likeMaximum shift
NEUT5.3.6
NEUT5.4.0
GENIE2.12.6

Ecorr2 rate

Ecorr2 rate

These criteria are introduced in the simulation software, and compared to the reconstructedenergy. Fig. 8.19 shows the reconstructed energy ratios of the events based on the normal MC
1.1

1.08
1.06

1.04

1.02

1.02

1

1

0.98

0.98

0.96

0.96

0.94

0.94

0.92
0.9

0.92
10 E

1.06

1.04

[MeV]

mring e-likeMaximum shift
NEUT5.3.6
NEUT5.4.0
GENIE2.12.6

1.1

1.08
1.06

1.04

1.02

[MeV]

mring µ-likeMaximum shift
NEUT5.3.6
NEUT5.4.0
GENIE2.12.6

1.02

1

1

0.98

0.98

0.96

0.96

0.94

0.94

0.92
0.9

104E

vis,ID

Ecorr2 rate

Ecorr2 rate

1.08

0.9

4

vis,ID

1.1

1 ring µ-likeMaximum shift
NEUT5.3.6
NEUT5.4.0
GENIE2.12.6

0.92
0.9

4

10 E

vis,ID

[MeV]

104E

vis,ID[MeV]

Figure 8.19: Reconstructed energy ratios of the events based on the normal MC to the events
based on the MC which has the neutron emission process from the π − capture. The red line
definition is the same as Fig. 8.12.
to the events based on the MC which has the neutron emission process from the π − capture. The
systematic error of the reconstructed energy in the resonance energy region is 1.8% at maximum.

Chapter 9

Oscillation Analysis
The eﬀect of the neutrino oscillation can be detected as the change of neutrino flux. The neutrino
oscillation depends on some parameters, which are mixing angles, CP phase, and mass squared
diﬀerences. Because the normal and inverted mass hierarchies are the diﬀerent patterns of themass squared diﬀerences, they can be identified based on the oscillation probability. To search
for the suitable parameters from the detected atmospheric neutrino flux and its likelihood,the two-dimensional distributions of the detected particle’s energy and incident direction arecompared between MC output and data. In this chapter, the current analysis method[136], the
new method, and their sensitivities are described.

9.1

Analysis Method

9.1.1

χ2 Calculation

To search for the suitable parameters, χ2 is calculated for each oscillation pattern. In this study,
χ2 is calculated based on the log likelihood-ratio test as follows:χ2 ≡ −2 lnL(E, O)
L(O, O)(9.1)

where E and O are the expected and observed fluxes, respectively. In this analysis, the likelihoods are calculated assuming the number of events in each bin following the probability basedon Poisson distribution. When the systematic error is not considered, χ2 is calculated as follows:(( ()/ ( ())∏ e−En E On )∏ e−Om OOm )

)
∑(
On
χ
=2
En − On + On ln
On !
Om !
En
n
m
n
(9.2)
where En and On are the expected and observed numbersofevents
in
n-th
bin.
The
systematic∑∑
errors are introduced by converting En → En (1 + i fni ϵi ) and adding χ2sys = i (ϵi /σi )2 as
follows:()∑∑∑ ( ϵi )2
O
n
i
2
2
2
∑
En (1 +
fn ϵi ) − On + On lnχ = χobs + χsys = 2+(9.3)
En (1 + i fni )σin2= χ2obs = −2 ln

n

m

i

i

127

CHAPTER 9. OSCILLATION ANALYSIS

128

where fni is the fractional change of the n-th bin’s event rate due to the i-th systematic error. ϵiis set to have the minimum χ2 value. In this equation, the systematic error is assumed to follow
Gaussian distribution.

9.1.2

FC Event Type Criteria

Based on the selection described in Sec. 7.1, the candidates of neutrino interactions are selected.
Because neutrino interaction have large types, for example charged and neutral current interaction, and their signals are largely diﬀerent, the oscillation analysis is done after the event
type separation. The final FC samples are separated to some event types, which are listed inTable 9.1. As mentioned in Sec. 8.2.1, the multi-GeV single ring samples and multi-ring samples use the corrected energy for the new oscillation analysis. The other samples use the same
method as the current analysis.
Event type(dcy=decay)

Sub-GeV
e-like 0 dcy
e-like 1 dcy
1-Ring π 0 -like
2-Ring π 0 -likeµ-like 0 dcyµ-like 1 dcyµ-like 2 dcy
Multi-GeV
single-ring νe likesingle-ring ν̄e likesingle-ring µ like
Multi-ring νe like
Multi-ring ν̄e like
Multi-ring otherMulti-ring µ likedecay electrons

The Number of
Cherenkov rings andparticle ID of
the most energetic ring0
≥1

1,electron
1,electron
1, electron
2, electron
1,muon1,muon1,muon

0
1
≥2
≥1
0

1,electron
1,electron
1, muon≥2, electron≥2, electron≥2≥2, muon

Other criteria

Evis ≤ 1.33 GeV
pe ≥ 0.1 GeV
85 MeV< mπ <215 MeV
Evis ≤ 1.33 GeV
pµ > 0.2 GeVEvis > 1.33 GeV
Evis > 1.33 GeV
pass νe /ν̄e separation
Fail νe /ν̄e separation
Evis > 1.33 GeV
or pµ > 0.6 GeV/c

Table 9.1: The event types and their separation criteria.

9.1.3

Separation Method of Multi-GeV Multi-ring Electron Like Events

Based on the likelihood analysis, the multi-GeV multi-ring electron like events are separated tothree types, which are Multi-ring νe like, Multi-ring ν̄e like, and Multi-ring other. The likelihood
is evaluated based on the variables below:
1. Likelihood about the particle identification

CHAPTER 9. OSCILLATION ANALYSIS

129

2. Momentum fraction of the most energetic ring
3. The number of decay electrons
4. The distance between the vertex of the neutrino interaction and muon decay
Multi-ring other samples are the events which are failed to be identified if the interacted particle
is a neutrino or an anti-neutrino. The detailed description about this separation is written in[136].

9.1.4

PC Event Types

PC events are separated only to two types, which are stopped and through-going PC events,
by checking if the produced particles stop in or go out the outer detector OD. The PC events
stopped in OD have better energy resolution than through-going events. This separation is doneby using the charge observed in the OD. If the most or second energetic ring is identified to beµ-like ring, which has long track, and the charge of OD PMTs is less than the expected charge
from through-going muon, the event is selected as a stopped PC event.

9.1.5

Binning Setting

To measure the diﬀerence between MC and data, the detected events are separated based on thereconstructed zenith angle of its incident direction and energy, which are reconstructed basedon Sec. 6 and Sec. 8 methods.
Fig. 9.1 shows the schema of the energy binning for all event types. Originally, the events inthe multi-GeV and multi-ring types are separated to 5 energy bins at maximum. To detect theresonant enhancement of atmospheric neutrino oscillation probability eﬃciently, the FC samples
with 2.5 - 12.5 GeV are separated to 8 bins with narrower width except for the Multi-ring Others
events.

9.2

Systematic Errors

As shown in Eq. (9.3), this χ2 calculation takes the systematic errors into account. In thissection, the systematic errors are described.

9.2.1

Systematic Errors Related to Neutrino Flux

9.2.1.1

Absolute Normalization

The error of the absolute normalized flux is taken as systematic error based on [137], taking theatmospheric muon data into account. Fig. 9.3 shows the uncertainties of the pion production(δπ ), Kaon production (δK ), hadronic interaction cross section (δσ ), and the atmospheric density
profile (δair ). In these uncertainties, the total uncertainty from δπ and δσ is taken as the absolutenormalization uncertainty, and δK and δair are not included in this error. The uncertainty of
the kaon production and the air density profile are included in the other errors. For Eν >10 GeV events, the uncertainty from the diﬀerence between the neutrino flux models, which are
Bartol[67], FLUKA[68] and Honda flux[50], has to be also considered. Therefore, 5% is assigned
as the additional normalization uncertainty for multi-GeV FC, PC, and upward-going stopping
muon samples.

CHAPTER 9. OSCILLATION ANALYSIS

130

Figure 9.1: Schema of the energy binning for oscillation analysis. The energy region in the whiteboxes are separated to 10 bins based on the zenith angle. On the other hand, the gray boxes
have only one bin.

CHAPTER 9. OSCILLATION ANALYSIS

131

Figure 9.2: Binning for 1-ring Multi-GeV events and Multi-ring e-like events and the oscillation
probability from νµ to νe based on normal hierarchy. The oscillation probability table is takenfrom [35].

CHAPTER 9. OSCILLATION ANALYSIS

132

Figure 9.3: The upper figure shows relative uncertainties of the error sources for the atmospheric
neutrino flux[137]. The total uncertainty (δtot ) includes the eﬀects of all the error sources in thisfigure, which are the pion and kaon production (δπ and δK ), hadronic interaction cross section(δσ ) and the atmospheric density profile (δair ). The lower figure shows the uncertainty of theabsolute normalization which is calculated based on the sum of δπ and δσ [136].

CHAPTER 9. OSCILLATION ANALYSIS
9.2.1.2

133

Flux Ratios

The uncertainty of the neutrino flux ratios between the diﬀerent neutrinos are considered. The
introduced uncertainties of the ratios are listed as follows:• (νµ + ν̄µ )/(νe + ν̄e ) Ratio• ν̄e /νe Ratio• ν̄µ /νµ Ratio• Up/Down Ratio• Horizontal/Vertical Ratio
The uncertainties are evaluated based on the comparison between Honda flux[137] and some
diﬀerent models, which are FLUKA[68] and the Bartol flux models[67].
9.2.1.3

K/π Ratio

When the primary cosmic ray has high energy, the K production probability is increased. theuncertainty of the production ratio between K and π is evaluated by using the data from
SPY experiment[138]. The errors are 5% for Eν < 100 GeV and it increases with energy
logarithmically until 20% at 1 TeV.
9.2.1.4

Solar Activity

Solar activity cycles with 11 years period, and it aﬀects to the primary cosmic-ray flux. Asa systematic error of this analysis, the neutrino flux variation caused by ±1 year variation is
introduced.
9.2.1.5

Neutrino Path Length

The neutrino path length depends on not only the neutrino’s incident direction, but also the
production height in the atmosphere. The atmospheric neutrino production probability depends
on the atmosphere density, and its uncertainty finally aﬀect to the neutrino flux at the detector
by changing the path length. The systematic error of the atmosphere density is set to 10% basedon the diﬀerence between US-standard 76, which is used for oscillation analysis, and MSISE-90
models[139].
9.2.1.6

Matter Eﬀect

As shown in Sec. 2.5, the neutrino oscillation probability in the matter depends on the electron
density. The matter density in the Earth has been already measured by [140]. However, theelectron density depends on the chemical composition of the Earth, and it still has uncertainty.
The Earth core is considered to consist of heavy elements such as iron. The expected flux
is estimated assuming the Earth consists of light elements, which have the same number ofneutrons, protons and electrons. The electron density of iron is 6.8% smaller than the lightelements. This diﬀerence is introduced for the systematic error of the matter eﬀect.

CHAPTER 9. OSCILLATION ANALYSIS

9.2.2

Systematic Errors Related to Neutrino Interaction

9.2.2.1

Charged Current Quasi-Elastic Scattering134

For the uncertainty on the charged current quasi-elastic scattering (CCQE) cross sections, somekinds of systematic errors, which are the cross sections of (νe + ν¯e ) and (νµ + ν¯µ ), and the crosssection ratios of (νµ + ν̄µ )/(νe + ν̄e ) and ν̄/ν, are introduced. Their values are assigned to thediﬀerence between Smith and Monitz Fermi-gas model[71] and the model of Nieves et al.[109].
9.2.2.2

Axial Mass

Assuming the uncertainty of axial mass based on the results of K2K[72] (MA = 1.20 ± 0.12 GeV)
and MiniBooNE[73] (MA = 1.23 ± 0.2 GeV), the diﬀerence of the four momentum transfer
distribution between MA = 1.2 GeV and MA = 1.1 GeV are introduced as the systematic error.
9.2.2.3

Meson Exchange Current

Meson exchange current has large uncertainty, and its cross section has not been measured yet.
Therefore, the uncertainty of its cross section is set to 100%.
9.2.2.4

Single Meson Production

The systematic error for the single meson production cross section is obtained by comparing theestimation with the experimental data of νµ p → µ− pπ + . Because there are few experimental
data of νµ p → µ− pπ + , additional 40% uncertainty is introduced based on the cross section
comparison with the model of Hernandez[141]. The uncertainty of the single pion production
cross section ratio between ν and ν̄ is also estimated based on the comparison with the Hernandez
model.
9.2.2.5

Deep Inelastic Scattering

When the neutrino energy is more than 10 GeV, the diﬀerence of the cross sections between theexperiments and calculation based on NEUT is within 5.0%. Because the uncertainty increases
at the low energy, additional systematic error is estimated based on the diﬀerence from CKMT
model[142]. As another systematic error, the uncertainty of the four momentum transfer distribution in deep inelastic scattering is also introduced based on the diﬀerence between GRV98
with and without the Bodek-Yang correction.
9.2.2.6

Coherent Pion Production

Based on the diﬀerence between the measurement by SciBooNE[143] and Rein and Sehgal model,
which is used for NEUT, the systematic errors of coherent pion production cross sections in NC
and CC νe interaction are estimated to 50%. On the other hand, the systematic error of the crosssection in CC νµ interaction is conservatively set to 100% because of non-existing experimentresult.
9.2.2.7

NC/CC Ratio

20% systematic error is assigned to the cross section ratio of NC interaction to CC interaction.

CHAPTER 9. OSCILLATION ANALYSIS
9.2.2.8135ντ CC Cross Section

25% systematic error is assigned to ντ CC total cross section based on the comparison with the
Hagiwara model[144].

9.2.3

Systematic Errors Related to Event Selection and Reconstruction

9.2.3.1

FC Event Selection

By comparing the cut variables for FC event selection between the observed data and MC, theuncertainties of the selection eﬃciencies are evaluated to 0.2%, 0.2%, 0.8%, and 1.3% for SK-I to
SK-IV, respectively. The uncertainty of the background contamination rate in the final sample
is changed by the period and the particle. The errors in all periods are estimated to be less than
0.5%.
9.2.3.2

PC Event Selection

By comparing the cut variables for PC event selection between the observed data and MC, theuncertainties of the selection eﬃciencies are evaluated to 2.4%, 4.8%, 0.5% and 1.0% for SK-I to
SK-IV, respectively. The uncertainty of the background contamination rate in the final sample
is changed by the period. They are evaluated to 0.2%, 0.7%, 1.8% and 0.49% for SK-I to SK-IV.
9.2.3.3

FC/PC Separation

As shown in Sec. 7.1.6 and Sec. 7.2.6, FC and PC events are separated based on NHITAC. Basedon the diﬀerence of NHITAC between the observed data and MC, the uncertainty of the FC/PC
separation eﬃciency is evaluated to 0.6%, 0.5%, 0.9%, and 0.02% for SK-I to SK-IV.
9.2.3.4

UPMU Event Selection

The systematic errors of the upward-going muon reduction eﬃciencies are evaluated by comparing the cut variables for UPMU reduction between the observed data and MC. These errors
are diﬀerent between the stopping and through-going muon. The errors for the stopping muons
are evaluated to 0.7%, 0.7%, 0.7%, and 0.5% for SK-I to SK-IV. On the other hand, the errors
for through-going muons are evaluated to 0.5%, 0.5%, 0.5%, and 0.3% for SK-I to SK-IV. The
uncertainties of the separation ratios of stopping muons to through-going muons and the ratioof showering muons to non-showering muons are also evaluated by comparing the data and MC.
The uncertainties of the former ratios are evaluated to 0.4%, 0.6%, 0.4%, and 0.6% for SK-I
to SK-IV, and those of the latter ratios are 3.4%, 4.4%, 2.4%, and 3.0% for SK-I to SK-IV. In
addition, the uncertainties of the number of the subtracted background events from the final
samples are also introduced.
9.2.3.5

π 0 Like Sample Selection

This uncertainties of 1-ring and 2-ring π 0 like sample selection eﬃciencies are estimated bycomparing the data with MC. The uncertainty of 1-ring selection eﬃciency has the momentum
dependence, and it was estimated for each momentum region. The uncertainty of two ring π 0like samples selection is estimated based on the uncertainty of the ring-counting and particle
identification.

CHAPTER 9. OSCILLATION ANALYSIS
9.2.3.6

136

Errors on Decay Electron Tagging

The number of tagged decay electrons depends on the π + decay probability and decay electron
tagging eﬃciency. The uncertainty of π + decay probability is estimated to 10% based on thediﬀerence between data and MC. This is estimated based on the uncertainty of π + − O crosssection, which is evaluated by comparing the cross section based on the data and MC. The
systematic error in the tagging eﬃciency is estimated to 1.5% for SK-I to III and 0.8% for
SK-IV based on the data of stopping muons produced by the cosmic rays.
9.2.3.7

Energy Scale

As mentioned in Sec. 4.3.4, the systematic errors for the energy scale of charged leptons isevaluated by using the data of time variation and the diﬀerence between MC and data, and theerrors are 3.26%, 2.75%, 2.43%, and 2.13% for SK-I to SK-IV. As described in the next section,
the systematic errors on the energy are introduced
9.2.3.8

Errors on Energy Correction Based on Secondary Interactions

As mentioned in Sec. 8, the uncertainties of the corrected energy based on the secondary interactions depend on visible energy. Therefore, the systematic errors of the corrected energy are
introduced as a function of visible energy. The maximum shifts which are the black lines inFig. 8.12, Fig. 8.13, Fig. 8.17, Fig. 8.19, are introduced as the systematic errors. By shifting thecorrected energy based on the systematic errors and comparing the energy distributions withand without this shift, the eﬀects to the oscillation analysis are evaluated.
9.2.3.9

Errors on Event Reconstruction

The systematic errors on the event reconstruction, which are ring separation and particle identification, are also considered. They are estimated based on the diﬀerence between the observed
data and MC. They depend on the visible energy, the number of ring, and event type. The
systematic errors of particle identification and ring counting distribute from 0.2% to 9.5% andfrom 0.4% to 8.2%, respectively.
9.2.3.10

Multi-GeV Multi-Ring νe and ν¯e Separation

As described in Sec. 9.1.3, the Multi-GeV Multi-Ring e-like samples are separated to νe and ν̄elike samples. The variables for the νe and ν¯e separation depend on the cross section of the final
state interactions of pions and pion multiplicity. These uncertainties are estimated based on
the diﬀerence between MC and experimental data[145][146]. Based on these uncertainties, thesystematic error for the number of νe and ν̄e in Multi-Ring Multi-GeV νe like events is estimated
to 7.16%, 7.91%, 7.68%, and 6.82% for SK-I to SK-IV.

9.3

Sensitivity Study

The oscillation analysis evaluates χ2 values for the normal and inverted hierarchies, separately.
The mass hierarchy is identified by selecting the hierarchy with the smaller χ2 value. If their χ2values have large diﬀerence, the mass hierarchy can be identified with high reliability. Therefore,
the sensitivity study uses the diﬀerence of χ2 between the hierarchies. The sensitivity study is

CHAPTER 9. OSCILLATION ANALYSIS

137

done based the live time until 2017 April (SK-I: 1489.2 days, SK-II: 798.6 days, SK-III: 518.1
days, SK-IV: 2867.1 days). The sensitivity study assumes the distributions based on ”true”
parameters as the observed distribution. The ”true” parameters are summarized in Table 9.2.
In this table, the parameters are separated to ”Fitted” and ”Fixed” parameters. The ”Fixed”
parameters of the candidates to compare are set to the true parameters. The ”Fitted” parameters
can be changed from the true parameters. These parameters are set to the same as the previous
study[2].
Parameters|∆m223 |[eV2 ]sin2 θ23
δCP
sin2 θ13∆m212 [eV2 ]sin2 θ12

True Parameters
2.5 × 10−3
0.575
4.19
2.19 × 10−2
7.65 × 10−5
0.309

Fixed/Fitted
FittedFittedFitted
Fixed
Fixed
Fixed

Table 9.2: The table of ”true” parameters and Fixed/Fitted types for the oscillation sensitivity
study.
Fig. 9.4 and Fig. 9.5 show the distributions of ∆χ2 (= χ2 −χ2min ) as a function of ∆m223 , sin2 θ23
and δCP when the true hierarchies are the normal and inverted hierarchies, respectively. χ2minis the minimum χ value of all the normal and inverted hierarchies candidates. The diﬀerence ofthe minimum χ2 values of the normal and inverted hierarchies candidates (∆χ2min ) are definedas the sensitivity for the mass hierarchy. The sensitivities for each condition are summarized inTable 9.3. As shown in the table, the new method improves the neutrino oscillation sensitivity by

Current Method
New Method

True Normal Hierarchy
3.04
4.05

True Inverted Hierarchy
1.68
2.05

Table 9.3: Oscillation sensitivities (in χ2 diﬀerence) when the true oscillation parameters follows
Table 9.2
33% and 22% for NH and IH, respectively. This improvement comes from the small bin width andbetter energy resolution. Fig. 9.6 shows the comparison of the reconstructed energy spectrum of
upward-going (−1 < cos Θ < −0.9) 1-ring Multi-GeV elike events between normal and inverted
hierarchies based on the current and new energy binning. Their oscillation parameters are set
to the true parameters in Table 9.2. These two spectra are compared based on χ2 test whichis shown in Eq. (9.3). En and On in Eq. (9.3) are assigned to the number of events in n-thenergy bin of the energy spectrum in Fig. 9.6 for normal and inverted hierarchies, respectively.
As shown in Fig. 9.7, the narrow binning and better energy resolution enables to detects theenergy region where the diﬀerence between the hierarchies are significantly shown. It enables toimprove the mass hierarchy sensitivity.
This new analysis method also improved the precision of mixing parameters. The summaryof 68% and 90% C.L. regions are shown in Table 9.4, Table 9.5, and Table 9.6.

∆χ 2

10

138

∆χ 2

∆χ 2

CHAPTER 9. OSCILLATION ANALYSIS

10

9

9

8

8

7

7

6

6

5

5

4

4

3

3

2

2

1

1

0
0.001

0

0.00150.002

0.00250.003

0.0035

0.004

4

5

0.0045|∆ m232 | [MeV2 ]

0.3

0.35

0.4

0.45

0.5

0.55

0.6

0.65

0.72
sin θ23

10
9
8
7
6
5
4
3
2
1
0

0

1

2

3

6

δCP

Figure 9.4: ∆χ2 (≡ χ2 − χ2min ) distributions as a function of ∆m223 , sin θ23 and δCP whose truehierarchy is the normal hierarchy. The red and blue lines are assigned to the sensitivities of thecurrent and new analysis method, respectively. The solid and broken lines are assigned to thesensitivities of the normal and inverted hierarchy candidates, respectively.

10

10

9

9

8

8

7

7

6

6

5

5

4

4

3

3

2

2

1

1

0
0.001

∆χ 2

139

∆χ 2

∆χ 2

CHAPTER 9. OSCILLATION ANALYSIS

0.00150.002

0.00250.003

0.0035

0.004

4

5

0.0045|∆ m232 | [MeV2 ]

0
0.3

0.35

0.4

0.45

0.5

0.55

0.6

0.65

0.72
sin θ23

10
9
8
7
6
5
4
3
2
1
0

0

1

2

3

6

δCP

Figure 9.5: ∆χ2 (≡ χ2 − χ2min ) distributions as a function of ∆m223 , sin θ23 and δCP whose truehierarchy is the inverted hierarchy. The red and blue lines are assigned to the sensitivities of
the current and new analysis method, respectively. The solid and broken lines are assigned tothe sensitivities of the normal and inverted hierarchy candidates, respectively.

CHAPTER 9. OSCILLATION ANALYSIS

140

Figure 9.6: Reconstructed Energy spectrum of upward-going (−1 < cos Θ < −0.9) 1-ring MultiGeV νe like samples. The left figure shows the spectrum based on the current binning setting
and reconstruction method, and the right figure shows the spectrum based on the new binningand method.

∆m232 (×10−3 eV)NH
IH

Current Method
New Method
Current Method
New Method

68% C.L.
2.31 < ∆m232 < 2.72
2.35 < ∆m232 < 2.68
2.29 < ∆m232 < 2.73
2.34 < ∆m232 < 2.69

90% C.L.
2.18 < ∆m232 < 2.88
2.24 < ∆m232 < 2.82
2.16 < ∆m232 < 2.91
2.23 < ∆m232 < 2.84

Table 9.4: Allowed region of ∆m232 at 68% and 90% C.L., where 68% (90%) C.L. is defined as
χ2 − χ2min = 1.0(2.7).

sin2 θ23
NH
IH

Current Method
New Method
Current Method
New Method

68% C.L.
0.525 < sin2 θ23 < 0.617
0.531 < sin2 θ23 < 0.616
0.477 < sin2 θ23 < 0.616
0.491 < sin2 θ23 < 0.616

90% C.L.
0.468 < sin2 θ23 < 0.633
0.478 < sin2 θ23 < 0.632
0.430 < sin2 θ23 < 0.632
0.433 < sin2 θ23 < 0.632

Table 9.5: Allowed region of sin2 θ23 at 68% and 90% C.L., where 68% (90%) C.L. is defined as
χ2 − χ2min = 1.0(2.7).

CHAPTER 9. OSCILLATION ANALYSIS

141

Figure 9.7: χ2 distribution of upward-going(−1 < cos Θ < −0.9) 1-ring Multi-GeV νe likesamples. χ2 in this figure is calculated by setting En and On following to the parameters in
Table 9.2 with normal and inverted hierarchies, respectively. The blue and red lines are assigned
to the distributions based on current and new reconstruction method and binning, respectively

δCP
NH
IH

Current Method
New Method
Current Method
New Method

68% C.L.
1.95 < δCP < 5.59
1.99 < δCP < 5.57
1.77 < δCP < 5.80
1.79 < δCP < 5.79

90% C.L.
All values
All values
All values
All values

Table 9.6: Allowed region of δCP at 68% and 90% C.L., where 68% (90%) C.L. is defined as
χ2 − χ2min = 1.0(2.7).

Chapter 10

Conclusion and Outlook
New energy reconstruction algorithm has been developed for neutrino mass hierarchy determination. This algorithm corrects the missing neutrino energy taken by the produced hadrons,
and improves the reconstructed energy resolution. Until SK-III, in which any neutrons could
not be detected, the energy resolution is improved in the range from 4% to 26%. For SK-IV, inwhich neutrons are tagged with 21.7% eﬃciency, the resolution is improved in the range from
12% to 39%. The degree of improvement depends on the event type. Better energy resolution
will enable to detect the resonance enhancement of atmospheric neutrino oscillation probability
more obviously. To realize this, some modifications of the oscillation analysis were required. Theenergy binning was changed and the resonance energy region (103.4 − 104.1 MeV) was separated
by 8 bins. In addition, the new systematic errors related to this energy correction, which are
the shifts of the corrected energy from the uncertainties of the interaction model and neutron
tagging eﬃciency and the shift from the neutrons produced by π − capture are introduced. Theseerrors depend on the neutrino energy, and it was evaluated to 5% at maximum in the resonanceenergy region. Based on this new analysis method, the mass hierarchy sensitivity was estimatedby using the live time until 2017 April (SK-I: 1489.2 days, SK-II: 798.6 days, SK-III: 518.1
days, SK-IV: 2867.1 days). As a result, the sensitivity is improved from ∆χ2 = 3.04 (1.68) to∆χ2 = 4.05 (2.05) in the true NH (IH) case.
New method shows large improvement of the mass hierarchy sensitivity. However, thissensitivity is not suﬃcient to identify the mass hierarchy with 5σ(∆χ2 = 25). In parallel with thisstudy, various studies for improving the sensitivity have been developed by some collaborators.
For example, the event reconstruction tool was changed from APFit to new tool, which is calledfiTQun. fiTQun enabled to expand the fiducial volume, and the inverted hierarchy has been
disfavored at 74%[3]. For another example, the identification method of neutrino and antineutrino based on the number of tagged neutrons was also developed.[1]. As a next step, we are
required to combine these development and introduce into the oscillation study. Currently, theworks for the introduction of Gadolinium into SK water tank is in progress. Because neutron
capture on Gadolinium produces γ-rays with higher energy and its capture time is much shorter
than the capture on hydrogen, the neutron tagging eﬃciency is expected to be improved. If theneutron tagging eﬃciency becomes higher, the uncertainty of the number of produced neutrons
as a function of the number of tagged neutrons will be smaller. Because the missing energy
largely depends on the number of produced neutrons, this improvement will make the energy
resolution better, and the sensitivity to the mass hierarchy will be improved.

142

Bibliography[1] T. J. Irvine, PhD thesis, University of Tokyo (2014).
http://www-sk.icrr.u-tokyo.ac.jp/sk/_pdf/articles/TristanJamesIrvine_
Thesis.pdf[2] P. Fernandez, PhD Thesis, University Autonomous of Madrid (2017)
http://www-sk.icrr.u-tokyo.ac.jp/doc/sk/_pdf/articles/2017/thesis_
PabloFernandez_March2017.pdf.
[3] M. Jiang et al., arXiv:1901.03230[hep-ex] (2019).
[4] M. Fukugita and T. Yanagida, Phys. Lett. B, 174, 45-47 (1986).
[5] LEP, ALEPH, DELPHI, L3, OPAL Collaborations and Line Shape Sub-Group of the LEP
Electroweak Working Group, Phys. Rept. 427,257-454 (2006).
[6] M. Nakagawa, Z. Maki, and S. Sakata, Prog. Theor. Phys. 28:870-880, (1962).
[7] Y. Fukuda et al., Phys. Rev. Lett. 81, 1662 (1998).
[8] Q.R. Ahmad et al., Phys. Rev. Lett. 87, 071301 (2001).
[9] Q.R. Ahmad et al., Phys. Rev. Lett. 89, 011301 (2002).
[10] Q.R. Ahmad et al., Phys. Rev. Lett. 92, 181301 (2004).
[11] L. Wolfenstein, Phys. Rev. D17 2369-2374 (1978).
[12] S. P. Mikheev and A. Yu. Smirnov, Sov. J. Nucl. Phys., 42:913-917 (1985).
[13] J.N. Bahcall, Phys. Rev. Lett., 12, 300 (1964)[14] K. S. Hirata et al., Phys. Rev. Lett. 63, 16 (1989).
[15] P. Anselmann et al., Phys. Lett. B, 285, 376-389 (1992).
[16] A.I. Abazov et al., Phys. Rev. Lett. 67, 3332 (1991).
[17] S. Fukuda et al., Phys. Rev. Lett. 86, 5651 (2001).
[18] J.N. Bahcall, M. H. Pinsonneault, Sarbani Basu, Astrophys. J. 555, 990 (2001).
[19] K. Hirata et al., Phys. Lett. B, 280, 146-152 (1992).
143

BIBLIOGRAPHY

144[20] Y. Fukuda et al., Phys. Lett. B, 335, 237-245 (1994).
[21] Y. Fukuda et al., Phys. Rev. Lett. 81, 1562 (1998).
[22] Z. Li et al., Phys. Rev. D98, 052006 (2018).
[23] C. E. Ortiz et al., Phys. Rev. Lett. 85, 2909 (2000).
[24] K. Eguchi et al., Phys. Rev. Lett. 90, 021802 (2003).
[25] M. Tanabashi et al., Phys. Rev. D 98, 030001 (2018).
[26] D.G. Michael et al., Phys. Rev. Lett. 97, 191801 (2006).
[27] M.H.Ahn et al., Phys. Rev. D 74, 072003 (2006).
[28] P. Adamson et al., Phys. Rev. Lett. 107, 021801 (2011).
[29] Y. Abe et al., Phys. Rev. D 86, 052008 (2012).
[30] F. P. An et al., Chinese Phys. C, 37, 011001 (2013).
[31] J. K. Ahn et al., Phys. Rev. Lett. 108, 191802 (2012).
[32] K. Abe et al., Phys. Rev. Lett. 112, 061802 (2014).
[33] K. H. Ackermann et al., Eur. Phys. J. C 73, 2330 (2013).
[34] A.Gando et al., Phys. Rev. Lett. 117, 082503 (2016).
[35] K. Abe et al. Phys. Rev. D 97, 072001 (2018).
[36] C. Alduino et al., Phys. Rev. Lett. 120, 132501 (2018).
[37] K. Hirata et al., Phys. Rev. Lett. 58, 1490 (1987).
[38] Y. Fukuda et al., Phys. Lett. B 335, 237-245 (1994).
[39] I. Frank and I. Tamm, C.R. Acad. Sci. USSR 14, 109 (1937).
[40] K. Fujita et al., Nucl. Instrum. Meth. A 501, 418-462 (2003).
[41] Y. Koshio, PhD thesis, University of Tokyo (1998).
http://www-sk.icrr.u-tokyo.ac.jp/sk/_pdf/articles/koshio-dthesis.pdf[42] Y. Nakano, PhD thesis, University of Tokyo (2016).
http://www-sk.icrr.u-tokyo.ac.jp/sk/_pdf/articles/2016/doc_thesis_naknao.
pdf[43] R. Becker-Szendy et al., Nucl. Instrum. Meth. A 324, 363-382 (1993).
[44] K. Abe et al., Nucl. Instrum. Meth. A 737, 253-272 (2014).
[45] C. Ishihara, PhD Thesis, University of Tokyo (2010).
http://www-sk.icrr.u-tokyo.ac.jp/sk/_pdf/articles/ishihara-phdthesis.pdf

BIBLIOGRAPHY

145[46] Nucl. Instrum. Meth. A 501, 418-462 (2003).
[47] E. Pope and R. Fry, Appl. Opt. 33 (1997) 8710-8723.
[48] V. S. Langford et al., J. Phys. Chem. A 105 (2001) 8916-8921.
[49] N. Sakurai, PhD thesis, University of Tokyo (2000).
http://www-sk.icrr.u-tokyo.ac.jp/sk/_pdf/articles/sakurai.pdf[50] M. Honda et al., Phys. Rev. D 83, 123001 (2011).
[51] J. Alcaraz et al., Phys. Lett. B 490, 27 (2000).
[52] J. Alcaraz et al., Phys. Lett. B 494, 193 (2000).
[53] S. Haino et al., Phys. Lett. B 594, 35 (2004).
[54] T. Sanuki et al., Astrophys. J. 545, 1135 (2000).
[55] T. K. Gaisser et al., Proc. of the 27th Int. Cosmic Ray Conf., Hamburg 5, 1643 (2001).
[56] T. K. Gaisser and M. Honda, Ann. Rev. Nucl. Part. Sci. 52, 153 (2002).
[57] H. Collaboration, Astropart. Phys. 30, 124 (2008).
[58] International Geomagnetic Reference Field https://www.ngdc.noaa.gov/IAGA/vmod/
igrf.html.
[59] United States Committee on Extension to the standard Atmosphere, U.S.Standard Atmosphere 1976.
https://ccmc.gsfc.nasa.gov/modelweb/atmos/us_standard.html[60] S.Roesler, R. Engel, and J. Ranft, arXiv:hep-ph/0012252v1 (2000).
[61] A. Capella et al., Phys. Rept. 236, 225 (1994).
[62] K. Niita et al., Radiat. Meas. 41, 1080 (2006).
[63] T. Sanuki, M. Honda, T. Kajita, K. Kasahara, and S. Midorikawa, Phys. Rev. D75, 043005(2007).
[64] S. Haino et al., Phys. Lett. B 594, 35 (2004).
[65] T. Sanuki et al.,Phys. Lett. B 541, 234 (2002).
[66] K. Abe et al.,Phys. Lett. B 564, 8 (2003).
[67] G. D. Barr et al., Phys. Rev. D 70, 023006 (2004).
[68] G. Battistoni, et al., Astropart. Phys. 19, 269 (2003).
[69] M. Nakahata et al., J. Phys. Soc. Jpn. 55, 3786 (1986).
[70] C.H. Llewellyn-Smith, Phys. Rep. C 3, 261 (1972).

BIBLIOGRAPHY[71] R.A. Smith, E.J. Moniz, Nucl. Phys. B 43, 605 (1972); B101, 547(E) (1975).
[72] R. Gran et al., Phys. Rev. D 74, 052002 (2006).
[73] A.A. Aguilar-Arevalo, Phys. Rev. D81, 092005 (2008).
[74] Y. Hayato, Acta Phys. Polon. B 40, 2477 (2009).
[75] D. Rein, L. M. Sehgal, Ann. Phys. 133, 79 (1981).
[76] K.M. Graczyk and J.T. Sobczyk, Phys.Rev. D 77,053001 (2008).
[77] C.H. Albright and C. Jarlskog, Nucl. Phys. B 84, 467 (1975).
[78] M. Guck et al. Eur. Phys. J C, 5:461, (1998).
[79] A. Bodek, U.K. Yang, arXiv:hep-ex/0308007 (2003).
[80] P. Musset, J.-P. Viale, Phys. Rep. 39, 1 (1978).
[81] J.E. Kim et al., Rev. Mod. Phys. 53, 211 (1981).
[82] S. J. Barish et al., Phys. Rev. D 17,1 (1978).
[83] T. Sjostrand et al. CERN-TH-7112-93 (1994).
[84] K.M. Graczyk and J.T. Sobczyk, Phys.Rev. D 77,053001 (2008).
[85] D. Rein, L.M. Sehgal, Nucl. Phys. B 223, 29 (1983).
[86] M. Hasegawa et al., Phys. Rev. C 95, 252301 (2005).
[87] K. Hiraide et al., Phys. Rev. D 78, 112004 (2008).
[88] J.L. Raaf, Ph.D. thesis, University of Cincinnati (2005).
[89] H. Faissner et al., Phys. Lett. B 125, 230 (1983).
[90] E. Isiksal, D. Rein, J.G. Morfin, Phys. Rev. Lett. 52, 1096 (1984).
[91] A. Kartavtsev et al., Phys. Rev. D 74, 054007 (2006).
[92] L. Alvarez-Ruso et al., Phys. Rev. C 75, 055501 (2007).
[93] J. Nieves et al., Phys. Rev. C 83, 045501 (2011).
[94] R. Woods, D. Saxon, Phys. Rev. 95, 577 (1954).
[95] L.L. Salcedo, et al., Nucl. Phys. A 484, 557 (1988).
[96] P. de Perio, PhD thesis, University of Toronto (2014).
https://tspace.library.utoronto.ca/handle/1807/68591[97] G. Rowe et al., Phys. Rev. C 18, 584 (1978).
[98] R. Seki, K. Masutani, Phys. Rev. C 27, 2799 (1983).

146

BIBLIOGRAPHY

147[99] C. Oh et al., Phys. Rev. C 56, 635 (1997).
[100] R. Arndt et al., Phys. Rev. C 69, 035213 (2004).
[101] R. Workman et al., SAID Partial Wave Analysis Program.
http://gwdac.phys.gwu.edu.
[102] B.R. Martin, M.K. Pidcock, Nucl. Phys. B 126, 266 (1977).
[103] B.R. Martin, M.K. Pidcock, Nucl. Phys. B 126, 285 (1977).
[104] J. S. Hyslop et al., Phys. Rev. D 46, 961 (1992).
[105] D.A. Sparrow, Proc. of the Conf. on the Intersection Between Particle and Nuclear Physics,
1019 (1984).
[106] GEANT, CERN Program Library Long Writeup W5013 (1994).
[107] C. Zeitnitz, T.A. Gabriel. Nucl. Instrum. Meth., A 349, 106-111 (1994).
[108] H.W. Bertini, Phys. Rev. C6, 631 (1972).
[109] J. Nieves, J.E. Amaro, M. Valverde, Phys. Rev. C 70, 055503 (2004).c[110] A. Bodek, U.K. Yang, arXiv:hep-ph/0508007 (2005).
[111] C. Bronner, talk at NuSTEC Workshop on Shallow- and Deep-Inelastic Scatteringhttps://indico.cern.ch/event/727283/contributions/3102160/[112] Ch. Berger, L.M. Sehgal, Phys. Rev. D 79, 053003 (2009).
[113] Y. Hayato, NuInt18 talkhttps://indico.cern.ch/event/703880/contributions/3157411/[114] E. S. Pinson Guerra, et al., Phys. Rev. C 95, 045203 (2017).
[115] K. Ieki, et al., Phys. Rev. C92, 035205 (2015).
[116] D. Ashery et al., Phys. Rev. C, 30(3):946 (1984).
[117] M. K. Jones et al., Phys. Rev. C, 48(6):2800 (1993).
[118] E. Bellotti, et al., II Nuovo Cimento A (1965-1970) 18, 75 (1973).
[119] E. Bellotti, et al., II Nuovo Cimento A (1965-1970) 14, 567 (1973).
[120] C. Andreopoulos et al., Nucl. Instrum. Meth. A 614, 87 (2010).
[121] Y. Takenaga, PhD thesis, University of Tokyo (2008).
http://www-sk.icrr.u-tokyo.ac.jp/sk/_pdf/articles/Takenaga-d.pdf[122] F. Dufour PhD thesis, Boston University (2009).
http://www-sk.icrr.u-tokyo.ac.jp/sk/_pdf/articles/fdufour-thesis.pdf[123] D. E.Groom, N. V. Mokhov, S. I. Striganov, Volume 78, Issue 2, July 2001, Pages 183-356.

BIBLIOGRAPHY

148[124] Z. Li, PhD Thesis, Duke University (2017).
http://www-sk.icrr.u-tokyo.ac.jp/sk/_pdf/articles/2018/zepeng-dthesis.pdf[125] M. R. Dziomba, PhD Thesis, University of Washington , Aug. 2012.
http://www-sk.icrr.u-tokyo.ac.jp/sk/_pdf/articles/Michael R. Dziomba_
thesis Aug. 2012.pdf[126] S. Desai, PhD thesis, Boston University (2004).
http://www-sk.icrr.u-tokyo.ac.jp/sk/_pdf/articles/desai.pdf[127] K. Nitta, PhD thesis, Osaka University (2003).
http://www-sk.icrr.u-tokyo.ac.jp/sk/_pdf/articles/nitta.pdf[128] J. Kameda, PhD thesis, University of Tokyo (2002).
http://www-sk.icrr.u-tokyo.ac.jp/sk/_pdf/articles/kameda_d.pdf[129] V. F. Sears, Neutron News 3, 26 (1992).
[130] D. Cokinos and E. Melkonian, Phys. Rev. C 15, 1636 (1977).
[131] C. Ingram et al., Phys. Rev. C 27, 1578 (1983).
[132] C. Patrignani et al.,(Particle Data Group), Chin. Phys. C 40, 100001 (2016).
[133] C. Andreopoulos et al., Nucl. Instrum. Meth. A614 87-104 (2010).
[134] C. Andreopoulos et al., arXiv:1510.05494 [hep-ph] (2015).
[135] J. Nieves and J. E. Sobczyk, Annals of Physics 383 (2017) 455-496.
[136] Lee Ka Pik, PhD Thesis, University of Tokyo (2012).
http://www-sk.icrr.u-tokyo.ac.jp/sk/_pdf/articles/Lee Ka Pik_Thesis_print_
v10.pdf[137] M. Honda et al., Phys. Rev. D 75, 043006 (2007).
[138] G. Ambrosini et al., Phys. Lett. B 420, 225 (1998).
[139] A.E. Hedin, J. Geophys. Res., 96, 11591172 (1991).
[140] A. M. Dziewonski, D. L. Anderson. Phys. Earth Planet. Inter., 25, 297-356 (1981).
[141] E. Hernandez et al., Phys. Rev. D 76, 033005 (2007).
[142] A. Capella et al., Phys. Lett. B, 337, 358 (1994).
[143] Y. Nakajima et al., Phys. Rev. D 83, 012005 (2011).
[144] K. Hagiwara et al., Nucl. Phys. B 668, 364 (2003).
[145] P. de Perio, AIP Conf. Proc. 1405, pp. 223-228 (2011).
[146] The CHORUS Collaboration, Eur. Phys. J. C 51, 775-785 (2007).

Acknowledgement
This thesis was completed owing to many people contribution. First, I am deeply grateful tomy supervisor, Masahiro Kuze. He had given me the opportunity to study neutrino physics.
He had devoted himself to realize my hope to study in Super-Kamiokande. Even after my leftof university, he had encouraged me to write this thesis. Without his contribution, this thesis
would not be completed.
I am grateful to Masato Shiozawa for his support. He had kindly accepted the researchguidance in Kamioka for me. He regarded me as a student of his oﬃce, and he gave me a lot ofadvice for my study.
I would like to thank to all Super-Kamiokande collaborators for their support. Especially, I
would like to thank to Roger Wendell, Pablo Fernandez, Christophe Bronner, Ryosuke Akutsu,
and Yutaro Sonoda. The crucial element of this thesis, the new energy reconstruction algorithm
and oscillation analysis method based on the new reconstructed energy, were developed under
their support. I also appreciate to Shigetaka Moriyama, Jun Kameda, Shinichi Mine, and all theother SK calibration group members. They advised me how to introduce the new calibration
trigger. Especially, Mine-san kindly told me many things, for example, how the SK detector iscalibrated. I would also like to thank to Masaki Ishitsuka, Yoshinari Hayato, Makoto Miura,
Hidekazu Tanaka, Shoei Nakayama, Yusuke Koshio, Motoyasu Ikeda, Takaomi Yano. They told
many things related to my study. Especially, Ishitsuka-san had supported me as an assistantprofessor in Kuze laboratory until 2016. He told me the basic things of the study.
I also appreciate to Hyper-Kamiokande collaborators. I had developed new MC simulator
and PMT for it in my master course. Owing to this experience, I understood how to simulate
the neutrino event and what is required to improve the detector performance. Especially, I
appreciate to Yasuhiro Nishimura, Seiko Hirota, Yusuke Suda. I had worked with them asa member of the photo-detector working group. They supported both my study and life in
Kamioka.
I also thank to Yuki Nakano, Masatoshi Kobayashi, Asato Orii, Jiang Miao, Daisuke Fukuda,
Akira Takenaka, Takumi Suzuki, Fumimasa Muto, and Tianmeng Lou and all the other staﬀs
and students who studied together in Kamioka. Owing to their kind support, my study in
Kamioka was very fun.
I appreciate for all of the current and previous members of Kuze laboratory, Yohei Yamaguchi, Michiru Kaneda, Tomoyuki Konno, Takayuki Kanno, Yosuke Abe, Takuya Nobe,
Ayumu Ito, Shinnosuke Yokoyama, Dai Kobayashi, Kota Higuchi, Ralitsa Sharankova, Masahiro
Tanaka, Kento Hosokawa, Yu Nishihara, Kohei Shimizu, Tomoyo Yoshida, Dai Hamabe, Satoshi
Kawaguchi, Yukitaka Yazawa, Lukas Berns, Tatsuyuki Sekine, and Yunjian He.
Finally, I deeply appreciate for my family. They encouraged me to do everything what I
want to do. Without their support, my study will not be completed.
149