33 RD I NTERNATIONAL C OSMIC R AY C ONFERENCE , R IO DE JANEIRO 2013
T HE A STROPARTICLE P HYSICS C ONFERENCE

Robust Statistics in IceCube Initial Muon Reconstruction
T HE I CE C UBE C OLLABORATION1 , B. R ECHT 2 , C. R É2
1 See special section in these proceedings

2 Dept. of Computer Sciences, University of Wisconsin, Madison, WI 53706, USA.

mwellons@cs.wisc.edu
Abstract: In the IceCube Neutrino Detector, muon tracks are reconstructed from the muon’s light emission. The
initial track “linefit” reconstruction serves as a starting point for more sophisticated track fitting, using detailed
knowledge of the ice and the detector. The new approach described here leads to a substantial improvement of theaccuracy in the initial track reconstruction for muons. Our approach is to couple simple physical models with robust
statistical techniques. Using the metric of median angular accuracy, a standard metric for track reconstruction, thissolution improves the accuracy in the reconstructed direction by 13%.
Corresponding authors: M. Wellons1

1 Computer Sciences Department, University of Wisconsin, Madison

Keywords: Icecube, Muons, Track Reconstruction.

Signal [VEM]

H. KOLANOSKI I CE T OP OVERVIEW

reconstruction algorithm used to generate
the initial track
S = 57.21 ± 2.11
PRELIMINARY
β = 2.93 ± 0.05
position and directionof detected muonsin the IceCube
10
detector. We achieve this improvement in accuracy with theaddition robust statistical techniques to the reconstruction
10
algorithm.
103

125

2

1

2

10-1

Background
High Gain

Low Gain

100
200detector
300 is composed
400
500of 5,160600
The IceCube
opticalPerp. distance to shower axis [m]detectors, each composed of a photomultiplier tube (PMT)
FigureReconstructionof [3].
showerthe over 86
and2:onboard
digitizer
Theparameters
PMTs arefrom
spreadlateraldistribution.
verticalstrings arranged in a hexagonal shape, with a total
instrumented volume of approximately one cubic kilometer.
on a given
are separated
bywideThebinsPMTs
for a durationof 422string
ns. DOMs,
electronicsverticallyandreadoutschemeare string-to-string
the same as for thein-ice detector.
17 m,
and theseparationis roughly 125 m.
the telescope,
it sometimes
The twoWhen
DOMsainneutrinoeach tank enters
are operated
at different PMT
withand generates
muon.
As the muon
gainsinteracts(1 · 105 and5 · the106 )ice
to covera dynamic arangeof morethoughthe charges
detector,radiates in
light
[4],
thantravels
104 . The
measuredareitexpressed
unitsof which is
Figure 1: The IceCube Observatory with its components ‘verticalobserved
by thePMTs(VEM)
and divided
intobydiscrete
equivalentmuons’
determinedcalibrat-hits [5]. A
Fig.DeepCore
1: The IceCube
neutrino detector in the Antarctic ice. ing each DOM with muons (see ref. [1]).
and IceTop.
collection of hits is called an event, and when the number
A picture of the Eiffel Tower is shown for scale.
of hitsthe
in readoutan event
sufficiently
large, the muontrack
To initiate
of is
DOMs,
a local coincidenceofreconstruction
algorithm
is triggered.
the twohigh gain DOMsof a stationis required. This re3of 1 km at a depth between 1450 m and 2450 m (Fig. 1). In
in a station trigger rate of about 30 Hz compared to
the lower part of the detector a section called DeepCore is sults2.1
Ray Muons
about 1600Cosmic
Hz of a single
high gain DOM at a threshold
moredensely
instrumented.
The
main
purposeof
IceCube1 IntroductionIn
additiontoneutrinos,
muons
can
be generated by
of
about
0.1
VEM.
The
data
are
writtento alsoa permanentis the detection of high energy neutrinos from astrophysicalstoragemedium,
and
are
thus
availablefor
analysis,
if thecosmic
rays.
One
of
the
simplertechniques
used
to separate
Thesources
IceCubeneutrino
detector
searchesforneutrinos
that
via the Cherenkov light of charged particles generof sixmuons
or morefrom
DOMscosmic
are launchedby a local
coinneutrino
ray muonsis reconstructing
areatedgeneratedby theuniverse’s
mostviolent
astrophysicalin neutrino
interactionsin theice orthe rockbelow the readouts
This leads
to aandtriggerthreshold ofwhether
about 300theTeV.
the muontrack
determining
muon wasevents:ice. exploding stars, gamma ray bursts, and cataclysmic cidence.
is always
read
in casea triggerout of thetravelingIceTop
downwards
into
theoutEarthor ofupwards
phenomena involving black holes and neutron stars [1]. Additionally,by anothercomponent(and vice
Earth.
Sincedetector
neutrinos
can penetrate
theversa).
EarthForbut cosmic
These neutrinos are detected by the charged particles, often issued
IceTop: The IceTop air shower array is located above each single tank above threshold, even without a local coray muons cannot, it follows that a muon traveling out ofa muon, produced in their interaction with the rock or ice
IceCube at a height of 2832 m above sea level, correspond- incidence, condensed data containing integrated charge andthe Earth must have been generated by a neutrino. Thus, bynear the detector.
ing to an atmospheric depth of about 680 g/cm2 . It consists time stamp are transmitted. These so-called SLC hits (SLC
selecting only the muons that are reconstructed as up-going,The
detector,roughlyone
cubic
kilometer
in
size,
is
of 162 ice Cherenkov tanks, placed at 81 stations and dis- = ‘soft local coincidence’) are useful for detecting single
the neutrino muons can, in principle, be identified.
locatednear
Poleis spacing
buried attributed
overthean geographic
area of 1 km2Southon a gridwithandmean
muons in showers where the electromagnetic component
While separation is possible in principle, the number
depths
1.5-2.5kminthe
Antarctic
ice
[2].
The
detectorisof 125 m (Fig. 1). In the center of the array, three stations has beenabsorbed (low energies, outer region of showers,ofillustratedin
Figure1andamore
complete
descriptionishave been installed at intermediate positions. Together inclinedobserved
showers). cosmic ray muons exceeds the number ofobserved neutrino muons by more than five orders ofgiven
Section2.
withinthe
neighbouring
stations they form an in-fill array for For monitoring
events via rate variations, the timeThis
describes
an improvementin themagnitudetransient[6]. Thus,
high-accuracy reconstructions aredensermanuscript
shower sampling.
Each station
comprises two cylinof single hits in different tanks with various thresholds aredrical tanks, 10 m apart from each other, with a diameter of histogrammed.
1.86 m and filled with 90 cm ice. The tanks are embedded into the snow so that their top surface is level with theICRC 2013 Template
33 RD I NTERNATIONAL C OSMIC R AY C ONFERENCE , R IO DE JANEIRO 2013

critical for preventing erroneously reconstructed cosmic
ray muons from dominating the neutrino analysis. While
this technique does not recover neutrinos from the top halfof the sky, there are alternative techniques that attempt torecover down-going neutrinos [7].

2.2Challenges in Neutrino Detection

There are several challenges for the reconstruction
algorithms used in the detector.
Modeling Difficulties The underlying physics of thesystem are nontrivial to model. The muon’s light is scatteredby the dust impurities and air bubbles in the ice medium.
This scattering cannot be analytically well-approximated,
and the scattering properties of the ice vary with depth [8].
These challenges make it difficult to design a complete
model of the muon’s light scattering.
Noise The outliers inherent in the data present anadditional challenge. The PMTs are so sensitive to lightthat they can record hits from the radioactive decay in thesurrounding glass [9].

3Improvements to Muon Track
Reconstruction

We now discuss the improvements we have made to thereconstruction algorithm. By augmenting the reconstruction
algorithm using robust data analysis techniques, we improve
the reconstruction algorithm’s accuracy.

3.1

Algorithm Improvement

The accuracy of the SPE reconstruction is dependent on
the accuracy of the seed. Given a seed that is inaccurate by
6◦ or more, SPE often cannot recover, and can produce areconstruction result that is inaccurate by 6◦ or more. In
addition, the likelihood space for the SPE reconstruction cancontain multiple local maxima, so improving the accuracyof a seed already near the true solution will improve theaccuracy of the SPE reconstruction . Thus, we focused ourwork on improving the quality of the seed.
As the muon travels through the detector, it generates hits.
As indicated in Equation 1, linefit fits a line to these hits,
weighting each hit quadratically in its distance from this line.
This quadratic weighting makes the reconstruction result
sensitive to outliers. There are two reasons why outliers
may appear far from the muon track:

Computational Constraints Reconstruction algorithms
need to be efficient enough to process about 3,000 muons
per second with the computing resources available at theSouth Pole. Thus, algorithms with excessive computational
demands are disfavored.

1. Some of the photons can scatter in the ice and get
delayed by more than a microsecond. When thesescattered photons are recorded by a PMT, the muon
will be over 300 m away, so these photons are nolonger useful indicators of the muon’s position.

2.3

2. While the noise reduction filters remove most of theoutlier noise, the noise hits that survive are unrelated
to the muon.

Prior IceCube Software

Starting with the positions and times of each hit, the detector
reconstructs the muon track. Once the data is collected, it ispassed through a series of filters that removes hits isolated
in space and time [10].
After removing outliers, the data is processed using asimple reconstruction algorithm, linefit, which finds the
track that minimizes the sum of the squares of the distances
between the track and the hits. More formally, assume
there are N hits; denote the position and time of the ith
hit as (~xi ,ti ) ∈ R3 × R. Let the muon have a reconstructedvelocity of ~v, and let (~x0 ,t0 ) be a point on the reconstructed
track. The linefit reconstruction solves the least-squares
optimization problemNmin ∑ ρi (t0 ,~x0 ,~v)2 ,(1)

ρi (t0 ,~x0 ,~v) = k~v(ti − t0 ) +~x0 −~xi k2 ,(2)t0 ,~x0 ,~v i=1

where

where k · k2 is the 2-norm.
The linefit reconstruction is primarily used to generate
an initial track to be used as a seed to a more sophisticatedreconstruction.
The reconstruction algorithm for the sophisticatedreconstruction is Single-Photo-Electron-Reconstruction(SPE) [6]. SPE takes the result of the least-squares
reconstruction and event data, and uses a likelihood
maximization algorithm to reconstruct the muon track. The
SPE reconstruction typically takes about two orders ofmagnitude more time to compute than linefit. The complete
reconstruction process is outlined in Figure 2.

To solve the outlier problem we made two changes:
improve the modeling of the scattering and replace theleast-squares optimization problem with a robust line-fitting
algorithm.
3.1.1 Improving the Scattering Model
The least-squares model does not model the scattering. Thus,
hits generated by photons that scattered for a significant
length of time are not useful predictors of the muon’s
position within this model. We found that a filter could
identify these scattered hits, and improve accuracy byalmost a factor of two by removing them prior to performing
the fit.
A hit (~xi ,ti ) is considered a scattered hit if there exists aneighboring hit (~x j ,t j ) that is within a distance of r and has
a time coordinate, t j , that is earlier than ti by an amount oftime given by t. If (~xi ,ti ) is a scattered hit, it is filtered out.
More formally, let H be the set of all hits for a particular
event. Then, we define the scattered hits as(~xi ,ti ) | ∃(~x j ,t j ) ∈ H : ~xi −~x j 2 ≤ r and ti − t j ≥ t .
(3)
Optimal values of r and t were found to be 156 m and
778 ns, respectively.
3.1.2 Adding Robustness to the Model
As described in Section 2.3, the least-squares model gives
all hits quadratic weight, whereas we would like to limit theweight of the outliers. Some models in classical statistics
marginalize the weight of outliers. We find that replacing
the least-squares model with a Huber reconstruction [11]improves the reconstruction accuracy.

ICRC 2013 Template
33 RD I NTERNATIONAL C OSMIC R AY C ONFERENCE , R IO DE JANEIRO 2013
























Fig. 2: The reconstruction pipeline used to process data in the IceCube detector. Each point indicates a PMT that recorded
a photon (PMTs that recorded nothing are omitted). After initial data is collected and passed though some noise filters, thedata is processed by a linefit (solid line), which is used as the seed for the SPE (dashed line). The SPE reconstruction isthen evaluated as a potential neutrino. Our work on the reconstruction problem makes changes to the linefit reconstruction
algorithm (indicated by the dashed box).
For the Huber reconstruction, we replace Equation 1 withthe optimization problem:Nmin ∑ φ (ρi (t0 ,~x0 ,~v)),

t0 ,~x0 ,~v i=1

where the Huber penalty function φ (ρ) is defined as 2
ρ
if ρ < µφ (ρ) ≡.
µ(2ρ − µ) if ρ ≥ µ

(4)

(5)

Here, ρi (t0 ,~x,~v) is defined in Equation 2 and µ is a constantcalibrated to the data (for this application, the optimal valueof µ is 153 m).
The Huber penalty function has two regimes. In thenear-hit regime (ρ < µ), hits are assumed to be strongly
correlated with the muon’s track, and the Huber penalty
function behaves like least squares, giving point quadratic
weight. In the far-hit regime (ρ ≥ µ), the Huber penalty
function gives points a weaker linear weight, as they are
more likely to be noise.
In addition to its attractive robustness properties, the
Huber reconstruction’s weight assignment also has theadded benefit that it inherently labels points as outliers(those with ρ ≥ µ). Thus, once the Huber reconstruction iscomputed, we can go one step further and simply remove thelabeled outliers from the dataset. A better reconstruction isthen obtained by computing the least-squares reconstruction
on the data with the outliers removed.
3.1.3 Implementation
Our scattering filter has a worst-case complexity that isquadratic in the number of PMTs that recorded a hit, but thisis typically only between 10 and 100 PMTs. Unlike linefit,
the Huber regression does not have a closed form solution,and thus must be solved iteratively. We use an alternating
direction method of multipliers [12] to implement the Huber
regression.

3.2

Results

We now present our empirical results, which validate
our changes to the linefit. We also present our runtime
performance results.
3.2.1 Accuracy Improvement
Our goal is to improve the accuracy of the reconstructionin order, and to better separate neutrinos from cosmic rays.

Thus we present three measurements: (1) the accuracychange between linefit and the new algorithm, (2) the
accuracy change when SPE is seeded with the newalgorithm, and (3) the improvement in separation between
neutrinos and cosmic rays.
To measure the accuracy improvement, we use the metricof median angular resolution θmed , which is the arc distance
between the reconstruction and the simulated true track.
Our dataset is simulated neutrino data designed to be similar
to that observed at the Pole. We find that we can improvethe median angular resolution of the simple reconstruction
by 57.6%, as shown in Table 1.
Table 1: Median angular resolution (degrees) forreconstruction improvements. The first line is the accuracyof the prior least-squares model, and the subsequent lines
are the accuracy measurements from cumulatively adding
improvements into the simple reconstruction algorithm.
Algorithm
θmed (◦ )Linefit Reconstruction (Least-Squares)
9.917
With Addition of Scattering Filter
5.205
With Addition of Huber Regression
4.672
With Addition of Outlier Removal
4.211
We also find that seeding the SPE reconstruction
with the improved simple reconstruction generates animprovement in the angular resolution of 12.9%, and thatthese improvements in the reconstruction algorithm resultin 10% fewer atmospheric muons erroneously reconstructed
as up-going, and 1% more muons correctly reconstructed
as up-going.
3.2.2 Runtime Performance
We now report the runtime of our implementation, whichis written in C++. The individual mean runtime of eachcomponent of the new algorithm is presented in Table 2.
As shown, our new algorithm is more computationallydemanding than linefit, but only by approximately a factorof six.

4

Conclusions

Muon track detection is a challenging problem in the
IceCube detector. We achieve a 13% improvement inICRC 2013 Template
33 RD I NTERNATIONAL C OSMIC R AY C ONFERENCE , R IO DE JANEIRO 2013

Table 2: The mean runtime for each component of the newsimple reconstruction, contrasted with the mean runtimeof the original linefit. As shown, the total runtime isapproximately six times that of the original linefit.
Algorithm
Runtime (µs)
Linefit Reconstruction (Least-Squares)
24.2
Scattering Filter
56.6
Huber Regression
47.5
Outlier Removal
51.8
reconstruction accuracy with the addition of a scatteringfilter, and a more robust line-fitting algorithm. We achievethese results with a reconstruction algorithm that is
only 6 times slower than the previous algorithm. Our
reconstruction software runs at the South Pole in thedetector, and is included in all IceCube analyses.

References[1] IceCube Collaboration, IceCube webpage
http://icecube.wisc.edu/.
[2] A. Achterberg et al, Astroparticle Physics 26 (2006)
155173.
[3] R. Abbasi et al, Nuclear Instruments and Methods inPhysics Research Section A 618 (2010) 139-152.
[4] R. Abbasi et al, Nuclear Instruments and Methods inPhysics Research Section A 703 (2012) 190-198.
[5] R. Abbasi et al, Nuclear Instruments and Methods inPhysics Research Section A 601 (2009) 294-316.
[6] J. Ahrens et al, Nuclear Instruments and Methods inPhysics Research Section A 524 (2004) 169-194.
[7] IceCube Collaboration, paper 0565 these proceedings.
[8] M. Wolf, Verification of South Pole glacial ice
simulations in IceCube and its relation to conventionaland new, accelerated photon tracking techniques.
Masters thesis University Heidelberg September 2010.
[9] R Abbasi et al, Astronomy & Astrophysics 535 (2011)18.
[10] M. Ackermann, Searches for signals from cosmic
point-like sources of high energy neutrinos in 5 years of
AMANDA-II data PhD thesis Humboldt-Universita?t zu
Berlin (2006).
[11] S. Boyd, and L. Vandenberghe, (2009) Convex
Optimization, Cambridge University Press.
[12] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J.
Eckstein, Foundations and Trends in Machine Learning
3 (2011) 1-122.