c Copyright 2012
Michael R. Dziomba

A Study of Neutrino Oscillation Models with Super-Kamiokande
Atmospheric Neutrino Data

Michael R. Dziomba

A dissertation
submitted in partial fulfillment of therequirements for the degree ofDoctor of Philosophy

University of Washington
2012

Reading Committee:
R. Jeffrey Wilkes, Chair
Ann E. Nelson
Thompson Burnett

Program Authorized to Offer Degree:
Physics

University of Washington
Abstract
A Study of Neutrino Oscillation Models with Super-Kamiokande Atmospheric Neutrino
Data
Michael R. Dziomba
Chair of the Supervisory Committee:
Professor R. Jeffrey Wilkes
Physics

The rapidly evolving field of neutrino oscillations continues to offer exciting experimental
results as it enters the era of precision measurements. The mixing angles and mass-squared
splittings governing neutrino oscillations have all been measured experimentally and only
one parameter, the CP-violating phase δ, remains unmeasured. Both atmospheric and longbaseline accelerator experiments are sensitive to the parameters θ23 and ∆m223 . In particular
the wide range of neutrino pathlengths and energies accessible to the Super-Kamiokande(SK) detector with atmospheric neutrinos offer a unique opportunity to measure theseparameters with high-statistics data. SK is also suitable for constraining exotic oscillation
models which include sterile neutrinos and non-standard matter interactions, due to thelarge variations in density through which atmospheric neutrinos can traverse. In this thesis
we use an expanded atmospheric neutrino dataset combined with improved Monte Carlosimulations to update the SK 2-flavor muon neutrino disappearance zenith angle analysis.
−3 eV2 and sin2 2θ
We measure values of ∆m223 = 2.30+0.17
23 > 0.95 (90% C.L.), in−0.22 × 10

excellent agreement with recent worldwide results. We also fit the SK data to a recent exotic
model which incorporates both sterile-active mixing and a new non-standard interaction inmatter. We find that SK rules out most of the model’s parameter space and restricts it to
the region where standard 2-flavor mixing is recovered.

TABLE OF CONTENTS

Page
List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv
Glossary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi
Chapter 1:

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1

1.1

Neutrinos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1

1.2

Atmospheric Neutrinos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2

1.3

Neutrino Oscillations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3

1.4

Sterile Neutrinos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

1.5

Oscillations in Matter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

9

1.6

Summary Of This Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

Chapter 2:

The Super-Kamiokande Detector . . . . . . . . . . . . . . . . . . . . . 11

2.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

2.2

Cerenkov Radiation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

2.3

The SK Detector . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

2.4

PMTs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

2.5

Data Acquisition System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

2.6

Air Purification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

2.7

Water Purification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

Chapter 3:

Calibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

3.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

3.2

Water Calibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

3.3

Charge calibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

3.4

Timing calibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

3.5

Energy scale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44

3.6

OD calibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

3.7

Supernova burst test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52i

3.8

GPS Timing System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

Chapter 4:

Atmospheric Neutrino Monte Carlo . . . . . . . . . . . . . . . . . . . . 55

4.1

Atmospheric neutrino flux . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55

4.2

Neutrino interaction simulation . . . . . . . . . . . . . . . . . . . . . . . . . . 59

4.3

Detector Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

Chapter 5:

Atmospheric Neutrino Data Reduction . . . . . . . . . . . . . . . . . . 71

5.1

Fully Contained Reduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72

5.2

Partially Contained Reduction . . . . . . . . . . . . . . . . . . . . . . . . . . 82

5.3

Upward-going Muon Reduction . . . . . . . . . . . . . . . . . . . . . . . . . . 88

Chapter 6:

Event Reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107

6.1

Vertex Fitting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107

6.2

Ring Counting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110

6.3

Precise Vertex Fitting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116

6.4

Momentum Determination . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120

6.5Reconstruction Improvements . . . . . . . . . . . . . . . . . . . . . . . . . . . 120

Chapter 7:

Two-Flavor Oscillation Analysis . . . . . . . . . . . . . . . . . . . . . . 122

7.1

Analysis Space . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123

7.2

Analysis Binning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123

7.3

Monte Carlo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126

7.4

Dataset . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126

7.5

Chi-Squared Analysis

7.6

Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139

Chapter 8:. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129

Apparent CPT Violation Model . . . . . . . . . . . . . . . . . . . . . . 147

8.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147

8.2

Mathematics of Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148

8.3

Fit to MINOS Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150

8.4

Effects on Atmospheric Neutrinos at SK . . . . . . . . . . . . . . . . . . . . . 150

8.5

Analysis Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155

8.6

Analysis Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156

Chapter 9:Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 164ii

Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 165

iii

LIST OF FIGURES

Figure Number
1.1

Page

The LSND allowed region[69] overlayed with the exclusion regions from KARMEN[68]and MiniBooNE[70] for short baseline oscillations with large ∆m2 . Imagetaken from [70]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

2.1
2.2

Location of the Super-Kamiokande experiment in Japan. Taken from Ref. [12] 12
An overview of the SK detector and its location within Mt. Ikenoyama.
Taken from Ref. [12] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.3 A schematic of the Cerenkov radiation emission angle relative to the charged
particle’s track. Taken from Ref. [112] . . . . . . . . . . . . . . . . . . . . . . 13
2.4 The scaffolding structure which divides the ID and OD, along with themounted PMTs. Taken from Ref. [12] . . . . . . . . . . . . . . . . . . . . . . 15
2.5 Schematic view of an ID 50 cm PMT. Taken from Ref. [12] . . . . . . . . . . 16
2.6 Quantum efficiency of the ID PMTs. The peak efficiency coincides with UVradiation wavelengths. Taken from Ref. [12] . . . . . . . . . . . . . . . . . . . 17
2.7 A schematic of the new DAQ system. Taken from Ref. [33] . . . . . . . . . . 20
2.8 A schematic of the new periodic software trigger. Taken from Ref. [33] . . . . 21
2.9 A schematic of the timing resolution as a function of input charge for a QBEE
and ID PMT. Taken from Ref. [33] . . . . . . . . . . . . . . . . . . . . . . . . 21
2.10 Timing chart for QTC processing of raw PMT signal. Taken from Ref. [33] . 22
2.11 The high-speed processing rate of the QBEEs. Frequency of the input signal
is on the horizontal axis. Taken from Ref. [33] . . . . . . . . . . . . . . . . . 22
2.12 Schematic view the radon-free air system. Taken from Ref. [5] . . . . . . . . 24
2.13 Schematic view the water purification system. . . . . . . . . . . . . . . . . . . 25
3.1
3.2
3.3
3.4
3.5

The index of refraction of the SK tank water modeled as a function of lightwavelength. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Angular distribution of scattered light from symmetric (left) and asymmetric(right) scattering. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
The positions of the eight laser light injectors in Super-Kamiokande. “Old”
refers to the sole original injection point used in SK-I. . . . . . . . . . . . . .
A typical laser light injection event (injected from the top of the tank). . . .
Measurement of the water parameters for scattering and absorption. . . . . .
iv

27
28
29
29
30

3.6

Cosmic ray muon water transparency measurement schematic. . . . . . . . . 31

3.7

Thru-mu measurement of attenuation length for SK-IV. . . . . . . . . . . . . 31

3.8

The installed locations of the standard PMTs (shown in red) within the ID. . 34

3.9

The scintillator ball setup for the PMT HV determination. A sample grouping of PMTs calibrated together is shown in blue. . . . . . . . . . . . . . . . . 34

3.10 The ratio of each PMTs corrected charge to the target charge (mean correctedcharge of standard PMTs) before (left) and after (right) HV adjustments. . . 35
3.11 The single p.e. distributions for SK-III (black) and SK-IV (red, adjusted byscaling factor) shown for two different horizontal scales (left and right). . . . 36
3.12 A schematic of the setup used to measure the relative PMT gains. . . . . . . 37
3.13 The distribution of relative gains for all ID PMTs after normalizing to mean 1. 38
3.14 The distribution of PMT hit rates for the barrel (top plot), top (bottom left)
and bottom (bottom right) PMTs. The blue is the calibration data and the
red the MC simulation data. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
3.15 The corrected charge divided by the QE and gain for PMTs divided geometrically. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
3.16 The nonlinearity as a function of expected charge for ID PMTs before (left)
and after (right) multi-hit correction. The red points show the average of theblack points at the corresponding horizontal charge bin. . . . . . . . . . . . . 42
3.17 The setup used for the timing calibration measurement. . . . . . . . . . . . . 43
3.18 A typical TQ distribution before correcting with TQ Map table. . . . . . . . 44
3.19 The data and MC energy distributions for stopping muon decay electron events. 45
3.20 The time variation of the decay electron momentum mean from SK-IV data.

46

3.21 The ratio of MC to data of the ratio of measured momentum to inferredmomentum from stopping muon opening angle. . . . . . . . . . . . . . . . . . 46
3.22 The ratio of momentum loss to range for stopping muons for both data andMC in initial SK-I calibration. . . . . . . . . . . . . . . . . . . . . . . . . . . 47
3.23 The time variation of the ratio of momentum to range for high-energy stopping muons for SK-IV. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
3.24 The invariant mass distribution of data and MC for 2-ring e-like events.
Taken from Ref. [12]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
3.25 A summary plot showing the data and MC agreement for all energy scalecalibration measurements as a function of energy. . . . . . . . . . . . . . . . . 49
3.26 Example 1-p.e. distributions for old (top) and new (bottom) OD PMTs
before (left) and after (right) SK-IV gain adjustment. . . . . . . . . . . . . . 50
3.27 Measured versus input charge for old (red) and new (black) OD PMTs fromlaser measurements. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51v

3.28 A photograph of the NIM modules used in the SN Burst test which the authorof this work arranged to simulate a SN signal. . . . . . . . . . . . . . . . . . . 52
3.29 A sample of the triggering profile during the SN Burst test taken during
SK-IV test installation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
4.1

The zenith angle dependence of the Honda 2011 flux for atmospheric neutrino
flux at Kamioka averaged over azimuthal angle, taken from [49]. . . . . . . . 56

4.2

The energy dependence of flux for each neutrino type (left) and important
flux ratios (right). HHKMS06 is the Honda 2006 flux[48], and Bartol[3] andFluka[4] are existing alternative flux calculations. Taken from [49]. . . . . . . 57

4.3

The ratios of the direction-averaged flux for Kamioka in the Honda 2011(“Thismodel”), Bartol[3] and Fluka[4] models relative to the Honda 2006 flux[48],Taken from [49]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

4.4

The νµ and ν¯µ CCQE cross-sections for NEUT (line) and various experimental measurements. See [83] for details. . . . . . . . . . . . . . . . . . . . . . . 60

4.5

The cross-sections for CC single-π production for (a) νµ p → µ− pπ + , (b)νµ n → µ− pπ 0 and (c) νµ n → µ− nπ + . The solid line shows the NEUT
calculation and the dots show the data points summarized in the bottom
right. Taken from [83]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62

4.6

The cross-sections from old and new pion interaction models in NEUT compared to data as a function of pion momentum. . . . . . . . . . . . . . . . . . 64

4.7

The QE differential cross-section for pion scattering on a thin-water target.
The new SGNUCEFF model agrees with scattering data much better thanthe previous SGPI model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66

4.8

Number of decay-electrons for a selected SK-4 event sample. The black isdata, the gray the MC used in previous SK analyses, and the blue and red
histograms are the predictions from the new MC with old and new pion
interaction cross-sections, respectively. . . . . . . . . . . . . . . . . . . . . . . 67

4.9

Distributions of SK-IV OD variables for cosmic ray data (black), MC with
old tuning (purple) and MC with the most recent tuning (red). . . . . . . . . 69

5.1

Typical neutrino (dashed) and resulting charged lepton (solid) tracks for the
fully contained (FC), partially contained (PC), upward-going stopping (UpStop) and upward through-going (Up-thru) event types, shown left to right. . 72

5.2

Location of the cable holes in the top of the tank. Red circles indicate holes
with veto counters and pink circles those without. Taken from [92]. . . . . . . 75

5.3

Typical flasher (i) and FC neutrino (ii) event timing distributions. Taken
from [87]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76

5.4

Typical NMIN100 distributions for data (a), MC (b) and post-FC3 data andnormalized MC (c) samples. Taken from [87]. . . . . . . . . . . . . . . . . . . 77
vi

5.5

A typical pair of flasher events caused by the same PMT separated by 5
minutes. Taken from [88]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79

5.6

The criterion for removing flasher events in FC4. Vertical axis shows correlation (r) and horizontal axis shows number of matching events. Taken from[88]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79

5.7

The removal of stopping muons in FC5. NHITAin is shown for (a) data, (b)MC and (c) final FC sample with normalized MC. Taken from [88]. . . . . . . 80

5.8

Number of PMT hits in top and bottom versus number in barrel for PC data.
Events above the red line are rejected. Taken from [19]. . . . . . . . . . . . . 84

5.9

Distributions of NHITAin for data (a) and MC (b) for the PC3 cut, and thefinal data and MC samples (c). Taken from [88]. . . . . . . . . . . . . . . . . 85

5.10 The effect of the PC4 cut which is based on the pfit direction and the direction
from the vertex to the earliest saturated PMT. Shown for SK-1 data(a),MC(b) and their final samples(c). Taken from [88]. . . . . . . . . . . . . . . . 86
5.11 Charged-current weak interaction which produces can produce an upmu.
Taken from [111]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
5.12 Neutrino (top) and muon (bottom) energies for the three upmu samples as
generated in SK-III MC (from [111]. . . . . . . . . . . . . . . . . . . . . . . . 91
5.13 A flow chart describing the umred1st selection and rejection steps. Taken
from [111]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
5.14 The dE/dx histogram for stopping muons used to estimate track length bymuboy[18]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
5.15 The distributions of precise fit direction and true direction for a 40-year SK-I
MC sample. Taken from [18]. . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
5.16 Muon energy losses in water. Curve (a) shows losses due to ionization, (b)shows radiative losses, and (c) shows the sum of both[91]. . . . . . . . . . . . 98
5.17 A 2-dimensional scatter plot of ∆ (vertical axis) and χ2 (horizontal axis) for
showering event selection for SK-IV data. Magenta events above the red line
are saved as showering, and blue ones below as non-showering. . . . . . . . . 101
5.18 The distributions of fitted tracklength (left) for through-going events andfitted momentum (right) for stopping events from precisefit for SK-IV data
and MC. Final data cuts are applied to the data, but not the MC. . . . . . . 102
5.19 The event rate output of umred1st (black), umred2nd (red), precisefit(green), umred3rd (blue) and umred4th (magenta) for SK-IV from the automated upmu reduction software. . . . . . . . . . . . . . . . . . . . . . . . . . 104
5.20 The background subtraction method for near-horizontal upward-going muons,shown for SK-IV through-going muons. The LHS shows the azimuthal distribution of events summed over all zenith angles, and the RHS shows thezenith angle distribution of events with the BG fit function overlaid. . . . . . 105
vii

6.1

Typical angular distribution of the observed, corrected charge in number of
p.e. (top) and its second derivative with respect to angle (bottom), which
are used in the ring edge search algorithm. Taken from [108]. . . . . . . . . . 109

6.2

Difference between true and reconstructed vertices for SK-IV FC multi-GeV
multi-ring MC. The shaded region covers 68% of all events. Taken from [92]. 111

6.3

The Hough transformation method for finding Cerenkov rings. Ring candidates are guessed by drawing circles with 42◦ opening angles surrounding all
hit PMTs and intersection points are assumed to be the center of the ring. . . 111

6.4

Ring-counting separation for sub-GeV events for SK-I,-II,-III and -IV. Data
are in black and oscillated MC (sin2 (2θ23 )=1.0,∆m223 =2.5×10−3 ) are the colored histograms. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113

6.5

Ring-counting separation for multi-GeV events for SK-I,-II,-III and -IV. Data
are in black and oscillated MC (sin2 (2θ23 )=1.0,∆m223 =2.5×10−3 ) are the colored histograms. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114

6.6

Event displays for typical (1 GeV) e-like (left) and µ-like (right) MC events. . 115

6.7

PID separation for SK-I,-II,-III and -IV FC sub-GeV single-ring samples ofdata (black) and MC (histograms). . . . . . . . . . . . . . . . . . . . . . . . . 117

6.8

PID separation for SK-I,-II,-III and -IV FC sub-GeV multi-ring samples ofdata (black) and MC (histograms). . . . . . . . . . . . . . . . . . . . . . . . . 118

6.9

Distance between true and reconstructed vertex for SK-IV FC MC e-like(left) and µ-like (right) samples. Sigma (shaded region) represents 68% of
the full area. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119

6.10 Angle between true and reconstructed direction for SK-IV FC MC e-like (left)
and µ-like (right) samples. Sigma (shaded region) represents 68% of the fullarea. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
6.11 Improvement to the reconstructed invariant π 0 mass with an updated fitting
algorithm shown in red compared to the original algorithm in blue. . . . . . . 121
7.1

The layout of bins used in the zenith-angle analyses. There are 370 for FC,
60 for PC, and 50 for Upmu events, giving 480 in total. White boxes indicate
those for which 10 zenith bins are defined and shaded indicates single zenithangle bins. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124

7.2

A typical CCνe interaction which produces an e-like multi-ring event in SK. . 125

7.3

Contours of the probability of νµ survival as a function of pathlength and energy assuming two-flavor oscillations with (∆m223 = 2.11 ×10−3 eV2 , sin2 (2θ23 )= 1.00). Red areas indicate no oscillations and green areas indicate maximum
oscillations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129

7.4

The components entering the absolute normalization uncertainty (top) andthe overall uncertainty (bottom) as functions of energy[48]. . . . . . . . . . . 132
viii

7.5

A plot of the ratios of the FLUKA[4] (dashed red) and Bartol[3] (solid black)flux models to the Honda[49] model. . . . . . . . . . . . . . . . . . . . . . . . 132

7.6

Predicted cross-sections (left) and ν̄/ν cross-section ratio for νµ p → µ− pπ +(solid), νµ n → µ− pπ 0 (dashed) and νµ n → µ− pπ + (dotted) for the Hernandez
model[85] normalized to the Rein-Sehgal model[101]. . . . . . . . . . . . . . . 135

7.7

Comparisons between data and NEUT simulations for π + quasi-elastic (left)
and absorption (right) interactions demonstrating the weighting method usedin the nuclear effects uncertainty term. . . . . . . . . . . . . . . . . . . . . . . 136

7.8

Pion multiplicity as a function of W2 for NEUT predictions (left) and CHORUS experimental measurements (right). . . . . . . . . . . . . . . . . . . . . . 138

7.9

The 2-d allowed regions for ∆m223 versus sin2 (2θ23 ) at 68%, 90% and 99%
C.L. for SK-I+II+III+IV data. The black star indicates the best-fit point. . . 140

7.10 The 1-dimensional χ2 distributions for ∆m223 (left) and sin2 (2θ23 (right, physical region only) fit to the SK-I+II+III+IV data. . . . . . . . . . . . . . . . . 141
7.11 The 90% C.L. allowed regions of the most recent 2012 results for this analysis,
the SK L/E analysis, the SK full three-parameter zenith-angle analysis, T2K
and MINOS. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142
7.12 The 2-dimensional 68%(blue), 90%(red) and 99%(green) C.L. allowed regions
for the SK-I data (left) and SK-II data (right) with the updated MC andreconstruction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142
7.13 The 2-dimensional 68%(blue), 90%(red) and 99%(green) C.L. allowed regions
for the SK-III data (left) and SK-IV data (right) with the updated MC andreconstruction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143
7.14 The zenith-angle distributions of the SK-I+II+III+IV data (black points),
the SK-I+II+III+IV MC oscillated with the best-fit parameters from thisanalysis (∆m223 , sin2 (2θ23 )) = (2.30 ×10−3 eV2 , 0.99) (red histogram) andthe MC with no oscillations (blue histogram). . . . . . . . . . . . . . . . . . . 144
7.15 The ratios i /σi for systematic error terms from the 2-flavor fit to the SKI+II+III+IV data. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145
7.16 The ratios i /σi for systematic error terms from the 2-flavor fit to the SKI+II+III+IV data. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146
8.1

Allowed regions from original fit to MINOS anti-neutrino data within 68.3%(blue), 99.5% (red), and 99.7% confidence limits for each pair of parameters,with the best fit point indicated by the star. . . . . . . . . . . . . . . . . . . . 151

8.2

Earth density profile as a function of radius according to PREM[20] on theleft, and the average line-integral density as a function of cosine zenith angle
on the right. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152
ix

8.3

Effective mass squared splittings in matter as a function of energy in the apparent CPT-violation model. The 2-flavor atmospheric mass squared splitg
ting is shown in the left-hand figure for comparison. A density of ρ = 4.0 cm
3for the Earth is assumed. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153
8.4 Effective mixing angle in matter as a function of energy. A density of ρ =g
4.0 cm3 for the Earth is assumed. . . . . . . . . . . . . . . . . . . . . . . . . . 154
8.5 Theoretical survival probabilities for muon neutrinos and anti-neutrinos in
the apparent CPT-violation model. . . . . . . . . . . . . . . . . . . . . . . . . 154
8.6 Difference of oscillation survival probabilities for muon neutrinos and antineutrinos between the apparent CPT-violation model and the standard 2flavor model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155
8.7 The 68% (blue), 90% (red) and 99% (green) allowed regions in the g versus
m space for the fit to the SK-I+II+III+IV data, with M minimized away. . . 157
8.8 The 68% (blue), 90% (red) and 99% (green) allowed regions in the g versus
M space for the fit to the SK-I+II+III+IV data, with m minimized away. . . 158
8.9 The 68% (blue), 90% (red) and 99% (green) allowed regions in the M versus
m space for the fit to the SK-I+II+III+IV data, with g minimized away. . . . 159
8.10 The 1-d χ2 distribution and 68% (blue), 90% (red) and 99% (green) C.L. forthe active neutrino mass m after fitting to the SK-I+II+III+IV data. . . . . 160
8.11 The 1-d χ2 distribution and 68% (blue), 90% (red) and 99% (green) C.L. forthe sterile neutrino mass M after fitting to the SK-I+II+III+IV data. . . . . 161
8.12 The 1-d χ2 distribution and 68% (blue), 90% (red) and 99% (green) C.L. forthe new coupling interaction g after fitting to the SK-I+II+III+IV data. . . . 162
8.13 The 1-d χ2 distribution and 68% (blue), 90% (red) and 99% (green) C.L. forthe new coupling interaction g in a sensitivity study for livetimes 0f 7.9, 15.8
and 50 years with the SK-I+II+III MC. . . . . . . . . . . . . . . . . . . . . . 163

x

GLOSSARYSK:Super-Kamiokande

ID:

Inner detector

OD:

Outer detector

FV:

Fiducial volume

MC:Monte Carlo

CC:

Charged current

NC:

Neutral current

FC:

Fully contained

PC:

Partially contained

FCFV: Fully contained fiducial volume

QBEE: QTC-based electronics with ethernet

QTC: Charge(Q) to time(T) converter(C)PMT: Photomultiplier tube

DAQ: Data acquisition system

TDC: Time to digital converterxi

ACKNOWLEDGMENTS

The work in this thesis would not have been possible without the help of many people.
First and foremost I would like to thank Professor Jeff Wilkes for advising me throughout
my graduate school career. His mentoring style suited me perfectly and allowed me toenjoy my experience as a graduate student. The introductions to the Super-Kamiokande
Collaboration and Japan that he gave me will leave a lasting mark on me.
I am indebted to Roger Wendell for personally guiding me through the analyses in thiswork. His endless efforts for both Super-Kamiokande as a whole and my personal progress
have been invaluable and his friendship has made Japan a memorable experience. I am also
very grateful for the wonderful and patient instruction I received from Jen Raaf, who wasalways willing to take the time to help me.
I would also like to extend thanks to Ed Kearns, Shiozawa-san, Chris Walter and Kate
Scholberg for their excellent guidance within the AtmPD group, as well as to all Japanesecolleagues whose tireless efforts have allowed Super-Kamiokande to be so successful.
I have thoroughly enjoyed my time as a graduate student at the University of Washington
because of the friends I have been lucky to surround myself with. Kevin Connolly has been agreat friend throughout and I hope our discussions about physics, language and life continue
for a long time. The rest of my graduate school comrades, including Laura Bodine, Eric
Wallace, Jason Dexter and many others have made Seattle a wonderful experience.
I have also shared many wonderful experiences with my fellow students in Japan, including but not limited to Kota, Piotr, Maggie, Joshua, T, Tanaka-san and Ikeda-san. I am also
very thankful for my many friends in Buffalo, especially Zack Long and Dave Stockburger,
who have not forgotten me as my dreams have taken me far away from home.
The love, support and encouragement of my family, especially Mom, Dad, and Grandma
Jane, has been invaluable and for that I cannot thank them enough.

xii

DEDICATION

For my family.

xiii

1

Chapter 1
INTRODUCTION

The study of neutrinos began in 1930 when Wolfgang Pauli postulated the existenceof a neutral particle with a very small mass in order to explain the continuous electronenergy spectrum found in β decays[99]. Since then, the field has been driven by vast
theoretical literature and exciting experimental results. The neutrino has been the targetof experimental studies for approximately 60 years and as of this writing new results arestill being released and future experiments are in their early design phases. The current andfuture experiments have taken neutrino physics beyond the discovery stage and into the eraof precision measurements. The present work builds upon the vast experimental framework
set up by previous results and also hints at the shift to precision measurements in the fieldof atmospheric neutrino oscillations.
1.1

Neutrinos

Neutrinos are electrically-neutral spin-1/2 leptons which only undergo weak interactions.
In the standard model (SM) of particle physics, neutrinos are massless and come in threeflavors, electron (νe ), muon (νµ ) and tau (ντ ), analogous to the charged leptons. They
are also experimentally known to exist only as left-handed neutrinos or right-handed antineutrinos[80].
Though existing copiously throughout the universe, they are experimentally very difficult
to detect as they can pass through tremendous amounts of matter without interacting.
Soon after their initial theoretical proposal by Pauli in 1930, Fermi generated a point-liketheory for neutrino interactions[74]. The currently verified unified electroweak theory which
describes neutrino interactions would not be fully developed until the 1960s and the weakinteraction mediating bosons, the W± and Z0 were eventually discovered in 1983. The large
masses of these vector bosons (MW = 80.6 GeV/c2 , MZ = 91.2 GeV/c2 ) result in a weak2interaction range of only about 10−18 m. The neutrino itself was first definitively detected
in 1956 by Cowan and Reines using the inverse β decay reaction ν̄ + p → e+ + n taking
place in a liquid scintillator surrounded by photomultiplier tubes (sensitive photon detectors
which will be described in more detail later)[90]. The source of the electron antineutrinosin this experiment was a nearby nuclear reactor facility.
Neutrinos have since been detected from natural sources such as the Sun[16], Earth’s
atmosphere[67], the core of the Earth[10], man-made accelerator beamlines[52] and a Supernova event in 1987[41]. In addition to electron flavor neutrinos, muon[29] and tau flavor
neutrinos[77] have also been detected, and experiments studying weak boson decays (inparticular, the Z0 boson which mediates NC weak interactions) currently limit the numberof light species which participate in weak interactions to these three[9].
1.2Atmospheric Neutrinos

This work focuses on the behavior of atmospheric neutrinos. Cosmic rays which are predominantly high-energy protons arrive from space and interact with nitrogen and oxygen
molecules in the Earth’s atmosphere. Mesons (quark-antiquark pairs) such as pions (π) andkaons (K) are produced in large numbers in these reactions. These mesons are not stable
and their decays produce neutrinos. A typical decay chain beginning with a pion has theform
π + → µ+ + ν̄µ(1.1)

which produces a muon antineutrino. The muon then decays by the processµ+ → e+ + νe + ν̄µ(1.2)

producing an additional muon antineutrino and a single electron neutrino. This process was
assumed to begin with a π + particle, but the charge conjugate of this process which produces
two muon neutrinos and one electron antineutrino is also possible. With these decay chainsproducing the atmospheric neutrino flux (the fraction of neutrinos produced by kaon decays
becomes significant at really high energies), we expect the flux ratio (νµ + ν¯µ /νe + ν¯e ) to beroughly two. In reality this ratio is energy-dependent, but this assumption for the ratio isgenerally valid for neutrino energies up to around 1 GeV. More energetic muons may reach3

the surface of the Earth before decaying and so a ratio smaller than two would be expected
for a terrestrial experiment at higher neutrino energies. Though the uncertainties on theabsolute neutrino and antineutrino fluxes are relatively large (they can be on the order of
10-25% depending on energy), this flavor ratio is estimated with an uncertainty of only
2% [48]. Significant deviations of this ratio from 2 offered the first evidence that neutrinos
exhibited behavior beyond their description in the standard model.

1.3

Neutrino Oscillations

Contrary to their masslessness in the SM, evidence for massive neutrinos has been found via
flavor transformations, or “oscillations”, in neutrinos produced from the sun[61], hadronic
decays in the atmosphere[65], accelerator beamlines[52, 36] and nuclear reactors[59, 27, 43].
Neutrino oscillations are sensitive to the differences in the squared masses of the states,
defined by ∆m2ij ≡ m2j − m2i . Thus, measurements of this parameter by neutrino oscillation
experiments places a lower bound on at least as many neutrino mass states as the numberof separate mass-squared splittings.
Independent of the oscillation experiments, upper constraints have been placed on theabsolute neutrino masses. Tritium beta-decay end-point measurements have restricted theeffective electron-neutrino mass to be mνe < 2.0 eV [63, 25, 34]. Relatively high limits foref ff
the effective muon and tau neutrino masses of mef
νµ < 0.19 MeV and mντ < 18.2 MeV

also exist [34]. These limits are taken from pion decay (π + → µ+ + νµ ) and tauon decay
spectrum shape analyses and are limited by the short lifetimes of the decaying particles
involved. Cosmology also places very rough bounds on the upper limit of the sum ofP
neutrino masses. Current limits are in the range ν ≤ (0.75 − 5.5) eV[34][86][2].
Neutrino oscillations occur when the neutrino flavor states as defined by charged-current
weak interactions and the mass eigenstates are not identical, and neutrinos also have nonzero masses. These were first proposed theoretically in the 1960s by Maki, Nakagawa and
Sakata[93] and Pontecorvo[100]. The basic assumption is that the flavor states (indexed byletter) can be written as a linear superposition of three mass states (traditionally indexedby number):4|να i =

3
X

∗
Uαi
|νi i

(1.3)

i=1

The transformation can be fully written out in terms of a unitary matrix transforming
between the two sets of basis states. It has the general form


νe





Ue1

Ue2

  
  
νµ  =  Uµ1 Uµ2
  Uτ 1 Uτ 2
ντ

 
ν
  1 
Uµ3  ν2  .
 
ν3
Uτ 3
Ue3

(1.4)

It then takes the form



νe





1

0

0



cosθ130 sinθ13 e−iδ

  

  

νµ  =  0 cosθ23 sinθ23  
0
1
  

ντ
0 −sinθ23 cosθ23−sinθ13 eiδ 00cosθ13cosθ12  −sinθ12

0

 
ν
  1 
cosθ12 0  ν2 
 
0
1
ν3
(1.5)
sinθ120where the so-called mixing angles θ12 , θ13 , θ23 and the CP-violating phase δ are the fourparameters of the matrix which can be experimentally constrained. The separation of the
sub-matrices roughly represent different areas of experimental research in neutrino oscillations. In particular, atmospheric and long baseline accelerator experiments such as SK havemeasured θ23 , reactor and accelerator experiments have recently measured θ13 , and solar
experiments (including SK) have measured θ12 . There has been no measurement of δ as ofyet. The current best-fit values for the mixing angles from a 2011 global analysis are shownin Table 1.1.

The experiments involve transformations of one neutrino flavor να into another, νβ , after
propagation of a baseline L. It is useful to express this mathematically via the oscillation
probability P(να → νβ ). In the case where α = β the oscillation probability becomesa survival probability. In the following calculations the letters a,b and c represent flavor
eigenstates, and the letters i,j, and k represent mass eigenstates. We can calculate the5

Parameter

Best-Fit (±1σ)

∆m223

2
−3
2.35+12−.09 × 10 eV

sin2 θ23

0.42+.08−.03
2−5
7.58+0.22−0.26 × 10 eV

δm212
sin2 θ12+.018
0.306−.015

sin2 θ13+.007
0.025−.007

Table 1.1: Global analysis constraints on the three mixing angles and two mass-squared
splittings governing neutrino oscillations[76].

∗ |ν i. Applying time
probability by assuming the initial state at t = 0 to be |νa = Uajjevolution, the state at time t is then∗ −iEj t|ν(t)i = Uaj
e
|νj i

(1.6)

We can then calculate the probability amplitude of detecting the neutrino in flavor state|νb i at time t:∗ −iEj t∗ −iEj t∗
A(νa → νb ; t) = hνb |ν(t)i = Uajehνb |νj i = Ubi Uajehνi |νj i = Ubj e−iEj t Uaj
.

(1.7)

Here summation over repeated indices is implied. To convert this amplitude into a probability, we take the product of this with its complex conjugate:∗ 2
P (νa → νb ; t) = |A(νa → νb ; t)|2 = |Ubj e−iEj t Uaj
| .

(1.8)

These oscillation probabilities are in general complicated to calculate due to the entries in
the matrix U , especially when mixing of all three flavors is considered. A simplified caserelevant for the analysis in this thesis will be shown in detail in the following subsection.

1.3.1

Two Flavor Oscillations

A good demonstration of a tractable closed-form oscillation probability is the two flavor
oscillation case. We will assume oscillations only of the form νµ ↔ ντ , which will be used in6both analyses in this thesis. The 2x2 mixing matrix then takes the form of the lower right
2x2 section of the first sub-matrix in Eq. 1.5, as follows:
 
  
ν
cosθ23 sinθ23ν  2 µ = ν3−sinθ23 cosθ23
ντ

(1.9)

We insert into Eq. 1.8 the matrix elements from this 2x2 matrix and the energies Ei ,which can be given byqm2
m2
Ei = p2 + m2i ≈ p + i ≈ p + i .
2p2E(1.10)

The approximations made here are valid for relativistic neutrinos for which p  mi , which
is true for atmospheric neutrino energies because of the upper limits placed on the absolute
masses. Inserting this into 1.8, we get22

P (νµ → ντ ; t) = sin 2θ sin




∆m2
t
4E

(1.11)

where ∆ m2 is equal to m23 - m22 . For relativistic neutrinos we can replace t with L, and if
we rewrite the probability in units useful for atmospheric oscillations with L in km and E
in GeV, we arrive at2LP (νµ → ντ ; L) = sin 2θ sin 1.27∆m
.
E
2

2

(1.12)

The first sine term in this probability is often referred to as the amplitude of the oscillations
and is solely dependent on the mixing angle θ23 . The second sine term contains an L/E
dependence which in combination with the mass-squared splitting ∆m223 leads to a signature
zenith-angle dependent deficit of muon neutrinos in atmospheric neutrino data characteristic of oscillations. The wide range of pathlengths and energies spanned by atmospheric
neutrinos provide are suitable for probing large regions of oscillation parameter space.
The two flavor oscillation case is only a simplification of full three-flavor mixing, but isa valid approximation when two conditions are met. The first is that oscillations driven by
the other mass-splitting do not contribute significantly, which occur when∆m212 L
 1.
E

(1.13)

7

This condition allows the simplification ∆m212 → 0 to be made. The second simplificationis the matrix element |Ue3 | = 0. This element was recently confirmed to be non-zero, but
the condition |Ue3 |  1 is still a respectable approximation given the global best-fit value
shown in 1.1. It is these two assumptions which allow the full three-flavor oscillation caseto be reduced to the simple two-flavor form of Eq. 1.12.
1.4

Sterile Neutrinos

In addition to the three mass eigenstates driving the standard oscillations, there is also
the possibility of additional neutrino mass eigenstates. With the addition of new mass
eigenstates, a new linear combination of the states can be constructed as in|νS i =NX∗
Uαi
|νi i,(1.14)

i=1

where N is now equal to the total number of mass eigenstates. Neutrino flavor states
described by this linear combination are called sterile neutrinos because they do not participate in Weak interactions. With sterile neutrinos, the unitary transformation matrix
expands to NxN, where N is the total number of active plus sterile species. In general theoscillation probabilities will change, usually becoming more complicated, and effects suchas extra active disappearance and short baseline appearance can be generated. There has
been no experimental confirmation of sterile neutrinos yet but their existence would explainsome controversial experimental results.
In 1995 the LSND experiment searched for oscillations of the form ν̄µ → ν̄e with neutrinos from µ+ decay at rest over very short baselines approximately 30 m with neutrinoenergies of a few tens of MeV[69]. The experiment saw an excess of ν̄e above the expectedbackground with 3.8σ significance. If interpreted as oscillations, the necessary oscillation
frequency to drive oscillations with this ratio of L/E would be on the order of 1 eV2 . This
is not compatible with the two measured ∆m2 values shown in Tab.1.1. The additionof a fourth mass eigenstate and by extension a sterile neutrino could explain the shortbaseline oscillations. However, experiments which sought to confirm the LSND result, suchas KARMEN[13, 68] and MiniBooNE[70], did not see evidence for the same type of oscillations and therefore ruled out most of LSND’s allowed regions, as shown in Fig. 1.1.

8

Figure 1.1: The LSND allowed region[69] overlayed with the exclusion regions fromKARMEN[68] and MiniBooNE[70] for short baseline oscillations with large ∆m2 . Imagetaken from [70].

However, a recent combined analysis done by the MiniBooNE collaboration claims to haveseen an excess consistent with LSND[71]. For atmospheric oscillations, Super-Kamiokande
conducted an analysis of their atmospheric neutrino dataset and concluded that for νµdisappearance the ντ mode was favored over any sterile neutrino mode[60].

Theorists have constructed many exotic models based on sterile neutrinos to explain thisoscillation anomaly while maintaining compatibility with other experimental constraints(see [96, 79] for two such models). These models in general employ 1 or 2 sterile neutrinos
in 3+1,3+1+1 or 3+2 frameworks, among others. In addition to terrestrial oscillation
experiments, cosmologists have also recently used radiation studies related to the cosmic
microwave background to suggest the existence of one or two sterile neutrinos[81]. Studies
of sterile neutrinos remain an open and interesting topic for experimental searches.

9

1.5

Oscillations in Matter

Neutrino oscillations in matter are in general not equivalent to those in vacuum. In particular the presence of electrons in ordinary matter further complicates mixing. In addition
to scattering and absorption interactions, neutrinos in matter experience effective potentialfrom coherent forward scattering in CC and NC Weak interactions. All flavors experience
an identical potential from NC scattering interactions given by
GF
VN C = ± √ Nn ,2

(1.15)

where the sign is positive for neutrinos and negative for antineutrinos. Here Nn represents
the number density of neutrons and GF is the Fermi constant. Ordinary matter is assumed
to be electrically neutral so the number of electrons and protons cancel and do not contributeto this potential. Additionally, electron neutrinos can undergo charged-current scatteringby exchange of a W boson due to the presence of electrons. This term can be written as√
VCC = ± 2GF Ne ,(1.16)

where Ne is now the electron number density and the same sign definitions as above hold.
The fact that these potentials are different for electron neutrinos than they are for muon andtau neutrinos can create interesting effects in oscillations. As an example, in the two flavor
case where νe only mixes with νµ , closed-form oscillation probabilities can be calculated.
The oscillation probability has the same form as in the vacuum two flavor case:L.
P (νe → νµ ; L) = sin2 2θ sin2 πlm(1.17)

In matter the “oscillation length” lm is equal tolm = r2π∆m2
2E cos2θ0 −√2GF Ne

2

+



∆m2
2E

2

.

(1.18)

2

sin 2θ0

Likewise, the expression for the mixing angle in matter also has a similar resonance form: 2 2∆m
sin2 2θ02E
2.
(1.19)
sin 2θ = 
2 
2
√
∆m2
∆m2
2
cos2θ
−
2G
N
+
sin
2θ
0
0
F e
2E
2E

10

The maximum value of the mixing angle in matter occurs when the condition√

2GF Ne =∆m2
cos2θ02E(1.20)

is true. Thus, large or even maximal mixing can occur in matter even when the vacuum
mixing angle is small. This particular condition is called the MSW resonance condition
and is used to explain the oscillation of νe states produced in the sun. The oscillations
assumed here are for νe ↔ νµ . In the two-flavor approximation, there is no effect frommatter in νµ ↔ ντ oscillations because both neutrinos experience only VN C (i.e. there
are no CC interactions because µ and τ are not present in ordinary matter). However, if
full three-flavor mixing is considered, matter effects apply because the two flavors can mix
with the electron neutrino. Additionally, mixing νµ and ντ with a sterile neutrino will also
introduce matter effects, because the sterile neutrino does not even undergo NC interactions(i.e. VN C = 0) and so a relative potential difference exists between the active and sterile
species, just as in the case when mixing with νe was assumed.
1.6

Summary Of This Work

This dissertation uses the entire Super-Kamiokande atmospheric neutrino dataset to placeconstraints on two neutrino oscillation models. The first model is the standard two-flavor
νµ disappearance analysis. This model has becoming very well constrained by both atmospheric and long-baseline accelerator experiments. This analysis differs from previous
Super-Kamiokande 2-flavor zenith-angle analyses by using new SK data (1,097 days livetime), an updated detector simulation, updated neutrino interaction simulations and a new
analysis binning technique (increasing the overall number of bins). These updates are sharedwith the 2-flavor L/E analysis and full 3-flavor analysis which are also being studied in otherSK dissertations this year[92, 94]. The second model which is studied in this work is an
exotic model which includes both sterile neutrino and matter oscillation effects. It was firstpublished in 2011 and has not been constrained beyond the authors’ initial constraints. The
model is analyzed with the same dataset and simulations as the 2-flavor model.

11

Chapter 2
THE SUPER-KAMIOKANDE DETECTOR
2.1

Introduction

Super-Kamiokande (SK) is a water Cerenkov detector at the Kamioka Observatory of theInstitute for Cosmic Ray Research, University of Tokyo. It is located in the town of Kamioka
in Gifu prefecture, Japan, west-northwest of Tokyo (see Fig.2.1). The detector sits inside
the Mozumi Mine, which is operated by the Kamioka Mining and Smelting Company. It is
accessible via a 1.8 km-long mine tunnel which is nearly horizontal.
The geographic coordinates of the detector are 36◦ 250 32.600 N, 137◦ 180 37.100 E[12]. The
detector cavity lies within Mt. Ikenoyama as shown in Fig. 2.2, covered by 1,000 meters ofrock overburden, or 2,7000 meters-water-equivalent(m.w.e). This overburden reduces therate of downward-going cosmic ray muons reaching the detector to ∼ 3 Hz.
It currently operates both as a standalone detector and also as the far detector for theT2K Experiment[39]. Construction of the detector began in 1991 and data-taking began on
April 1, 1996. The detector has taken data in four distinct phases, labeled SK-I(1996-2001),
SK-II(2002-2005), SK-III(2006-2008), and SK-IV(2008-present). In addition to neutrino
oscillation physics, the detector also conducts research on neutrino astrophysics[44], nucleon
decay searches[32], dark matter searches[62], and various other particle physics topics.
2.2

Cerenkov Radiation

SK relies on the phenomenon of Cerenkov radiation to detect neutrino interactions. Coherent electromagnetic radiation is generated when a charged particle moves through a medium
faster than the local speed of light, analogous to the ’sonic boom’ shock wave generated byan object moving faster than the speed of sound in air. The radiation is conical and allows
us to determine the incoming direction of neutrinos that have interacted in the detector.
The local speed of light in a given medium is given by

12

Figure 2.1: Location of the Super-Kamiokande experiment in Japan. Taken from Ref. [12]Figure 2.2: An overview of the SK detector and its location within Mt. Ikenoyama. Taken
from Ref. [12]

13Figure 2.3: A schematic of the Cerenkov radiation emission angle relative to the charged
particle’s track. Taken from Ref. [112]

v=

c
n

(2.1)

where n is the index of refraction of the medium. Charged particles traveling at speeds
equal to greater than this will emit Cerenkov radiation. This speed corresponds to Cerenkov
emission energy thresholds for electrons of about 0.77 MeV and for muons of about 160 MeV.
The angle of the Cerenkov radiation cone relative to the charged particle’s track also
depends on n and is equal to [89]cos(θc ) =

1
nβ

(2.2)

.
A schematic showing the Cerenkov radiation emission angle for a particle moving withspeed v/c is shown in Fig. 2.3.
In the ultra-relativistic limit where β → 1, the angle only depends on the index ofrefraction. In this case, pure water (n ∼ 1.33) yields an angle of 42◦ .
The amount of energy lost per unit track length via this radiation is given by theexpression [89]
dE
(ze)2
= 2
dx
c

Z
n(ω)>(1/β)

1
1− 2 2β n (ω)


ω dω,

(2.3)

14

which equals approximately 500 eV/cm for the radiation wavelengths relevant to SK(300-600 nm). Note that this loss is small compared to the overall energy loss of charged
particles in water due to scattering and ionization, which for minimum ionizing muons isaround 2 MeV/cm. The number of photons emitted per unit track length is given by ze 2 1
d2 N
=
dxdω
c
~


1−

1
2
β n2


.

(2.4)

For particles such as electrons and muons with z = 1, e = 1 and β ∼ 1, there areapproximately 200 photons per cm emitted along the particle’s track in pure water. As will
be described in following sections, accurate detection of the timing and quantity of photons
will allow reconstruction of Cerenkov rings and thus the direction and energy informationof charged particles in SK.
2.3

The SK Detector

The detector consists of a steel cylindrical tank 39.3 meters in diameter and 41.4 meterstall, with a nominal water capacity of 50 ktons. On top of the tank there is a dome area
where five huts containing the electronics for the experiment are located. There is a also acontrol room connected to the dome area by a tunnel that is used for monitoring shifts, as
shown in Fig. 2.2. The daytime shift is done form inside the control room, while evening
shifts are conducted remotely.
Within the tank is a steel scaffolding structure approximately 2 meters from the outerwall which houses photomultiplier tubes (PMTs) and divides the detector volume into two
optically separated concentric cylinders, an inner detector (ID) and outer detector (OD).
The structure and PMTs are shown in Fig. 2.4. The OD acts as a veto and is defined by
the area extending 2.05 meters inward from the tank walls and 2.2 meters from the tank top
and bottom, containing approximately 18 ktons of water. The ID is defined as the insideof the scaffolding and contains 32 kilotons of water. The scaffolding structure itself is 0.55
meters in thickness, which is dead space through which light is not able to pass.
The ID side of the scaffolding is covered in a black low-density polyethylene sheet which
absorbs light. This prevents light from leaking from the ID into the OD (and vice versa),

15

Figure 2.4: The scaffolding structure which divides the ID and OD, along with the mounted
PMTs. Taken from Ref. [12]and also prevents reflected or scattered light from bouncing around the ID. The OD sideof the scaffolding and the inside of the outer detector wall are covered with white reflective
Tyvek R (manufactured by DuPont) sheets. This increases the light collection efficiency of
the OD PMTs. Tyvek has a reflectivity of ∼ 90% for a wavelength of 400 nm.
2.4

PMTs

The photomultiplier tubes housed on the scaffolding act as light-collection detector elements
viewing the ID and OD water volumes. PMTs operate by detecting photoelectrons emitted
via the photoelectric effect when a photon strikes their photocathode surface. This initial
photoelectron is multiplied into a large cascade of electrons by a series of interactions withdynodes at increasingly higher potentials. The gain factor for the PMTs used in SK isabout 107 , which means that 1 photoelectron becomes an output signal of roughly 107
photoelectrons after the final dynode (so that 1 input p.e. would produce a charge outputof ∼ 2 pC).

16

Figure 2.5: Schematic view of an ID 50 cm PMT. Taken from Ref. [12]The ID contains 11,129 inward-facing Hamamatsu Type R3600-006 PMTs 50 cm indiameter. A schematic of an ID PMT is shown in Fig. 2.5.
The PMTs have a bi-alkali photocathode, with total photocathode coverage of the ID
surface of effectively 40%. The number of PMTs and their effective photocathode coveragein the detector was significantly reduced in SK-II due to a major accident in 2001[11].
Reconstruction work between the SK-II and SK-III phases brought the detector back to
its original PMT configuration. Additionally, the ID PMTs were placed inside fiberglass
enclosures with acrylic hemispheres covering the photocathode surface, to prevent similar
accidents in the future.
The ID PMTs have a peak quantum efficiency of about 20% within the range of wavelengths 360-400nm, as shown in Fig. 2.6, and are operated at voltages between 1,700 and
2,000 V. Each PMT has a dynamic range from 1 pe to ∼ 1200 pe (to be discussed in moredetail in the section about the DAQ, which sets these limits), allowing SK to detect eventswith energies between 4.5 MeV and over 1 TeV, a wide range of energies spanning fromlow-energy solar neutrinos to high-energy atmospheric neutrinos. The timing resolution of
the ID PMTs is 2-3 ns.
The OD contains 1,885 outward-facing PMTs each 20 cm in diameter. They are a
combination of old PMTs recycled from the IMB experiment [56], and new Hamamatsu

17

Figure 2.6: Quantum efficiency of the ID PMTs. The peak efficiency coincides with UVradiation wavelengths. Taken from Ref. [12]R1408 PMTs. The IMB PMTs are concentrated near the top and upper part of the barrel
while the majority of the barrel and the bottom of the tank have the new PMTs. Squareacrylic wavelength shifting (WS) panels (60 cm per side) are attached to each PMT toenhance the light collection efficiency, achieving an enhancement ratio of 1.5 relative to abare PMT. The plates absorb UV light and emit blue-green light closer to the sensitivity
peak of the PMTs. The timing resolution of the OD PMTs is approximately 15 ns, inpart due to the presence of the WS panels. While this is much poorer than the ID PMTs,
the extra light collection efficiency it creates is valuable for use as a calorimeter and vetocounter.

A magnetic field of 450mG exists naturally in the detector cavity volume. The timing
and photoelectric trajectories of the PMTs are sensitive to such a magnetic field. To reduce
the field strength, 26 sets of horizontal and vertical Helmholtz coils are placed inside thetank. Their operation reduces the average field inside the detector cavity to about 50mG.
A ∼ 70 m coaxial cable is connected to each PMT for the high voltage supply and datatransfer. The PMT signals are sent up through these to one of four electronics huts locatedin the dome area on top of the tank.

18

2.5

Data Acquisition System

During the first three phases (SK-I, SK-II, and SK-III; 1996-2008) of the Super-Kamiokande
experiment, the data acquisition (DAQ) system did not undergo any major upgrades. In
2008, in anticipation of the SK-IV running period and the T2K Experiment[39], the DAQ
front-end electronics underwent a complete upgrade. The data analyzed in this work spans
all four phases of the experiment. We present a brief overview of the original DAQ system,
and a more detailed overview of the new SK-IV system. The reader can refer to Ref. [12]for a more detailed description of the original DAQ system, and Ref. [33] for further details
on the specifics of the new SK-IV electronics.

2.5.1

Original DAQ

The original DAQ featured independent PMT readout electronics for the ID and OD PMTs.
The ID PMT signals were read out by Analog Timing Modules (ATMs), each of which
managed 12 PMTs, located in the electronics huts. Each ATM featured two channels per
PMT and switched between the channels to reduce readout dead time. Square wave pulses
were created by the ATM upon receiving a PMT signal, and then sent to a global hardware
analog-sum “hitsum” trigger. This trigger activated at a threshold of 29 PMT hits. In theevent of a global trigger the digitized charge and timing information for each ID PMT was
saved for a time window from 400 ns prior to 900 ns after the trigger. The event information
was finally sent to offline computers by custom built “TKO” crate interfaces via an optical
network. The ATMs were saturated at charge inputs of around 600 pC (a feature which
was improved with the new SK-IV system).
The OD PMT signals were readout by “paddle cards”, which each handled signals from
12 PMTs. Charge-to-time conversion (QTC) modules produced a square wave with thetiming information defined by the leading edge and the overall charge defined by the width.
The QTC output signals were digitized by LeCroy Module 1877 TDCs with 0.4 ns resolution,with a FASTBUS interface to the OD workstation, a Sun SPARCstation (later replaced bya Linux PC). For global OD triggers, information was saved between 10 µs prior to and 6µs after the trigger.

19

In both cases, individual PMT signal thresholds were set to 0.25 p.e. This DAQ electronicsystem was used until the end of the SK-III period in 2008.

2.5.2

New DAQ

The new DAQ system provides a stable setup for the next 10 to 20 years and many direct
benefits to SK physics analyses. To lower the energy threshold for solar neutrino studies, ahigh-speed DAQ is necessary to combat the increase in background rates. Also, a supernova
burst in the galactic center could produce tens of thousands of events in SK in a 10s time
window. The original DAQ was capable of recording approximately 60,000 events in thistime span, but the new DAQ increases that capacity by a factor of 100. The new deadtime free system also allows all PMT hits to be recorded. In the original system only
PMT hit information within the defined global trigger windows was recorded. This feature
also has benefits for supernova relic neutrinos studies. These events are distinguished fromatmospheric neutrino background events by a 2.2 MeV gamma ray emitted with a 200 µs
delay from neutron capture following the initial neutrino interaction. It becomes possible
to apply software triggers to select these special events because all PMT hit information isrecorded in the new system. A much larger dynamic charge range also raises the saturationlevel, which benefits higher-energy neutrino studies.
The ATM boards in the ID are replaced by QBEE modules, which stands for “QTC
(charge-to-time converter) Based Electronics with Ethernet”. Each QBEE board collects
analog input signals from 24 PMTs. The QTC outputs a digital pulse with the leadingedge defining the timing information and the pulse width defining the total charge. The
pulse information is read out to the on-line computers by Ethernet cables using the TCP/IP
protocol. This greatly improves the transfer time for data and eliminates bottlenecks whichwere sometimes caused by the (custom-built) TKO readouts in the original DAQ. The same
QBEE electronics were used for the OD PMTs as well. Although the OD boards were new,the QTC readout principle was carried over from the previous OD setup.
The analog sum hardware trigger has been replaced by a software trigger that records all
PMT hits. Event windows are defined by a 60MHz master clock which continuously applies

20

Figure 2.7: A schematic of the new DAQ system. Taken from Ref. [33]a periodic 17µs trigger to the QBEEs. The data readout rate which is necessary to save
all hits is 1.5MB/s/QBEE, which is easily achieved by this system. A schematic diagramof the new system with the QBEEs and master clock is shown in Fig. 2.7. A schematic
showing the new periodic software trigger system is shown in Fig. 2.8.
The timing resolution of the QBEEs at the 1 p.e. signal level is 0.3 ns. This timing
resolution is much better than the 2.8ns timing resolution for the 50cm (20 inch) ID PMTs
at the same charge level. A plot of the timing resolution for the QBEE and the PMT as afunction of input charge is shown in Fig. 2.9. The QBEEs also feature a much improved
dynamic charge range. There are three input channels per chip, each with three gain ratios1
of 1 : 17 : 49
, combining to give a dynamic range from 0.2 pC to 2500 pC. A timing chart

for the process of converting a raw PMT current signal to a digital output signal by the
QTC is shown in Fig. 2.10. As shown the built-in discriminator is used to produce a QTC
out signal with the hit timing (leading edge) and charge (pulse width) information encoded.
The process time is about 1 µs/hit.
The QBEEs also feature the high-speed processing necessary for recording all the data
in a supernova burst event. A high-rate test was performed by feeding the same signal into
the QBEE and no data loss was seen up to frequencies of 130 kHz, as shown in Fig. 2.11.
This represents a factor of 100 improvement over the original DAQ readout.

21

Figure 2.8: A schematic of the new periodic software trigger. Taken from Ref. [33]Figure 2.9: A schematic of the timing resolution as a function of input charge for a QBEE
and ID PMT. Taken from Ref. [33]22

Figure 2.10: Timing chart for QTC processing of raw PMT signal. Taken from Ref. [33]Figure 2.11: The high-speed processing rate of the QBEEs. Frequency of the input signal
is on the horizontal axis. Taken from Ref. [33]23

Trigger

SK-I

SK-II

SK-III

SK-IV

SK-IV Average Rates

SLE-260 mV-110 mV-186 mV-196 mV

3260 HzLE-320 mV-152 mV-302 mV-324 mV

36 Hz

HE-340 mV-180 mV-320 mV-300 mV

11 Hz

OD

19 hits

19 hits

19 hits

19 hits

3 Hz

Table 2.1: A summary of each analysis trigger threshold for each detector period, alongwith the average trigger rates for SK-IV.

2.5.3

Trigger

The SK-I,-II, and -III detector configurations used a hardware hitsum for the analyses
triggers. With the new SK-IV electronics, a new software trigger is used that selects trigger
events after saving all hit information. In both setups, four analyses triggers were defined:
three ID triggers (super-low energy [SLE], low energy [LE], and high energy [HE]) and an
OD trigger. The SLE and LE triggers are used mainly for solar neutrino events and SN
relic events, while the HE trigger is mainly used for atmospheric neutrino and high-energy
cosmic ray muon events. The OD trigger threshold is defined within a 200 ns coincidence
window. The SLE threshold is limited from below by backgrounds such as rock wall and
PMT glass radioactivity which increase as the threshold decreases. All trigger thresholds
were decreased during SK-II to reflect the reduced photocathode coverage. A summary of
the thresholds and approximate trigger rates for each analysis trigger in each run period
are shown in Table 2.5.3.

2.6

Air Purification

Radon in the mine air (from the U/Th decay chain) can decay and the products can become alow-energy background for the experiment. It can also be a health hazard to people working
in the mine, at more extreme levels. Natural levels of radon vary throughout the year,
ranging from about 200 Bq/m3 to 3,000 Bq/m3 depending on the season. Multiple actions

24

Figure 2.12: Schematic view the radon-free air system. Taken from Ref. [5]are taken to minimize this. The mine cavity surface is coated with MineGuard R (made by
Urylon), a radon-tight paint. Additionally, a radon-free air system is housed just outside theentrance to the mine. An air blower takes radon-free air from outside the mine and pumps it
into the mine at a slight over-pressure at a rate of about 10 m3 /min. This system maintains
a Radon level of about 40 Bq/m3 in the dome and control room areas. A schematic of the
air purification system is shown in Fig. 2.12.
2.7

Water Purification

Maintaining very pure water inside the tank is crucial to the successful operation of awater Cerenkov detector. A water purification system is used which continuously recycles
the water in the tank. It is a closed cycle system which processes water at a rate of
approximately 60 metric tons per hour. Water is passed through the same purification
system when refilling the detector as well, but with raw mine water used as the source.
Water is supplied to the bottom of the ID, and removed at the top and bottom of the OD,
as well as the top of the ID. This flow is chosen to minimize the stagnation of water in thetank. A schematic of the water system is shown in Fig. 2.13.
Water from the SK tank is initially passed through 1 µm mesh filters to remove particulate matter. This prevents reduction of transparency for Cerenkov photons. The water

25

Figure 2.13: Schematic view the water purification system.

also passes through two heat exchanges, which cool the water to about 13 ◦ C. This reduces
the dark noise rate of the PMTs and also suppresses growth of bacteria which can further
reduce water transparency. A UV sterilizer is used to kill any remaining bacteria.
Rn-reduced air is dissolved into the water, and a vacuum degasifier is used to remove
dissolved radon and oxygen gas. Radon gas is removed with ∼ 96% efficiency at this stage.
An ultra filter removes particles down to about 10 nm in diameter. Finally a membranedegasifier removes remaining radon gas dissolved in the water with more than 80% efficiency.
The water system underwent an upgrade in 2012 (the upgraded system is shown in the figure)to supply colder water to the SK tank. The water temperature exhibits an anti-correlation
with the transparency length, so colder water is better for reducing absorption of Cerenkov
light. A new He-4 heat exchanger and membrane degasifier were installed in the system,
shown in the red box in the schematic. The temperature of the water in the tank is slightly
above 13 ◦ C.

26

Chapter 3
CALIBRATION
3.1

Introduction

It is essential for the detector response to be well known so that the data can be accurately
analyzed. Many types of calibration are performed in the detector to measure its response.
In these calibrations a known source of light is used to simulate events in the detector andthe output is compared to expectations. The following are some of the main calibration
tests.
3.2

Water Calibration

3.2.1

Laser absorption and scattering

The ability of SK to operate as a water Cerenkov detector relies on having pure waterwith well-understood transparency, absorption and scattering properties. First, the water’s
index of refraction was modeled by a complicated function[110] which calculated the effective
index as a function of wavelength as shown in Fig. 3.1. As seen in the figure the index ofrefraction varies between 1.33 and 1.36 for the Cerenkov radiation wavelengths of interest.
Additionally, a propagation model for light in the tank water including the effects ofabsorption and scattering was developed. These effects reduce the intensity of light as it
travels in the water. Intensity as a function of distance traveled is given byL

I = I0 e − W T .

(3.1)

In this expression L is the distance traveled by the light, W T is the water transparency
parameter (the inverse of the attenuation length) including scattering and absorption effects,
I0 is the original light intensity, and I is the intensity after traveling a distance L. The goalof the calibration described here is to estimate the water transparency parameter W T . This

27

Figure 3.1: The index of refraction of the SK tank water modeled as a function of lightwavelength.

parameter is defined as

W T (λ) =1.
αsctt (λ) + αabs (λ)(3.2)

The parameters αsctt and αabs in this expression have units of 1/m. Light scattering
is modeled to include both Rayleigh and Mie scattering. Rayleigh scattering applies towavelengths of light that are much larger than the size of the scattering object (a watermolecule, for example), whereas Mie scattering applies when the two length scales arecomparable to each other (for example, light scattering from a large impurity in the water).
The scattering parameter can then be divided into these two components, as given by

Latten = W T (λ) =1αRay (λ) + αM ie (λ) + αabs (λ)

.

(3.3)

However, for actual calibration measurements the scattering parameters are most easily separated by their angular distributions. Fig. 3.2 shows the scattered light angular
distribution from symmetric (left) and asymmetric (right) scattering, respectively. Light
from Rayleigh scattering obeys the left distribution, and that from Mie scattering is approximately described by the sum of both distributions. Thus, the water transparency

28

Figure 3.2: Angular distribution of scattered light from symmetric (left) and asymmetric(right) scattering.

parameter can be rewritten as

W T (λ) =1
αsym (λ) + αasy (λ) + αabs (λ)

.

(3.4)

Each of these α coefficients are found to have λ−4 dependence to leading order, and theabsorption parameter has an additional term which dominates at longer wavelengths (it has
a complicated functional form that is plotted in 3.5).
For this calibration measurement, a set of eight laser light injectors which send laser
light of tunable wavelengths into the SK tank via an optical fiber were used. A schematic
showing the positions of these injectors in the detector is shown in Fig. 3.3.
Wavelengths of 337, 365, 400 and 420 nm produced by a nitrogen laser (337 nm wavelength) combined with wavelength-shifting laser dyes for were used in this system for SK-I,-II and -III. For SK-IV, the system was updated, replacing the nitrogen/dye laser systemwith a laser diode (LD) system. The new wavelengths in the SK-IV LD system were 337,
375, 405 and 440 nm. An event display for a sample laser injection event is shown in Fig.
3.4 As shown in the figure the laser light is not passed through a diffuser, so the majority of
the light causes PMT hits at the bottom of the detector. The other PMT hits in the eventare caused by either backgrounds (PMT dark noise, for example) or scattered light.
The hit and timing distributions of the laser light data were used to tune the waterparameters until sufficient agreement was seen with MC simulation events. The results ofthe measurements of the scattering and absorption parameters are shown in Fig. 3.5. The
inverse of the overall water transparency parameter is shown by the black line in both plots.

29

Figure 3.3: The positions of the eight laser light injectors in Super-Kamiokande. “Old”
refers to the sole original injection point used in SK-I.

Figure 3.4: A typical laser light injection event (injected from the top of the tank).

30

Figure 3.5: Measurement of the water parameters for scattering and absorption.

An attenuation length of approximately 100m was achieved for SK-IV.

3.2.2

Cosmic ray muons

Downward-going cosmic ray muons also serve as a calibration tool for measuring watertransparency. These particles are the majority of HE triggers in the SK data and occur at
a rate of a few Hz. They are typically minimum ionizing particles which deposit energy(mostly through ionization) steadily at a rate of about 2 MeV/cm, independent of theirenergy. One advantage of using cosmic ray muons over laser events is that the measurementwith cosmic ray muons is naturally averaged over the spectrum of Cerenkov radiation,
instead of a monochromatic distribution. The detected number of photo-electrons is givenby

Q = Q0f (θ) − L le atten .
l(3.5)

where Q0 is a constant, f (θ) is the PMT acceptance (geometric photosensitive area), lis the photon path length, and Latten is the attenuation length. A schematic showing thepath length of the photons and their arrival angle θ at the PMT is shown in Fig. 3.6.

31

Figure 3.6: Cosmic ray muon water transparency measurement schematic.

The detected charge is measured over the angle θ and plotted as a function of the path
length that the light took to reach the PMT. This plot is then fit for the normalization factor
Q0 and attenuation length Latten . The attenuation length measured from these studies isfound to vary between 80 and 90m for SK-IV. The plot for SK-IV is shown in Fig. 3.7. Thisis comparable to the measured values found for SK-I,-II, and -III (see Ref. [109] for moredetails).

Figure 3.7: Thru-mu measurement of attenuation length for SK-IV.

32

3.2.3

Water quality position and time dependence

The SK water quality has a position and time dependence because of the flow from the water
circulation system. For the position dependence, a top-bottom asymmetry in the water ismeasured by a monthly calibration with a Ni-Cf source positioned at three heights (Z =+12m, 0m and -12m) in the detector. The PMT hit rates are found to be a few percent
lower near the top of the tank than near the bottom [38]. A xenon lamp connected to a
scintillator ball is also used as a complementary measurement of the position dependence.
These light sources are regularly used to monitor the time dependence of the water quality
as well. More details can be found in Ref. [12].
3.3Charge calibration

In addition to the water properties, a thorough understanding of the response of the PMTs
to the Cerenkov radiation photons is necessary for precise event energy reconstruction andMC simulation. As described in the previous chapter, PMTs amplify photo-electron signals(the amplification factor is the “gain” of the PMT) and send the amplified signal to thereadout electronics. The following sections describe calibration studies done to understandthe charge and timing properties of the PMT signals.
3.3.1

HV determination

The charge response of the PMT for a given amount of light must we well understood forsufficiently accurate energy reconstruction. The response can be quantified by the gain andquantum efficiency (QE) values. The gain of each PMT is controlled by the amount of HV
supplied to the dynodes. Generally, the relationship between supplied HV and gain for a
PMT is given by

Gain ∝ αHVβ(3.6)

where α and β are unique to each PMT. The gain calibration includes measurements ofaverage and individual PMT gains, respectively. PMT gains can in general be determinedfrom the single photo-electron response to the Nickel γ-ray calibration source. The general

33

principle is to first adjust PMT HV values so that the output charge, Qobs , from each PMT(after corrections for the number of p.e.s detected) is uniform across all PMTs. The output
charge can be calculated by the expression

Qobs ∝ Nphoton (i) × QE(i) × Gain

(3.7)

where Nphoton represents the number of photons striking the photo-cathode, QE represents the quantum efficiency, and i represents the ith PMT. After adjusting the HV toachieve sufficient uniformity, additional independent measurements were made of the PMT
gains and quantum efficiencies. These serve as necessary inputs for the MC simulation of
the detector.
A group of “standard PMTs” were calibrated so that their output charges matched eachother to within 1.3% before installation in the ID. This so-called “pre-calibration” was donein 2005 before the start of SK-III data-taking. In total 420 PMTs were calibrated by placing
them in a dark box with a Xe lamp source passed through a UV filter (passing only lightwith λ < 400nm) and then a scintillator ball. The HV of each of these PMTs was adjusted
so that the charge output matched that of a reference avalanche photodiode (APD).
The locations of the standard PMTs in the ID are shown in Fig. 3.8. All ID PMTs weregrouped together (25 groups in total, such that 9-12 standard PMTs were included in eachgroup) and the charge responses to the calibration light sources were matched to the averagecharge response of the standard PMTs within that same group. This reduced uncertainties
due to water quality and photon propagation for PMTs located in different areas in the ID.
Light flashed from the Xe lamp was passed through the UV filter and scintillator ball. The
scintillator ball was positioned in the center of the talk and rotated around its vertical axis
to account for any spherical asymmetry. Each PMT received roughly 50 p.e.s per flash oflight. A schematic of the setup showing the scintillator ball and a sample group of PMTs
is shown in Fig. 3.9.
The corrected charge, Qcorr , observed at each PMT was then calculated from the observed charge, Qobs with the expression

34

Figure 3.8: The installed locations of the standard PMTs (shown in red) within the ID.

Figure 3.9: The scintillator ball setup for the PMT HV determination. A sample groupingof PMTs calibrated together is shown in blue.

35

Figure 3.10: The ratio of each PMTs corrected charge to the target charge (mean correctedcharge of standard PMTs) before (left) and after (right) HV adjustments.

Qcorr =

Qobs × r2
faccept (θ)(3.8)

where r is the distance from the scintillator ball to the PMT and θ the angle from the
PMT to the scintillator ball. The HV for each PMT was adjusted so that the correctedcharge equaled that of the average of the standard PMTs within the same group. The
results of this calibration can be seen in the contrast of the left and right plots in Fig. 3.10.
The figure demonstrates that after the tuning the variation of observed charge relative
to the standard charge decreased significantly. The observed charges for all PMTs wereequal to each other within 1.3% after this calibration procedure, compared to almost 10%
beforehand.
3.3.2

Absolute gain calibration

An overall “absolute gain” calibration was also updated for SK-IV and compared to previous detector configurations. A photon source called the “nickel” source was used for thismeasurement. It was made up of a source of Californium (252 Cf) surrounded by a sphericalmixture of NiO (6.5 kg) and polyethylene (3.5 kg). The spontaneous fission of the 252 Cf
source produces neutrons which are then thermalized (brought to thermal equilibrium) by

36

Figure 3.11: The single p.e. distributions for SK-III (black) and SK-IV (red, adjusted byscaling factor) shown for two different horizontal scales (left and right).

the polyethylene. After thermalization the neutrons can be captured by the nickel, producing gamma rays of about 9 MeV through the process Ni(n,γ)Ni (where the initial and final
Ni elements are different isotopes). Each of these gamma rays subsequently produces about
100 photons detected by the PMTs. The energy of these events were low enough that theprobability of PMTs receiving detecting only a single p.e. was very high. The nickel source
was positioned near the center of the SK tank for the measurement.
The resulting 1 p.e. distribution for SK-IV was similar to that of SK-III, with somedifferences in the low charge tail region. The QTC chips on the new SK-IV QBEE boards
had a poorer resolution in the low charge per p.e. region than that of SK-I through SK-III,
so the SK-IV single p.e. distribution was adjusted to overcome this. The SK-IV single p.e.
data were multiplied by a factor of 0.844, which gave agreement between the SK-III and
SK-IV peaks. The SK-III and corrected SK-IV distributions are shown in Fig. 3.11.
From this measurement the conversion factor between p.e. input and final PMT charge
output can be determined. This is the “gain” factor that was the goal of the measurement.
This scaling factor is given the name “pC2pe”. It represents the amount of charge in unitsof pC generated as the output signal from 1 input p.e., or in other words the gain (notethat this value is consistent with the rough estimate of the gain of the PMTs of ∼ 107

37

Figure 3.12: A schematic of the setup used to measure the relative PMT gains.

quoted in the previous chapter). For SK-IV the global pC2pe value was found to be 2.658.
This differed from the corresponding SK-III value, 2.243, by the scaling factor which wasmentioned above.
3.3.3

Relative gain calibration

The relative gain calibration describes the individual PMT response relative to the average
response described in the previous section. The setup for this measurement is shown in Fig.
3.12. A dye/nitrogen laser connected to a diffuser ball was used as the light source. As withall other calibration measurements the light source was partially redirected to a dedicatedtrigger PMT which indicated exactly when the source was flashed.
The relative gain for each PMT can be determined from the ratio

Gain(i) ∝Qobs (i)Nhit (i)(3.9)

where i is the i-th PMT, Qobs is that PMTs observed charge, and Nhit is the numberof photons reaching the PMT cathode. Thus, the calibration measurement should give
both values in the RHS of this equation. This value was then normalized by the meanvalue over all PMTs, giving the relative gain. The measurement was taken twice, with the38

Figure 3.13: The distribution of relative gains for all ID PMTs after normalizing to mean
1.

light set to low and high intensities, respectively. This was done to get around the fact that
individual PMT hits in a low-intensity setup could be from background and therefore cannot
be guaranteed to be coming from the light source (this was not a problem in the absolute
gain measurement described previously because the data were not studied separately for
individual PMTs.) The low intensity produced 0.3 [hit/event/PMT] and the high intensity
produced 50 [pe/event/PMT], on average. The charge observed by each PMT for the highintensity run was the numerator for 3.9, while the number of hits in the low intensity run
was the denominator. The distribution of this quantity divided by the mean of the value
from all PMTs is shown in Fig. 3.13. This information served as input to the detector MC
simulation.

3.3.4

QE measurement

The quantum efficiency, or QE, was defined as the fraction of photons incident on the PMT
cathode which produced a p.e. detected by the dynodes and subsequently read out by theelectronics. The nickel gamma-ray source introduced previously was used to measure thehit rate of each PMT for this measurement. As mentioned above, each gamma-ray emission
resulted in about 100 detectable photons, so the probability of PMTs measuring multiple

39

hits was very low. In this case the QE is proportional to the measured hit rate. The PMT
hit rate was adjusted by MC predictions of the number of photons expected to reach the
PMT, to account for the remaining uncertainty in the number of photons reaching each
PMT.
Several selection cuts applied to the data from these events reduced unwanted backgrounds. Only events with reconstructed vertices near the nickel source (¡ 400cm away,
because the gamma ray is expected to produce photons near its emission point), time differences > 5µs after the previous event (to remove noise), and passing other various selection
cuts were selected. The resulting distribution of PMT hit rates from this calibration isshown in Fig. 3.14. In this figure the blue dots represent the calibration data and the red
dots the detector MC simulation. The top plot shows the hit rate distribution of barrel
PMTs with the vertical position (z) along the horizontal axis. The bottom left and rightplots show the distribution for top and bottom PMTs, respectively, as a function of squareof the distance from the center (cm2 ). As seen in these plots the hit rate differs from 1
for some PMTs near the edges of the tank. This is from photon reflections off of neighboring PMTs and water scattering. However, the figure also shows that the MC simulation
correctly reproduces these tank edge PMT effects.

3.3.5

PMT calibration consistency check

So far calibration of the PMT HV values and independent measurements of the gain and
QE for individual PMTs have been described. A further measurement was made to test
the consistency of these PMT calibration steps. The “diffuser ball” source was placed at
the center of the tank and used as the light source. The charge observed by each PMTfrom this source was corrected for the distance and acceptance angle relative to the source.
This corrected charge was then divided by the QE and relative gain factors as measured
in previous sections. The resulting distribution of normalized, corrected charge is shownin Fig. 3.15 for the PMTs divided by their locations. The plot was mostly flat around
1.0, confirming the consistency of the calibration, while indicating a very slight top-bottom
asymmetry of less than 1%.

40

Figure 3.14: The distribution of PMT hit rates for the barrel (top plot), top (bottom left)
and bottom (bottom right) PMTs. The blue is the calibration data and the red the MC
simulation data.

Figure 3.15: The corrected charge divided by the QE and gain for PMTs divided geometrically.

41

3.3.6

QBEE Charge non-linearity

The linearity of the charge response of the PMTs with the new SK-IV QBEE electronics was
also measured. The linearity was defined as the ratio of the observed charge to the expectedcharge across as a function of increasing charge. The new QBEE electronics had a wide
dynamic range, capable of processing an input range of 1-1000 p.e.s. The response linearity
in this entire range was investigated. In the SK-I,-II and -III detector configuration, theoriginal electronic readout boards (ATMs) exhibited a ∼ 6% systematic upward shift in the
output charge introduced in the analog-to-digital conversion process. This was corrected
for by applying a scaling function to the output charge. The measurement described here
examined the existence of a similar effect in SK-IV. The nitrogen laser connected to the
diffuser ball was used as a light source.
First special monitor PMTs, whose outputs were connected directly to ADCs with wellknown response characteristics, were chosen as references. The monitor PMTs were illuminated by the diffuser ball at various gain settings and various distances from the lightsource (to get data at a range of charge outputs), and their linearity as a function of charge
was confirmed. After confirming their linearity, their response to the light was used in eachcalibration event to determine the expected charge for all other PMTs.
Next, normalizing runs where the diffuser ball light inside the tank was adjusted so thatthe output of each PMT was around the 10 p.e. level were completed. This establisheda base output charge ratio between each PMT and the monitoring PMT for that run, toaccount for geometric acceptance and location variation. To check each ID PMT/QBEE’slinearity, the laser source was flashed many times, and the mean output charge per event
Qmean was recorded for each PMT. The ratio of this mean output charge to the expectedcharge should be a constant value of 1.0 independent of the charge level for a perfectly
linear PMT/QBEE response. The expected charge Qexp was determined by the monitor
PMT Qmon and the normalizing run data and was given by the expression(i,n)

m
Q(i,m)
exp = Qmon ×

Qmean
(n)

(3.10)

Qmon

where m is the calibration run index, i the PMT index, and n the normalizing run

42

Figure 3.16: The nonlinearity as a function of expected charge for ID PMTs before (left)
and after (right) multi-hit correction. The red points show the average of the black points
at the corresponding horizontal charge bin.

index. The result of the ratio as a function of expected charge for all PMTs is shown in theleft-hand plot of Fig. 3.16. The linearity value varies from 1.0, especially at higher expected
p.e. levels.
This variation is caused by two separate effects. The first comes from the QBEE andoccurs when the input signal is read in by more than one channel (of the three possible channels) on the QBEE board. The second comes from the PMT itself, for which nonlinearityof about 2% of greater appears at high expected charge levels. The latter effect cannot beeasily corrected, but for the QBEE multi-hit effect a correcting function was applied andthe corrected result is shown in the right-hand side of Fig. 3.16. As shown the linearity is
maintained to within ±1%, the desired level for data-taking. The data in the high charge
region in these plots is bounded from above by QBEE saturation, which sets in at around
1200 p.e.
3.4

Timing calibration

The relative timing of the PMTs is also important for accurate reconstruction of event
vertices and directions. The time elapsed between Cerenkov radiation emission and thecomplete readout of event information includes the Cerenkov photon time of flight (the43

Figure 3.17: The setup used for the timing calibration measurement.

time for the photons to travel from their emission point on the charged particle’s track
to the PMT surface), the time for the photoelectron avalanche to pass through the PMT,
travel the length of the cable from the PMT to the QBEE, and finally be collected and readout by the QBEE boards. The time interval between the photon hitting the PMT and theglobal trigger not only differs for each PMT, but also exhibits a dependence on the chargeof the event (termed the “time walk”) due to the fact that events with more charge pass
the QBEE discriminator threshold level sooner than events with less charge.
The goal of this calibration is to measure these time intervals for each PMT as a functionof charge (from 0 p.e. up to QBEE saturation) with a calibration light source, model the
distribution in the time-charge plane (the “TQ Distribution”) by a polynomial function (the“TQ Map”), and then make a correction table which applies to each PMT so that the hit
timing is a constant value independent of event charge. The N2 /dye laser is pulsed in shortbursts, split and redirected to a monitor PMT for event triggering, passed through a filter
which varies the light intensity and then spread out in the tank via the diffuser ball. A
schematic is shown in Fig. 3.17.
A sample TQ distribution for a PMT channel before corrections from the TQ Map is44

Figure 3.18: A typical TQ distribution before correcting with TQ Map table.

shown in Fig. 3.18. Note here that smaller values correspond to hits later in time because
the axis represents time elapsed between the hit and the eventual global trigger. For eachcharge bin, Qbin , on the horizontal access, scattering and reflection causes a vertical spreadof time intervals. To make the TQ Map, the hit times in each bin are fit to a Gaussian tofind the peak and standard deviation. The values of these peaks times for each charge bin
are used to fit a seventh-order polynomial to all data points.

3.5

Energy scale

3.5.1

LINAC

Super-Kamiokande operates a LINAC on top of the tank to create electrons of known
energies and sends them into the tank via a pipe whose endpoint position can be adjusted.
This calibrates the absolute energy scale for lower energies (solar neutrinos, etc.) because
the LINAC electrons have energies in the range 5-16 MeV. The accuracy of the absoluteenergy scale is found to be better than 1% with this measurement. This section will focus
on calibration studies done at higher energies, so for more details on this the reader canrefer to Ref. [51].

45

Figure 3.19: The data and MC energy distributions for stopping muon decay electron events.

3.5.2

Decay electrons

Decay electrons for stopping muons exist in large numbers in the SK data and are usefulfor checking the energy scale in the tens of MeV range. Fig. 3.19 shows the comparison
between the measured energy spectrum and the Monte Carlo prediction for these events inunits of MeV for SK-I data. The data and MC agree well, with the means differing by only
about 2%. These data are also used to monitor the stability of the energy scale over time,showing stability within ± 1%. The SK-IV time variation of the mean is shown in Fig. 3.20.

3.5.3

Stopping muonsParticles emitting Cherenkov radiation have an opening radiation angle which is a known
function of particle momentum. Estimating this opening angle allows one to calculate particle momentum independent of the PMT charge measurement and check for consistency.
Stopping muons are used in this calibration. The opening angle function approaches a limiting value as the particle momentum increases, so only stopping muons with momentum lessthan 400 MeV/c are useful here. The ratio of the measured momentum to the momentum
calculated via Cherenkov angle for MC and data are both taken, and the ratio of theseratios is calculated (a double ratio). This double ratio is plotted for SK-I in Fig. 3.23. As

46

Figure 3.20: The time variation of the decay electron momentum mean from SK-IV data.

Figure 3.21: The ratio of MC to data of the ratio of measured momentum to inferredmomentum from stopping muon opening angle.

seen in the figure the double ratio is consistent to within ± 1% for the momentum range
between 0.2 and 0.4 GeV/c.
For higher energy stopping muons, the momentum loss is proportional to the track
length, or range (this is because high energy muons lose energy at an approximately constantvalue per unit track length). The range is calculated by assuming the end point to be
the location of the vertex of the decay electron event. The ratio of the momentum loss
to the range as a function of range for a sample of stopping muons from an initial SK-I
measurement is shown in Fig. 3.22. The data and MC have good general agreement, with

47

Figure 3.22: The ratio of momentum loss to range for stopping muons for both data andMC in initial SK-I calibration.

the ratio between the two staying almost constant over the measured values.
The SK-IV time variation of the momentum to range ratio is plotted in Fig. and is seen
to be constant within ± 1%.

3.5.4

Neutral pions

Atmospheric neutrinos which interact in the detector via the neutral current often produce
a single π 0 in the final state. These neutral pions are easily detected by looking for thesignature of their decay to two gamma rays (∼ 99% branching ratio.) These gamma rays
create two showering (e-like) ring patterns. The invariant mass of 2-ring e-like events isreconstructed and the distribution is compared to the expected peak from π 0 decay. The
resulting distributions of data and MC are shown in Fig. 3.24
Though the peak at 140 MeV/c is shifted relative to the value of the π 0 mass (135
MeV/c), the SK data and MC agree well with each other. The shift is caused by vertex
reconstruction bias and other nuclear effects [12]. Most events in this sample have visibleenergy less than 1 GeV, because efficiency of cleanly detecting the two rings decreases
with increasing π 0 energy from overlapping rings (decays with small opening angles) andasymmetric decays producing only one Cherenkov ring.

48

Figure 3.23: The time variation of the ratio of momentum to range for high-energy stoppingmuons for SK-IV.

Figure 3.24: The invariant mass distribution of data and MC for 2-ring e-like events. Taken
from Ref. [12].

49(MC-DATA)/DATA (%)

SK-IV
8
6
4
2
0
-2
-4
decay electron
neutral pion
stop muon (sub-GeV)
stop muon (multi-GeV)

-6
-8
10

2

10

3

4

10
10
momentum (MeV/c)

5

Figure 3.25: A summary plot showing the data and MC agreement for all energy scalecalibration measurements as a function of energy.

3.5.5

Energy scale calibration summary

A summary of the agreement between data and MC for all of the above-mentioned energy
scale calibrations in SK-IV is shown in Fig. 3.25. As shown in the figure the agreement
between data and MC is within ± 2% across a wide range of energies including the energies
of atmospheric neutrinos analyzed in this work.
3.6

OD calibration

The outer detector is mainly used as a veto for incoming cosmic ray muons and to identify
the exit points of outgoing charged particles. Therefore, the timing and charge precision
requirements are not as strict as in the ID. Calibration of the OD was most recently performed at the start of SK-IV upon installation of the new electronics, and that calibrationis briefly described here. At the time of the SK-IV upgrade, OD PMT gains were lowered(200 V for new PMTs and 400 V for old PMTs) to reduce crosstalk with ID electronics.
Light from the same N2 /Dye laser used in ID calibration is split and redirected to the OD

50

Figure 3.26: Example 1-p.e. distributions for old (top) and new (bottom) OD PMTs before(left) and after (right) SK-IV gain adjustment.

via optical fibers as one calibration light source. The fiber ends were mounted to various
positions along the top, barrel, and bottom of the outer wall of the detector.
The single-p.e. distributions are in general different for old and new OD PMTs. Fig.
3.26 shows an example 1-p.e. distribution for an old (top) and new (bottom) PMT for both
before (left side) and after (right side) the SK-IV gain adjustment.
The OD PMT gain measurement was performed by both the laser light and dark ratemeasurements. The laser light was flashed at low light levels to determine charge per p.e.
for each PMT. The mean number of output ADC counts from this measurement was used as
the single p.e. level. In addition to this the dark rate noise hits outside the trigger window,
which have a high probability of being sinle p.e. hits, were also used to measure the same
ADC count value after pedestal subtraction. The ADC measurements from this were foundto agree to the laser measurements to within 10%, so this method was primarily used. For
SK-IV, typical counts-per-p.e. ranged from 2-6 for OD PMTs.
The measured charge as a function of input charge for a wider charge range was alsomeasured with the laser. The curve from this measurement is shown in Fig. 3.27. A
saturation effect is seen at high p.e. levels (hundred’s of p.e.s upward) and this effect is also

51

Figure 3.27: Measured versus input charge for old (red) and new (black) OD PMTs fromlaser measurements.

modeled in the MC simulation code.
The goal of the OD timing check was to confirm that the timing offsets relative to the
ID PMTs and to the other OD PMTs was not too large (≤ few ns). The cable length foreach PMT, ranging from 70 to 78 m, was taken into account for this study. The offsets were
determined by flashing the laser in both the ID and OD with identical cable and diffuserball, and this confirmed that the offset was sufficiently small. Cosmic ray muons were also
used to check the agreement by comparing the timing difference for each pair of ID and OD
PMTs that was closest to the muon track at any given point along the track.
Finally, a new calibration study was introduced in SK-IV to check the relative efficienciesof the old and new OD PMTs. A laser diode and diffuser ball was inserted into the barrel
portion of the OD and flashed at 100 Hz at various heights. Old and new PMTs weregrouped into pairs with each other such that the distance and relative position (above or
below) for each PMT in the pair was identical (though not necessarily from the same diffuserball height). With these conditions 15 pairs of PMTs were chosen. For each PMT and eachrun the occupancy, defined as

52

Figure 3.28: A photograph of the NIM modules used in the SN Burst test which the authorof this work arranged to simulate a SN signal.

occupancy =ontime − N of f time
Nhits
hits
Ntriggers(3.11)

was measured and then the ratio taken of old and new PMTs for comparison. The
average of the ratio for all 15 pairs was used as the tuning input parameter for the detector
simulation, which is described in more detail in the following chapter.
3.7

Supernova burst test

With the installation of the new SK-IV front-end electronics the SK detector’s ability toprocess high-rate data from nearby supernova events was improved. A new Supernova
Burst test system was constructed and tested prior to SK-IV QBEE installation. The
system consisted of NIM logic modules arranged to simulate a supernova burst signal. Thiswas done by using a NIM pulse generator in combination with gate and delay modules togenerate sequences with as many as ∼ 28 million pulses over a 10-second time window whosetime profile simulated that of neutrino events in SK from a galactic supernova. A pictureof the NIM module setup is shown in Fig.3.28.

53

Figure 3.29: A sample of the triggering profile during the SN Burst test taken during SK-IV
test installation.

These pulses were fed into an LED light system whose 375 nm output was directed into
the SK tank from above via an optical fiber. A diffuser ball was attached to the end of the
optical fiber to spread the light throught the tank. An example of the triggering pattern that
occurred in the detector is shown in the time profile in Fig. 3.29. This profile closely mimics
the expected supernova neutrino curve. The tesitng system was first installed during the
SK-IV test installation in summer of 2008 and has been in place throughout the durationof SK-IV.

3.8

GPS Timing System

SK incorporates a GPS time system that timestamps each trigger with absolute accuracy on
the order of 50 ns. The GPS time system is critical for the long baseline experiments K2K
and T2K, where SK event timestamps are compared to neutrino beam spill timestamps to
identify beam-induced events. In addition, when a supernova occurs, the GPS time data
will be critical for coordination with other observatories and analysis of the detailed time
structure of early neutrino production. The GPS system was built at UW and consisted of

54

two independent GPS receivers, a rubidium atomic clock, and a Local Time Clock (custom
electronics board) that integrated receiver and atomic clock data and interfaced to the OD
DAQ system through a VME crate.

55

Chapter 4
ATMOSPHERIC NEUTRINO MONTE CARLO

Experimental particle physics analyses generally proceed by comparing final data samples to predictions from Monte Carlo which simulates the models being analyzed as wellas the response of the detector to the interactions under study. Atmospheric neutrino oscillation studies at SK use this same procedure, and the goal of this chapter is to describe
in detail the Monte Carlo samples produced for these studies. The three main pieces ofthe Monte Carlo which will be described are the atmospheric neutrino flux, the neutrino
interactions, and the detector response. Monte Carlo sample are generated for each uniquedetector phase (SK-I, SK-II, SK-III and SK-IV). The analyses in this work feature updated
MC samples for the first three SK running periods, as well as SK-IV samples which are used
for the first time in physics analyses.

4.1

Atmospheric neutrino flux

The Monte Carlo begins with a simulation of the atmospheric neutrino flux starting withprimary cosmic rays. At low energies (. 1 GeV) the primary flux is well-known, but theearth’s magnetic field and solar activity have large effects on cosmic rays in this region. The
Earth’s magnetic field creates a rigidity (momentum/charge) cutoff for incoming charged
particles. Additionally, the solar wind reduces cosmic ray flux via its magnetic field. Periodsof peak solar activity create flux reductions by as large as a factor of two. At energies ∼10 GeV and upwards, the flux becomes less well-known, but the influence of geomagnetic
and solar effects diminishes. The accuracy of cosmic ray flux measurements below 100 GeV
primary energy is about 5% and gets worse with increasing energy value[65]. The flux
model used for the numbers and figures in this work is the Honda 2011 model[49], which
is an update to the model used in previous SK analyses [48, 47]. The monte carlo eventswere initially generated with the Honda 2006 flux[48] (at the time of production the most

56

Figure 4.1: The zenith angle dependence of the Honda 2011 flux for atmospheric neutrino
flux at Kamioka averaged over azimuthal angle, taken from [49].

recent available model) and then reweighted afterwards to match the Honda 2011 model.
The solar wind and geomagnetic effects are accounted for in the flux calculations. The local
magnetic field strength near Kamioka is larger than the global average, resulting in a slightup-down asymmetry in the neutrino flux.
The Honda flux assumes the primary cosmic rays are dominantly H nuclei with small
contributions from He and other nuclei. The subsequent hadronic interactions are modeled
to produce an atmospheric muon flux. A combination of two hadronic interaction models isused in the updated 2011 model. At low energies (atmospheric muon Energy < 32 GeV/c)a modified version of the JAM model[40], found to have good agreement with atmosphericmuon flux data, is used for the first time. A modified DMPJET-III model[103], used in theprevious versions of this calculation, is used for primaries that produce atmospheric muons
above 32 GeV/c. The main difference as a result of the new low-energy model is a small
increase in the atmospheric muon flux below 1 GeV/c. The model is three-dimensional,
meaning that the secondary particles are not confined to travel in the same direction as
their parent particle. The zenith-angle dependent aziumuth-averaged neutrino flux for threedifferent energies is shown in Fig. 4.1.

57

Figure 4.2: The energy dependence of flux for each neutrino type (left) and important flux
ratios (right). HHKMS06 is the Honda 2006 flux[48], and Bartol[3] and Fluka[4] are existingalternative flux calculations. Taken from [49].

The zenith-angle dependence shows a peak in the near-horizontal direction. Secondary
particles travel a longer distance in the atmosphere when coming from these directions andtherefore have a greater probability to decay and produce a neutrino than those from otherzenith angles. The deflection of low energy secondary particles by the geomagnetic fieldenhances the peak at lower energies. There is an additional geometric enhancement fromthe 3D nature of the simulation (relative to previous 1D models) because the effective areafrom which a neutrino can be produced is larger near the horizon when the neutrino is
not constrained to the same dimension as the parent cosmic ray. The resulting calculated
neutrino flux as a function of energy is shown for each neutrino type in the LHS of Fig.
4.2. The ratios of µ-type to e-type and ν to ν̄ are shown on the RHS of the figure. Overallν +ν¯

there are more neutrinos than anti-neutrinos expected in the flux and the ratio νµe +ν¯µe is
approximately two at low energies and increases at higher energies.

58

Figure 4.3: The ratios of the direction-averaged flux for Kamioka in the Honda 2011(“Thismodel”), Bartol[3] and Fluka[4] models relative to the Honda 2006 flux[48], Taken from [49].

59

4.2

Neutrino interaction simulation

Neutrino interactions on nucleons and nuclei are simulated by the NEUT program[82, 83], acustom program for SK written in Fortran. Interactions of neutrinos on protons, neutrons,
oxygen and other nuclei are modeled and interactions can take place in the detector wateror in the rock surrounding the tank. Interactions within the detector are assumed to occur
on the nucleons and nuclei of water molecules, while those in the rock surrounding thedetector assume silicon dioxide (SiO2 ). NEUT covers a wide range of energies spanningfrom several tens of MeV up to hundreds of TeV. Charged and neutral current interactionsof the following modes are included in the program:

CC/NC quasi-elastic: νN → lN

00

CC/NC single-π production: νN → lN π0

CC/NC single-γ production: νN → lN γ
CC/NC single-K production: νN → lΛK
0

CC/NC single-η production: νN → lN η
0

CC/NC deep inelastic scattering: νN → lN hadrons
CC/NC coherent-π production: ν 16 O → lπX
0where N and N are differing nucleons (protons or neutrons), l is a lepton and X is
the post-interaction nucleus. NEUT provides the cross-sections of the interactions listed
above and selects an interaction mode for each event. Interactions of neutrinos on electrons
are ignored because their cross-sections are smaller than nucleon cross-sections by a few
about 103 at energies around 1 GeV. NEUT also simulates meson interactions in the target(nucleus) as well as nucleon re-scattering.

4.2.1

Elastic and Quasi-elastic Scattering

For quasi-elastic scattering NEUT considers both free and bound nuclei. Quasi-elastic
scattering on a free nucleon is simulated by the model of Llewellyn-Smith[106]. For bound60

Figure 4.4: The νµ and ν¯µ CCQE cross-sections for NEUT (line) and various experimental
measurements. See [83] for details.

nucleons in a nucleus the model of Smith and Moniz[107] is used. The latter is a relativistic
Fermi gas model which assumes a flat momentum distribution with a Fermi surface (which
follows from the Pauli exclusion principle, allowing only one fermion per quantum state)of pF = 225 MeV/c. The axial-vector form factor (the component of the weak interaction
current in which the non-conservation of parity is manifested) is assumed to be a dipole
and the axial-vector mass, MA , is set to a value of 1.21 GeV/c, which is in agreement withrecent experimental results[45, 24]. The predicted charged-current quasi-elastic (CCQE)
cross-sections from NEUT as well as experimental measurements for νµ and ν¯µ are shownin Fig. 4.4.
The calculation of neutral current quasi-elastic cross-sections was recently updated
within NEUT. The NCQE interaction between a neutrino and oxygen can produce γ-ray
decays with energies of a few to 20 MeV. These gamma rays can be a serious background forsupernova relic neutrino (SRN) detection, so estimating the interaction’s cross-section accurately can improve background estimates. Previously, NCQE cross sections were roughly
estimated as ratios of known CCQE cross-sections. After the update, the correct form factors for these interactions were used to directly calculate the cross-sections. An additional
nuclear medium correction, the spectral function, for ν −16 O interactions was included. This

61

function calculates the probability of removing a nucleon from a nucleus, given a nucleon
momentum of p and a residual nucleus state of energy E, and is used in calculating thedifferential cross-section. After this addition the ν + p cross-sections were nearly identical,
but the new ν + n cross-sections differed by ∼ 20 - 40 percent as a function of energy, so
this update was limited to only oxygen cross-sections.
4.2.2

Meson and Photon Production

Production of mesons such as the π, K and η (hadronic particles each made of a quark
and anti-quark) is simulated by the Rein and Sehgal model[101]. The model assumes aninteraction via a nuclear resonance of the form
ν + N → l + N∗
N∗ → m + N

0where N ∗ is the baryonic resonance and m is a meson such as a pion. The Rein and Sehgal
model provides the probabilities for resonance production. NEUT then multiplies that
number by the probability for each resonance to decay to one meson (π, K, η) or γ, and one
nucleon to calculate the cross-section for their respective productions. In these calculations
a total of 18 resonances below 2 GeV/c2 are considered. As in quasi-elastic scattering, theaxial-mass is set to 1.2 GeV/c2 . For pion production, the angular pion distribution was
calculated as in the Rein and Sehgal model for the ∆(1232) resonance, and was assumed
to be isotropic for the other resonances. NEUT’s predicted π + angular distribution agrees
well with the measured experimental results for ν + p → µ− + p + π from The Pauli blocking
effect is also included. This effect is based on the Pauli exclusion principle and suppresses
cross-sections by a few percent by forbidding scattering states which produce fermions of
momentum less than the Fermi momentum. Delta decays which produce a lepton andnucleon but no pions are also included and this accounts for ∼ 20% of all decays.
Coherent pion production is a special subset of meson production which is also included
in NEUT. In this process the nucleus is left in an unchanged state after an interaction which
produces a single pion (π ± in CC and neutral π in NC). The cross-section is calculated analytically and relies on pion-nucleon cross-sections taken from Rein and Sehgal’s model[101].

62

Figure 4.5: The cross-sections for CC single-π production for (a) νµ p → µ− pπ + , (b) νµ n →µ− pπ 0 and (c) νµ n → µ− nπ + . The solid line shows the NEUT calculation and the dotsshow the data points summarized in the bottom right. Taken from [83].

63

The NC process is especially important for backgrounds in accelerator beamline experiments
searching for νe appearance candidates because the π 0 decay γ’s can mimic e-like rings inthe νe signal energy region.

4.2.3

Deep Inelastic Scattering

Deep inelastic scattering (DIS) becomes the most dominant interaction neutrino crosssection at higher energies. This is the interaction which was used to probe the inside ofnucleons and helped establish the first evidence of quark structure. In DIS interactions thetarget nucleon is not kept in tact and there are many hadronic particles in the final state, so
the invariant mass of the hadronic system is larger than that of the original nucleon. NEUT
directly calculates the CC DIS cross-section by integrating the differential cross-section
above hadronic invariant masses of 1.3 GeV/c2 . The parton distribution functions, which
describe the momentum distributions of quarks inside the proton, are taken from GRV98 [78]with corrections from Bodek and Yang [6]. A pion multiplicity function is used to describe
multiple pion production. Single pion production was handled as described in the previous
section to avoid double counting. This function took the form nπ = 0.09 + 1.83lnW 2 , where
nπ is the mean number of pions produced and W is the hadronic system invariant mass.
The NC multiple pion production rates were treated simply by multiplying CC interaction
rates by coefficients, whose value was different depending on energy region[83]. Hadronic
final states are simulated by a combination of custom NEUT code[50] used up to masses of
2.0 GeV/c2 , and Pythia/JETSET[105] for higher energy interactions.

4.2.4

Nuclear Effects

The particles which are generated in the primary neutrino interaction can undergo interactions in the nucleus before escaping. Therefore, the particles eventually observed by thedetector outside the nucleus may not be the same as those initially generated by the interaction. The nuclear effects are therefore very important and NEUT simulates these effects forpions, kaons, etas and nucleons. These final state interactions (FSI) for pions were recently
updated. Pions in particular are important because the cross-sections of pion production

64

Figure 4.6: The cross-sections from old and new pion interaction models in NEUT compared
to data as a function of pion momentum.

in neutrino interactions and subsequent pion-nucleon interactions are relatively large. The
propagation of the π through the nucleus is done by stepping it through the nucleus andcalculating the probability of interaction at each point using the mean free path (MFP)
information. The interactions of pions within the nucleus that occur with high probabilities include charge exchange (π ± → π 0 ), quasi-elastic scattering (affects kinematics), andabsorption (loss of a pion, typically produces nucleons).
For low energy pions, pπ < 500 MeV/c, the interaction probability per step is calculatedby the model of Salcedo and Oset[104]. Their model includes Fermi motion in the nucleus
and Pauli blocking effects, as well as a density-dependent Fermi surface momentum. As aresult the pion interaction probabilities depend on the momentum as well as the position
within the nucleus. The LE probability table was recently tuned to achieve better agreementwith π + + C scattering data. This was done by altering the momentum dependence of theinteraction probability per step.
The scattering kinematics are determined by a phase-shift analysis of free nucleon scattering data[30, 58] which was unchanged from previous analyses. For pπ > 500 MeV/c, thenucleons are assumed to be quasifree particles. The π ± P cross-sections are used directly to
calculate interaction probabilities. As part of the recent updates, the elastic and chargeexchange cross-sections, σel and σCX were scaled by a factor of 1.8, motivated by the most

65

recent high energy π + -C data. The scattering kinematic calculations were also changed
to use the same phase-shift method as the low energy scattering. In the case of pion absorption, NEUT previously did not generate any additional particles. This was updated sothat now ejection of protons and neutrons after such pion interactions (absorption and also
scattering) is modeled. This is important for vertex activity studies in neutrino interactions
in current and future detectors. This implementation referenced both π-nucleus absorption
data [26] and π-d data[102] for consistency.
4.3

Detector Simulation

After the neutrino interactions modeled by NEUT, the outgoing particles are fed into asimulation which tracks them in the detector. This detector simulation is responsible fortracking the charged particles path in water, their cherenkov radiation, any secondary interactions they may undergo, reflection from Tyvek in the detector scaffolding and finally
the PMT and electronics responses. The detector simulation is a customized program called“skdetsim” based on CERN’s GEANT package[7] (in particular, GEANT3). The simulation
includes many physics effects within the tracking. Photons can undergo pair production andcompton scattering. The photo-electric effect at the photocathode surface is also modeled.
All charged particles (e± , µ± , Hadrons) can emit Cherenkov radiation, undergo multiple
scattering, ionzation, Bremsstrahlung, as well as in decay in flight (µ and Hadrons), among
other processes (see [88] for more details).
The photon propagation in the simulation models the Cherenkov radation with a wavelengthdependent group velocity. The calibration procedures described in the previous chapter
determine the parameters which go into the effective index of refraction, which is used to
calculate the velocity. Photons can be absorbed or scattered by water molecules and theirintensity is assumed to decrease according to equation 3.1. As mentioned in the Calibration
chapter, both Mie and Rayleigh scattering are considered and are tuned for agreement withthe calibration data. The reflection and absorption of photons on the detector materials
such as PMTs and black scaffolding sheets are calculated from probability functions of thephoton incident angle that are estimated from direct calibration measurements. The PMT
response is modeled by assuming the measured quantum efficiency for each tube. The to-

66

Figure 4.7: The QE differential cross-section for pion scattering on a thin-water target. The
new SGNUCEFF model agrees with scattering data much better than the previous SGPI
model.

tal charged detected by each PMT is taken to be the summation of all p.e.s weighted by
the single-p.e. distribution. The hit timing is also smeared by the measured PMT timing
resolution. The simulation also accounts for the reduced number of PMTs in the SK-II run
period.
The hadronic interaction model for pions in water was recently updated in SKDetsim to
improve agreement with recent data. Pion interactions are modeled by two separate codes
depending on the pion momentum. In both cases the pion hadronic interaction length isdetermined by the total cross-section given by

σT = σOxygen + 2σH .

(4.1)

The GCALOR package[114] is used for modeling high-energy interactions (pπ gt 500
MeV/c) and shows good agreement with the data. For low-energy pions (pπ ≤ 500 MeV/c)custom codes[50] are used and these have been recently updated. The previous code[50]called SGPI relied on older data which was tuned to 27 Al pion scattering. The new model,
called SGNUCEFF, uses an improved measurement of Eq. 4.1 tuned to 16 O data with

67

Figure 4.8: Number of decay-electrons for a selected SK-4 event sample. The black is data,
the gray the MC used in previous SK analyses, and the blue and red histograms are thepredictions from the new MC with old and new pion interaction cross-sections, respectively.

the interaction channels determined by the NEUT cascade model. The 1 H cross-sections
are modeled as fractions of scattering cross-sections on H2 O. The improvements to the QE
differential cross-section of pion scattering (pπ = 213 MeV/c) from a thin-target water
simulation are shown in Fig. 4.7. The SPNUCEFF model shows good agreement withdata at low pion momentum, unlike the previous SGPI model which showed a forward-peakbias. These updated pion interactions were included in the monte carlo detector simulation
used in the analyses described in this work. These changes in general decreased the pion
interaction cross-sections, leading to more pion decays producing decay-e, an example ofwhich is shown in Fig. 4.8 (the data sample shown in this plot is explained in the following
chapter). As seen in the plot the mean number of decay electrons predicted in the montecarlo (red histogram) is increased by the updated pion cross-sections relative to the old MC(gray histogram).

68

4.3.1

OD Tuning

As part of the detector simulation process there are a set of parameters specific to the OD
which its response and are specific to each phase of the detector. These parameters describe the PMT response (dark rate, collection efficiency, single-p.e. response, wavelengthshifting plate characteristics) as well as the PMT support structure materials (reflectivity/transmission rate of the Tyvek sheets) independently for the top, barrel and bottom
portions of the detector. The tuning method was updated recently during the SK-IV run
period. Brief descriptions of the original methods and the SK-IV improvements are described here. A more detailed description of the previous method can be found in [112].
Originally the tuning was done manually by looking at distributions of all relevant
variables until sufficient achievement is found. This was done for detector phases SK-I,
SK-II, SK-III, and initially in SK-IV. Cosmic-ray muon samples were used as the data fortuning. Through-going muons which enter from the top and stopping muons which enterfrom the side were chosen for the sample. MC muons with the same directions and energies
as the reconstructed muons from the data are generated and the detector response between
the two is compared and tuned until agreement is satisfactory. The light reflection model of
the Tyvek includes both Gaussian and Lambertian (isotropic diffuse reflection) components.
The fraction of each component is based on the measurements of [8].
The single-p.e. response is measured by taking hits in the muon data sample which
occur before the event itself because these dark noise mostly represent one photoelectron.
A PMT response function made from the combination of a plateau response at low chargeand exponential response at higher charges was constructed for each detector phase. PMT
responses in the simulation randomly sampled this distribution. PMT saturation is relevant
for hit PMTs close to the entrance or exit points of energetic muons and two parameters in
the tuning account for this.
Recently the tuning process has undergone an upgrade which introduces some automation. In SK-IV, the simulated charge in the OD top was determined to be too low relative
to the data samples after initial tuning. This is because the cosmic-ray muons samples
mentioned above do not provide a good sample for tuning the top of the OD. To improve

69

Figure 4.9: Distributions of SK-IV OD variables for cosmic ray data (black), MC with oldtuning (purple) and MC with the most recent tuning (red).

70

this, a monochromatic MC sample of muons which exit the top of the detector was generated and added to the set of data histograms used in the tuning. Additionally, the process
was partially automated by defining a grid for all the OD tuning parameters and generating distributions (comparing updated data to MC samples) on all points of this grid. A
Kolmogorov-Smirnov test is performed for each set of generated distributions and the tuning parameters which provide the best result from this test are chosen. The distributions
used in this tuning included nhitac (number of hits in the OD cluster with the most charge),
nhitac2 (number of hits in the second-largest OD charge cluster), Odrmean (average distance among all pairs of OD hits), total OD charge and top OD charge. The resulting
distributions of the tuned MC for the original and updated SK-IV tunings are shown in Fig.
4.9. In general the agreement is improved by the new tuning.
The old (from IMB) and new (made for SK) OD PMTs have different relative efficiencies
which were also included in the tuning. The new calibration study described in the previous
chapter was performed and the average relative efficiency for old versus new tubes wasmeasured to be 0.96 ± 0.03.

71

Chapter 5
ATMOSPHERIC NEUTRINO DATA REDUCTION

The SK data stream processes about 106 trigger events per day. The large majority arebackground events, such as cosmic ray muons, gamma rays from radioactive atoms suchas Radon, and PMT internal discharges triggering a flasher event. These are filtered out
by the reduction process, thus producing the SK data set requires a thorough and efficientsequence of reduction steps. There are three main reduction chains categorized broadly
by event geometries: Fully Contained (FC), Partially Contained (PC) and Upward-Going
Muons (UPMU), listed in order of increasing average energy. For atmospheric neutrino
analyses only events with a visible energy above 30 MeV are passed on to reduction.
The first two categories, FC and PC, comprise events with an interaction vertex locatedwithin the fiducial volume of the ID, which is defined to be greater than 2m away fromthe IDs surface. FC events deposit all of their lepton’s Cherenkov energy in the ID (inother words, the Cherenkov radiation emission track does not make it out of the ID. PC
events have higher energies on average and are characterized by the charged lepton exiting
the detector and creating an exit point cluster in the OD in the process. PC events are
assumed to be from muons because electrons generally do not travel far enough to exit thedetector.
The UPMU sample is made up of high-energy muon events which have an interaction vertex in the rock beneath and surrounding the detector. Only upward-going muon
events are chosen because neutrino-induced muons cannot be isolated from the high rate ofdownward-going cosmic ray muons penetrating the detector. These are further separated
into two categories, stopping and through-going, based on whether or not the reconstructedmuon track exits the detector. Schematics of the main event categories are shown in Fig.
5.1. Each of these three broad event categories contains more detailed event sub-categories
whose characteristics will be described in the following sections. All data samples are di-

72

Figure 5.1: Typical neutrino (dashed) and resulting charged lepton (solid) tracks for the fully
contained (FC), partially contained (PC), upward-going stopping (Up-Stop) and upward
through-going (Up-thru) event types, shown left to right.

vided by the four detector phases: SK-I (1996-2001), SK-II (2002-2005), SK-III (2006-2008),
and SK-IV (2008-present). The analyses in this work uses SK-IV data for the first time.
The SK data is divided into run files which are limited to 24 hours at most. Runs canbe less than 24 hours especially in cases where calibration or other detector maintenance istaking place. A typical 24-hour run corresponds to about 150 GB of data, the bulk of which
is just low energy background events. Each run is divided into many subruns which areeach stopped at either ten minutes of data or ten MB of data (whichever limit is reachedfirst). Each of the three reduction streams described here have become automated within
the past few years so that the most recent data is run through the reduction daily.
5.1

Fully Contained Reduction

FC events are first separated from PC events by requiring the number of hits in the largest
OD cluster to be fewer than 10. The following are the five main steps in the reduction
process:
5.1.1

First Reduction

An event passes the first reduction step if it satisfies the following criteria:
1. PEtotal > 200 p.e.s (100 p.e.s for SK-II), where P.E. is the total ID charge collected

73

in a 300 nsec time window.
2. NHITA800 ≤ 50 where NHITA is the number of OD hits in an 800 nsec window
surrounding the trigger timing.
The first criterion eliminates low-energy background events from radioactivity. A 200p.e. event corresponds to an electron momentum of 22 MeV/c. Because the analysis event
energy threshold is 30 MeV, this cut does not eliminate analysis-worthy events. The second
cut removes cosmic ray muons, which deposit energy in the OD as they enter the detector.
The output rate for the FC1 reduction is 1500 events/day on average.
5.1.2

Second Reduction

1. NHITA800 ≤ 25 for events with PEtotal < 100, 000 p.e.s (50,000 p.e.s for SK-II)
2. PEmax /PE300 < 0.5, where PEmax is the maximum number of p.e.s observed by asingle PMT within the event and PE300 is the total number of p.e.s observed in the
ID in a 300 nsec sliding time window
The first criterion performs a stricter cut on cosmic ray muon events. The second
criterion helps remove events caused by electronics noise, which tend to have a single PMT
discharge with a signal much larger than an average PMT signal and make up the majorityof charge for an event. After FC2 reduction there are about 200 events/day remaining.
5.1.3

Third Reduction

The FC3 reduction cuts reduce the event rate to about 30 events/day. The majority ofremaining backgrounds which are put into the third reduction are low-energy radioactivity
events and cosmic ray muons which do not produce much OD activity. The following are
specific categories of background which are removed by targeted reduction cuts.
Through-going muonsThrough-going muons deposit large amounts of energy in the ID and if their OD activity is
minimal they survive the first two reduction cuts. For events with at least 1000 ID PMTs

74

each detecting more than 230 p.e.s a special through-going muon fitter is applied. With theassumption that such a muon caused the ID light pattern, the fitter sets the entrance point
to be the earliest hit PMT among the entrance cluster. The exit point is defined to be thecenter of all PMTs which are saturated around the muon’s exit location. A goodness-of-fit
parameter is constructed by the function

goodness =1
Σi σ12i×(ti − Ti )2
exp
−
σi2
2(1.5 × σi2 )

X 1
i




(5.1)

where i is the PMT index, ti and σi are the i-th PMTs observed hit timing and timing
resolution, and Ti is the estimated hit timing based on the muon’s assumed entering point
and track. After this is calculated, the cut criteria are defined as
1. goodness > 0.75
2. NHITAin ≥ 10 or NHITAout ≥ 10
where in (out) represents the number of hit OD PMTs located within 8 m of the assumed
entrance (exit) point within an 800 nsec time window. Events which pass these cuts are
rejected as through-going muons.
Stopping Muon Cut
Cosmic ray muons which stop in the detector are also eliminated in a similar manner. The
fitter estimates the entrance point via the same method as through-going events, and relieson the same goodness of fit as in Eq. 5.1. All events which have NHITAin ≥ 10 OR eventswith NHITAin ≥ 5 and goodness > 0.5 are rejected in this cut. The muon direction isreconstructed so that the number of p.e.s within the cone half-opening angle of 42◦ (therelativistic limit of the Cherenkov radiation opening angle) is maximized.
Cable Hole Muons
The top of the detector includes twelve holes through which the cables supplying the highvoltage and return signals of the PMTs pass through. Because of the cable bundles there
are no OD PMTs to act as vetos installed at these positions. As a result, cosmic ray muons

75

Figure 5.2: Location of the cable holes in the top of the tank. Red circles indicate holes
with veto counters and pink circles those without. Taken from [92].

which pass through these holes can enter the ID without depositing appreciable energy in
the OD. Special veto counters and a dedicated set of cuts were used to remove these eventsfrom the reduction. These cuts have an inefficiency (incorrectly eliminating true neutrino
events) of less than 0.01%. The veto counters, which are 2 × 2.5 m plastic scintillation
counters, were installed in four of the twelve cable holes in 1997. A schematic of the twelve
holes and the four veto counters is shown in Fig. 5.2.
Events with at least one veto counter hit and a reconstructed vertex within 4 m of anyof the veto counters were rejected. A new cut was introduced in SK-IV to reject events
which pass through cable holes that do not have veto counters. The criteria for this cut are:1. goodness ≥ 0.4 (Eq. 5.1
2. PEtotal > 4000 p.e.s (total number of ID p.e.s)
3. reconstructed as downward-going (cosθ < −0.6)
4. Lveto < 2.5m, where Lveto is the distance from reconstructed vertex to cable hole

76

Figure 5.3: Typical flasher (i) and FC neutrino (ii) event timing distributions. Taken from[87].

Flasher Event Cut
Flasher events are caused by repetitive internal discharges of PMTs and have a broader event
timing than actual neutrino-induced events. They are removed as events which satisfy thefollowing criteria:
1. SK-I: NMIN100 ≥ 10 OR NMIN100 ≥ 10 if the number of ID PMT hits ≤ 800
2. SK-II,-III and -IV: NMIN100 ≥ 20
where NMIN100 is the minimum number of ID hits found in a sliding 100-nsec window that
follows the event trigger. The timing distributions of typical flasher and FC neutrino events
are shown in Fig. 5.3. The distributions of the variables used in these selection cuts for FC
data and MC are shown in Fig. 5.4. As shown in the figure the cuts remove a majority of
flasher events in the data which tend to have larger NHIT100 values.
Other Cuts
Accidental coincidence events occur when a low energy background event causes a trigger
to be applied while a cosmic ray muon coincidentally passes through in the same trigger

77

Figure 5.4: Typical NMIN100 distributions for data (a), MC (b) and post-FC3 data andnormalized MC (c) samples. Taken from [87].

window. These events can survive the first FC reduction steps because their OD activity
surrounding the initial low-energy event timing is zero or negligible while the muon event
deposits a significant amount of energy in the ID. These events can be removed by requiringthat NHITAof f ≥ 20, where NHITAof f is the number of OD hits within +400 nsec to+900 nsec after the event timing, and PEof f > 5000 p.e.s (2000 for SK-II), where PEof f
is the number of total p.e.s detected in the ID within the same 500 nsec timing window.
Remaining low energy events from radioactive backgrounds or electronics noise are removed
as events which have NHIT50 < 50 (25 for SK-II) where NHIT50 is the number of hits ina sliding 50 nsec window. The PMT timings are corrected to include photon TOF after
assuming a single vertex source for all photons which is determined by the peak of the timing
residual distribution. After all cuts in FC3 reduction are applied, there are approximately
35 events/day remaining.

78

5.1.4

Fourth Reduction

The FC fourth reduction is dedicated to remove remaining flasher events. Periodic flasher
events tend to have a physical charge pattern which is highly correlated to all other events
caused by the same flashing PMT. An example of two flasher events separated by about 5
minutes is shown in Fig. 5.5. An algorithm designed to identify flashers is applied to for
the fourth reduction. The algorithm includes:

1. Divide the ID surface into 1450 patches, each including between 6 and 9 PMTs
2. Calculate a correlation r for each patch based on its charge for pairs of events, where
r is
r=

A
B
B
1 X (QAi − hQ i) × ((Qi − hQ i)),NσA × σB(5.2)

iwhere N is the number of patches, i is the PMT index, A and B represent the two
events in consideration, hQi is the averaged charge and σ is the charge distribution’s
standard deviation.
3. Calculate the distance DISTmax between the PMTs which have maximum charge ineach event and apply an offset r = r + 0.15 if DISTmax < 0.75 m.
4. Mark A and B as matching events if the r > rth where rth is a threshold quantity thatis a function of the total number of ID p.e.s in each event.

The same procedure is applied for 10,000 events before and after the target event eachmaking a pair with the target event, and the number of matching events are tallied. Events
are then removed by a 2-dimensional cut using r and Nmatch as the selection variables. A
sample plot of this selection for SK-I is shown in Fig. 5.6.

5.1.5

Fifth Reduction

The fifth and final stage of FC reduction aims to remove all remaining background events.
The specific backgrounds which are targeted are as follows.

79

Figure 5.5: A typical pair of flasher events caused by the same PMT separated by 5 minutes.
Taken from [88].

Figure 5.6: The criterion for removing flasher events in FC4. Vertical axis shows correlation(r) and horizontal axis shows number of matching events. Taken from [88].

80

Figure 5.7: The removal of stopping muons in FC5. NHITAin is shown for (a) data, (b)MC and (c) final FC sample with normalized MC. Taken from [88].

Stopping Muons
Stopping muons which were not eliminated in FC3 are now removed with a similar cut thatreconstructs the muon entrance point in a more precise manner. Events with NHITAin ≥ 5
where NHITAin is the number of OD PMT hits within 8 m of the reconstructed entrance
point in a sliding 200 nsec window that covers ± 400 nsec around the event time. The
entrance point is estimated by constructing the distribution of PMT residual times for eachpossible candidate vertex and selecting the point which yields the sharpest peak in the residual time distribution. The entrance point is then estimated as the backward extrapolationof the vertex. Fig. 5.7 shows the effect of this cut.

Invisible Muons
Muons which are below Cherenkov threshold can escape detection yet produce a decay
electron which is reconstructed as an FC neutrino event. These events can be identified by
an OD signal (from the muon before it goes under Cherenkov threshold) detected ahead of

81

a low-energy signal in the ID from the electron. Events satisfying the following criteria areremoved:1. PEtot inID < 1000 p.e.s
2. NHITACearly > 4, where NHITACearly is the maximum number of OD hits in a cluster
in a large timing window (-8800 nsec to -100 nsec) preceding the ID event.
3. NHITACearly + NHITAC500 > 9 if DISTclust < 500 cm OR NHITACearly > 9 if
DISTclust ≥ 500 cm, where NHITAC500 is the number of OD PMT hits in the cluster ina (-100 nsec, +400 nsec) window surrounding the event, and DISTclust is the distance
between the two OD clusters described by the other variables.

Other Cuts
Accidental coincidences of a low-energy event and a muon are further removed by selecting
events with more than 20 late hits in a 200 nsec sliding window which starts searching 400
nsec after the event timing. This cut is found to be effective in reducing background without
significantly removing any data. Finally, flasher events which have a long timing distribution
tail are removed by a stricter version of the corresponding FC3 cut. Events found to havea goodness 5.1 < 0.4 and at least five hit ID PMTs in a 100 nsec sliding window covering
the interval [+300,+800] nsec relative to event timing are selected for removal.

5.1.6

Final FC Sample

The official event criterion for FC events remaining after all reduction stages are (1) vertex
inside ID fiducial volume, (2) NHITAC (number of hits in the largest OD cluster) < 9,(3) visible energy (assuming electron event) > 30 MeV. The efficiency (calculated as percentage of true FC events that survive the reduction) of FC reduction is estimated by theatmospheric FC MC and is found to be greater than 99%, with an estimated systematic uncertainty of less than 1%. The fiducial volume cut in addition to the cuts mentioned above
greatly reduce the background from radioactive sources and cosmic ray muons. These background contaminations are each estimated to be less than half a percent for all FC event

82

samples[88]. The final fiducial volume FC sample (FCFV) recorded approximately 8.1 ±0.1 events/day.
5.2

Partially Contained Reduction

The PC sample represents events with higher average energy than the FC sample. PC
event tracks penetrate the wall separating the ID and OD and cause at least 10 OD hits.
The additional OD activity is the only main difference from FC events but this causes
the reduction steps to differ drastically because many of the FC cuts are OD-based. PC
events are predominantly muons because electrons rarely travel far enough above Cherenkov
threshold to exit the OD. The main goal of the PC reduction process is to eliminate cosmic
ray muons background events, which can either be through-going or stopping. An overview
of the reduction steps is written here, while more thorough descriptions can be found in thetheses [88, 19].

5.2.1

First Reduction

The first step in the PC reduction, which is run online, removes easily identifiable low-energy
events and cosmic ray muons. Events satisfying the criteria

1. Total number of ID p.e.s ≥ 1000
2. Width of hit OD PMTs’ timing distribution ≤ 260 ns
3. Number of OD hit clusters ≤ 1

pass through the first reduction. The first cut corresponds to a minimum momentum of 310
MeV/c for muons, safely below the event threshold of 500 MeV/c for exiting PC events which
is used. The second and third cuts are effective at removing cosmic ray muon backgrounds.
Through-going muon backgrounds pass through the OD twice, so they typically produce
multiple OD hit clusters with a broad overall timing distribution. Here a cluster is definedas all PMTs with hits within 8 m of a PMT that detects at least 8 p.e.s. Approximately
14,000 events/day remain after PC1 reduction.

83

5.2.2

Second Reduction

PC2 reduction uses a different cluster definition algorithm to remove backgrounds. The
ID and OD walls are divided into a grid of isolated patches, 21x21 and 11x11 in number,respectively. The charge in each grid (sum of charges of all PMTs in the grid box) is
calculated, and the charge vector gradient is taken for each grid to separate all grid boxes
into clusters. After constructing the clusters events meeting the following conditions are
saved:
1. Number of OD hit clusters with more than 6 hit PMTs ≤ 1
2. Endcaps vs wall cut: events which a small number of hits in the endcaps (top andbottom of the detector) relative to their wall (barrel) hits are saved. The 2-d function
which determines the threshold is shown in Fig. 5.8. Events above the red line arerejected as corner clipper events.
3. Number of p.e.s within 200 cm of the PMT with the highest charge in the ID hit
cluster closest to the OD cluster is above a threshold which is a function of the OD
cluster hits (see [19] for details).
The first cut is designed to remove more through-going muons. The second cut is a topologybased cut which identifies corner clipper events that deposit most of their energy on the edgeof the detector between the top and barrel or the bottom and barrel. The third cut here
identifies and removes stopping muons which deposit charge in the OD but not significantlyin the ID. These cuts reduce the PC event rate to approximately 2000 events/day.

5.2.3

Third Reduction

A test for flasher events is applied in the third reduction by cutting events with broad timing
distributions, as in FC reduction. This is done by removing events for which N M IN100(the minimum number of hit ID PMTs in a sliding 100nsec window from +300 to +800
nsec after the event) > 14 or N M IN100 > 9 if the number of hit ID PMTs is < 800. A
fitter called pfit is used to reconstruct the particle’s interaction vertex and direction. The
entrance point is estimated by extrapolating the vertex backward with this direction, and84

Figure 5.8: Number of PMT hits in top and bottom versus number in barrel for PC data.
Events above the red line are rejected. Taken from [19].

events for which the number of hits in the OD within 8 m of this entrance point is fewer thanor equal to ten survive the reduction. Events are required to meet both of these criteria to
be saved. Plots showing the removal of events based on the last cut for data, MC and thefinal samples for both are shown in 5.9

5.2.4

Fourth Reduction

Fourth reduction relies on the muon fitters pfit and mfmufit to remove mostly cosmic raymuons with little OD activity. The latter fitter determines the entrance point of the particle
to be the location of the earliest ID cluster. Events with a total ID charge (Qismsk) greaterthan 1000 p.e.s are saved. Events with longer track lengths are typically cosmic muonsand those with a length greater than 3000 cm are removed. A goodness of fit for the muon
fitters is calculated and events for which this value is above 0.85 (a good fit indicating a highprobability of a through-going muon) are rejected. Events in which the distance between
the vertex as calculated by the point fitter and the corner of the tank is greater less than150 cm are also removed. These events would eventually be removed by the fiducial volume
cut but their removal during PC4 saves computation time.

85

Figure 5.9: Distributions of NHITAin for data (a) and MC (b) for the PC3 cut, and thefinal data and MC samples (c). Taken from [88].

Finally only events satisfying d~pf it · d~P M T > −0.85 are saved. These vectors are the track
direction as found by the point fitter (d~pf it ) and the direction from the vertex to the earliest
saturated ID PMT (d~P M T ) and this cut removes cosmic muons which are likely to have
opposite directions for these two vectors. A sample plot for SK-1 data and MC showing
the effect of the vector cut is shown in Fig. 5.10. The output of PC4 is approximately 20
events/day.
5.2.5

Fifth Reduction

The fifth and final stage of the PC reduction includes many specific cuts targeting theremaining backgrounds. In this step more precise fitters (called apfit, stmfit, pfdofit
and msfit) are used to fit the event. Events with a total (ID) charge greater than 3000
p.e.s are saved, removing low energy events. A detection of 3000 p.e.s indicates a muon of
approximate momentum 500 MeV/c and because muons which exit the ID and reach the OD
typically have at least 700 MeV/c, this cut is safe. The number of OD hits, nhitac, is used
to distinguish stopping (< 16 PMTs) and through-going (≥ 16 PMTs) PC events. Next a86

Figure 5.10: The effect of the PC4 cut which is based on the pfit direction and the direction
from the vertex to the earliest saturated PMT. Shown for SK-1 data(a), MC(b) and their
final samples(c). Taken from [88].

87

combination of hits based on OD cluster information are applied. Multiple OD clusters are
assumed and they are ranked according to their charge. Events which have fewer than 2 hits
in their third-highest charge cluster, a distance between the two highest-charged clusters ofless than 20 m (through-going muons are likely to fail this cut) or which have fewer than 10
p.e.s in the second-highest charge cluster are all saved (only one of these three conditions
needs to be satisfied). Additionally there are cuts made on events which fall into calibration
time windows and runs which exhibited problems with the electronics huts.
Junk cuts include those events for which activity spread across two triggers or in which
the OD had a problem are also applied. Cable hole muon cuts are applied similarly to the
FC reduction. Events with one veto counter hit and d~ring · d~veto−vertex > −0.8, were d~ringis the reconstructed ring direction and d~veto−vertex is the direction from the reconstructedvertex to the hit veto counter, are considered cable hole muons and removed.
A set of cuts for removing remaining stopping and through-going muons is also added
here. To eliminate stopping muons, events which have more than ten entry hits in theentrance cluster found by msfit are removed. The charge in a cone of radius 8 m aroundthe ID entrance location from the reverse extrapolation of the pfdofit direction is required
to make up 60% of the ID total and contain fewer than six hits. The number of hits within
8 m of the OD entrance point (Umehit8m, same as above) as calculated by apfit is required
to be fewer than ten. The angle between the fitted direction and the direction from thevertex to the highest charge OD cluster (all calculated with the precise fitters) is required
to be no more than 90◦ .
Through-going muons which travel vertically downward in the region between the ID and
OD leave a track which resembles a millipede are removed by the so-called geji1 cut[19].
A separate clustering algorithm called grad cluster is used to identify OD clusters andthe two with the most hits are each required to have fewer than 10 hits. If the entrance
and exit hit clusters as found by msfit have at least five hits each, then events for which
the time difference between the average top and bottom hit timing divided by the distanceseparating them is between 0.75 and 1.50 are removed as through-going candidates. The
1“Geji” is the word for millipede in Japanese.

88

visible energy as calculated by apfit (a different version of the variable Evis) along with theestimated track length between vertex and OD exit position are used to remove remaining
corner clipper events. Finally low energy events with Evis < 350 MeV within the fiducial
volume are removed.

5.2.6

PC Reduction Summary

The PC reduction was found to have efficiencies of 88% (SK-I), 83% (SK-II), 96% (SK-III)and 96% (SK-IV) for the data sets used in the analyses described in this work. These
numbers are taken from collaboration meeting slides. The final background contamination
comes mainly from cosmic ray muons and is accounted for by manual eye scanning of the
final event sample. The majority of this background is outside the FV and the contamination within the FV is found to be less than 1.0% for each detector period. A systematic
uncertainty is estimated for the PC reduction by comparing data and MC (for ID cuts)and generating toy MC sets with different OD tuning parameters and seeing how they differed (for OD cuts). The uncertainties from each stage of reduction were then added inquadrature to calculate the overall uncertainty. The efficiencies for each reduction stage
and run period along with the estimated background contaminations are shown in Table
5.2.6. The selection efficiencies shown in the table are for events with real vertices in thefiducial volume, number of OD hits > 15 (9 for SK-II) and total observed charge > 3000
p.e.s. The final fiducial volume PC sample recorded 0.66 ± 0.02 events/day on average.

5.3

Upward-going Muon Reduction

The upmu sample represents muon-like events with the highest average energy among all
atmospheric neutrino data samples. Upmus are created when a muon neutrino scattersfrom a nuclei via a charged-current weak interaction in the rock either beneath or next tothe detector cavity. An example Feynman diagram for this interaction is shown in Fig.
5.11. There is also a small contribution to the upmu sample from neutrinos scattering fromoxygen nuclei in the OD water and then entering the ID with upward momentum. The

89

PC
reduction step

SK-I

SK-II

SK-III

SK-IV

selection efficiency %

1st

98.98

98.58

99.09

99.63

2nd

96.74

93.43

98.52

98.73

3rd

95.69

92.32

98.51

98.68

4th

89.86

84.60

97.87

97.42

5th

88.66

74.80

96.61

96.15

5th (FV)

80.98

74.80

88.80

86.30

Cosmic Ray BG Contamination

0.2%

0.7%

1.8%

0.6%

Table 5.1: A summary of the PC reduction efficiencies for each detector phase (top) andthe estimates of the corresponding cosmic ray background contamination percentages.

Figure 5.11: Charged-current weak interaction which produces can produce an upmu. Taken
from [111].

90

produced muon enters SK with net upward momentum and deposits energy in both the ID
and OD. The fundamental difference between the FC/PC samples and the Upmu sample
is that the vertex of the neutrino interaction is outside the ID. This topology implies thatthe measured energy represents at best a lower bound on the muon energy. The upmu
sample is used in zenith angle oscillation analyses such as the ones described here, and isalso effective for astrophysical neutrino studies because astrophysical neutrino signals are
expected to have a harder energy spectrum than atmospheric neutrinos (see [111],[109] for
SK analyses on such topics). The final upmu data samples contain between 1.4 and 1.5
events/day.
Upmu events are divided into the three sub-categories of stopping (muons which stop in
the detector), through-going non-showering (more energetic events which pass through theentire detector) and through-going showering (the highest energy events). A χ2 measureis used to separate the through-going sample into showering and non-showering events.
Showering events are characterized by the drastic (less uniform than typical ionizationenergy loss) mechanisms of energy loss which occur such as Bremsstrahlung radiation, pair
production and photo-nuclear processes. Typical parent νµ energies are tens of GeV for
stopping, few hundreds of GeV for thru, and many hundreds of GeV to one thousand GeV
for showering. The muons created from their interaction are roughly a few GeV, tens ofGeV, and a few hundred GeVs for the same categories, respectively [111]. True energy
distributions for SK-III MC for these particles and event categories are shown in Fig. 5.12.
The upmu reduction is designed to reject background events while saving true upmu
events and nearly-horizontal events. Much like FC and PC reductions, multiple fitters areapplied during the reduction and nearly-horizontal events which were initially classified as
downward-going may become upward-going after the application of more precise fitters.
Additionally nearly-horizontal events are used to estimate the background of cosmic ray
muon events which up-scatter into the detector and are mistakenly categorized as an upmu.
As detailed below eye scanning will constitute the final check of the upmu sample andthis helps to remove backgrounds such as corner clippers, multi-cluster events (multiple
muons typically) and very high energy downward-going showering events which saturate
the detector and are mis-reconstructed as upward-going.

91

Figure 5.12: Neutrino (top) and muon (bottom) energies for the three upmu samples as
generated in SK-III MC (from [111].

5.3.1

Upmu First Reduction

The first step in this reduction chain is the umred1st program which runs continuously
online. The subruns (ten minutes of data or ten MB of data, whichever occurs first) are
merged together into larger subruns) which are easier for processing. A flow chart of theumred1st selection criteria is shown in Fig. 5.13. Events from the HE trigger stream arepassed to the first step saves ultra-high energy (UHE) events, characterized by having acharge greater than 1.75 Mp.e.s, for astrophysical studies. UHE upmu events are very rare
and a past analysis with SK-I data found only one event, against an expected backgroundof 0.46 ± 0.23 events[46].
Next two routines called stopmu1st and muboy are applied. The former uses mostly
OD information to apply one of eight flags to each event, as shown in Table 5.2 without
estimating any directional information. The OD cluster-finding algorithm is based on theidea that for typical muon events passing through a PMT result in nearly the same timing
and charge information for all nearby surrounding PMTs. Clusters are found by centering
the search on PMTs which have a local maximum charge value and requiring at least seven
nearby neighbors to have at least 8 pe within a specified 800 nsec window. One check forthis algorithm is the number of clusters found; events which are found to have exactly two

92

Figure 5.13: A flow chart describing the umred1st selection and rejection steps. Taken from[111].

clusters are classified as thru-going while those with two or more clusters separated by atleast 8 m are tagged as multi-muon events. More details can be found in Reference [18].
The latter fitter muboy is very important in the upmu reduction process because it mostbackground events are rejected by this fitter. It began as a fitter for muons which produce
spallation products that are backgrounds in low-energy analyses, and was adapted for use
in high-energy studies[14]. It is similar to stopmu1st but in addition estimates a unit vector
for the direction of the muon (muboy dir[i], where i can correspond to x(1), y(2) or z(3)directions), a track length and an associated goodness of fit quantity (muboy goodness). It
is important to note here that the index i represents the cosine of the unit direction vector
with respect to each axis, so that upward-going muons will have a dir[3] component which
is negative because this means the z-component of the direction vector is greater than 90◦
down from the upward-pointing local z-axis.
The track length is calculated by constructing a dE/dx histogram of the muon’s energy
loss. The point along the track where the number of photons emitted in a 1-meter intervalgoes below less than 40% of the average from the first 1.5 m of the track is defined as
the stopping point. A typical dE/dx histogram for this calculation is shown in Fig. 5.14.
The table of possible muboy flags is shown in 5.3. Muboy mainly uses ID information for

93

Figure 5.14: The dE/dx histogram for stopping muons used to estimate track length bymuboy[18].

94

Flag

Meaning0

thru-mu candidate

1

ID off / calibration2

decay electron3

Q < 8,000 p.e.s (3,000 p.e.s for SK-II)

4

stopping candidate

5

multi-muon

6

no result7

pedestal event

Table 5.2: A summary of all possible flags for each event that passes through stopmu1st.

assigning these flags.
Events which have muboy goodness > 0.3, muboy class (flag number) 6= 0, and 0≥ muboy dir[3] ≥ 1 are saved as upmu candidates. Near-horizontal events for whichmuboy dir[3] indicates a nearly horizontal direction are also saved with a near-horizontal
flag. Events which are identified as being downward-going, having multiple muon rings,
very short tracklengths or large forward/backward charge ratios were rejected here. The
upmu first reduction output contains about 9,000 events/day on average.
The output of first reduction is passed to a merger program which merges the many
subruns into a smaller number. After merging each run typically contains 3 or 4 merged
subruns.

5.3.2

Upmu Second Reduction

The second upmu reduction, umred2nd, has more time-consuming fitting algorithms designed to more carefully separate upmus from downward-going muons and other backgrounds. The first fitter is called stopmu2nd and its responsibility is to fit the directionof events which were classified as stopping in the first reduction. The algorithm first makes
guesses for the entrance and exit points and then calculates goodness parameters which are

95

Flag

Meaning0

invalid fit

1

thru-mu, good fit

2stop-mu, good fit

3

multi-muon, type 1

4

multi-muon, type 2

5corner clipper

Table 5.3: A summary of all possible flags for each event that passes through muboy.

based on how many PMT hits fill the center of the Cherenkov cone ring and the sharpness/fuzziness of the ring edges (more details in [18]). The entrance point is then fixed andthe muon track direction varied while another goodness of fit parameter (similar to 5.1)
based on residual PMT timing is maximized.

The event’s directional information puts it into one of six distinct categories: failed
fit, downward-going (cos(θ) < -0.12), near-horizontal (-0.12 ≤ cos(θ) < 0), regular upmu,
downward-going (top entry) and nearly-horizontal (top entry). Two fitters designed forfitting through-going events, called thrumu1st and thrumu2nd, are called next. The former
assigns essentially the same flags as stopmu2nd except for the definition of the thresholdfor downward-going, which is set to -0.1 instead. It is designed to be a quick fitter andas a result is more likely to return misfits than the other fitters. The latter specializes inmessy fits such as showering muons and multi-muon events and has the same set of flags for
assignment. The flags from the two fitters are combined and events which are categorized as
downward-going are rejected, whereas horizontal and upmu events are saved along with theirstopping/thru-going decision. The output of upmu second reduction is typically around 250
events/day.

96

5.3.3

Precise Fitter

The next step applies the most precise fitting in the reduction chain with a program called
precisefit. The fit is not applied to events with fewer than 8000 p.e.s (amount of light
too small to find good fit) and events with more than 175000 p.e.s (saturation of the ID
prevented the fitter from working well) (for SK collaborators this information can be foundin the header file umparm.h). It uses three fitters upmufit, OD-fit, TDCfit. upmufit was
taken from the FC/PC msfit algorithm and adapted for upmu fitting. It works by taking agiven entry point and direction and compares the expected number of p.e. to the observed
number in each PMT, then calculates a χ2 -statistic to quantify this. The calculation of the
expected number of p.e.s from Cherenkov light is given by the expression

Niexp =

!


1
li
2
δ
αu
sin
θ
+
N
(θ
)
exp
−
f (Θ) + SiiiidθLattenli (sinθi + li dx
)

(5.3)

where i is the PMT index, αu is a normalization factor, θi is the angle of incidence
relative to normal, Niδ is the number of p.e.s expected from interactions caused by the ionizedelectrons (sometimes called δ-rays, hence the index), La tten is the water attenuation lengthfrom 3.3 and f (Θ) is the angular dependence of the PMT acceptance. The term following
the normalization factor represents the change in the Cherenkov opening angle due to themuon losing energy. A log-likelihood is calculated from this number and the number ofobserved p.e.s and this is used to create a χ2 -statistic. This statistic is calculated for all
possible entry and exit points and the points which minimize this value are taken for thefit.
The second fitter OD-fit uses only OD clustering information to fit a vertex and direction. This is useful for events which undergo energy loss via Bremsstrahlung radiation,
which creates a fuzzier ring edge than typical ionization events. Clusters are defined to haveat least 8 nearby neighbor PMTs surrounding a PMT with hits between 800 ns and 1300
ns timing relative to the trigger. The cluster with the earliest average timing is chosen as
the entrance cluster and the cluster with the maximum number of hit PMTs is assumed tobe the exit cluster.
The third and final fitter in this step is TDCfit, which is also adapted from a contained

97

Figure 5.15: The distributions of precise fit direction and true direction for a 40-year SK-I
MC sample. Taken from [18].

event fitter. The first step to this fitter is to find a rough vertex and direction by minimizing
the differences of all timing residuals. Next, a charge statistic is constructed to test thequality of the Cherenkov ring edge location for varying track directions (details in ??). The
track which yields the maximum effective charge in the ring for this statistic is temporarily
stored. Finally a finer vertex fit is achieved by taking into account light scattering effects
and testing the goodness of fit based on the residual timing information again. These steps
are iterated until successive vertices are found which differ by no more than 50cm. For
all events the three above-mentioned fitters are applied and the eventual final fit is chosenfrom one of them. The number of events which met the charge requirements outlined at the
beginning of this section and were subsequently passed through the precisefit algorithms
was approximately 100 events/day.
The precise fitter is found to have good agreement between the true direction and its
fitted direction. This was tested with a 40-year MC sample when it was first implemented in
SK-I. A plot of the angle between the precisefit and true directions is shown for thru-going
muons and stopping muons in Fig. 5.15. As shown in the figure the resolutions are foundto be approximately 1◦ .

98

Figure 5.16: Muon energy losses in water. Curve (a) shows losses due to ionization, (b)shows radiative losses, and (c) shows the sum of both[91].

Showering Algorithm
A showering algorithm used to identify very high-energy muons which lose energy throughprocesses such as pair production, Bremsstrahlung and other electromagnetic shower cascades as opposed to ionization (the most dominant form of energy loss for muons up to many
hundreds of GeV) is applied as a subroutine within the precise fitter. The energy losses fromthese showering mechanisms become equal to those from ionization for muons with energies
around 1 TeV. A diagram of muon energy losses in water is shown in Fig. 5.16[91]. Events
which are tagged as through-going pass through this algorithm. A χ2 -statistic is defined
which adds contributions from fluctuations in charge deposition above the average value.
Ionizing muons lose energy at the roughly steady value of 2.2 MeV/cm, whereas showeringmuons undergo more sudden, drastic energy loss.

99.
The track of the muon is broken up into 50 cm intervals and the corrected charge
deposited on the PMTs from each segment (Qicorr ) is compared to the average from all
segments (hQcorr i). Additionally this average corrected charge is compared to the averagevalue for non-showering muons, which is given as a function of the track length. The χ2 statistic takes the form

2

χ =

N
−2
X
i=3

(

[Qicorr − hQicorr i]2 [hQicorr i − Q(l)]2
+
2
2
σQ
σQ(l)
i

)
(5.4)

corr

where the σ variables are the statistical errors of their corresponding charges. The
bounds on the sum ignore the first and last meter of the track. Since showering events
are by definition through-going, this corresponds to ignoring the contributions near thewalls of the detector, where the corrected charge is harder to calculate. This is calculated
for all track segments which have Qicorr − Q(l) ≥ −2.0. This latter condition prevents
the calculation from including events which have large energy loss fluctuations below theaverage which are not considered contributions to a showering-type event. This occurs forapproximately 0.5% of events in the MC - usually through-going muons with poor fits or
stopping muons which were incorrectly classified as through-going. A second variable, ∆,is also calculated to identify these types of events. It is defined by∆ = hQcorr i − Q(l) − δQ

(5.5)

where δQ is a tuning parameter that changes for each SK run period and is set so thatthe peak of the ∆ histogram for MC is at zero. A table of the tuning values for each SK
run period is shown in 5.4. After this tuning the distribution of ∆ is not perfectly peakedat zero and a systematic uncertainty which will be described in detail later in this section is
calculated to account for this. The variables ∆ and χ2 are calculated for each through-going
event (and saved in the variables sh delta and sh chi1p, for those using the SK data).

100

Phase

δQ

SK-I

0.10

SK-II

2.50

SK-III

-0.30

SK-IV

-0.30

Table 5.4: A table of the δ Q tuning parameters for the showering upmu selection.

The showering events are then officially selected by the following cuts:∆ > 0.5,

χ2 > 50.0∆ > 4.5 − 0.08χ2 ,χ2 ≤ 50.0

Scatter plots of these two variables for SK-I,-II,-III and -IV MC samples (taken from 200years samples) are shown in Fig.5.17. The SK-II scatter plot exhibits a different topology
due to the significantly reduced photocathode coverage. A study done with SK-III MC
found that this cut selects about 70% of true showering events properly (for details pleasesee [109]), which was deemed to be sufficient for studies of astrophysical neutrinos. In thecase of atmospheric neutrino oscillations, the showering sample does not undergo significant
oscillations due to its very high energies.

5.3.4

Upmu Third Reduction

The third reduction applies the selection cuts to the upmu data after the events have been
given precise fits. The main goal here is to separate the stopping and through-going events.
The criteria for through-going events are the following:

1. Number of entrance hits within 8m of entrance point (Um ehit8m) > 10 (16 for SK-II)
2. Number of exit hits within 8m of entrance point (Um ohit8m) ≥ 10 (16 for SK-II)
3. precisefit fitted track length ≥ 700 cm
4. precisefit fitted z-direction (Fit dir[3]) ≤ 0

101

Figure 5.17: A 2-dimensional scatter plot of ∆ (vertical axis) and χ2 (horizontal axis) for
showering event selection for SK-IV data. Magenta events above the red line are saved asshowering, and blue ones below as non-showering.

and the criteria for stopping events are:
1. Number of entrance hits within 8m of entrance point (Um ehit8m) > 10 (16 for SK-II)
2. Number of exit hits within 8m of entrance point (Um ohit8m) < 10 (16 for SK-II)
3. precisefit fitted momentum ≥ 1.6 GeV/c
4. precisefit fitted z-direction (Fit dir[3]) ≤ 0
The precisefit fitted momentum is estimated from the total number of p.e.s in the
Cherenkov ring detected in the ID. Note that in these cuts the direction is taken such thatcos(z) = -1 corresponds to upward-going, and cos(z) = 0 is downward-going. For throughgoing events, the track length is easily estimated by the combination of entry and exit
points. Thus, the track length cut is used to select those events. However, for stoppingmuons, there is no exit point and therefore the track length cannot be estimated this way.
Therefore a cut is made on the momentum instead because it is known that a muon with areconstructed momentum of 1.6 GeV/c has a range of approximately 7 m in water (see [109]102

Figure 5.18: The distributions of fitted tracklength (left) for through-going events and fitted
momentum (right) for stopping events from precisefit for SK-IV data and MC. Final data
cuts are applied to the data, but not the MC.

for details). The distributions of MC and data for SK-IV are shown for the fitted tracklengthand momentum in Fig. 5.18. After all third reduction selection cuts are applied, there areapproximately 2.5 events/day remaining including all three event types.

5.3.5

Upmu Fourth Reduction

This is the final step in automated portion of upmu reduction and was introduced for the
first time in SK-III. Most of the cosmic ray background is removed by the previous reduction
steps. The majority of the remaining background are corner clipper events which deposit
charge near one of the detector edges and have a very short track length. These events are
not fit well by the precise fitter in many cases. Before the beginning of SK-III, the OD was
optically segmented so that the top/bottom and barrel regions became distinct. After thisseparation identification of corner clipper was made easier[109]. Additionally in SK-III an
extra step in the reduction process was added to reduce these events even further. This
uses a precise OD-cluster fitting algorithm to remove events which had at least two OD
clusters and whose fits from muboy and precisefit differed considerably. The total upmu

103

data rate after this stage was approximately 1.75 events/day. From MC studies it is foundthat the efficiency of all automated upmu reduction steps was 98% or greater (dependingon the zenith angle; upward-going events had nearly 100% efficiency). The purity of theevent samples before the fourth reduction is only ∼ 55%, but the fourth reduction removes
a majority of background events and leaves the upmu sample with about ∼ 85% purity inthe end.

5.3.6

Automated Reduction

The author of this work has recently completed a software package which automates the
upmu reduction steps mentioned above. Until this was implemented, each reduction stage
after umred1st was run manually. In the new automated framework, a daily cron job
manages a set of shell scripts which determine which new umred1st files are made each day,and then run them through all reduction steps. The reduction jobs are submitted to theparallel-computing cluster (nqs) and the number of events saved in each reduction step istallied and shown in plots on a web page. The web page also includes a shift check so that
AtmPD group members can monitor the quality of the new data as it is processed. A plot
showing the event rates for all stages of upmu reduction which was made by this software
package is shown in Fig. 5.19. The event rate for all stages is steady throughout the SK-IV
period. The scripts and organization for this software are based upon those written for the
FC automated reduction.

5.3.7

Eye Scanning

The final stage in upmu reduction is a manual scanning of all events which make it throughall prior reduction steps. An expert trained to distinguish good upmu events from miscategorized background scans the event displays for each event individually and decides
whether to accept or reject each event. This is mainly used to remove corner clipper eventswith poor fits and very high-energy down-going muons which receive an upward-going fit
because they have nearly saturated the detector. Since the addition of upmu fourth reduction in SK-III the purity of the eye scanning sample is very high to begin with. Of the104

Figure 5.19: The event rate output of umred1st (black), umred2nd (red), precisefit(green), umred3rd (blue) and umred4th (magenta) for SK-IV from the automated upmu
reduction software.

approximately 1.75 events/day that survive all automated reduction cuts approximately 1.5
events/day survive the eye-scanning stage as good upmu events.
The eye-scanners look at event displays of each event and use the quality of the fits (frommuboy and precisefit), the timing and charge information and the OD cluster information
to save good upmus and reject backgrounds. The charge deposited on the PMTs can becolor-coded by quantity (charge mode) or by time of signal (timing mode). The latter isuseful in reconstructing the direction of events for which the fit is ambiguous. The eyescanner can also choose to apply a manual fit to the event, and if the direction from thismanual fit differs from the computer’s fit by at least 5◦ then the event is rejected, and byless than 5◦ the event is saved (even in the case of slightly downward-going events; this will
be rejected by the final official event cut which requires an upward-going muon).

5.3.8

Background Subtraction

The upmu reduction process combined with eye-scanning produces a data set which ismostly free from backgrounds. The signals of high-energy upward-going muons in SK are

105

Figure 5.20: The background subtraction method for near-horizontal upward-going muons,shown for SK-IV through-going muons. The LHS shows the azimuthal distribution of eventssummed over all zenith angles, and the RHS shows the zenith angle distribution of eventswith the BG fit function overlaid.

mostly well-fit by the fitters mentioned previously and easily verified by their OD clustering,
charge and timing information. The only significant background sources for upmu events arecosmic ray muons which up-scatter just before entering the detector. These events produce
a muon-like track which is upward going and can easily be accidentally included in the final
data set.
The SK detector is located under an overburden of 2,7000 m.w.e. which greatly reduces
the cosmic ray background. However, the rock surrounding the detector is not uniformly
thick in all azimuthal directions. There are thin regions of rock which allow more cosmic
ray muons to enter than others. A map showing through-going muon data as a function of
azimuthal angle integrated over all zenith angles is shown in Fig. 5.20. The background isestimated by fitting the three parameters p0, p1 and p2 as a function of the cosine of thezenith angle with the dependence
f [cos(θ)] = p0 + ep1+p2cos(θ) .

(5.6)

The parameter p0 is first quantified as the average event rate per day per steradianfrom region 1, which does not contain either of the thin regions. The number of events

106

in the thin regions of the rock, which can clearly be seen as the places where the numberof through-going muons spikes within region 2, is then fit with this function and fixed p0
to determine p1 and p2. The best-fit parameters after the fit to region 2 are then used toestimate the number of BG events entering the two most-horizontal bins by integrating over
the relevant zenith angles:
Z 0
NBG =

ep1+p2cos(θ) dcos(θ)

(5.7)

cosθx

The RHS of Fig. 5.20 shows the zenith-angle distribution of through-going muons in SK-IV
with the result of the fit function overlaid. The number of BG events estimated in each of
the two most horizontal bins is taken from this fit and subtracted from the actual datasets
before their use in analyses. The uncertainty of the BG subtraction is estimated by the 1σ
deviations in the fitted curve itself, as shown by the dashed lines surrounding the fit. The
upmu reduction recorded 1.5 ± 0.1 events/day on average in the final sample.

107

Chapter 6
EVENT RECONSTRUCTION

As part of the reduction processes described in the previous chapter, events undergo
reconstruction which converts their PMT timing and charge information into kinematics
variables suitable for physics analyses. The reconstruction processes are identical for data
and MC samples. The information which the reconstruction calculates includes vertex fitting, ring finding an counting, particle identification, direction and momentum reconstruction among others. The analyses described in this work eventually bin all events according
to their direction (cosine of the zenith angle), their e-like or µ-like nature, and their lepton
momentum. The reconstruction methods described below mainly apply to the contained
events in the data sample. The PC and Upmu events are all treated as single-ring µ-likeevents so some, but not all, of the following methods apply to their reconstruction. The
following are the main steps in the reconstruction process (for a more detailed description,please refer to [108].
6.1

Vertex Fitting

The reconstruction begins by vertex fitting, which is done in three steps.
6.1.1

Point Source Fit

The first step in vertex fitting is simple and assumes a single point source for all ID photons.
The time of flight is then subtracted for all photons detected at the PMTs according to theirrelative position to the assumed point source. After this subtraction a distribution of PMT
hit residual times is constructed. The vertex is then estimated to be the point at which amaximal goodness-of-fit is achieved assuming a Gaussian distribution of the form

GP =



(ti − t0 )2
1 X
exp −
N
2(1.5 × σ)2
i

(6.1)

108

where N is the total number of PMTs registering hits, i is the i-th hit PMT, ti is theresidual time, t0 is a free parameter chosen to maximize GP and σ is the PMT timing
resolution which is taken to be a constant 2.5 nsec. The point which returns the maximum
value of GP is taken to be the vertex. The particle direction d~o is then estimated by taking
a charge-weighted vector sum over all hit PMTs, with the form:

d~0 =

X
i

qi ×

~0
P~i − O~ 0|
|P~i − O(6.2)~ 0 is the vertex position found by the above method, P~i is the vector from that
Here O
point to the i-th PMT, and qi is the charge on each PMT.

6.1.2

Ring Edge Search

Assuming the point-fit vertex position the edge of the brightest ring in the Cerenkov cone
is measured next. The angle between the axis of the particle’s direction and the edge of the
Cerenkov cone is taken to be the opening angle. The edge of the ring is determined by an
estimator, which is a function of opening angle, given by
R θedgeQ(θedge ) =0

P E(θ)dθ×
sinθedge

dP E(θ)dθ2(θedge − θexp )2
× exp −
2σθ2

(6.3)

where P E(θ) is the angular distribution of the observed charge as a function of the
opening angle θ, θedge is the opening angle at the current tested ring edge, θexp is the
expected opening angle from the charge in the cone and σθ is the angle’s resolution. The
observed charge in this calculation is corrected for PMT angular acceptance, attenuation
length and scattering. The ring edge position θedge is required to be greater than the angle atwhich the most charge is deposited (θpeak and the second derivative of the corrected charge
distribution with respect to angle is required to vanish. The corrected charge distribution
P E(θ) and its second derivative are shown for a typical event in Fig. 6.1.
The estimator in 6.3 is calculated while varying the both the particle direction d~o foundin 6.2 and the opening angle. The pair which maximize the estimator are used to determine
the ring edge.

109

Figure 6.1: Typical angular distribution of the observed, corrected charge in number of p.e.
(top) and its second derivative with respect to angle (bottom), which are used in the ringedge search algorithm. Taken from [108].

110

6.1.3

Precise Vertex Fit

After locating the edge of the Cerenkov ring a more precise vertex fitter called TDC-fit isapplied. Casting away the previous point-source vertex assumption, this fitter now takes
into account the track length of the charged lepton and also scattering effects for the photons
emitted. The fact that photons are now generated along the length of the track is reflectedin the definition of the PMT timing residuals, defined as:1~ i − O|~ − n × |P~i − X~ i |, for PMTs inside Cerenkov ring
× |X
c
c
n
~ for PMTs outside Cerenkov ringti = t0i − × |P~i − O|,cti = t0i −~ is the vertex position, X~ i is emission point along the track
Here i is the PMT index, O
for photons reaching PMT i, n is the refractive index of water, and P~i and ti describe theposition and hit timing of the i-th PMT.
The goodness-of-fit is now calculated in separate ways for PMTs inside and outsideof the Cerenkov ring. For those inside the ring the goodness-of-fit has the same form as
6.1. For PMTs outside the main ring a modified goodness is calculated which includes
scattering effects for those hits which occur later than t0 (for details please see [108]. The
final goodness-of-fit is estimated by combining the in-ring and out-ring parameters and thevertex which maximizes this is chosen. The vertex position from this calculation is used forevents which will be found to have multiple rings. A distribution of the difference between
true and reconstructed vertex position for one such sample, the SK-IV FC multi-GeV multiring µ-like sample, is shown in Fig. 6.2. The resolution for this sample is found to be 87.8
cm. Resolutions for older samples are detailed in [108] and in all cases are found to be near
100 cm or less.
6.2

Ring Counting

After determining the vertex and the main Cherenkov ring, a search for other ring candidates
based on the Hough transformation method[15] is performed. A schematic showing thebasic principle of this method is shown in Fig. 6.3. This image depicts hit PMTs and their
proposed Cherenkov ring viewed from a plane perpendicular to the ring’s plane. The hit

111

Figure 6.2: Difference between true and reconstructed vertices for SK-IV FC multi-GeV
multi-ring MC. The shaded region covers 68% of all events. Taken from [92].

Figure 6.3: The Hough transformation method for finding Cerenkov rings. Ring candidates are guessed by drawing circles with 42◦ opening angles surrounding all hit PMTs andintersection points are assumed to be the center of the ring.

112

PMTs are surrounded by dashed circles, which are ring candidates with opening angles of42◦ centered on each hit PMT. The intersection point of many of these candidate rings
is taken to be the center of the final proposed Cherenkov ring. In practice the expecteddistribution of charge around hit PMTs is not a circle with a fixed radius, so instead a more
complicated charge distribution function f(θ) is drawn on a spherical coordinates plane(Φ, Θ) for each hit PMT. The distributions from all initial candidate rings overlap so thatpeaks of charge, rather than intersecting lines, represent likely true ring centers.
The hypothesis of the existence of rings is quantitatively tested with a Likelihood
method. Given N rings previously found, the likelihood for the N+1-th ring is definedby:

LN +1 =

Xlog prob qiobs ,

i

N
+1
X

!!
exp
αn × qi,n

(6.4)

n=1

where the hits from all N+1 PMTs are included in the sum. In this expression qiobs and qiexp
are the observed and expected charge in the i-th PMT and the index n represents the ringindex, so that the sum inside the probability function represents the sum of charges expectedon each PMT from each contributing ring candidate. The coefficients αn are varied so thatthe expression is maximized. The probability function is described by the convulotion of
the single-p.e. response and a Poisson distribution for low charge (q e xp < 20 p.e.) and a
Gaussian distribution at higher charge.
If this algorithm finds a valid second ring, then it will continue searching for a thirdring, and so on, until it reaches a maximum of five rings. N + 1 ring candidates where
LN +1 ≥ LN are considered and if they pass quality cuts based on the p.e. distributions
under both assumptions then they are kept. Efficiencies for ring-counting with this method
are generally found to be above 90% for almost all relevant event samples. Differences
between data and MC are taken as sources of systematic uncertainty in the event selection
for analyses. Plots showing the ring-counting seperation for single-ring and multi-ring FC
events for data and MC are shown in 6.5

113

Figure 6.4: Ring-counting separation for sub-GeV events for SK-I,-II,-III and -IV. Data arein black and oscillated MC (sin2 (2θ23 )=1.0,∆m223 =2.5×10−3 ) are the colored histograms.

114

Figure 6.5: Ring-counting separation for multi-GeV events for SK-I,-II,-III and -IV. Data
are in black and oscillated MC (sin2 (2θ23 )=1.0,∆m223 =2.5×10−3 ) are the colored histograms.

115

Figure 6.6: Event displays for typical (1 GeV) e-like (left) and µ-like (right) MC events.

6.2.1

Particle Identification

A particle identification algorithm called PID is next applied to classify the rings as showering (e-like) or non-showering (µ-like). Showering rings can be produced by electrons as wellas γ-rays, which have diffuse, “fuzzy” ring edges due to electromagnetic shower productionand electron multiple scattering effects. This is because electrons are likely to scatter multiple times along their track and also emit hard photons upon scattering, which may cause
pair-production and further electron scattering. Non-showering rings are produced by more
massive particles such as muons and pions and have a very distinct, sharp ring edge. Event
displays of e-like and µ-like rings of 1 GeV monte carlo events for comparison are shown in
Fig. 6.6.
The PID algorithm first calculates the expected number of p.e.s in each PMT fromeither an electron or muon. In the electron case, a MC simulation of test events with various
momenta is generated with a test vertex at a distance 16.9 m from the wall. These momenta
distributions are used in the calculation of the expected light distribution for e-like events.
The expected number of p.e.s in this calculation accounts for scattering, PMT angular
acceptance, light attenuation and geometric factors such as distance dependence of the116light intensity. For the muon case the expected light distribution is calculated analyticallyusing muon energy loss and Cherenkov radiation models. This calculation includes effects
from knock-on electrons as well as the decrease in the Cherenkov opening angle as themuon loses energy (note that this effect is not included in the electron simulation, where isis assumed electrons are always highly-relativistic and thus have opening angles of 42◦ .
A likelihood function is then calculated from the expected charge distributions of bothe-like and µ-like assumptions for each ring in an event. From this likelihood function aprobability function is calculated, and the overall PID probability distribution is defined as

PP ID =

pp−logPµ − −logPe .

(6.5)

This function is defined so that PP ID < 0 is e-like and PP ID ≥ 0 is µ-like. Distributions
for this likelihood variable for SK-I,-II,-III and -IV data and MC are shown in Fig. 6.7 fora single-ring sample and 6.8 for a multi-ring sample. The separation is generally better for
single-ring events.

6.3

Precise Vertex Fitting

A more precise vertex fitter is applied next. The previously-mentioned fitter, TDC-fit used
only the PMT timing information. This results in a marginal performance for µ-like events
which have longer tracks because the point-source emission assumption for these events
does not accurately reflect time-of-flight corrections. This precise fitter, called MS-fit, also
incorporates the PID and Cherenkov opening angle information. The proposed direction
and vertex are varied while holding the Cherenkov ring angle fixed and the expected lightpattern is compared to the observed pattern. The process is iterated until the goodness-of-fit
along with the PID likelihood are both maximized. The fitter performs well and examplesof the vertex and direction resolutions are shown in Figs. 6.9 and 6.10 for both e-like andµ-like events.

117

Figure 6.7: PID separation for SK-I,-II,-III and -IV FC sub-GeV single-ring samples of data
(black) and MC (histograms).

118

Figure 6.8: PID separation for SK-I,-II,-III and -IV FC sub-GeV multi-ring samples of data
(black) and MC (histograms).

119

Figure 6.9: Distance between true and reconstructed vertex for SK-IV FC MC e-like (left)
and µ-like (right) samples. Sigma (shaded region) represents 68% of the full area.

Figure 6.10: Angle between true and reconstructed direction for SK-IV FC MC e-like (left)
and µ-like (right) samples. Sigma (shaded region) represents 68% of the full area.

120

6.4

Momentum Determination

Precise momenta are also assigned to events after ring-separation and particle identification
have been performed. The momentum of each particle is estimated by the number ofobserved p.e.s found within a cone of half-angle 70◦ surrounding the reconstructed direction.
The observed p.e.s are separated according to their contributions in each ring in order tofind individual ring momenta. The total number of p.e.s is calculated after correcting forlight attenuation, scattering, and PMT acceptance effects, and their sum is taken. The sum
is restricted to the time window from 50 nsec prior until 250 nsec following the peak of thehit timing distribution after subtracting for time-of-flight. These limits exclude effects frommuon decay electrons. The reconstructed momenta are estimated to have resolutions of 1.7
pp+ 0.7/ p GeV/c% and 0.6 + 2.6/ p GeV/c% for electrons and muons, respectively.
6.5

Reconstruction Improvements

Improvements were made to the reconstruction of events recently and these changes wereincorporated into the datasets used in this analysis. A very brief summary of the changes
will be described here. The msfit algorithm was updated after not having been changed forall of SK-I,-II and -III. The photon density table, which is used to estimate the number ofcorrected photons as a function of distance from the charged particle’s track, was limited to
distances of less than 15 m and angles of less than 90 degrees from the reconstructed vertex.
This was improved to include distances up to 40 m and angles within 180 degrees. This
improved the agreement between MC and data for corrected photons at large distances fromthe vertex. This resulted in significant improvements in reconstructed vertex resolutions.
In particular, for Multi-GeV electrons, the 68% vertex resolution improved from 44.5 cm
to 26.7 cm. Additionally, the POLfit program, which is used in reconstructing π 0 events
which decay to two gammas and produce two e-like rings, was also improved recently. The
RNGSEP algorithm, used in reconstruction and identification of rings, was not applied to
the reconstruction of the two gammas in these events. Applying RNGSEP results in a better
resolution of the reconstructed invariant π 0 mass peak. The effect is shown in Fig.6.11, where
the old fitter is shown in blue and the improved fitter in red. The peak and resolution of the121

Figure 6.11: Improvement to the reconstructed invariant π 0 mass with an updated fitting
algorithm shown in red compared to the original algorithm in blue.

π 0 mass were 134.9 and 15.7 MeV/c after the update, compared to 133.0 and 17.7 MeV/c
before. This in turn improved the e/π 0 separation efficiency by reducing the number of π 0
backgrounds leaking into the νe samples while having no effect on the νe efficiency itself.

122

Chapter 7
TWO-FLAVOR OSCILLATION ANALYSIS

There are multiple atmospheric neutrino oscillation analyses that have recently been
updated this year for the first time with the addition of the new SK-IV dataset onto theSK-I+II+III data. They can be divided into two main categories according to their binningschemes - zenith angle analyses and L/E analyses. The L/E analysis focuses on a smaller
subset of the atmospheric neutrino data which has a well-defined L/E (for example, someof the FC and PC but no Upmu data are used). Details of this analysis can be foundin [64]. Unlike the L/E analysis the zenith-angle analysis does not try to define a precise
pathlength and energy to the neutrinos on an event-by-event basis. The L/E analysis isdesigned to produce a better constraint on the mass-squared splitting, while it is found thatthe zenith-angle analysis produces a slightly more restrictive bound on the mixing angle.
The analysis which is the subject of this chapter is a zenith angle analysis which assumes
two flavor oscillations of the channel νµ → ντ . The oscillation probability for this channel
is expressed as
P (νµ → ντ ) = sin2 (2θ23 )sin2



1.27∆m223 L
Eν


(7.1)

where ∆m2 is in units of eV2 , L units of km, and E units of GeV. This analysis is thereforea νµ disappearance search as the majority of the ντ s do not reach the energy threshold of
3.4 GeV necessary for a tau particle to be produced in a weak charge-current interaction.
The two parameters which are constrained by this analysis are the mass-squared splitting∆m223 and the mixing angle θ23 . SK first published evidence pointing to νµ disappearance
consistent with oscillations in 1998[66] and then published an updated zenith-angle analysis
with the full SK-I dataset in 2005[65]. This analysis constrained the mass-squared splitting
to 1.5 ×10−3 < ∆m2 < 3.4 × 10−3 eV2 at 90% C.L. and the mixing angle to sin2 2θ > 0.92.
Enhanced analyses searching for sub-leading effects from θ13 [57] and differences in oscillations between neutrinos and anti-neutrinos[37] using the same zenith-angle binning and123

the SK-I+II+III datasets have since been published. Additionally a full three-flavor analysis
using the same dataset and binning scheme as this work has also been completed recently
and is described in [92]. The search for differences in the signatures of normal hierarchy
versus inverted hierarchy was the motivation for the update to the binning scheme that will
be described in this chapter.
7.1

Analysis SpaceThe range of parameters tested in this zenith-angle analysis is motivated by the current constraints on the two atmospheric parameters by previous SK analyses and other experiments
such as K2K and MINOS[34]. In particular, the most recent published result from the νµ dis−3 eV2 and sin2 2θ
appearance analysis by MINOS[54] found |∆m223 |= 2.32+0.12
23 > 0.90−0.08 ×10

at 90% C.L. In the present analysis we use the parameter ranges 0.001 < ∆m223 < 0.01 eV2and 0.7 < sin2 2θ23 < 1.2. A grid is defined in which each unique pair of (∆m223 ,sin2 2θ23 )is assumed and fit against the data. There were 81 values of ∆m223 tested and they were
distributed on a log scale, while 71 values of the mixing angle were tested on a linear scale
resulting in a total of 5,751 grid points tested.
7.2

Analysis Binning

The data and MC are binned according to their zenith angle direction and event classification based on the information from the reduction and reconstruction methods describedin previous chapters. The binning scheme has been updated from previous SK zenith angle
analyses. The number of bins has been increased from 420 to 480. Figure 7.1 is a schematic
showing all the bins in the new binning scheme. The vertical axis represents increasing
momentum with increasing height. All samples, except for the shaded ones in figure, areseparated into 10 zenith-angle bins for each momentum bin. The binning is chosen suchthat there at least five expected MC events in each bin after scaling to data livetime. Thisguarantees the χ2 function, which will be introduced later in this chapter, can be properly
applied in this analysis.
The additional 60 bins come from a new event classification scheme in the FC multiGeV single-ring and multi-ring samples. These events are now binned according to their

124

Figure 7.1: The layout of bins used in the zenith-angle analyses. There are 370 for FC, 60
for PC, and 50 for Upmu events, giving 480 in total. White boxes indicate those for which
10 zenith bins are defined and shaded indicates single zenith-angle bins.

νe -like and ν̄e -like characteristics. Though SK cannot distinguish events as being induced byneutrinos or anti-neutrinos on an event-by-event basis, there are certain characteristics ofthe FC multi-GeV sample that allow it to be separated with reasonable purity. This binningscheme is motivated by the full three-flavor analysis which seeks to improve sensitivity to
the neutrino mass hierarchy[92]. An enhancement of electron neutrino appearance from thematter effect occurs for neutrinos in the normal hierarchy and ant-neutrinos in the inverted
hierarchy. The multi-GeV νe -like and ν̄e -like event separation criteria are as follows. For
single-ring events, a simple separation where all events with 0 decay electrons are consideredν̄e -like and events with 1 or more decay electrons are considered νe -like is used. The CC
interactions of νe occur as:νe + N → e− + N 0 + π +
π + → µ+ + νµ

125Figure 7.2: A typical CCνe interaction which produces an e-like multi-ring event in SK.

The µ+ from the pion decay then decays via µ+ → e+ + νe + ν¯µ and the e+ can bedetected as a decay electron. The splitting of these samples generates reasonable purities of
62.8% and 36.7%[92] when tested with SK-IV MC. For multi-ring events, a more complicatedselection using a likelihood function is applied. A typical interaction producing a multi-ring
event is shown in Fig. 7.2. For these interactions the energy transfer is typically larger
for CCνe than it is for CCν¯e and as a consequence, the following facts can be used inconstructing a likelihood function to separate the sample:
1. CCνe events have a smaller momentum fraction for the most energetic ring
2. CCνe events have a greater number of Cherenkov rings (from more outgoing pions)
3. CCνe events have larger transverse momentum
4. CCνe events have more decay electrons
The likelihood separation is found to produce νe -like and ν̄e -like samples with 59.4% and21.0% purities, which are modest values but are increases of a few percent above the original
purities found in the non-separated sample[92].
The two flavor analysis assumes no oscillations involving electron-flavor neutrinos andanti-neutrinos, so there is no matter effect or sensitivity to the hierarchies in this work. However, the analyses are otherwise similar (they both use the same dataset, MC, zenith-angle
binning and chi-squared statistic) and we use the same binning scheme here for consistency.
The effect of splitting the bins is at most very small and would only come from muon-flavor

126

impurities (which can undergo oscillations) in the FC multi-GeV.
7.3

Monte Carlo

The monte carlo samples used in this work are 500-year samples of the FC, PC and tau,and 200-years of upmu event simulation. The MC samples used here have been newly
regenerated from previous MC samples and have not been used in official analyses until now.
They are produced with the Honda 2011 flux model[49] as well as the other modifications
made to NEUT, skdetsim, and event reconstruction described previously. Oscillations are notapplied to the MC during initial production. The monte carlo events are each individuallybinned according to their zenith angle, momentum and event classification, and after thisbinning each event is reweighted by the oscillation probability calculated from the energy
and pathlength (which is determined from the true zenith angle direction) truth information.
The pathlength calculation also includes an averaging estimation of the neutrino production height in the atmosphere. For each event 20 oscillation probabilities from a sample
of production heights are calculated and the average of these is taken as the final oscillation
weight. The events are also reweighted to account for the solar wind’s modulation of neutrino fluxes. Primary cosmic ray flux is at a minimum (maximum) when the solar wind ismost (least) active and this affects the production of atmospheric neutrinos. After binning
the weighted MC each combination of oscillation parameters is tested in this analysis sothat for every pair of parameters, a uniquely weighted set of monte carlo including all event
samples is binned with all the data and compared.
7.4

DatasetThe datasets used in this analysis include SK-I,-II,-III, and -IV. Up to and including the
SK-III phase, the livetimes for the FC and PC data were not equal to that of the Upmudata. The FC and PC data are much lower energies than the Upmu data and so are more
susceptible to problems in the detector. In SK-IV (and thereafter) the livetimes for all
three samples are equal. Table 7.4 shows the livetimes for all datasets used in this analysis.
The number of data and oscillation-weighted,livetime-normalized MC events in each event
sample is listed in table 7.4.1.

127

Data sample

SK-I(days)

SK-II(days)

SK-III(days)

SK-IV(days)FC

1489.2

798.6

518.1

1096.66PC

1489.2

798.6

518.1

1096.66

Upmu

1645.9

827.7

635.6

1096.66

Table 7.1: A summary of the livetimes (units of days) for the data used in this analysis.

7.4.1

Event Categories

The official event categories for the analysis dataset were as follows. FC events were divided
into Sub-GeV (visible energy < 1.33 GeV) and Multi-GeV (visible energy > 1.33 GeV). They
were also subdivided according to various combinations of their ring PID (e-like, µ-like or
π 0 -like), the number of rings found (single or multi), the number of decay electrons, andtheir νe -like or ν̄e -like probability. PC events were divided into stopping and through-going.
PC stopping events were put into two momentum bins, while PC through-going were put
into four momentum bins. The upmu sample was divided into stopping (3 momentum bins),
through non-showering and showering (no momentum binning for through-going events). In
most cases there 10 zenith divisions used, with FC and PC events binned into the range-1.0 < cosθz < 1.0 and upmu events in the range -1.0 < cosθz < 0.0. The samples whichwere put into one large zenith bin are shown as the shaded bins in Fig. 7.1.

7.4.2

Two-Flavor Oscillation Effects

The two-flavor oscillation channel νµ → ντ exhibits muon-type neutrino disappearnce andassumes no νe appearance. Assuming only mixing between the νµ and ντ implies that theµ-like samples will undergo oscillations while the e-like samples should be largely unaffected.
A 2-dimensional plot of the predicted oscillation effects as a function of neutrino energy andpathlength assuming the SK-I+II+III best-fit values is shown in Fig. 7.3.

128

SK-I

SK-II

SK-III

SK-IV

Data

MC

Data

MC

Data

MC

Data

MC

0-decay e

2992

2705.41573

1445.4

1092

945.3

2098

1934.9

1-decay e

301

248.1

172

138.9

118

85.3

243

198.4

π 0 -like

176

160.0

111

96.3

58

53.8

116

96.2

0-decay e

1025

893.7

561

501.9

336

311.8

405

366.3

1-decay e

2012

1883.01037

1006.7

742

664.1

1833

1654.1

2-decay e

147

130.4

86

71.3

61

46.6

174

132.2

2-ring π 0 -like

524

492.8

266

259.8

182

172.2

380

355.9

νe -like

191

152.8

79

78.4

68

54.9

156

135.9

ν¯e -like

665

656.2

317

349.5

206

231.6

423

432.8

µ-like

712

775.3

400

415.7

238

266.4

420

554.8

νe -like

216

224.7

143

121.9

65

81.8

175

161.9

ν¯e -like

227

219.7

134

121.1

80

72.4

212

179.1

µ-like

603

640.1

337

337.0

228

231.4

479

499.0

Stop

143

141.4

77

71.3

54

52.0

109

115.0

Thru

759

772.4

350

377.3

290

304.0

626

629.9

Stop

429

425.5

210

211.6

193

166.3

284

276.4

Non-showering1564

1340.1

725.3

689.6

612.9

497.8

1048

863.1

Showering

271.6

294.7

108.1

109.2

110

127.7

282.9

217.0

FC Sub-GeV
1-ring e-likeµ-likeFC multi-GeV
1-ringmulti-ring

PC

Upmu

Table 7.2: A summary of the data and MC samples (with their oscillation weights applied)
used in this analysis. The MC samples are normalized to the data livetimes.

129

Figure 7.3: Contours of the probability of νµ survival as a function of pathlength and energy
assuming two-flavor oscillations with (∆m223 = 2.11 ×10−3 eV2 , sin2 (2θ23 ) = 1.00). Red
areas indicate no oscillations and green areas indicate maximum oscillations.

7.5

Chi-Squared Analysis

The likelihood used in this analysis is based on Poisson statistics in order to reflect the finebinning. The Poisson distribution for the number of expected events E given a number ofobserved events O is given by
P (E, O) =E O −E
e .
O!

(7.2)

There is in general an expected contribution of 1 unit from each bin to the reduced χ2 inthe case where a small number (∼ 1) events are expected in each bin. Assuming N total
bins with index n a likelihood of the form is constructed:
Y EO
n −En
L(E, O) =
e
.
O
n!
n

(7.3)

Taking the log likelihood ratio of this expression (the ratio of the quantities expected given
observed and observed given observed) the expression becomesX
On
L(E, O)=2(En − On + On ln)
χ2 ≡ −2ln
L(O, O)
E
n
n

(7.4)

To incorporate the effects of independent systematics into the analysis, the so-called “pull”
method for the χ2 minimization is used[75]. This is equivalent to introducing a shift in the

130

number of expected events in each bin that is linearly proportional to a systematic error
parameter. Mathematically we have
En → En (1 +X

fni i )(7.5)

iwhere i is the parameter for systematic error i and fni is the fractional change (which isdependent upon individual bin contents and calculated with the MC ahead of time assuming
the best-fit oscillations from [65]; For more details see [112]) that i effects on the numberof expected events. For each bin the effects of all systematic errors are summed (hence thesum over i) while allowing the i parameters to vary. The i parameters are varied so thatthe χ2 is minimized during the fit and the distribution of i /σi , where σi is the expected
1-σ variation of each systematic uncertainty, is expected to be Gaussian. A penalty term
is added to Eq. 7.4 to constrain the size of the effect that the systematic error parameters
have on the fit, so that the final χ2 statistic takes the form:
X
On
χ =2
(En − On + On ln
)+
En
n
2

X

i

2

χ =2

X
n

En → En (1 +

X
i



i
oi

2

On
P
fni i ) − On + On ln
En → En (1 + i fni i )

(7.6)

!
+

X  i 2
i

oi
(7.7)

With this equation in hand, the minimization of the χ2 becomes equivalent to solving the2equation ∂χm = 0 for all m. This minimization is analogous to solving a system of linear
equations, one for each systematic error parameter (for details see [65, 112]). After theminimization is done and the systematic error parameters are evaluated the data from theindividual SK-I,-II,-III and -IV bins are merged together.

7.5.1

Systematic Error Terms

There are 150 i systematic terms in total in this analysis, including terms common to allrun periods (physics-related errors) and those unique for each run period (detector systematics). There are 19 terms each for neutrino flux uncertainties and neutrino interaction
uncertainties, 68 related to event reduction (17 for each run period × 4 run periods) and 44

131(11 × 4) related to event reconstruction and selection. The resulting pull values from thefit will be shown later in this chapter.

7.5.2

Neutrino Flux Errors

Absolute Normalization
The absolute normalization of the neutrino flux is based on the work of Honda et al[48]which is calibrated with atmospheric muon data. In that work they calculate an energydepedent normalization uncertainty which is around 25% for low energies, decreases to to
below 10% for energies between 1 and 10 GeV, and then increases again steadily for higher
energies. A plot of this uncertainty is shown in 7.4. Two uncertainty terms for this enter theχ2 minimization, one for E < 1 GeV and one for E > 1 GeV, and the energy dependenceof these terms reflects the Honda calculation. The individual uncertainties of the pion
production in hadronic interactions (δπ ), kaon production (δK ), hadronic interaction crosssections (δσ ) and the density profile of the atmosphere (δair ) entering the calculation are
shown in the top plot in Fig. 7.4.

Flux Ratios
The value of the flavor ratio (νµ + ν¯µ )/(νe + ν¯e ) is known much better than the absolutenormalization itself. The uncertainties used here are calculated from comparisons of the
Honda[48], FLUKA[4] and Bartol[3] models. A plot of the ratios of these models to the
Honda model is shown in Fig. 7.5. The uncertainty is estimated here to be 2% for E < 1
GeV, 3% for 1 GeV < E < 10 GeV, 5% for 10 GeV < E < 30 GeV, and then increasing loglinearly in E up to 30% at 1 TeV. The convention here is that a positive value for the error
parameter  generates an increase the fraction of muon neutrinos. The ratio of neutrinos toanti-neutrinos νµ /ν¯µ and νe /ν¯e depends on the production ratios π + /π − and K+ /K− fromthe interactions in the flux calculations. This uncertainty is also estimated by comparing
the flux models. The electron ratio is taken to be 5% below 10 GeV, 8% between 10 and
100 GeV, and a log-linearly increasing function which reaches 30% at 1 TeV. The muon
ratio is smaller (2%) at low energies up to 1 GeV, 6% between 1 and 50 GeV, and then an

132

Figure 7.4: The components entering the absolute normalization uncertainty (top) and theoverall uncertainty (bottom) as functions of energy[48].

Figure 7.5: A plot of the ratios of the FLUKA[4] (dashed red) and Bartol[3] (solid black)flux models to the Honda[49] model.

133

increasing value reaching 60% at 1 TeV.
The up/down flux ratio also has an energy-dependent uncertainty from the geomagnetic
field. For lower energies this uncertainty has less of an effect because the angular correlation of the charged lepton with the parent neutrino gets worse with decreasing energy. A
comparison of the zenith angle distributions of the three flux models is used to estimate
the uncertainty, which is taken to be 1%. The ratio of horizontal/vertical flux is also assigned at 1% uncertainty. The uncertainty comes from the differences of the 3-dimensional(or 1-dimensional) calculations of the flux in each of the models at lower energies, but isdominated by the production ratio of K/π at higher energies. This error is also estimatedby comparing the zenith-angle distributions of the three flux models.

Production Ratios
Atmospheric neutrinos come primarily from pion decays at low energies. However at primarycosmic ray energies of a few tens of GeV and greater the component of K decay becomesincreasingly large so the K/π ratio is taken as a systematic error. Recent measurementsof the ratio[28] indicate an uncertainty of 3%. From these results and the dependence ofneutrino energy on the hadron momentum the uncertainty in the K/π ratio is estimated to
be 5% below 100 GeV and linearly increasing to 20% as a function of energy from 100 GeV
to 1 TeV.
The ν production height is used in the calculation of the pathlength and this is a large
percentage of the pathlength for downward and horizontal neutrinos, so a systematic error
term is added for those events. The error is estimated by assuming a 10% variation in thedensity structure of the atmosphere (determined from comparisons of US-Standard ’76 [98]and MSISE90[84]). The difference in the neutrino flux from path length variations in the
different density structures is taken as the systematic error and its value is estimated to be
10%.
The primary cosmic ray energy spectrum is well modeled by the form E −2.74 with alarger uncertainty for larger energies (i.e. above 100 GeV). Differences in the energy spectrum among the flux models considered are not sufficient to account for the different fluxes

134

calculated in the respective models, and so a sample-by-sample normalization error is applied to each of the FC, PC and Upmu samples to account for these differences. The error
is set at 5% for each of the samples.
Solar flux activity is also taken into account when modifying the flux predictions in themonte carlo. The primary cosmic ray flux is modified on an 11-year cycle by the solar windwith an uncertainty of ±1 year. This corresponds to errors on the solar activity term of
20%, 50%, 20% and 10% for SK-I,-II,-III and -IV, respectively.

Neutrino Interaction Systematics
Uncertainties in cross-sections and other terms related to the interactions are common
across SK-I,-II,-III and -IV. The axial mass parameter, relevant for CCQE and single meson
production interactions, is assumed to be 1.21 GeV in the MC production described in thiswork, with an uncertainty of 10%. The 10% comes from the difference in the spectra ofsquared-momentum transfer due to between this axial mass value and a value of 1.1 GeVfrom previous experimental measurements and used in previous SK analyses.
The CCQE interaction cross-section uncertainty for bound neuclei is estimated from thedifference in the q2 distribution of the Fermi gas model of Smith-Moniz model used in thiswork and that of Nieves et al[35]. It is taken to have a 1-σ width of 1%. The uncertainties
for the CCQE cross-section ν̄/ν ratio and (νµ + ν̄µ )/(νe + ν̄e ) flavor ratio are estimated in
the same manner and taken to be 1.0% each. Uncertainty of single meson production isestimated by comparing expected to experimental results of νµ p → µ− pπ + and is taken to be
20%. The production of single π 0 via the interaction νµ n → µ− pπ 0 is less well understood.
An estimate of 40% for the ratio π 0 /π ± uncertainty is used. This is estimated by comparing
the models used in this work to the Hernandez model[85]. Comparisons to the Hernandez
model are also used to estimate the ratio ν/ν̄ in single meson production, which is set to
1% for both e and µ.
Deep inelastic scattering uncertainties are divided into two regions: above and below 10
GeV. Above 10 GeV predictions agree with experiment to within 5% and this is used as
the uncertainty. Below 10 GeV comparisons with the CKMT model are used and the 1-σ135

Figure 7.6: Predicted cross-sections (left) and ν̄/ν cross-section ratio for νµ p → µ− pπ +(solid), νµ n → µ− pπ 0 (dashed) and νµ n → µ− pπ + (dotted) for the Hernandez model[85]normalized to the Rein-Sehgal model[101].

width is set to 1.0%. The uncertainty of the squared-momentum transfater Q2 distribution
is estimated by comparing the present model to the original GRV98 model. Coherent pion
production is not well-known. The error for NC and CC νe interactions is calculated to be
50% from the difference between the Rein-Sehgal[101] and experiment. The CCνµ coherent
pion production channel is completely unmeasured and so the uncertainty is set to 100%.
The overall NC/CC ratio uncertainty is included and is set to be 20%. Nuclear effects inoxygen were updated in the MC sample as described in a previous chapter and the methodof calculating their uncertainty was updated as well. The cross-section of the pion final
state interactions in the updated NEUT model is compared to experiment to estimate theuncertainty[17]. A new weighting method is applied which weights all possible interaction
channels until a 1-σ difference between the simulation and data are found, as shown in 7.7.

7.5.3

Event Selection and Reconstruction Uncertainties

FC/PC Reduction
Errors are calculated for the event reductions by looking at the differences between data
and MC for distributions of the variables used in the selection cuts. The FC reduction
uncertainties are 0.2%, 0.2%, 0.8% and 0.3% for SK-I,-II,-III and -IV, respectively. For PC,
the same reduction uncertainties were 2.5%, 4.8%, 0.5% and 1.0%. The FC/PC separation

136

Figure 7.7: Comparisons between data and NEUT simulations for π + quasi-elastic (left)
and absorption (right) interactions demonstrating the weighting method used in the nuclear
effects uncertainty term.

is based on the number of hit OD PMTs in a cluster and the uncertainty is based oncomparing the MC and data distributions of the NHITAC variable. These terms are assumed
to be anti-correlated because a change in the number of FC events will simultaneously resultin a change in the opposite direction for the number of PC events. This uncertainty for eachrun period is estimated to be 0.6%, 0.5%, 0.9% and 0.02%. Similarly, the PC stop/thru
separation relies on OD charge distributions (through the ratio PEanti /PEexp ) and terms
for the top, barrel and bottom are separately included in the fit.

Upmu Reduction
Uncertainties for the upmu reduction are classified according to the three event categories:
stop, thru and showering. Reduction uncertainties are calculated by fitting the reduction
efficiency for MC events as a function of zenith bin to a constant and its uncertainty fromthe fitting is taken as the error. The errors for stopping are 0.7%, 0.7%, 0.7% and 0.5%,
and through-going 0.5%, 0.5%, 0.5% and 0.3% for SK-I,-II,-III and -IV, respectively. The
stop/thru separation is also treated as a separate error. This is calculated by varying the
cut value of the variable Um ohit8m (the number of hits within 8 meters of exit point in theOD) by ± 10% and calculating the percantage change in the stop/thru ratio for MC events.

137

These values for each run period are 0.5%, 0.8%, 0.7% and 1.0%.
Uncertainty from the Upmu stopping selection cut is accounted for by varying the reconstructed momentum cut value used in the selection (Fit mom) by ±∆E (the energy scale
uncertainty) and calculating the percentage change in the number of surviving MC events isthe uncertainty estimate. These magnitudes of the errors are 0.9%, 1.3%, 2.0% and 1.7% for
SK-I,-II,-III and -IV. For the through-going sample, the reconstructed pathlength (Fit len)cut is varied and the percentage change in the number of through-going events is used as
the uncertainty. These values are 1.5%, 2.3%, 2.8% and 1.5% for each detector phase.
Finally the separation of the showering and non-showering samples is also accounted forby varying the shifting the Sh delta cut for MC events an amount based on the difference
between the data and MC distributions for the same variable. The resulting percentagechange in the showering/non-showering ratio is used for the uncertainty. The estimatedvalues were 11.7%, 11.5%, 6.4% and 10.7% for SK-I,-II,-III and -IV, respectively.

Contamination Backgrounds
A handful of errors accounting for contamination backgrounds in the event samples are also
used in the analysis. Examples of these contaminations include propagation of hardons
which are reconstructed as µ-like for single-ring NC events, low-energy backgrounds fore-like events such as flasher PMTs and neutron-induced events, and contamination of allµ-like events from cosmic ray muons. Additionally the upmu sample undergoes backgroundsubtraction as described in the data reduction chapter and from this subtraction the uncertainty in the number of background events in the most-horizontal bins is accounted for inthe χ2 minimization.

7.5.4

Multi-GeV multi-ring separation

The new binning method includes the separation of the multi-GeV multi-ring sample into νe like and ν¯e -like categories (the single-ring sample is also separated but the selection criteria
is much simpler than for multi-ring, using only the number of decay electrons). A systematic
uncertainty is assigned to this separation based on the contents of the separation likelihood,

138

Figure 7.8: Pion multiplicity as a function of W2 for NEUT predictions (left) and CHORUS
experimental measurements (right).

which are the number of decay electrons, the number of rings and in turn pion effects (final
state interactions, multiplicity, etc.). Pion final state interaction uncertainties are included
as mentioned above by the newly-introduced MC weighting methods. Uncertainties in pion
multiplicity are estimated by comparing the NEUT predictions to data from the CHORUS
experiment[23]. Plots showing the pion multiplicity as a function of W2 (W is the invariant
mass of the final state particles) for NEUT predictions and CHORUS data are shown in
7.8.
A pion multiplicity weight is added to each MC event in addition to its FSI cross-section
weight. The difference in the CCνe purities of the original likelihood and the weighted
likelihood is than calculated and taken as the systematic uncertainty. The values are foundto be 7.16%, 7.91%, 7.68% and 6.82% for SK-I,-II,-III and -IV, respectively.

7.5.5

π 0 -like selection

The identification of π 0 -like events are dependent upon the π 0 mass reconstruction andlikelihood distributions. Uncertainties here are again estimated by comparing MC and data
distributions with separate error terms for each momentum region (because the π 0 likelihood
distribution depends on momentum). The uncertainty also depends on the energy scale
uncertainty and the calculation includes this contribution as well. Two-ring π 0 -like events
are also part of the final dataset and are dependent upon ring-counting and PID algorithms.

139

Their uncertainties are estimated to be 2.0% for SK-I and -II, 1.0% for SK-III and 1.8% for
SK-IV.

7.5.6

Decay Electron Tagging

A cosmic ray muon sample is used to estimate the efficiency of tagging decay electrons.
This study was completed in SK-I[67] and is assumed to be 1.5% for all detector periods.
A component from π + decay is also considered because positive pions have a probability
to decay before interacting in water. This decay probability depends on their momentum
but can be more than 10% for some momenta. This decay uncertainty is included in the
decay-electron tagging efficiency for the FC sub-GeV single-ring e-like sample.

7.5.7

Other Reconstruction Uncertainties

Ring separation and PID likelihood uncertainties are estimated by comparing their respective likelihood MC distributions to data. The absolute energy scale is also assigned an
uncertainty, and the values are 1.1% for SK-I, 1.7% for SK-II, 2.7% for SK-III and 2.3% for
SK-IV. This uncertainty also indirectly enters into the upmu stopping sample selection cut
mentioned in earlier in this chapter. The energy scale uncertainty parameter is defined suchthat a positive value will increase the amount of visible energy in the MC. Uncertainty inthe up/down asymmetry in the energy scale is also included and these values are 0.6% for
SK-I and -II, 1.3% for SK-III and 0.3% for SK-IV.
7.6

Results

The results of the fit to the SK-I+II+III+IV full dataset are shown in Figs. 7.9 and 7.10
for the physical region. The 2-dimensional contour in Fig. 7.9 shows the allowed regions
at 68% (blue), 90% (red) and 99% (green). The allowed regions represent regions in which∆χ2 (≡ χ2 − χ2min , or the difference between the χ2 from the fit and the minimum χ2 atthe best-fit point) satisfies the conditions ∆χ2 ≤ 2.3, 4.6 and 9.2 for 68%, 90% and 99%
C.L., respectively. The best-fit point is found at ∆m223 , sin2 (2θ23 ) = (2.30 ×10−3 eV2 ,0.99) with χ2 /d.o.f = 555/478 and is indicated by the black star in Fig. 7.9. The 1σ and140

Figure 7.9: The 2-d allowed regions for ∆m223 versus sin2 (2θ23 ) at 68%, 90% and 99% C.L.
for SK-I+II+III+IV data. The black star indicates the best-fit point.

141

Figure 7.10: The 1-dimensional χ2 distributions for ∆m223 (left) and sin2 (2θ23 (right, physicalregion only) fit to the SK-I+II+III+IV data.

2σ constraints on the individual parameters are determined by their 1-d χ2 distributions,which are shown in Fig.7.10. The individual parameter constraints are found to be ∆m223−3 eV2 (1σ) and sin2 (2θ ) > 0.95 (90% C.L.).
= 2.30+0.17
23−0.22 × 10

The allowed regions for this analysis are found to be similar to, but shifted relative to theprevious SK-I+II+III official result. In particular the allowed region has moved upward in
the mass-squared splitting space. This effect is mostly from the addition of the SK-IV data
to the dataset, which prefers higher values for ∆m223 . This analysis gives excellent agreementwith the most recent results from the other active analyses of the same parameters. An
overlay of all 2012 results is shown in Fig. 7.11. The physical allowed regions for eachindividual SK dataset (SK-I,-II,-III and -IV) fit to 2-flavor mixing are shown in 7.12 and7.13. As the MC was completely regenerated to incorporate the updated NEUT flux models,pion FSI interaction models, and reconstruction updates, the fits to the older datasets (SKI,-II and -III) are also updated.
The zenith-angle distributions of all the data samples (with SK-I+II+III+IV binned
together) are shown in Fig. 7.14. In this figure the red histogram shows the MC oscillatedwith the best-fit parameters from this analysis and the blue histogram represents the no142

Figure 7.11: The 90% C.L. allowed regions of the most recent 2012 results for this analysis,
the SK L/E analysis, the SK full three-parameter zenith-angle analysis, T2K and MINOS.

Figure 7.12: The 2-dimensional 68%(blue), 90%(red) and 99%(green) C.L. allowed regions
for the SK-I data (left) and SK-II data (right) with the updated MC and reconstruction.

143

Figure 7.13: The 2-dimensional 68%(blue), 90%(red) and 99%(green) C.L. allowed regions
for the SK-III data (left) and SK-IV data (right) with the updated MC and reconstruction.

oscillations prediction.
7.6.1

Systematic Error Pulls

The results of all the systematic error pulls are shown in Figs.7.15 and 7.16. In each ofthese figures the vertical axis represents the value of i divided by the 1σ uncertainty forthe i-th error term. The majority of pull terms are well within their 1σ estimation and thelargest pull terms are found to be on the order of 1.5σ.

144

Figure 7.14: The zenith-angle distributions of the SK-I+II+III+IV data (black points),
the SK-I+II+III+IV MC oscillated with the best-fit parameters from this analysis (∆m223 ,sin2 (2θ23 )) = (2.30 ×10−3 eV2 , 0.99) (red histogram) and the MC with no oscillations (bluehistogram).

145

Figure 7.15: The ratios i /σi for systematic error terms from the 2-flavor fit to the SKI+II+III+IV data.

146

Figure 7.16: The ratios i /σi for systematic error terms from the 2-flavor fit to the SKI+II+III+IV data.

147

Chapter 8
APPARENT CPT VIOLATION MODEL
8.1

Introduction

The present atmospheric neutrino oscillation data from experiments worldwide are well-fit
by the simplified 2-flavor oscillation model, as discussed in detail in previous chapters. Thismodel assumes that only oscillations of the form νµ → ντ occur. However, motivated byresults from the LSND[69] and, more recently, MINOS[73, 72] experiments, theorists havesuggested more exotic models which reproduce the effects of standard 2-flavor atmosphericmixing as well as simultaneously incorporating additional effects these experiments claim
to observe. Common features of these models include the introduction of sterile neutrinos,
matter effects, and CPT-violation in neutrino mixing.
One such model, proposed at the University of Washington by Nelson, Engelhardt, andWalsh[21], attempts to explain both standard atmospheric mixing and recent anti-neutrino
oscillation results by introducing a sterile neutrino and a non-standard matter interaction potential without introducing CPT-violation. The model was primarily motivated byMINOS[53], a long-baseline accelerator experiment which has constrained both neutrino
and anti-neutrino oscillation parameters.
MINOS’ initial anti-neutrino constraints, taken from a low-statistics sample of antineutrino events found in a neutrino beam, did not agree well with their corresponding
neutrino parameters[73, 72]. The constraints were later updated with larger statistics froma dedicated anti-neutrino beam running period, and the discrepancy between neutrinos andanti-neutrinos has since disappeared[55]. The model which is the subject of this analysis
was motivated by the earlier results from MINOS. There currently exists no compelling
evidence for a difference between neutrino and anti-neutrino oscillations. Nonetheless, thismodel remains worthy of consideration as it predicts interesting effects in neutrino- and antineutrino oscillations in matter without invoking CPT-violation while claiming compatibility

148with existing experimental oscillation results.
8.2

Mathematics of Model

The main features of the model are the addition of a sterile neutrino which participatesin active-sterile mixing and an exotic matter interaction mediated by a newly-introduced
gauge vector boson. With the addition of the sterile neutrino state, the mixing matrix
now expands to a 4x4 matrix to include the new mixing angle θ34 . The mixing angles
θ12 and θ13 are assumed to be negligible because they do not contribute significantly tomixing at the MINOS neutrino energies and baselines. These assumptions also apply tothe atmospheric neutrinos detected at SK. Additionally, maximal atmospheric mixing isassumed with θ23 = π4 . The resulting simplified 4x4 matrix which transforms between themass and flavor bases takes the form


1 0

0

0







 0 1

0
0


Us = 

 0 0 cosθ34 sinθ34 
0 0 −sinθ34 cosθ34(8.1)

where the angle θ34 quantifies mixing between the heaviest active state and the sterile
state.
The neutrino mass hierarchy also expands to include the sterile neutrino as the fourthmass state. The model then simplifies this expanded mass hierarchy to

Mm = diag(0, 0, m, M )(8.2)

by assuming that m1 = m2 = 0. This assumption is analogous to neglecting solar terms
in the oscillations. This assumption is valid for long-baseline accelerator and atmospheric
neutrino oscillations. The remaining two non-zero mass states correspond to the largest
active neutrino mass, m, and the newly-introduced sterile neutrino of mass M.
The model also introduces a new matter interaction for neutrinos. From standard model
electroweak interactions, neutrinos traveling in matter experience charged current and neutral current potentials of the forms

149

√
GF
VCC = ± 2GF Ne and VNC = ± √ NN ,
2

(8.3)

respectively. Here GF represents the Fermi constant, Ne the local electron density, andNn the local neutron density. The model’s new B − L gauge interaction is mediated by a
new vector boson, previously proposed by Nelson and Walsh[97] to reconcile the results of
the MiniBooNE[70] and LSND[69] experiments with the standard oscillations picture. This
interaction causes neutrinos to experience an additional potential in matter,

VB−L = ±

g̃ 2
Nn
m2V

(8.4)

where g̃ represents the coupling constant for the new interaction and mv the mass of theinteraction’s gauge boson. This term is positive (negative) for active (anti-)neutrinos andnegative (positive) for sterile (anti-)neutrinos, and enhances the difference between neutrino
and anti-neutrino oscillations in matter originally coming from the MSW effect[95, 113].
Combining these potentials, the full potential matrix in the flavor basis can be writtenasVf = (Vcc − Vnc − VB−L , −Vnc − VB−L , −Vnc − VB−L , VB−L ).

(8.5)

The assumptions made so far in the model lead to mixing between the sterile neutrino state and a combination of νµ and ντ . For atmospheric oscillation studies at SuperKamiokande, the main oscillation channel is assumed to be νµ → ντ . This muon neutrino
disappearance is expressed mathematically by the muon neutrino survival probability, as
shown in Eqn. 7.1. The muon neutrino survival probability in this new model can be
expressed in the interaction basis, where the overall Hamiltonian is diagonal, as

P (νµ → νµ ) = 1−cos2 (θs )sin2




m̃2 L
−sin2 (θs )sin2
4E!
M̃ 2 L−cos2 (θs )sin2 (θs )sin2
4E(M̃ 2 − m̃2 )L
4E
(8.6)

The variables m̃2 and M̃ 2 in this equation are the effective mass-squared splittings in theinteraction basis. Likewise, θs is the effective sterile mixing angle in the interaction basis.

!
.

150

The reader can refer to [21] and the Appendix for more details. As seen in this equation thisνµ survival probability allows for more exotic oscillation effects to appear than in standard
2-flavor mixing, especially in the presence of matter which affects the values of M̃ and m̃.
8.3

Fit to MINOS Data

The creators of the model constrained it with MINOS’ early muon anti-neutrino data sample,
which came from anti-neutrino contamination of the muon neutrino beam[73]. For thisanalysis the parameter V = 12 VN C + VB−L was defined to quantify the matter interaction.
The best-fit values were found to bem = 0.0394 eV, M = 0.157 eV, and V = 2.01 × 10−4 neV.

(8.7)

g
For this analysis the local Earth density was assumed to be a constant ρ = 2.7 cm
3 , whichis the matter density in the crust. The best-fit potential V is larger from the contribution
expected from SM interactions alone, which is VSM = 0.25 × 10−4 neV. The best-fit value
for the active neutrino mass is smaller than the square root of the 2-flavor atmosphericpmass-squared splitting ( ∆m223 ≈ 0.049eV). These two values should become equal in thelimit where new exotic oscillation effects are removed from the model.
The allowed regions at 68.3% (blue), 99.5% (red), and 99.7% confidence limits for eachpair of parameters are shown in Fig. 8.1. There are two main allowed regions in this analysis.
The first includes the best-fit point and has a constant active neutrino mass near sqrtδm223 ,a sterile neutrino mass M which is large relative to m (minimal active-sterile mixing), andvalues of the coupling constant near the Fermi constant. The second region extends to muchsmaller values of m, smaller values of M (more sterile mixing), and larger values of the newinteraction’s coupling constant. In this allowed region with a great contribution from the
new interaction m and M are inversely related.
8.4

Effects on Atmospheric Neutrinos at SK

The atmospheric neutrino data collected at SK cover a wide range of energies and baselines,
unlike the MINOS data to which this model was originally fit, which cover a relatively narrow

151

10

10

5

1
0.1
V !neV"M !eV"

2
1
0.5

10!3
10!4

0.2
0.1
10!4

10!2

4"10!410!3

4"10!3
m !eV"

10!2

4"10!2

10!1

10!5
10!4

4"10!410!3

4"10!3
m !eV"

10!2

4"10!210!1

10
1

V !neV"

0.1
10!2
10!3
10!4
10!5
0.1

0.2

0.5

1
M !eV"

2

5

10

Figure 8.1: Allowed regions from original fit to MINOS anti-neutrino data within 68.3%(blue), 99.5% (red), and 99.7% confidence limits for each pair of parameters, with the best
fit point indicated by the star.

152

Figure 8.2: Earth density profile as a function of radius according to PREM[20] on the left,
and the average line-integral density as a function of cosine zenith angle on the right.

distribution of neutrino energies at a fixed baseline of 735km[53]. The model’s oscillations
are sensitive to both neutrino energy and baseline (see Eqn. 8.6), so atmospheric neutrinos
at SK offer a good opportunity for constraining the model. Because the density of the Earthvaries with radial distance, the zenith angle of the incoming neutrino, which determines its
pathlength, specifies the Earth matter density along the path. In order to calculate theaverage matter density traversed by each incoming neutrino, a radial density profile for the
Earth given by PREM[20] was assumed. The earth density profile from this model is shownin Fig. 8.2.
Average density as a function of zenith angle was calculated by taking the line-averaged
integral of the density along the path for all zenith angles. The average matter potential
for each neutrino was calculated by assuming this line-averaged density to be a constantdensity along the path. The matter potential as a function of zenith angle is shown on the
RHS in Fig. 8.2.
The model’s variables m̃2 , M̃ 2 and θs as a function of energy, assuming the parametersin Eqn. 8.7 are shown in Figs. 8.3 and 8.4, respectively. These plots assume an Earth153

Figure 8.3: Effective mass squared splittings in matter as a function of energy in the apparent
CPT-violation model. The 2-flavor atmospheric mass squared splitting is shown in the leftg
hand figure for comparison. A density of ρ = 4.0 cm
3 for the Earth is assumed.

gdensity of ρ = 4.0 cm
3 , since neutrinos detected at SK are not limited to passing throughonly the crust.
With the interaction potential distribution known, the oscillation survival probability
can be calculated for atmospheric muon neutrinos. A 2-dimensional contour plot of themodel’s predicted muon neutrino and anti-neutrino survival probabilities for a wide range
of energies and zenith angles, including those of the atmospheric neutrinos arriving at SK,
is shown in Fig. 8.5. The model input parameters for the neutrino masses and interaction
coupling are assumed to be the values from Eqn. 8.7.
In this figure it can be seen that the model predicts rapid oscillations as a function ofneutrino baseline and energy, especially in the anti-neutrino case. Contours showing thedifference between this model and the standard 2-flavor model (as shown in Fig. 7.3) are
shown in Fig. 8.6. The regions of large difference predict noticeable disappearance effects
in the SK atmospheric data.

154

gFigure 8.4: Effective mixing angle in matter as a function of energy. A density of ρ = 4.0 cm
3for the Earth is assumed.

Figure 8.5: Theoretical survival probabilities for muon neutrinos and anti-neutrinos in theapparent CPT-violation model.

155

Figure 8.6: Difference of oscillation survival probabilities for muon neutrinos and antineutrinos between the apparent CPT-violation model and the standard 2-flavor model.

8.5

Analysis Method

The model was fit to the SK-I+II+III+IV datasets using the same pull method [75] usedin the 2-flavor analysis from the previous chapter. The details of the datasets are shown in
Table 7.4.1. For the official monte carlo samples, 500 simulated livetime years of FC and PC
and 200 years of Upmu events were used for each generation of SK. In this analysis threeparameters were varied on a grid to find the best fit to the SK data. The active neutrinomass (m) and sterile neutrino mass (M), defined just as in the original model, were the firsttwo parameters fitted. The third and final fit parameter, g, describes the overall (weak plus
exotic) interaction coupling and is equal to√ g̃ 2
g = GF + 2 2 2 .
mV

(8.8)

This coupling constant reduces to the Fermi constant in the absence of the interaction.
For a given matter density g is directly proportional to the matter potential V. However, unlike the accelerator beam neutrinos at MINOS, the neutrinos detected at SK travel througha non-uniform matter density, so the potential V varies along the path even for the sameneutrino. Therefore, g was used in the fitting instead of V. The ranges of each of theparameters was

156

0.0001 eV < m < 0.1 eV(8.9)

0.1 eV < M < 100.0 eV(8.10)

1.167 x10−5 GeV−2 < g < 1.0 GeV−2 ,(8.11)

with a total of 96,400 points tested. These ranges completely covered the values tested
in the original model’s paper. At each analysis point, each MC event was assigned an
oscillation weight using Eqn. 8.6 with the model parameters taken at that point. Afterassigning oscillation weights, both the data and oscillated MC events were placed into oneof the 480 analysis bins shown in Fig. 7.1, based on their event type, energy and zenith-angle
values. For each analysis point a total chi-squared was calculated by summing the individual
chi-squared contributions of each bin using Eqn. 7.7. Each of the 150 systematic error pullterms was allowed to fluctuate to minimize the overall chi-squared for each analysis point
just as in the 2-flavor analysis.
8.6

Analysis Results

The 2-dimensional contours showing the 68%, 90% and 99% allowed regions of the modelfrom the fit to the SK-I+II+III+IV data are shown in Figs. 8.7, 8.8 and 8.9.
The SK dataset rules out a large portion of the original model’s allowed regions. The 1dimensional χ2 distributions for each of the three parameters is shown in Figs. 8.10,8.11,8.12.
The minimum chi-squared was found at the combination of parameters.
m = 0.046 eV, M = 10 eV, g = 1.44 × 10−5 GeV−2(8.12)

with χ2min = 554/478 d.o.f.
The best-fit value for the lightest non-zero neutrino mass in the model agrees well withthe square root of the atmospheric mass-splitting value predicted by the standard 2-flavorpmodel (m ∼ ∆m223 ) if the lighter mass states are assumed to be zero. The best-fit mass
for the heavier sterile neutrino state allows only a modest amount of sterile-active mixing
because of the result M  m (this limit suppresses active-sterile mixing). This is consistent
with a previous analysis at SK constraining sterile-active mixing[60].

157

Figure 8.7: The 68% (blue), 90% (red) and 99% (green) allowed regions in the g versus mspace for the fit to the SK-I+II+III+IV data, with M minimized away.

158

Figure 8.8: The 68% (blue), 90% (red) and 99% (green) allowed regions in the g versus M
space for the fit to the SK-I+II+III+IV data, with m minimized away.

159

Figure 8.9: The 68% (blue), 90% (red) and 99% (green) allowed regions in the M versus mspace for the fit to the SK-I+II+III+IV data, with g minimized away.

160

Figure 8.10: The 1-d χ2 distribution and 68% (blue), 90% (red) and 99% (green) C.L. forthe active neutrino mass m after fitting to the SK-I+II+III+IV data.

In conclusion the only allowed region at 90% C.L. corresponds to small sterile-active
mixing, active neutrino disappearance consistent with 2-flavor atmospheric mixing, and avery small contribution from the new gauge boson in the interaction coupling constant. A
lower bound is placed on the sterile neutrino mass in this model of M > 1.4 eV at 90% C.L.
An upper bound is placed on the component of the new interaction in the Fermi couplingof g < 8.0 × 10−4 GeV−2 at 90% C.L.
In addition to the largest allowed region, smaller disconnected allowed regions exist.
These regions correspond to a smaller amount of traditional muon disappearance combinedwith increased sterile-active mixing. These regions are allowed only at 99% and this isbecause the particular combination of parameters there almost reproduces 2-flavor disappearance for neutrinos, but diverges significantly from the 2-flavor model for anti-neutrinos,
thus ruling them out at 90% C.L.
Studies were done on this region to see if the SK data were capable of ruling them out

161

Figure 8.11: The 1-d χ2 distribution and 68% (blue), 90% (red) and 99% (green) C.L. forthe sterile neutrino mass M after fitting to the SK-I+II+III+IV data.

with increased livetime. One such study was a sensitivity study using SK-I+II+III MC
simulated at livetimes of 7.9 years (the actual SK-I+II+III livetime), 15.8 years (double the
livetime) and 50 years. The result of the expected sensitivity to the exotic coupling parameter g is shown in Fig. 8.13. This study shows that for large values of g, the χ2 becomes
independent of g and flattens off, even with 50 years livetime. Thus, the flattening of the
distribution is an inherent limitation in SK’s ability to constrain such physics. However,the degree at which these higher values of g are ruled out by SK increases correspondinglywith livetime.

162

Figure 8.12: The 1-d χ2 distribution and 68% (blue), 90% (red) and 99% (green) C.L. forthe new coupling interaction g after fitting to the SK-I+II+III+IV data.

163

Figure 8.13: The 1-d χ2 distribution and 68% (blue), 90% (red) and 99% (green) C.L. forthe new coupling interaction g in a sensitivity study for livetimes 0f 7.9, 15.8 and 50 years
with the SK-I+II+III MC.

164

Chapter 9
CONCLUSIONThe Super-Kamiokande experiment has been taking atmospheric neutrino data sinceApril 1, 1996 and has operated in four distinct phases. We have conducted an analysis of
two neutrino oscillation models using the atmospheric neutrino data set collected in all four
periods. The updates in this analysis include the addition of 1,096.66 days of additional
livetime collected during the SK-IV period, improvements to the Monte Carlo simulations,
including flux, nuclear interaction and detector response models, an upgraded data acquisition system with improved dynamic charge and high-rate performance, an updated event
binning scheme and fully reprocessed Monte Carlo samples for fully contained, partiallycontained and upward-going muon reduction.
We fit the SK-I+II+III+IV data to the 2-flavor muon neutrino disappearance modelusing the improved MC sample and a zenith-angle binning technique. We measured values−3 eV2 and sin2 2θ
of ∆m223 = 2.30+0.17
23 > 0.95, 90% C.L. for the parameters in−0.22 ×10

the model. We find that these results have excellent agreement with recent worldwide
constraints on the same parameters.
We also fit the same dataset to a non-standard interaction model which includes bothsterile neutrino mixing and a new matter interaction. We find that the SK data rule out
most of the model’s parameter space and limit the model to the region where standard
2-flavor mixing is recovered.

165

BIBLIOGRAPHY[1] E. Kh. Akhmedov. Neutrino physics. arXiv:hep-ph/0001264v2, 2000.
[2] V. Barger, Danny Marfatia, and Adam Tregre. Neutrino mass limits from sdss, 2dfgrs,
and wmap. Physics Letters B, 595:55–59, 2004.
[3] G. D. Barr, T.K. Gaisser, P.Lipari, S. Robbins, and T. Stanev. Three-dimensional
calculation of atmospheric neutrinos. Physical Review D, 70:023006, 2004.
[4] G. Battistoni, A. Ferrari, T. Montaruli, and P.R. Sala. The fluka atmospheric neutrino
flux calculation. Astroparticle Physics, 19:269–290, 2004.
[5] Kirk Bays. PhD thesis, University of California-Irvine, 2012.
[6] A. Bodek and U.K. Yang. Modeling neutrino and electron scattering inelastic crosssections. arXiv:hep-ex/0308007v1.
[7] R. Brun and F. Carminati. Geant: Detector description and simulation tool. CERN
Program Library Long Writeup, W5013, 1993.
[8] A. Chavarria. http://www.phy.duke.edu/ schol/superk/alvaro thesis.pdf. Senior Thesis, 2007.
[9] ALEPH Collaboration. Determination of the number of light neutrino species. Physics
Letters B, 231:519–529, 1989.
[10] Borexino Collaboration. Observation of geoneutrinos. Physics Letters B, 687:299–304,
2010.
[11] Super-Kamiokande Collaboration, 2002.
committee/index-e.html.

http://www-sk.icrr.u-tokyo.ac.jp/cause-[12] Super-Kamiokande Collaboration. The super-kamiokande detector. Nuclear Instruments and Methods in Physics Research, 501A:418–462, 2003.
[13] T. Jannakos (KARMEN Collaboration). The karmen limit on ν̄µ → ν̄e oscillationand its implication for the lsnd result. Nuclear Physics B (Proceedings Supplements),
85:84–90, 2000.

166[14] Zoa Conner. PhD thesis, University of Maryland, 1997.
[15] E. R. Davies. Machine Vision: Theory, Algorithms, Practicalities. Academic Press,
San Diego, 1997.
[16] Raymond Davis Jr., Don Harmer S., and Kenneth Hoffman C. Search for neutrinos
from the sun. Physical Review Letters, 20:1205–1209, 1968.
[17] Patrick de Perio. Neut pion fsi. AIP Conference Proceedings, 1405:223–228, 2011.
[18] Shantanu Desai. PhD thesis, Boston University, 2004.
[19] Fanny M. Dufour. PhD thesis, Boston University, 2009.
[20] Adam M. Dziewonski and Don L. Anderson. Physics of the Earth and Planetary
Interiors, 25(297), 1981.
[21] Netta Engelhardt, Ann E. Nelson, and Jonathan R. Walsh. Apparent cpt violation in
neutrino oscillation experiments. Physical Review D, D81(113001), 2010.
[22] A. A. Aguilar-Arevalo et al. A combined νµ → νe & ν̄µ → ν̄e oscillation analysis ofminiboone excesses. arXiv:hep-ex/1207.4809v1, 2012.
[23] A. Kayis-Topaksu et al. Charged-particle multiplicities in charged-current neutrinoand anti-neutrino-nucelus interactions. The European Physical Journal C, 51(4):775–
785, 2007.
[24] A.A. Aguilar-Arevalo et al. Measurement of muon neutrino quasi-elastic scatteringon carbon. Physical Review Letters, 100:032301, 2008.
[25] Ch. Kraus et al. Most recent results of the mainz neutrino mass experiment. Nuclear
Physics B (Proceedings Supplements), 118:482, 2003.
[26] D. Rowntree et al. π + absorption on n and ar. Physical Review C, 60:054610, 1999.
[27] F.P. An et al. Observation of electron-antineutrino disappearance at daya bay.
arXiv:1203.1669v2 [hep-ex].
[28] G. Ambrosini et al. Production ratios from 450 gev/c protons on beryllium. Physics
Letters B, 420:225, 1998.
[29] G. Danby et al. Observation of high-energy neutrino reactions and the existence of
two kinds of neutrinos. Physical Review Letters, 9:36–44, 1962.

167[30] Glenn Rowe et al. Energy-dependent phase shift analysis of pion-nulceon scattering
below 400 mev. Physical Review C, 18:584–589, 1978.
[31] H. Nishino et al. High-speed charge-to-time coverter asic for the super-kamiokande
detector. Nuclear Instruments and Methods in Physics Research, 610A:710–717, 2009.
[32] H. Nishino et al. Search for proton decay via p → e+ π 0 and p → µ+ π 0 in a large
water cherenkov detector. Physical Review Letters, 102, 2009.
[33] H. Nishino et al. The new front-end electronics for the super-kamiokande experiment.
Nuclear Science Symposium Conference Record, NSS ’07, IEEE, pages 127–132, October 26 2007–November 3 2007.
[34] J. Beringer et al. 2012 review of particle physics. Physical Review D, 86:010001, 2012.
[35] J. Nieves et al. Inclusive quasielastic charged-current neutrino-nucelus interactions.
Physical Review C, 70:055503, 2004.
[36] K. Abe et al. Indication of electron neutrino appearance from an accelerator-produced
off-axis muon neutrino beam. Physical Review Letters, 107:041801, 2011.
[37] K. Abe et al. Search for difference in oscillation parameters for atmospheric neutrinos
and anti-neutrinos at super-kamiokande. Physical Review Letters, 107:241801, 2011.
[38] K. Abe et al. Solar neutrino results in super-kamiokande-iii. Physical Review D,
83:052010, 2011.
[39] K. Abe et al. The t2k experiment. Nuclear Instruments and Methods in Physics
Research, 659A:106–135, 2011.
[40] K. Niita et al. Phits - a particle and heavy ion transport code system. Radiation
Measurements, 41:1080, 2006.
[41] K. S. Hirata et al. Observation in the kamiokande-ii detector of the neutrino burst
from supernova sn1987a. Physical Review D, 38(2):448–458, 1988.
[42] K. S. Hirata et al. Observation of 8 b solar neutrinos in the kamiokande-ii detector.
Physical Review Letters, 63(1):16–19, 1989.
[43] Kim Soo-Bong et al. Observation of reactor electron antineutrino disappearance inthe reno experiment. arXiv:1204.0626v2 [hep-ex].
[44] Kirk Bays et al. Supernova relic neutrino search at super-kamiokande. Physical Review
D, 85, 2012.

168[45] M. Ahn et al. Measurement of neutrino oscillation by the k2k experiment. PhysicalReview D, 74:072003, 2006.
[46] M. E. C. Swanson et al. Search for diffuse astrophysical neutrino flux using ultrahigh-energy upward-going muons in super-kamiokande-i. The Astrophysical Journal,652(1), 2006.
[47] M. Honda et al. New calculation of the atmospheric neutrino flux in a threedimensional scheme. Physical Review D, 70:043008, 2004.
[48] M. Honda et al. Calculation of atmospheric neutrino flux using the interaction model
calibrated with atmospheric muon data. Physical Review D, 75:043006, 2007.
[49] M. Honda et al. Improvement of low energy atmospheric neutrino flux calculation
using the jam nuclear interaction model. Physical Review D, 83:123001, 2011.
[50] M. Nakahata et al. Atmospheric neutrino background and pion nuclear effect for
kamioka nucleon decay experiment. Journal of the Physical Society of Japan, 55:3786,
1986.
[51] M. Nakahata et al. Calibration of super-kamiokande using an electron linac. Nuclear
Instruments and Methods in Physics Research, A421:113–129, 1999.
[52] M.H. Ahn et al. Measurement of neutrino oscillation by the k2k experiment. PhysicalReview D, 74:072003, 2006.
[53] P. Adamson et al. Measurement of neutrino oscillations with the minos detectors in
the numi beam. Phys. Rev. Lett., 101:131802, Sep 2008.
[54] P. Adamson et al. Measurement of the neutrino mass splitting and flavor mixing by
minos. Physical Review Letters, 106:181801, 2011.
[55] P. Adamson et al. An improved measurement of muon antineutrino disappearance inminos. Physical Review Letter, 108:191801, 2012.
[56] R. Becker-Szendy et al. Imb-3: a large water cherenkov detector for nucleon decay
and neutrino interactions. Nuclear Instrumentation & Methods in Physics Research,
A324:363, 1993.
[57] R. Wendell et al. Atmospheric neutrino oscillation analysis with subleading effects insuper-kamiokande i,ii and iii. Physical Review D, 81:092004, 2010.
[58] R.A. Arndt et al. Dispersion relation constrained partial wave analysis of πn elastic
and πn → ηn scattering data: The baryon spectrum. Physical Review C, 69:035213,
2004.

169

[59] S. Abe et al. Precision measurement of neutrino oscillation parameters. PhysicalReview Letters, 100:221803, 2008.
[60] S. Fukuda et al. Tau neutrinos favored over sterile neutrinos in atmospheric muon
neutrino oscillations. Physical Review Letters, 85:3999–4003, 2000.
[61] S.N. Ahmed et al. Measurement of the total active 8 b solar neutrino flux at the
sudbury neutrino observatory with enhanced neutral current sensitivity. PhysicalReview Letters, 92:181301, 2004.
[62] Tanaka T. et al. An indirect search for weakly interacting massive particles in the sunusing 3109.6 days of upward-going muons in super-kamiokande. The AstrophysicalJournal, 742, 2011.
[63] V.M. Lobashev et al. Direct search for neutrino mass and anomaly in the tritium
beta-spectrum: Status of “troitsk neutrino mass” experiment. Nuclear Physics B(Proceedings Supplements), 91:280–286, 2001.
[64] Y. Ashie et al. Evidence for an oscillatory signature in atmospheric neutrino oscillations. Physical Review Letters, 93:101801, 2004.
[65] Y. Ashie et al. A measurement of atmospheric neutrino oscillation parameters bysuper-kamiokande i. Physical Review D, 71:112005, 2005.
[66] Y. Fukuda et al. Evidence for oscillation of atmospheric neutrinos. Physical Review
Letters, 81:1562–1567, 1998.
[67] Y. Fukuda et al. Measurement of a small atmoshperic νµ /νe ratio. Physical Letters
B, 433:9–18, 1998.
[68] B. Ambruster et al. (KARMEN Collaboration). Upper limits for neutrino oscillations
muon-antineutrino to electron-antineutrino from muon decay at rest. Physical Review
D, 65:112001, 2002.
[69] A. Aguilar et al. (LSND). Evidence for neutrino oscillations from the observation ofν̄e appearance in a ν̄µ beam. Physical Review D, D64(112007), 2001.
[70] A. A. Aguilar-Arevalo et al. (MiniBooNE). Physical Review Letters, 98(231801), 2007.
[71] A. A. Aguilar-Arevalo et al. (MiniBooNE). A combined νµ to νe and ν̄µ to ν̄e oscillation
analyses of the mineboone excesses. arXiv:1207.4809v1 [hep-ex], 2012.
[72] P. Adamson et al. (MINOS Collaboration). First direct observation of muon antineutrino disappearance. Physical Review Letters, 107:021801, 2011.

170[73] J. J. Evans. PhD thesis, fERMILAB-THESIS-2009-14, 2009.
[74] Enrico Fermi. Zeitschrift fur Physik, 88:161, 1934.
[75] G. L. Fogli, E. Lisi, A. Marrone, D. Montanino, and A. Palazzo. Getting the most fromthe statistical analysis of solar neutrino oscillations. Physical Review D, D66(053010),
2002.
[76] G.L. Fogli, E. Lisi, A. Marrone, A. Palazzo, and A.M. Rotunno. Evidence oftheta(13)>0 from global neutrino data analysis. 2011. arXiv:1106.6028v2[hep-ph].
[77] Thomas Patzak for the DONUT Collaboration. First direct observation of the tau
neutrino. Europhysics news, 32(2):56–57, 2001.
[78] M. Gluck, E. Reya, and A. Vogt. Dynamical parton distributions revisited. The
European Physical Journal C, 5(3):461–470, 1998.
[79] S.N. Gninenko. Lsnd/miniboone excess events and heavy neutrino from µ and k
decays. Physical Review D, 83:093010, 2011.
[80] Maurice Goldhaber, L. Grodzins, and A.W. Sunyar. Helicity of neutrinos. PhysicalReview, 109:1015–1017, 1958.
[81] Jan Hamann, Steen Hannestad, Georg G. Raffelt, Irene Tamborra, and Yvone Y.Y.
Wong. Cosmology favoring extra radiation and sub-ev mass sterile neutrinos as an
option. Physical Review Letters, 105:181301, 2010.
[82] Y. Hayato. Neut. Nuclear Physics B (Proceedings Supplements), 112:171–176, 2002.
[83] Y. Hayato. A neutrino interaction simulation library neut. Acta Physica Polona B,
40:2477, 2009.
[84] A. E. Hedin. Msise model (1990). Planetary and Space Science, 40:556, 1992.
[85] E. Hernandez, J. Nieves, and M. Valverde. Weak pion production off the nucleon.
Physical Review D, 76:033005, 2007.
[86] Kazuhide Ichikawa, Masataka Fukugita, and Masahiro Kawasaki. Constraining neutrino masses by cmb experiments alone. Physical Review D, 71:043001, 2005.
[87] Chizue Ishihara. PhD thesis, University of Tokyo, 2010.
[88] Masaki Ishitsuka. PhD thesis, University of Tokyo, 2004.

171[89] J.D. Jackson. Classical Electrodynamics. John Wiley & Sons, Inc., 3rd edition, 1998.
[90] C.L. Cowan Jr., F. Reines, F.B. Harrison, H.W. Kruse, and A. D. McGuire. Detectionof the free neutrino: A confirmation. Science, 124(3212):103–104, 1956.
[91] S. I. Klimushin, E. V. Bugaev, and Igor A. Sokalski. Precise parameterizations ofmuon energy losses in water. 2001. arXiv:hep-ph/0106010v1.
[92] Maggie Ka-Pik Lee. PhD thesis, University of Tokyo, 2012.
[93] Ziro Maki, Masami Nakagawa, and Shoichi Sakata. Remarks on the unified model of
elementary particles. Progress of Theoretical Physics, 28(5):870–880, 1962.
[94] Thomas McLachlan. PhD thesis, University of Tokyo, 2012.
[95] S. P. Mikheyev and A. Y. Smirnov. Resonant enhancement of oscillations in matter.
Soviet Journal Nuclear Physics, 42:913–917, 1985.
[96] Ann E. Nelson. Effects of cp violation from neutral heavy fermions on neutrino
oscillations, and the lsnd/miniboone anomalies. arXiv:1010.3970v1, 2010.
[97] Ann E. Nelson and Jonathan R. Walsh. Physical Review D, D77(033001), 2008.
[98] NOAA, NASA, and US Air Force. U.S. Standard Atmosphere. U.S. Government
Printing Office, Washington D.C., 1976.
[99] Wolfgang Pauli. Open letter to radioactive ladies and gentlemen (1930). Physics
Today, 31:27, 1978.
[100] Bruno Pontecorvo. Neutrino experiments and the problem of conservation of leptonic
charge. Soviet Physics JETP, 26:984, 1968.
[101] Dieter Rein and Lalit M. Sehgal. Neutrino-excitation of baryon resonances and single
pion production. Annals of Physics, 133:79–151, 1979.
[102] B.G. Ritchie. Parameterization of total and differential cross sections for πd → pp
below 1 gev. Physical Review C, 44:533–536, 1991.
[103] S. Roesler, R. Engel, and J. Ranft. The monte carlo event generator dpmjet-iii.
arXiv:hep-ph/0012252.
[104] L.L. Salcedo and E. Oset et al. Computer simulation of inclusive pion nuclear reactions. Nuclear Physics A, 484:557–592, 1988.

172[105] T. Sjoestrand. High-energy-physics event generation with pythia 5.7 and jetset 7.4.
Computer Physics Communications, 82:74–89, 1994.
[106] C.H. Llewellyn Smith. Neutrino reactions at accelerator energies. Physics Reports,
3:261–379, 1972.
[107] R. A. Smith and E. J. Moniz. Neutrino reactions on nuclear targets. Nuclear Physics
B, 43:605–622, 1972.
[108] Yumiko Takenaga. PhD thesis, University of Tokyo, 2008.
[109] Takayuki Tanaka. PhD thesis, University of Nagoya, 2011.
[110] I. Thormahlen, K. Scheffler, and J. Straub. Water and Steam: Proceedings of the 9th
ICPS. Pergamon, Oxford, 1980.
[111] Eric Thrane. PhD thesis, University of Washington, 2008.
[112] Roger A. Wendell. PhD thesis, University of North Carolina, 2008.
[113] L. Wolfenstein. Neutrino oscillations in matter. Physical Review D, 17:2369–2374,
1978.
[114] C. Zeitnitz and T.A. Gabriel. The geant-calor interface and benchmark test calculations of zeus calorimeters. Nuclear Instruments and Methods in Physics Research A,
349:106–111, 1994.

173

VITA

Michael Dziomba was born in Buffalo, NY in 1984. He earned a BS in Physics, graduating Summa Cum Laude, from the State University of New York at Buffalo in 2006.
He earned his MS (2007) and PhD (2012) in Physics from the University of Washington.
His thesis focus was atmospheric neutrino oscillation studies with the Super-Kamiokande
detector, with advisor R.J. Wilkes.

.