Muonic event reconstuction in KM3NeT-ORCA

(1)

MSc Physics and Astronomy

Gravitation and Astroparticle Physics Amsterdam

Master Thesis

Muonic event reconstruction in KM3NeT-ORCA

written by Max Merlijn Briel

UvA ID: 10606513 September 2019

60 EC 2018-2019

Supervisor/Examiner: Second Examiner

(2)

Abstract

The KM3NeT collaboration aims to solve the unknown neutrino mass ordering by mea-suring the direction and energy flux of atmospheric neutrinos. At a depth of over 2 kilometres in the Mediterranean Sea, it uses photomultiplier tubes (PMTs) to record Cherenkov radiation from the high-energy neutrino reaction products.

A charged current interaction with a muon neutrino leaves a path, named a track, of photon hits through the detector. Using the hits and the properties of the PMTs, the direction of the track is approximated in a process called reconstruction and comprises two successive algorithms. The first generates a set of starting values for the second main algorithm, which uses a maximum likelihood method to find the best fitting track to the data. This research determines that the initial directional accuracy the main reconstruction has to be 10° or lower for the algorithm to achieve its full potential. Since the first algorithm only reaches this limit in 10% of the events, we implement new quality functions and likelihoods in the context of re-ranking and minimisation. They contain different parametrisations of the hit information, which have been developed for this work. We show that the ranking of the main algorithm is optimal in selecting the best track set, but a new likelihood improves their internal ranking. Minimisation with the new likelihoods increases the percentage below the limit up to 25%.

A more complete likelihood is introduced with hit and no-hit information. Its like-lihood space shows potential, and minimisation with the truth parameters leads to a significant improvement in the mean directional error from 1.92° to 0.81° compared to the current main reconstruction. 60% of the reconstructed tracks are now at a sub-degree accuracy. However, its likelihood contains many local minima and the pre-reconstruction is not optimised for this comprehensive likelihood putting it at a disadvantage. With better optimisation of the full reconstruction chain, the extended algorithm could provide better constraints on the direction and the neutrino mass ordering.

(3)

Introduction

The Standard Model is the best description of the subatomic world. It contains the 12 known particles, their antiparticles, three fundamental forces with their 4 gauge bosons, and generates mass using a single scalar boson [1–3]. The theory results from many years of theoretical and experimental research and has predicted the existence of new particles, such as the W±/Z bosons [4–7], gluon [8–11], and the Higgs [12, 13].

The model, however, is incomplete and contains only three of the four fundamental forces. Gravity has been difficult to merge with the microscopic scale of particles. At this moment, its extreme weakness allows it to be safely ignored in experimental predictions [1]. More pressing issues with the Standard Model are at the small scales, such as neutrino masses.

Early neutrino experiments, sensitive to only solar electron neutrinos, measured fewer interactions than predicted [14–20]. The most natural explanation was an in-flight flavour change of the neutrino [21]. For example, when an electron neutrino (νe) is created with

energy Eν and travels a distance L, it has a probability to be detected as a muon neutrino

(νµ). Through the transition νe → νµ, muon neutrinos appear in the detector, while

electron neutrinos disappear, resulting in less than expected νe. Experiments measuring

all three neutrino types were consistent with predictions and confirmed the neutrino mixing [22–30]. As a result, neutrinos are massive, which is in contrast with the massless Standard Model neutrinos. The mass states of the three neutrinos have to be unique and the flavour states are a linear combination of these mass eigenstates [31, 32].

The current neutrino oscillation framework comprises three flavours (νe, νµ, ντ) and

three mass eigenstates (ν1, ν2, ν3), although it is still an open question whether more

non-interacting neutrinos exist [33]. Eν and L, and a set of constants determine the mixing

probability. Most of them are contained within the Pontecorvo-Maki-Nakagawa-Sakata matrix (PMNS matrix) [31, 32], which parametrises the mixing using three mixing angles and one CP violating phase, similar to the quark sector [34, 35]. The other constants are the squared mass differences between the mass eigenstates. With three eigenstates, two mass differences are independent. So far, experiments are unable to measure the sign of one of the mass differences, which results in a degeneracy in the mass ordering. This is known as the Mass Hierarchy problem, also named the Mass Ordering problem.

A promising method for determining the neutrino mass ordering is to exploit the effect of matter on the oscillation probability of neutrinos and antineutrinos [36], which alters the oscillation probability depending on energy, distance travelled, and the electron

(6)

density of the matter traversed. The ORCA (Oscillation Research with Cosmics in the Abyss) telescope of the KM3NeT (Cubic Kilometer Neutrino telescope) infrastructures plans to use this method with neutrinos and antineutrinos produced by cosmic rays. Cosmic rays are highly energetic charged particles that hit the atmosphere and create a shower of particles including νe/µ and ¯νe/µ. These atmospheric neutrinos oscillate

through the Earth and create high-energy interaction products, such as muons, taus, and electrons.

ORCA’s detector structure is optimised for the 3 to 20 GeV range, where matter effects have a strong presence. It is, however, unable to distinguish between neutrino flavour of the interaction directly. Instead, the detector can only separate two interaction signatures: showers and tracks. Neutrinos interacting with hadrons or electrons through the neutral current cause the former, while charged current interactions from electron neutrinos result in a similar signal. Tracks, on the other hand, are caused by muons that travel far through the detector. They originate from charged current muon neutrino interactions, but the decay of a tau can also generate muons. This thesis focuses on muonic events because of the excellent directional information provided by the long path of the muon through the detector. The muon’s direction links back to the amount of matter traversed by the neutrino. Together with the neutrino’s energy, they are essential in determining the mixing angles, mass differences, and mass ordering.

To balance the extremely small neutrino interaction rate, the detector will encompass 5 megatons of sea water. At a depth of over 2 kilometres in the Mediterranean Sea, the ORCA detector uses photomultiplier tubes (PMTs) housed in pressure-resistant glass spheres called Digital Optical Modules (DOM) to record Cherenkov radiation from the high-energy reaction products. When a photon is registered on a PMT, it records a hit containing the time and duration of the measured signal. From a collection of such hits in combination with the direction and orientation of the PMTs (an event), we approximate the original path of the muon through the detector. This process is called reconstruction and consists of two successive algorithms: JPrefit and JGandalf. The former generates a set of 36 starting values for the latter. The JGandalf algorithm implements a maximum likelihood method to find the best fitting track to the data using these given initial guesses. However, it is sensitive to the directional accuracy of the tracks provided by JPrefit. Besides quantifying the sensitivity, this research implements new likelihood functions to improve the directional accuracy of the JPrefit tracks. The probability a PMT was hit during an event can provide additional information for a better directional reconstruction. In this research, different parametrisations of these probabilities are created from the analytical description of the emission, propagation, and detection of light from high energy charged particles. The new likelihoods with parametrisations are used to re-rank or minimise the JPrefit tracks to give better starting values for JGandalf.

The absence of light on a PMT also provides information on the muon’s path through the detector. If a PMT is not hit during an event, it is known as a no-hit. Multiple no-hits can indicate a region where the track cannot have passed through due to the lack of hits. Using the full probability functions, a new more complete likelihood with hit and no-hit information is implemented: JMerlin.

(7)

This thesis starts with a description of the neutrino oscillation framework and the KM3NeT detector in Chapter 2 and 3, respectively. Chapter 5 evaluates JGandalf and the current tracks provided by the JPrefit algorithm. The tracks are improved upon in Chapter 6 using the likelihood with parametrisations. JMerlin is introduced and tested in Chapter 7. Chapter 8 ends the thesis with a discussion and conclusion about the new likelihoods and their implementation.

(8)

Chapter 2

Neutrino Oscillations

2.1 The History of the Neutrino

2.1.1 β decay

The history of the neutrino starts at the end of the 19th century with the discovery of β− decay by Ernest Rutherford [37]. The nucleus changes charge from Z to Z + 1 alongside the emission of an electron depicted in 2.1. If we assume the interaction to be a two-body decay, as depicted, the electron’s energy spectrum is constant. Experiments have instead shown that the spectrum is continuous [38].

A(N, Z) → A(N − 1, Z + 1) + e− (2.1)

One of the many proposed solutions came from Wolfgang Pauli. He suggested a third light particle with a mass of 0.01mp, a spin 1₂, and named it ”neutron”. After the

discovery of the neutral nucleon, the neutron [39], the idea was forgotten until Fermi implemented the light particle in his theory for β-decay. He gave it the name we know today: neutrino; ”the light neutral one” [40].

2.1.2 Neutrino Discoveries: νe νµ ντ

When Pauli included the particle in his theory, he worried it would be extremely diffi-cult to detect [40]. Neutrino cross section calculations strengthened this idea, and the particle was expected to remain undetected [41]. With the rise of nuclear fission, intense neutrino and anti-neutrino sources became available, which vastly increased the chance of measuring a neutrino interaction. In the 40s and 50s, many experiments used nuclear reactors as their neutrino source [42–44].

In 1956, Reines and Cowan successfully measured a neutrino interaction using inverse beta decay (Eq. 2.2) for the first time [45]. An electron antineutrino from the nuclear reactor interacts with a proton and emits a positron and neutron. The positron immedi-ately annihilates with electrons in the medium to photons, while a nucleus captures the neutron and emits delayed light allowing for background reduction [46].

¯

(9)

W− νl l f f0 (a) W+ ¯ νl ¯l f f0 (b) Z0 νl(¯νl) νl(¯νl) f (f0) f (f0) (c)

Figure 2.1: Charged Current neutrino (a), Charged Current antineutrino (b), and Neutral Current (c) interactions, where l (¯l) is a (anti)lepton and f ( ¯f ) is a (anti)fermion

After the discovery of the electron neutrino, it only took six years to find the muon neutrino [47]. Muons were observed to originate from π± decays. Thus, a νµ had to be

involved for a conserved interaction. The discovery of the tau lepton suggested a third neutrino flavour [48], which CERN experiments, such as DELPHI, indirectly confirmed through the Z0 _{decay width [3]. A direct measurement of a tau neutrino interaction took}

till 2000, when the DONUT collaboration imaged the interaction [49].

2.2 Neutrino Interactions

Neutrino interactions are complex and depend heavily on the target particles. As the literature thoroughly discusses many aspects of the interactions [46, 50, 51], this thesis will summarise those relevant for the KM3NeT-ORCA detector.

The weak force mediates the interactions between the neutrino and ordinary matter [3], in which only left-handed flavour eigenstates of ν can partake. In the neutral current (NC) channel, the neutrino remains a neutrino, as seen in Figure 2.1c. Salam, Glashow and Weinberg postulated this channel [52–54] and the Gargamelle neutrino collaboration at CERN confirmed it in 1973 [55]. During the interaction, only some momentum of the neutrino is transferred to the other particle through a neutral Z boson.

ν detection experiments used the charged current (CC) channel to confirm the neu-trino’s existence. The neutrino transforms into a charged lepton of the same flavour (Figure 2.1a). It determines the flavour of the neutrino, while a W± boson mediates the charge. In ORCA’s energy region from 3 to 20 GeV, the CC interaction cross section is 3 times larger than the NC channel [56]. For each interaction, an opposite one exists for the antimatter counterpart with one such example being depicted in Figure 2.1b.

In the energy range of the KM3NeT-ORCA detector different types of interactions can take place between the neutrino and nucleons, such as (quasi-)elastic scattering, resonance production, and deep inelastic scattering [51]. Their combined cross section for neutrinos is twice as large as those for antineutrinos [56], as shown in Figure 2.2.

(10)

(GeV) E -1 10 1 10 ₁₀2 / GeV) 2 cm -38 (10 cross section / E 0 0.2 0.4 0.6 0.8 1 1.2 1.4 (GeV) E -1 10 1 10 ₁₀2 / GeV) 2 cm -38 (10 cross section / E 0 0.2 0.4 0.6 0.8 1 1.2 1.4 TOTAL QE DIS RES (GeV) E -1 10 1 10 ₁₀2 / GeV) 2 cm -38 (10 cross section / E 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 (GeV) E -1 10 1 10 ₁₀2 / GeV) 2 cm -38 (10 cross section / E 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 TOTAL QE _RES DIS

Figure 2.2: Nucleon interaction cross sections of the neutrino (left) and antineutrino (right) on an isoscalar target. The different scattering processes are marked with RES for RESonance production, QE for Quasi-Elastic scattering, and Deep Inelastic Scattering as DIS. Figures taken from [51].

The KM3NeT’s target particle, however, is not a single proton or neutron. Their bound state in the water molecules increases the interaction rate with 1% through coherent scattering. Other effects include the motion of the nucleon in the nucleus [57,58] and final state interactions [59]. In the high energy regime, the relative contribution of the motion is negligible, and the KM3NeT infrastructure cannot distinguish individual particles, except the muon.

2.3 Neutrino Oscillations

As seen in Section 2.2, a neutrino takes part in an interaction with a specific flavour: electron (νe), muon (νµ), or tau (ντ). These are eigenstates of the weak interaction, so

the ν is created in one of these three flavours. A neutrino, however, does not have to stay in its original flavour state. The Homestake experiment measured the νe flux from the

sun using a chlorine bath to capture neutrinos using inverse beta decay, which creates measurable radioactive argon, see Equation 2.3. Only around 1/3 of the expected electron neutrino flux was measured [14, 17] and this deficit became known as the Solar Neutrino Problem.

ν +37Cl → e−+37Ar (2.3)

The most natural solution was already proposed many years earlier as a reaction to a previous experiment from Davis using the same reaction 2.3 to measure reactor antineutrinos, which should be impossible. However, rumours spread that he had found such events and reached Bruno Pontecorvo, who came up with neutrino-antineutrino oscillations as an explanation [60, 61]. While the rumours turned out to be wrong, the idea of oscillations stayed. After the discovery of the second generation of neutrinos,

(11)

Maki, Nakagawa, and Sakata took the idea and proposed that the ”true neutrinos” are a linear combination of different neutrino flavours [32]. Pontecorvo extended this idea to include flavour oscillations between the different flavours, such as νe νµ [31]. These

oscillations were the most natural solution to the Solar Neutrino Problem [21] and have major consequences. The neutrinos need to have distinct masses, and the flavour states are a linear combination of these mass eigenstates.

|να(t = 0)i =

X

j

Uαj|νji with (α = e, µ, τ & j = 1, 2, 3), (2.4)

in which Uαj, the Pontecorvo, Maki, Nakagawa, and Sakata (PMNS matrix), describes

the mixing of the mass states in the flavour eigenstates. 3 mixing angles and 6 phases parametrise this unitary matrix [3]. The number of phases that are real depends on the type of neutrino, Dirac or Majorana [34, 62]. In neutrino oscillation experiments, the difference is unimportant and the simplest form, the Dirac case, can be considered [63]. Five of the six phases are absorbed and one charge-parity violating phase (δCP) remains.

Together, with the three mixing angles, they parametrise the PMNS matrix, similar to the CKM matrix in the quark sector [3]. This is shown in Equation 2.5, where cij and

sij are the cosine and sinus of the mixing angles: θ12, θ13, and θ23.

U =   c12c13 s12c13 s13e−iδ −s12c23− c12s13s23eiδ c12c23− s12s13s23eiδ c13s23 s12s23− c12s13c23eiδ −c12s23− s12s13c23eiδ c13c23   =   1 0 0 0 c23 s23 0 −s₂₃ c23     c13 0 s13e−iδCP 0 1 0 −s₁₃eiδCP ₀ _c 13     c12 s12 0 −s12 c12 0 0 0 1   (2.5)

Neutrinos are always created with a specific flavour depending on the associated charged lepton, but each mass eigenstate propagates with its own phase factor: e−iEjt_. The initial flavour state becomes a mixture of the other flavour states (Eq. 2.6).

The probability of measuring a specific state is given by the probability amplitude of the interaction. For measuring a να as flavour β with α 6= β the probability is:

P (να→ νβ) = |hνβ|να(t)i|2 = X j UβjUαj? e−iEjt 2 (2.7) In the ultra-relativistic limit (E ≈ p), where the neutrinos of many experiments reside, Ej can be rewritten and Taylor expanded. In natural units, the time and distance

(12)

Ej = q p2_{+ m}2 j = p 1 + m2_j p2 ! ≈ p 1 + m 2 j 2p2 ! = p + m 2 j 2p = E + m2_j 2E (2.8)

We can rewrite the oscillation probability to a real and imaginary part. The real part is charge-parity invariant, while the imaginary part is charge-parity violating and switches sign for the antineutrino.

P (να → νβ) = δαβ − 4 X i<j ReUαiUβj? Uαj? Uβj sin2 ∆m2_ji 4E L ! + 2X i<j

ImUαiUβi?Uαj? Uβj sin

∆m2_ji 2E L

! (2.9)

Equation 2.9 shows that the neutrino oscillation probability depends on elements of the PMNS matrix, the squared mass difference between mass eigenstates (∆m2

ji), the

neutrino’s energy (E), and the distance travelled (L). To show the dependence on the PMNS matrix elements explicitly, we consider the two-flavour ν case. In this limit, the oscillation probability depends only on a single mass difference and a single mixing angle, as shown in Equations 2.10 and 2.11.

P (να→ να) = 1 − sin2(2θ) sin2 ∆m2_L 4E , (2.10) P (να→ νβ) = sin2(2θ) sin2 ∆m2_L 4E (2.11) In the full neutrino flavour model, three mass differences are distinguished with two being independent. The convention defines these as ∆m2₂₁ and ∆m2₃₁. Together they can describe ∆m2₃₂ = ∆m2₃₁− ∆m2

21, but a degeneracy remains. Since the sign of ∆m231 is

unknown, it is impossible to determine the ordering of the neutrino mass eigenstates. ∆m2₂₁ is defined as positive with m1 < m2 from experimental data. m3 can either be

heavier than m2 or lighter than m1. This is the neutrino mass ordering or mass hierarchy

problem. The first case, m1 < m2 < m3, is referred to as the normal ordering with

∆m2₃₁> 0. In the inverted ordering case, m3 < m1 < m2, ∆m231 is negative.

Determining the exact order has become a focus of neutrino oscillation experiments, because it rules out or strengthens around half the proposed particle physics models [64]. It will help constraint result in leptonic CP violation, the absolute neutrino mass, and the Majorana-Dirac nature of the neutrino [65]

(13)

2.4 Matter Oscillations

W−

e− νe

νe e−

Figure 2.3: Electron and neutrino interaction called coherent forward scatter-ing resultscatter-ing in the MSW effect. 2 4 6 8 10 Energy [GeV] 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 P (νµ → νe ) Vac IO NO Vacuum

Matter (Normal Ordering) Matter (Inverted Ordering)

Figure 2.4: Electron neutrino appearance probability for a 1300 km baseline. Oscillation parameters are taken from [66] and calculated using [67]. Vacuum oscillations for NO (dashed blue), matter oscillation for NO (blue) and IO (red) are shown.

The main idea to determine the neutrino mass ordering is to use the influence of the matter on the neutrinos. When travelling through a medium, the neutrino feels a potential caused by coherent forward scattering. The neutrino interacts with surrounding nucleons and electrons through the NC and CC interactions without disappearing [68]. Neutral current interaction influences all neutrino types equally, but CC interactions only take place between the electron (anti)neutrino and electron. This changes the effective mass and the oscillation probability of the neutrinos. In the two-flavour limit, we can rewrite this probability like their vacuum counterparts:

Pm(να → να) = 1 − sin2(2θm) sin2 ∆m2 mL 4E (2.12) Pm(να → νβ) = sin2(2θm) sin2 ∆m2 mL 4E (2.13)

The matter alters the mixing angle θm and the mass difference ∆m2m. But they still

relate back to their vacuum values through Equations 2.14 and 2.15, where GF is Fermi’s

constant, ne the electron density in the matter, and the sign depends on the involvement

of a neutrino (+) or antineutrino (−).

∆m2_m=p(∆m2_{cos(2θ) − A)}2_{+ (∆m}2_{sin 2θ)}2_{, A = ±2}√_2EG

(14)

sin2θm= ∆m2sin 2θ ∆m2 m = 1 2 1 +A − ∆m 2_{cos 2θ} ∆m2 m (2.15) From Equation 2.15, the amount of mixing between the two flavours can be maximised when ∆m2cos 2θ = 2√2GFneE. This resonance enhancement allows for mixing angles

many times larger than in vacuum [68]. Whether the condition is satisfied depends on the neutrino mass ordering. With normal ordering, matter boosts the neutrino mixing, while with inverted ordering the antineutrino is boosted, as depicted in Figure 2.4.

This is known as the Mikheyev-Smirnov-Wolfenstein effect, which fully explains the solar neutrino problem [69, 70]. Moreover, with a varying electron density, the mixing parameters themselves change, introducing an extra phase to the oscillations. The litera-ture provides more details about this effect and other oscillation altering processes, such as parametric enhancement [71–73].

2.4.1 Usage within KM3NeT-ORCA

To achieve the full effectiveness of the resonance, a large mass is required. KM3NeT-ORCA employs the Earth for this purpose and measures atmospheric neutrino, which are products from cosmic charged particles interacting with the atmosphere. These so called cosmic rays are very energetic and create showers of particles including pions, which decay to νµand νe. The neutrinos travel through the Earth towards the detector

while oscillating. Since the matter boosts either the electron neutrinos or antineutrinos, their fluxes provide information on the mass ordering. Details can be found in [36, 74, 75]. For the Earth, the resonance is largest in the region around 3 and 6 GeV, which is within ORCA’s sensitivity range [65].

2.5 Current Parameter Constraints

Three different research groups have constraint the parameters describing the neutrino oscillations using data from several experiments [66]. Table 2.1 shows their current values and limits. The leading mixing angles, θ12 and θ13have been measured accurately, but

θ23still spreads a large part of the parameter space and δCP is almost unconstrained. At

this moment, no preference for either mass ordering exists.

The PMNS model has been shown to explain neutrino oscillations from experiments, but some anomalies remain. Gallium experiments have been measuring a lower rate from their radioactive sources than expected [76–83]. Furthermore, LSND [84, 85] and Micro-BooNE [86] have measured an excess in the low energy region of electron antineutrinos and is still unexplained [87]. Finally, a recalculation of the ¯νe flux showed that all short

baseline reactor experiments have measured fewer events than expected [88–90]. More complete models for nuclear reactor decays have been suggested as a suitable solution for this reactor anomaly [91, 92]. For a full discussion of the neutrino anomalies, see [93]. Since most are only at a 2-3σ level, they might be a statistical fluke or the result of a

(15)

more complex neutrino theory [94]. For now, the 3ν scheme is the best current model in explaining neutrino oscillations.

Parameter Normal Ordering Inverted Ordering θ12 [deg] 33.82+0.78−0.76 33.82+0.78−0.76 θ23 [deg] 49.6+1.0−1.2 49.8+1.0−1.1 θ13 [deg] 8.61+0.13−0.13 8.65+0.13−0.13 δCP [deg] 215+40−29 284+27−29 ∆m2 21[10−5eV2] 7.39+0.21−0.20 7.39+0.31−0.20 ∆m2_3l[10−3eV2] 2.525+0.034_−0.032 −2.512+0.034_−0.032

Table 2.1: Current neutrino oscillation parameters for both mass orderings from [66] without Super-Kamiokande data. ∆m2_3l is ∆m2₃₁ for normal ordering and ∆m2₃₂ for inverted ordering. Errors are 1σ deviations.

(16)

Chapter 3

KM3NeT Infrastructure

3.1 Neutrino Telescopes: ARCA & ORCA

For a better understanding of neutrinos and their astrophysical sources, two next genera-tion neutrino telescopes are being built in the Mediterranean Sea as part of the KM3NeT infrastructure [75]. The ARCA telescope, Astroparticle Research with Cosmics in the Abyss, will look at the universe to find sources of high-energy neutrinos. These travel undisturbed from their source to the Earth and provide information about the sources and the particle acceleration mechanisms. Recently, a similar ice-based Cherenkov neu-trino telescope at Antarctica, IceCube, has identified a blazar as the first high-energy astrophysical neutrino source [95]. The neutrino’s energy ranges from a few TeV to the PeV range. The second KM3NeT telescope, ORCA, will look at neutrinos with GeV energies to measure the neutrino oscillation parameters. Moreover, it will use the method described in Section 2.4 to determine the mass ordering of the neutrino mass states.

Both detectors use the same 3-inch Photomultiplier Tubes (PMTs) housed in pressure-resistant glass spheres called Digital Optical Modules (DOMs) to record Cherenkov radiation from high-energy neutrino interaction products, such as muons, taus, and electrons. Each DOM contains 31 PMTs and other electronics, such as an accelerometer, compass, and many more components. Detailed information about all the contained electronics in a DOM can be found in the KM3NeT 2.0 Letter of Intent and the technical report [75, 96].

18 DOMs are combined in a long vertical string, called the Detection Unit (DU). 115 DUs with all support infrastructure is called a building block. The separation between DUs and DOMs is optimised for cost and performance per detector based on the expected neutrino flux and the energy regime. ARCA will be build 100km offshore from Porto Palo di Capo Passero, Sicily, Italy at a depth of 3500m. 40km from Toulon, France in the Mediterranean Sea at a depth of 2450m the ORCA detector is under construction.

The ARCA detector has DUs with 36m between each DOM and an average of 90m between each DU. ORCA’s design is denser to measure GeV neutrinos. Its DOMs are only separated by 9 meters, and DUs by 20 meters. The full ARCA detector will consist of 2 building blocks, adding up to a total volume of around 1 km3, while the ORCA detector will only be 1 building block with a volume 3 orders of magnitude smaller, around 5.4 × 106 m3.

(17)

3.1.1 Photomultiplier Tubes

The PMTs are photosensitive devices with a 3-inch diameter capable of measuring single photons with wavelength from 280nm to 720nm [97]. When a photon hits a PMT, the photocathode layer releases an electron through the photoelectric effect. The efficiency of this process is the quantum efficiency and is wavelength dependent. Accelerated by an internal electric field, the released electron will hit several dynodes in quick succession, each releasing more electrons. If the charge pulse at the final dynode exceeds a set threshold, the charge pulse information is passed to the electronic within the DOM and a so called ”hit” is registered. Further steps are described in Section 3.4.

KM3NeT’s 31 PMTs design provides several advantages over non-multi-PMT designs, such as ANTARES’s, IceCube’s, and Baikal Gigaton Volume Detector’s (GVD). The smaller PMT size reduces cost per photosensitive area. Furthermore, physical signals are easily recognised by requiring simultaneous hits on multiple PMTs of the same DOM. During extremely bright events, when normally PMTs are over-saturated, away-pointing PMTs are still able to provide valuable information on the event. And finally, the direction of the individual PMTs provides essential information about the direction and location of the event [98].

3.2 Detection Principle

3.2.1 Cherenkov Emission

Neutrino telescopes do not measure the neutrino directly, but through the photons from reaction products. As described in Section 2.2, when a neutrino interacts through CC interactions, it can create neutral and charged particles with an very high energies. When the latter move through a dielectric medium, their surroundings are excited and emit electromagnetic waves spherically. With a sufficiently high velocity of the charge particle, the emitted waves in the direction of movement constructively interfere to a wave front. This cone-like emission is known as Cherenkov radiation. Figure 3.1 depicts how the spherical emission interferes to create the cone-like Cherenkov emission.

cos (θC) =

1

βn (3.1)

The light has a characteristic emission angle, θC given by Equation 3.1. For a

relativis-tic parrelativis-ticle (β = 1) in seawater (n = 1.35) the Cherenkov angle is approximately 42°. The amount of Cherenkov light emitted around this angle is described by the Frank-Tamm equation [3]: d2N dxdλ = 2π αλ2 1 − 1 n2_β2 (3.2) It describes the number of photons per distance travelled (dx) per wavelength (dλ) and depends on the wavelength λ of the emission. A charged particle will continue to emit Cherenkov radiation until its velocity β drops below 1/n.

(18)

ϑ

_c

v=βc

Figure 3.1: Visualisation of Cherenkov emission for charged particles moving through a dielectric medium to the right (red arrow). The particle’s movement causes the medium to emit electromagnetic waves spherically, which interfere to create a cone-like emission due to the high velocity of the charged particle. This is called Cherenkov radiation (blue arrows).

3.2.2 Interaction Signatures

Depending on the neutrino flavour and type of interaction between the neutrino and nucleons, the KM3NeT-ORCA detector distinguishes three main detection signatures. When the energy of the neutrino is large enough, it can recoil a quark and destroy the original nucleon, creating hadronic particles. They interact and decay to create more and more hadrons, resulting in a hadronic shower. This process is part of all neutrino interaction in the detector. Except for the NC channel, it is combined with other signatures. These come from charged current interactions with each neutrino flavour having its own unique signature. Electrons from electron neutrino interactions cause a cascade of photons and electron-positron pairs: an electromagnetic shower. A muon, on the other hand, propagates through the detector with an energy dependent straight path length, and is, thus, referred to as a track.

Besides emitting Cherenkov radiation, the charged particles lose energy through ioni-sation (δ-rays), Bremsstrahlung, and e−e+pair production. In the low GeV regime ionisa-tion is the most dominant energy loss process, but the contribuionisa-tion from Bremsstrahlung

(19)

and pair production increases linearly with energy and is dominant above 1 TeV, which Section 4.3.1 discusses in more detail. The emitted photons through these processes can cause small electromagnetic showers along the track of a muon through pair production or Compton scattering. Tracks, hadronic, and electromagnetic showers are also a signature of tau neutrino interactions, where the tau decays either into a muon, hadron, or electron. These have distinct signatures from the normal interaction due to the time delay in the tau decay, but can be indistinguishable because of the position resolution of the ORCA detector. All interactions and their signatures are summarised in Table 3.1.

Interaction Particle signature Detector signature

νµ CC hadronic shower and µ track

track-like hadronic shower and µ track

(τ±→ µ±_ν

µντ,∼ 17% BR)

ντ CC hadronic and EM shower_(τ_±

→ e±νeντ,∼ 18% BR) point-like or shower-like hadronic showers (τ±→ hadrons, ∼ 65% BR)

νe CC hadronic and EM shower

ν NC hadronic shower

Table 3.1: Different neutrino interactions with neutrinos (dashed black line), muons (orange), taus (green), hadronic showers (blue), and EM showers (red). The particle signature happens and the detector signature is how KM3NeT measured the event. Figure from [98].

(20)

3.3 Background Sources

The Mediterranean Sea is not a perfectly controlled experimental setup, and several background sources obscure the neutrino signal.40K radioactive decay, bioluminescent organisms in the water, and the dark count rate result in the PMT registering ”hits” that do not originate from a photon from a muon. The average rate of these background hits is 8 kHz for a single PMT. Signal-like background comes from atmospheric muons from above the detector.

• K-40 decay The salt water contains radioactive isotopes that emit an electron in its beta decay or electron capture. A dominant source in seawater is the isotope Potassium-40, whose resulting electron has a maximum energy above the Cherenkov threshold. Other40_{K decay products, such as} 40_{Ar, emit photons of 1.46 MeV that}

can Compton scatter to Cherenkov emitting electrons. These processes provide a steady background rate for the detector.

• Dark count rate Even in the absence of light, thermal noise can release an electron from the photocathode layer. These are registered as hits and provide a significant background.

• Bioluminescence Many organisms live within the sea, including a collection of luminescent creatures, such as the pyrosoma and siphonophores. Through a chemi-cal reaction they emit light in the optichemi-cal range, which is visible to the KM3NeT PMTs. It is a slowly changing seasonal background rate, but large disturbances excites the organisms into emitting a bright burst of well above the background rate [99].

• Atmospheric muons Atmospheric showers from cosmic ray interactions create ν, µ and other particles of which ORCA only wishes to measure the neutrinos. The Earth shields the detector from the other shower particles in the up-going direction. In the down-going direction, a few kilometres of water shield the detector. If muons from above are sufficiently energetic, they reach the detector. This atmospheric muon rate is a factor 105 larger than the neutrino rate from below. It is, there-fore, important to reconstruct the direction of the tracks properly to reduce this atmospheric muon background.

3.4 Data Acquisition

When a pulse charge on a PMT is over threshold, its analogue signal is digitised into a time (t) and a time-over-threshold (T oT ). The combination of these two data values is known as an ”L0 hit”, also often referred to as just a ”hit”. Every 100 ms each DOM sends an identically sized time window containing all L0 hits to shore. To maintain time consistency between each DOM, a fibre-optic network, an on-shore White Rabbit switch, and electronics embedded in the DOM work together to synchronise the complete detector

(21)

up to nanosecond precision [75, 100]. After the time slice arrives on-shore over Ethernet, the data stream has to be reduced. Without it, a full building block would result in a stream of 25 Gb/s; much of which is background noise. Therefore, many software triggers run over the timeslices in parallel to search for physics events and reduce the number of background events in the data stream. They look for ≥ N causally related hits from two different PMTs on the same DOM within a 10ns time window (L1 hits). If found, all the L0 hits at the same time as the L1 hits, and those in a time window of approximately 10µs before and after the causally related hits are combined into a single event. This is written to disk and passed on to the reconstruction algorithms. A full overview and description of triggers can be found in [101].

3.5 Monte Carlo Simulation

Monte Carlo simulations are the only method to understand the detector response and develop new algorithms for analysis. They mimic the ORCA detector and its response to neutrino events. For this thesis, the neutrinos and their interactions are generated by gSeaGen [102]. Their energy is between 10 to 100 GeV to allow for better analysis of the reconstruction algorithms, which results in muons between 0 to 100 GeV. For the propagation of the charged particles and photons, a full simulation, KM3Sim, or a tabulated response, JSirene, can be used. The latter only propagates the initial lepton and uses probability density functions to calculate the expected amount and arrival time of light on PMTs [103, 104]. KM3Sim also propagates the initial charge particle, but continues to fully simulate the creation and propagation of photons and new charged particles using a Geant4 based algorithm [105]. The full simulation is used in this thesis. The final steps are to add background hits to the Monte Carlo events and simulate the detector and trigger response, which is done using the program, JTriggerEfficiency [106].

Version 5 of the Monte Carlo production using the above chain for a full ORCA build-ing block with 20m horizontal and 9m vertical separation is used in the analysis of this thesis. Only a limited number of events can be analysed due to computational constraints, thus, in most cases the 1969 events from the file mcv5.0.gsg muon-CC 10-100GeV .km3sim.jte.1.root are used. When more statistics are required, additional Monte Carlo files are included.

(22)

Chapter 4

Current Muonic Event

Reconstruction

and Probability Density Functions

4.1 Track Parameters

After the trigger algorithms from Section 3.4 have found a potential event, the details of the ν-interaction have to be reconstructed. As discussed, the event contains L0 hits; each with the hit time and the time over threshold of the charge pulse. Together with the position and direction of the PMTs, this information is used to estimate the original lepton’s properties. This process is the event reconstruction. For which KM3NeT imple-ments it as a multistage system that distinguishes between shower-like and track-like events. This chapter covers the reconstruction of track-like event, which consists of a pre-reconstruction, JPrefit, and a main-reconstruction phase, JGandalf. Their details are discussed in Sections 4.4 and 4.5. For the shower reconstruction see [75] and [107].

For a track-like event, the position of the lepton is described by three parameters: x, y, and t. No z position is required due to a degeneracy with the time parameter, t. Similarly, the direction of the lepton can be described a unit vector with two independent angles because of the cylindrical symmetry of the track. Instead of the actual angles, two parameters of the unit vector are used in the reconstruction: dx and dy. They are retrieved from the angle using simple geometry, shown in Equations 4.1, where θ and φ are the angles from the z-axis and around the z-axis, respectively. The directional parameters dx and dy can, therefore, only take on values between -1 and 1. The final property that is fitted is the muon energy E. All the track parameters are encompassed in a vector ~θtrack for an easy of notation.

dz = cos θ dx = sin θ cos φ

(23)

4.2 Maximum Likelihood Method

The reconstruction chain uses a maximum likelihood method in which the probability that the data is generated by a specific statistical model is maximised [108]. Each data point (xi) has a probability (P (xi| ~θ)) that it was generated under an assumed model

with corresponding parameters (~θ). Their product gives a likelihood function: L(~θ|~x) =Y

i

P (xi|~θ) (4.2)

When maximised, it gives the most probable model parameters, the maximum likeli-hood estimate. The likelilikeli-hood can, of course, be comprised of more complex combinations of probabilities, as we will see in Section 4.5.

The product of many small probabilities often runs into numerical issues. Thus, the likelihood is represented as a − log(L). This also makes it compatible with many minimisation algorithms. The optimum searching step is often the most complex due to the non-trivial landscape of the likelihood function and requires many initial starting values to cover a large parameters space.

4.3 PMT’s Probability Density Functions

For KM3NeT, ~θtrack contains the model parameters, and the probability (P (x|~θ)) is a

description of the expected number of photo-electrons on a PMT as a function of light arrival time, also known as a probability density function (PDF). Besides depending on the details of the interaction, it is influenced by the water properties and geometry of the event with respect to the PMT.

The PDFs are generated either by transforming data from a full Monte Carlo sim-ulation into chance tables or by using a semi-analytical description. The latter only considers direct and single scattered photons, which is justified by the long scattering length in water. This thesis uses the analytical approach, because it shows the relation between the probability and physical parameters. Since they are an integral part of the reconstruction, we will discuss their origin and dependencies in three steps, based on an internal note by Maarten de Jong [109]. In this approach we will build a description of the expected number of photo-electrons per unit time on a PMT, which can be turned into a PDF by normalising the function with its integral.

1. Light Emission 2. Light Propagation 3. Light Detection

For a simpler PDF, the topology of an event can be rotated and transformed, such that the track is the z-axis, as depicted in Figure 4.1. The PMT is rotated to be in the x − y plane with a closest distance R, zenith angle θ, and azimuth angle φ. This allows for easier storage of the probability values and a clear dependence on three parameters.

(24)

6

z R u (0, 0, 0) u (R, 0, 0) u (0, 0, z) θ0 6 - * z0 x0 y0 J J J J J J J ] P P P P P u φ℘ θ℘

Figure 4.1: Topology of a track event along the z-axis and the PMT in the x-y plane with a shortest distance R. θ℘ and φ℘ are the zenith and azimuth angles of the PMT.

Image from [109].

4.3.1 Light Emission

Cherenkov Emission

Several of the light emission processes are already mentioned in Section 3.2. For minimum ionising particles, such as the low GeV muons, the main emission process is Cherenkov radiation. The amount of light emitted is described by Equation 3.2, but the number of photons at the PMT is different. Using the cone hypothesis, the flux of photons becomes Equation 4.3, with R the shortest distance between the PMT and track.

Φ0(R, λ) = d2N dxdλ 1 2πR sin θC (4.3)

(25)

Showers

Depending on the energy of the muon, different energy loss processes can cause small showers along the track. Due to their cascading nature, showers have a typical longitudinal profile of particles. After an initial steep increase in particles, their energy drops below the pair production limit and their total number decreases. Each particle above the Cherenkov limits emits their own radiation, which results in an effective emission angle of the shower around the normal Cherenkov angle in sea water. The longitudinal profile, (dP_dz) and angular emission profile (d2_{P/[d cos θ}

0dφ0]) are parametrisations from Monte

Carlo data. Appendix A provides more information on the the parametrisations. The longitudinal and angular profiles are combined into a photon flux at a PMT in Equation 4.4 and depend on the emission wavelength (λ) and angles (θ0, φ0). Moreover, dx/dE is

the distance propagated by the shower per unit energy, which is around 4m per GeV in sea water. Φ1(cos θ0, λ) = dx dE d2N dxdλ d2P d cos θ0dφ0 (4.4)

Muon Energy Loss

−dE

dx = a(E) + b(E)E (4.5)

The energy losses are classified ”constant” (a[E]) or linear with energy (b[E]). Ionisa-tion and δ-rays fall in the first category, while bremsstrahlung and e+e− pair production increase linearly and become the dominant energy loss process above 1 TeV, as depicted in Figure 4.2. Since they create small showers, their measurable flux becomes:

Φ2(cos θ0, E, λ) = b(E)EΦ1(cos θ0, λ) (4.6)

The contribution from ionisation is minimal, but δ-rays are measurable. δ-rays are knocked-on electrons with a kinetic energy T . Their energy loss (T d2N/dT dx) is depen-dent on the kinetic energy, which decreases while travelling. By performing an integral over the allowed energies, the number of photons from δ-rays is found. Tmin and Tmax

are constraint by kinematics of the interaction and can be found in [3]. Assuming that their emission is isotropic, the emitted number of photons is described by:

Φ3(cos θ0, E, λ) = d2N dEdλ 1 4π Z Tmax Tmin dT T d 2_N dT dx (4.7) 4.3.2 Light Propagation

Along the way to the PMT, the emitted photons are absorbed, scattered, or unaffected. This can be taken into account by using a full Monte Carlo simulation, a parametrisations, or a simplified analytical solution. The latter is used in this thesis to keep the relation between the PDFs and the physics clear.

(26)

Absorption is taken into account by adding an exponential factor (e−d/λabs_{) to the} expected number of photons, where d is the distance between the emission location and the PMT. λabs is the absorption length. The exact distance is wavelength dependent and

shown in Figure 4.3.

Only single scattering is taken into account for the semi-analytical solution. Other methods can go up to higher order, but all of them will reduce the number of photons hitting a PMT. For the case considered here, a similar exponential factor to the absorption is used with a scattering length (λscat), which is also wavelength dependent, as depicted

in Figure 4.3. At the same time, light that would originally miss a PMT can now scatter into it. We model this indirect light by an effective attenuation length (λatt) instead of

the scattering and absorption coefficients [109].

Figure 4.2: Showing the energy loss per unit track from a muon. The ionisation loss is almost constant over energy (dot-dashed), while the bremsstrahlung (thin solid line) and the pair production (dashed) increase linearly with energy. Nuclear in-teractions (dotted) will not be considered here. Image from [110].

300 350 400 450 500 550 600 650 [nm] λ 0 20 40 60 80 100 120 140 length [m] scattering absorption scattering absorption

Figure 4.3: The scattering (red) and absorption (black) distances as a function of wavelength for KM3NeT adapted from [109].

4.3.3 Light Detection

Two more obstacles stand in the way before a PMT measures a photon. First, the angle (θ) between the incoming photon and the direction of the PMT determines whether

it can actually see the photon. This is the angular acceptance () and is created from Monte Carlo simulations and measurements. It is tabulated to allow for easy access and interpolation between calculated values. Figure 4.4 shows the angular acceptance over

(27)

cos θ’s range. The second obstacle is the quantum efficiency of the PMT. As described in

Section 3.1.1, it is the efficiency per wavelength with which a photon releases an electron in the photocathode layer. Figure 4.5 shows that it further reduces the chance a PMT registers a hit. 1 − −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 ∅ θ cos 0 10 20 30 40 50 60 70 ∈

Figure 4.4: Angular acceptance of a PMT. The cosine of the incoming photon angle (cos θ) determines if the PMT registers

it. Image from [109] and originally from Monte Carlo data.

300 350 400 450 500 550 600 650 700 [nm] λ 0 0.05 0.1 0.15 0.2 0.25 0.3 QE

Figure 4.5: Quantum efficiency. The ratio of released electrons and incoming photons in the photocathode layer per unit wave-length. Image from [109].

4.3.4 Direct and Indirect Light as PDFs

The emission, propagation, and detection parametrisations are combined to form a function to describe the expected number of photo-electrons with a specific arrival time on a PMT. This can be divided by its integral to give an actual probability density function for the expected arrival time. However, we will refer to the non-normalised functions as PDFs because of their implementation within the KM3NeT software.

In the analytical case, a total of 6 function can be distinguished; two for each emission process. Their subtleties can be found in [109]. But here, we will consider a few of their important aspects. The direct light PDFs are the most straightforward with the three creation steps clearly visible; see Equations 4.8, 4.9, and 4.10. All three equations are dependent on the effective area of the PMT. The first has a direct dependence with the photocathode area (A), while the latter two are related to it by the solid angle of the PMT (dΩ). Furthermore, the first two direct light emission processes are wavelength dependent (∂t/∂λ). The arrival time from the muon energy loss, on the other hand, is dependent on the z position along the track (dt/dz).

Indirect light, on the other hand, requires integration over several angles and distances, as seen in Equations 4.11, 4.12, and 4.13. These include the wavelength (λ), z position

(28)

(not for 4.12), zenith emission angle to the track (θ0), and for the last two PDFs the

azimuth emission angle (φ0). Other important changes compared to the direct case are

the factors 1/λs, dN/dx, and 1/2π. These, respectively, are the probability of scattering

per unit length, the number of scatterings per track length, and a normalisation for the number of photons after the integration over angles. The latter is only required in the Cherenkov light PDF. As discussed in Section 4.3.2, the scattering and absorption have been replaced by a single attenuation factor (e−d/λatt_{) and a scattering probability} dPs/dΩs. Further details can be found in [109].

Direct Muon Cherenkov Light dnp.e.

dt = Φ0(R, λ) A ∂t

∂λ −1

(cos θ) QE(λ) e−d/λabs e−d/λs (4.8)

Direct Shower Light dnp.e.

dt = Φ1(cos θ0, λ) ∂t

∂λ −1

(cos θ) QE(λ)e−d/λabse−d/λsdΩ (4.9)

Direct Light from Muon Energy Loss dnp.e. dt = Z dλ X z=z1,z2 dt dz −1

Φ2,3(cos θ0, E, λ) dΩ (cos θ)QE(λ)e−e/λabse−d/λs

(4.10) Indirect Muon Cherenkov Light

dnp.e. dt = Z Z Z dλdzdφ0 1 2π dN dx 1 λs ∂t ∂u −1

(cos θ) QE(λ)e−d/λatt

dPs

dΩs

dΩ (4.11)

Indirect Shower Light dnp.e. dt = Z Z Z dλdφ0d cos θ0Φ1(cos θ0, λ) 1 λs ∂t ∂u −1

(cos θ) QE(λ)e−d/λatt

dPs

dΩs

dΩ (4.12) Indirect Light from Muon Energy Loss

dnp.e. dt = Z Z Z Z dλdzdφ0d cos θ0Φ2,3(cos θ0, E, λ) 1 λs ∂t ∂u −1 (4.13) (cos θ) QE(λ)e−d/λatt

dPs

dΩs

(29)

Tabulated PDFs

Each time the probability for a hit on a PMT is required, all the integrals must be done for each wavelength. This is time consuming and inefficient. To improve performance, the PDFs are integrated over wavelengths from 300 to 700 nm; the range where the PMT is most sensitive. The total probability is calculated for a grid of the following parameters:

• Distance between track and PMT (R) • Zenith angle of the PMT (θ)

• Azimuth angle of the PMT (φ)

• Time difference from the Cherenkov light arrival time (dt)

The calculated values are stored in a 4 dimensional table that allows for interpolation between the calculated grid values. These PDFs store the expected number of photo-electrons and calculate the probability using stored integrated values when calculating the likelihood.

4.4 The Pre-reconstruction: JPrefit

The pre-reconstruction phase, JPrefit, generates starting values for JGandalf, the main reconstruction algorithm. It does this by 3D clustering causally correlated L1 hits and finding their centre of weight. By performing a directional scan over the whole sky with 5° between guesses, around 800 tracks are created around the centre of hits. At this stage, the tracks length is assumed to be infinite. Per assumed direction a hit selection is performed using a 1D clustering algorithm that looks for causality between the track and the L0 hits within a 50m radius cylinder around the track. The radius of this cylinder is referred to as the road width. The algorithm uses the time information of the selected hits to create the χ2 function is Equation 4.14, where ti and σi are the time of the hit

and the time resolution of the PMT, while tC is the expected Cherenkov arrival time.

χ2 = hits X i (ti− tC)2 σ2 i (4.14) JPrefit searches for the optimal set of parameters that minimise the χ2. The set consist of the time parameter (t) and two positional parameters orthogonal to the direction of the track: x and y. The three parameters are combined with the direction to give a track guess. A quality is attributed to the track parameter, as defined in Equation 4.15, where NDF is the number of hits minus 3 for the fit parameters. The best 36 quality tracks are passed on to JGandalf. In some cases more tracks are selected based on whether downward pointing tracks are already included in the 36 selected tracks.

Q = NDF −1 4

χ2

(30)

4.5 The Main Reconstruction: JGandalf

The main reconstruction uses the maximum likelihood method to find the event hypoth-esis that best fits the observed data. It also performs a L0 hit selection within a road width of 50m. The time of the hits ( ~D) are compared to the probability of light to arrive at that time (P (Di|~θtrack)). Using these, the likelihood in Equation 4.16 is minimised for

all track parameters (~θtrack), expect energy. Since the arrival time probability is required,

the tabulated PDFs from Section 4.3 are used. JGandalf’s likelihood gives the quality of the track and its parameters. The process is repeated for all 36 input tracks from JPrefit and they are ranked according to their likelihood values. The highest ranked track should, in principle, provide the best possible track parameters, if the likelihood is properly defined. L(~θtrack| ~D) = hits Y i=1 P (Di|~θtrack) (4.16)

4.6 Other Algorithms: JStart & JEnergy

After JGandalf the main reconstruction is over, but two more steps exist. The first, JStart, removes the assumption that the track is infinite and searches for the begin and end point of the track. This allows for a selection of hits consisting mostly of signal hits. This is required for the next step: JEnergy, which tries to reconstruct the energy of the muon and set a lower limit on it. This work does not consider these steps, as they no longer influence the final direction of the track.

4.7 Intrinsic Limits of the Reconstruction

The true interest of the reconstruction is to find the original neutrino direction and energy, but even with the perfect reconstruction algorithm only the true muon direction and a lower limit for its energy can be found. Due to the scattering processes, the energy and direction of the muon will differ from the neutrino. This introduces a fundamental limit for the reconstruction of the direction and the energy. The angle between the true neutrino and the muon direction shows the intrinsic limit of the reconstruction in Figure 4.6 (left). While the size of the limit is not necessarily a problem, it is the spread that introduces uncertainty in the neutrino’s direction. Luckily, higher neutrino energies lead to more forward boosting and push the median and the spread to lower values, as expected. Furthermore, the combined fraction of neutrino and antineutrino energy given to the muon decreases with higher energies, as Figure 4.6 (right) shows. The change from quasi-elastic and resonance scattering to deep inelastic scattering opens up the phase space for the muons energy with a larger component going to the shower. Another important element for the energy is whether the muon is contained within the detector. If it travels even a bit outside the detector volume, it becomes impossible to estimate

(31)

the original muon energy and only a lower limit can be set. These aspects will all restrict the measurements for the mass hierarchy, but a good directional reconstruction is still necessary to reduce background events.

1 10 102 [GeV] true ν E 0 5 10 15 20 25 30 35 40 45 [deg]_µ , ν θ 1 10 102 [GeV] true ν E 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 true ν / E true µ E

Figure 4.6: The median of the intrinsic limit (blue) for the angle between the muon and (anti)neutrino (left) and the energy fraction from the (anti)neutrino going to the muon (right). Both neutrinos and antineutrinos are included. The 16% and 84% quantiles are

(32)

Chapter 5

Evaluation of JGandalf ’s Input

Within the ORCA detector, the energy and direction are the most important parameters to reconstruct because of their sensitivity to the mass hierarchy. Their reconstruction principles are quite different, and in this thesis we will only consider the directional reconstruction of the muon. The accuracy of the measured direction is influenced by many parameters of the reconstruction, as introduced in Chapter 4, for example, the hit selection, the number of track guesses, and assumed coincidence window. These parameters are constantly optimised for the best results. Therefore, we focus on improving the reconstruction algorithms themselves.

Section 6 and 7 discuss the improvements to JGandalf’s input and to JGandalf itself. It is important to first understand the current reconstruction framework and its limitations. Therefore, this section covers the input JGandalf requirements and whether JPrefit reaches these.

5.1 The Positions of JGandalf ’s Minima

JGandalf uses a Levenberg-Marquardt minimisation algorithm, which is robust and finds solutions far from its starting value [111]. However, the KM3NeT collaboration knows that JGandalf does not perform well when starting values are far away from the true values. No quantitative limits for this exist and their nature has not yet been explored. For an event selected based on its high directional accuracy after the JGandalf algorithm, the -log likelihood space for the x parameter is shown in Figure 5.1. The other fit parameters are fixed to their truth values. The space has a minimum at the correct true value, indicated in red. However, the global minimum is not at the true parameter. Similar results can be seen in Appendix B for the other fit parameters and a collection of events from good to badly reconstructed events.

The examples indicate a structural problem with the global minimum in the JGandalf likelihood. To explore this further, a scan over all parameters of 1969 events is performed. The global minimum is marked in the ranges -200m to 200m for the positional parameters and -1 and 1 for dx and dy. The global minimum positions from the X and Y parameters in Figure 5.2 show that most minima lie at the scan boundaries. Their combined count is an order of magnitude larger than at the true parameter. The same behaviour, although less extreme, shows in the directional fit parameters in Figure 5.3. The amount of tracks

(33)

at the boundaries is double those at the true parameter. JGandalf’s likelihood is not properly defined. Minimisation to the global truth will lead to a wrong set of parameters. It is, therefore, the most probable cause for the initial value problem.

150 − −100 −50 0 50 100 150 Distance [m] 100 200 300 400 500 600 700 800 900 -Log LL

Figure 5.1: A negative log likelihood space scan for the X parameters done by JGandalf for a properly reconstructed event, where the Monte Carlo truth value is marked in red.

200 − −150 −100 −50 0 50 100 150 200 X deviation 1 10 2 10 3 10 Counts (a) 200 − −150 −100 −50 0 50 100 150 200 Y deviation 1 10 2 10 3 10 Counts (b)

Figure 5.2: The location of the global minimum for the X (a) and the Y (b) parameters. The true value is at 0.

(34)

1 − −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 DX deviation 1 10 2 10 3 10 Counts (a) 1 − −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 DY deviation 10 2 10 3 10 Counts (b)

Figure 5.3: The location of the global minimum for the DX (a) and the DY (b) parameters. The true value is at 0.

5.2 Quantifying JGandalf ’s Directional Input

Let us quantify the influence of JGandalf’s initial value problem on the reconstruction. By searching till what deviation from the true direction JGandalf is still able to improve the directional information, a limit can be set. For this purpose, a directional accuracy is defined as the dot product between the fitted or initial track with the true direction. This is also known as the directional error.

cos θARes= ~dinitial/fit· ~dtrue (5.1)

5.2.1 Deviates from the True Direction

The initial and output directional accuracy is compared by generating 639 tracks with 15 different deviations between 0 and 40 degrees from the Monte Carlo truth. Other fit parameters have been set to their corresponding truth value. After JGandalf minimises the likelihood, the median of the final directional accuracy as a function of E_µtrue is calculated for each input track. The process is repeated for numerous events to reach enough statistics for Figure 5.4, where the relative change compared to the directional accuracy of the initial track is shown.

JGandalf is successful when the final angle is smaller than the initial angle. This is the case for initial angles up to 10° above 20 GeV and up to 15° in the 10-20 GeV range. Below this energy range, the reconstruction is unable to properly reconstruct events due to the low number of signal photons measured. Tracks starting closest to their true value minimise to a worse directional accuracy, as the red in the lower left in Figure 5.4 indicates.

Above 10° JGandalf is unable to reconstruct up to the same level of accuracy, but a slight improvement can still be seen for high energies and a large improvement around 10 GeV as can be seen in Figure 5.5. These effects decrease as the initial angle gets

(35)

larger. Above 90 GeV statistical effects are present, but with a sample of more high-energy neutrinos these should disappear. The best possible reconstruction achievable by JGandalf is indicated by the yellow line in Figure 5.5, because it starts at the true parameters. For initial directional deviations of 10° or less, the errors on the final tracks are similar to starting at the true values. For example, the 5.71429° deviate reconstructs, on average, to the same directional accuracy as the true tracks.

Commonly, the reconstructed angles are plotted without the energy dependence as a count plot in log scale to show the sub-degree nature of the reconstruction. JGandalf’s best achievable directional accuracy is shown in Figure 5.6. Brute forcing many tracks would also result in this distribution if the likelihood is correctly defined. This, however, is computationally intensive and smarter methods, such as in KM3NeT, are required. This does result in an acceptable loss in the directional accuracy.

0 10 20 30 40 50 60 70 80 90 100 true µ E 0 5 10 15 20 25 30

Initial Directional Accuracy [deg]

1 − 0.8 − 0.6 − 0.4 − 0.2 − 0 0.2 0.4 0.6 0.8 1

Relative Directional Accuracy Change

Median of Final ARes

Figure 5.4: The relative directional accuracy change by JGandalf as a function of initial directional accuracy and true muon energy. The initial tracks are directional deviates from the Monte Carlo truth. Above 90 GeV statistical effects play a role.

5.2.2 Deviates using JPrefit tracks

One drawback of the MC deviates is that it only takes into account change in directional parameters. In reality, the position also contains some deviations and will influence the reconstruction. For this purpose, we use JPrefit’s 36 track output to generate a sample of angles and positions. This does generate a biased set of tracks, since a selection and fit has been applied to the tracks.

The relative change in directional accuracy in Figure 5.7 shows similar patterns as the MC deviates. However, the energy threshold has dropped to 10 GeV and the angle limit has increased to 15°. Figure 5.8 also shows that the improvements around 10 GeV

(36)

0 10 20 30 40 50 60 70 80 90 100 µ true E 0 5 10 15 20 25 30

Median Final Directional Accuracy [deg]

Median of Final ARes

Initial Dir.Acc. [deg] MC Truth; 0 5.71429 14.2857 20 28,5714

Median of Final ARes

Figure 5.5: Median final directional accuracy for five different initial directions. The yellow line is the reconstruction started from the Monte Carlo truth and is the best achievable directional accuracy for JGandalf.

have been removed. Angles above the limit show much more improvements in the final median directional error than with the Monte Carlo deviates. The track selection by JPrefit is the most probable cause for the increase in these limits, because it only selects track from its collection of 800 tracks with a high quality factor.

(37)

3

−

10 10−2 10−1 1 10 102

Directional Accuracy [deg] 0 10 20 30 40 50 60 70 80 90 Counts

Figure 5.6: The best achievable directional accuracy for JGandalf by starting the reconstruction at the Monte Carlo truth parameters.

0 10 20 30 40 50 60 70 80 90 100 true µ E 0 5 10 15 20 25 30

Initial Directional Accuracy [deg]

1 − 0.8 − 0.6 − 0.4 − 0.2 − 0 0.2 0.4 0.6 0.8 1

Relative Directional Accuracy Change

fmedianARes

Figure 5.7: The relative directional accuracy change as a function of initial directional accuracy and true muon energy. The initial tracks are the 36 selected track from JPrefit. These have been run through JGandalf to give the fitted tracks.

(38)

0 10 20 30 40 50 60 70 80 90 100 µ true E 0 2 4 6 8 10 12 14 16 18 20 22

Median Final Directional Accuracy [deg]

Median of Final ARes

Initial Dir.Acc. [deg]

3 2 0-2 3 1 - 5 3 2 2 3 2 8 - 10 3 2 16 - 18 - 24 3 1 21

Median of Final ARes

Figure 5.8: Median final directional accuracy for five different initial directional accuracy ranges. Initial tracks are JPrefit tracks and are run through JGandalf to get a fitted track.

5.3 JGandalf ’s Input Limits

We conclude that the initial directional accuracy for the JGandalf algorithm has to be below 10° to 15° for energies above 20 GeV. The influence of the errors in the positional parameters on the directional reconstruction is still unclear, but the JPrefit tracks increase the range in initial directional accuracy up to 20° and the energy range goes down to 10 GeV. Within these limits, the JGandalf algorithm can still, on average, reconstruct the events to a similar accuracy as when starting at the true parameters. The increase in parameter space for JPrefit tracks is most likely caused by the pre-selection of tracks.

5.4 Analysis of JPrefit’s Output

The JPrefit algorithm increases the directional error JGandalf is able to handle, but the question remains whether the selected 36 tracks are within the directional input limits of JGandalf.

JPrefit calculates 800 tracks, but only passes on the best 36. The normalised direc-tional accuracy distributions for a large collection of events are differ significantly, as shown in Figure 5.9a, which are not corrected for their phase space. The best 36 tracks

Muonic event reconstuction in KM3NeT-ORCA

MSc Physics and Astronomy

Master Thesis

Muonic event reconstruction in KM3NeT-ORCA

Contents

Chapter 1

Introduction

Chapter 2

Neutrino Oscillations

2.1

The History of the Neutrino

2.2

Neutrino Interactions

2.3

Neutrino Oscillations

2.4

Matter Oscillations

2.5

Current Parameter Constraints

Chapter 3

KM3NeT Infrastructure

3.1

Neutrino Telescopes: ARCA & ORCA

3.2

Detection Principle

ϑ

v=βc

3.3

Background Sources

3.4

Data Acquisition

3.5

Monte Carlo Simulation

Chapter 4

Current Muonic Event

Reconstruction

and Probability Density Functions

4.1

Track Parameters

4.2

Maximum Likelihood Method

4.3

PMT’s Probability Density Functions

6

4.4

The Pre-reconstruction: JPrefit

4.5

The Main Reconstruction: JGandalf

4.6

Other Algorithms: JStart & JEnergy

4.7

Intrinsic Limits of the Reconstruction

Chapter 5

Evaluation of JGandalf ’s Input

5.1

The Positions of JGandalf ’s Minima

5.2

Quantifying JGandalf ’s Directional Input

Median of Final ARes

Median of Final ARes

Median of Final ARes

Median of Final ARes

5.3

JGandalf ’s Input Limits

5.4

Analysis of JPrefit’s Output