# Reconstruction of Events Recorded with the Surface Detector of the Pierre
Auger Observatory

The Pierre Auger Collaboration

###### Abstract

Cosmic rays arriving at Earth collide with the upper parts of the atmosphere,
thereby inducing extensive air showers. When secondary particles from the
cascade arrive at the ground, they are measured by surface detector arrays. We
describe the methods applied to the measurements of the surface detector of
the Pierre Auger Observatory to reconstruct events with zenith angles less
than $60^{\circ}$ using the timing and signal information recorded using the
water-Cherenkov detector stations. In addition, we assess the accuracy of
these methods in reconstructing the arrival directions of the primary cosmic
ray particles and the sizes of the induced showers.

##  1 Introduction

The intensity of cosmic rays with energies above ${\sim}10^{14}$ eV is only a
few particles per square meter per day and thus too low to allow for direct
measurement with satisfactory statistical precision. The phenomenon of
extensive air-showers must therefore be exploited to study cosmic rays at
higher energies. In extensive air-showers, the interaction of a primary
cosmic-ray with Earth’s atmosphere induces a cascade of particles, where the
particles reaching the ground may be sampled by large arrays of detectors.

The higher the energies of primary cosmic rays, the lower is their arrival
rate such that increasingly larger detection areas are needed to counteract
the lower intensity. Meanwhile, the higher the primary energy, the larger are
the dimensions of the resulting air shower such that an increased spacing
between detectors is sufficient to ensure shower detection. For example, the
footprint at sea level of a shower generated by a primary cosmic-ray with an
energy of $10^{19}$ eV is ${\sim}10$ km2, and it has been found empirically
over the years that having several detectors inside such an area is sufficient
to reconstruct the properties of such showers. Above $10^{19}$ eV, the
intensity is only ${\sim}1$ per km2 per year, which therefore requires huge
collecting areas to gather adequate statistics.

A shower can be thought of as a thin, radially extended, and slightly curved
disk of particles propagating longitudinally at the speed of light along the
initial direction of the primary cosmic-ray. At every stage of development,
the particle density in the disk is largest at the center, or core, of the
shower and decreases with radial distance from the shower axis. At sea level,
about 90% of the energy flow lies within less than 100 m of the shower axis;
however, gathering a reasonable number of events here would require a
prohibitively dense array and a large dynamic range of the particle detectors.
A sparse array of detectors at the ground can only sample the disk at a fixed
stage of shower development and only beyond a certain distance from the core.
The arrival direction and a proxy of the shower size can nevertheless be
reconstructed and interpreted in terms of the properties of the initial cosmic
ray without the need for deeper insight into the complex physics of extensive
air-showers, where the latter quantity can be used to infer the energy of the
primary particle.

The arrival direction of the primary cosmic-ray can be determined from the
relative times of arrival of the shower disc at the dispersed detectors. This
was first demonstrated in the 1950s by the researchers in the MIT group [1],
who used their newly-developed liquid scintillation detectors to show, for the
first time, that the shower disc was sufficiently thin to enable determination
of the direction of the shower axis to within a few degrees. Since then, the
detectors of choice have either been scintillators (liquid or plastic) or
water-Cherenkov detectors, the use of which was pioneered a few years later
[2].

The number of particles in the shower (also referred to as the shower size) at
the observation level is related to the primary energy, but because of
fluctuations in shower development, the association is not a simple one. A
lower-energy primary particle that interacts deep in the atmosphere can create
a cascade which has the same size at the ground as the one produced by a more
energetic particle interacting higher in the atmosphere. Measuring the shower
size at the depth of the shower maximum (or better, the integral of the
deposited energy along the shower profile) reduces the impact of such
fluctuations; however, observations at the shower maximum have only been
possible at high altitudes and for relatively low-energy events (see e.g.
[3]), or by the observation of the profile with fluorescence detectors [4, 5],
albeit with an order of magnitude lower duty cycle.

When considering depths lower in the atmosphere than the position of the
shower maximum, and when using a widely-spaced array of detectors, the
situation is more complex. Determining the fall-off of the signal as a
function of distance from the axis is difficult when the spacing of the
detectors is large, and it cannot be determined with adequate precision on an
event-by-event basis. A technique to circumvent this problem was introduced by
Hillas [6], who proposed using the signal at a distance appropriate for the
separation of the detectors at which the uncertainty in the determination of
this signal is least susceptible to the lack of knowledge of the fall-off of
the signal with distance.

The Pierre Auger Observatory [7] is the world’s largest extensive air-shower
detector and makes use of these techniques to study cosmic rays at the highest
energies. The hybrid design of the Observatory consists of 27 telescopes of
the fluorescence detector (FD) [8] distributed between four sites overlooking
the huge array of the surface detector (SD). The SD, described in the
following section, provides the bulk of the events recorded by the
Observatory, which are principally used to measure the distribution of arrival
directions (e.g. [9]) and the energy spectrum (e.g. [10]) of the highest-
energy cosmic rays. The reconstruction of these SD events is the focus of this
paper.

In this document, we provide a detailed description of the methods used to
reconstruct the arrival direction and size of an extensive air-shower from the
magnitudes and temporal distributions of signals detected in each SD station.
A critical evaluation of the accuracy of the reconstruction procedure is also
an essential part of this work. The methodology described here is used only
for events with a zenith angle $\theta<60^{\circ}$, as the reconstruction of
more inclined showers requires different methods due to the substantial
asymmetry induced in the spatial distribution of the shower particles by the
geomagnetic field and due to geometrical effects [11].

Within the Pierre Auger Collaboration, two independent frameworks have been
developed for the reconstruction of these events with the corresponding
releases of reconstructed data dubbed _Herald_ and _Observer_ [12]. The
existence of these independent analyses facilitates easy access to cross-
checks, validations, and verification of systematic effects. The two
reconstructions are based on a similar strategy; however, they include slight
differences in the description of the temporal and lateral structure of
showers. Such differences are highlighted where relevant.

The paper is organized as follows. In Section 2, we give a brief introduction
to the features of the SD relevant for the reconstruction procedures. In
Section 3, we describe the treatment of the data from individual detectors,
including the derivation of the start-times and the magnitudes of the recorded
signals. In Section 4, we describe how shower events are built prior to the
reconstruction procedure. In turn, this reconstruction procedure is described
in detail in Section 5, beginning with the methods used to reconstruct the
geometry of the shower (arrival direction of the shower and the impact point
of the core) and ending with the estimation of the shower size. In Section 6
and in Section 7, the uncertainties associated with the reconstructed
parameters are discussed. In Section 8, we finalize the description of the
reconstruction procedures by briefly outlining the assignment of energy.

##  2 The surface detector

The SD consists of more than 1600 water-Cherenkov stations arranged on a
isometric triangular grid covering 3000 km2. Station elevations range between
1300 and 1600 m above sea level, and each station is separated from its six
nearest neighbors by approximately 1500 m. For example, a primary cosmic ray
with an energy of $10^{19}$ eV will trigger at least 5, but on average
${\sim}8$, stations. Within the boundary of the 1500 m array, there are two
denser regions: one extending over 28 km2 within which the stations are spaced
750 m apart from one another, and the other covering 2 km2 with 433 m spacing.
For events recorded in these parts of the SD, specific reconstructions are
developed based on the reconstruction given here, but appropriately scaled to
a smaller spacing. A detailed description of the SD can be found in [7], and
only the main features relevant to the work presented in this paper are
summarized here.

Each station is a cylindrical container with a base area of 10 m2. These
containers are filled to a depth of 1.2 m with highly purified water enclosed
by a diffusively-reflective liner. The volume of water is viewed from above by
three 9-inch photomultiplier tubes (PMTs), which detect Cherenkov light
emitted during the passage of charged particles. Each PMT provides two signals
which are tagged with the GPS time stamps to an absolute time accuracy of
${\sim}12$ ns [13] and are digitized by 40 MHz, 10-bit Flash Analog-to-Digital
Converters (FADCs). The _low-gain_ signal is taken directly from the anode of
the PMT, while the _high-gain_ signal is provided by the last dynode and
amplified to be nominally 32 times larger than the low gain, thus enlarging
the total dynamic range to span more than three orders of magnitude in
integrated signal.

The data acquisition is governed by a hierarchical system of triggers
implemented in hardware and software. The first two levels, T1 and T2, are
formed locally by each detector station [14]. The set of stations passing the
T2-level criteria are continuously examined remotely for temporal and spatial
coincidences in at least three stations. Upon identification of such a
coincidence, an array-level trigger T3 is issued by the central data
acquisition system (CDAS). After the T3 trigger is raised, FADC data from the
T1 and T2-triggered stations, as well as calibration and monitoring data, are
transmitted to CDAS by means of local-area-network radio-communication.

Figure 1: Example of a shower event detected by the SD, which is reconstructed
with $\theta=55.2^{\circ}$ and an energy of 38.7 EeV. _Left:_ Schematic
depiction of the surface detector array in which each dot represents a surface
detector station. The colored stations correspond to those participating in
this example event. Rays correspond to the field of view of the telescopes at
the four FD sites. _Right:_ A zoomed-in, top-down view of the example event
falling within the array. The projection of the shower axis on the ground is
represented by the black line, ending at the impact-point of the shower-core
(red square). The SD stations are colored according to their trigger time
(blue is early, green is late) with their area proportional to the logarithm
of the signal.

Additional higher-level triggers are then implemented off-line to select
shower events (physics trigger T4) and to ensure the shower is well contained
in the array. With respect to the latter, we employ a fiducial trigger $N$T5,
where $N$ is the number of functioning stations surrounding the station with
the highest signal at the time of the event. In this work we consider only the
most stringent condition, the 6T5, where all six neighboring stations must
have been present and functioning. This ensures an accurate reconstruction of
the impact point. For the data set used in this work, 6T5 events111 This
requirement not only ensures adequate sampling of the shower but also allows
evaluation of the aperture of the SD in a purely geometrical manner in the
regime where the array is fully efficient [14]. from the period of 1 January
2004 to 31 December 2018 have been selected, resulting in a set of more than
400 000 events. However, the results discussed in this paper also apply to the
events with a less conservative trigger 5T5, which have been used in
anisotropy studies where it is advantageous to counterbalance the slightly
reduced reconstruction accuracy with larger statistics [9, 15, 16]. The whole
acquisition chain, starting from the single-station trigger up to the off-line
event selection, reduces the counting rate of single stations from about
$3000/\text{s}$, due mainly to single uncorrelated atmospheric muons [17],
down to about $3{\times}10^{-5}/\text{s}$ for high quality physics events.
Consequently, a rate of events across the complete array of about
$0.05/\text{s}$ is measured, which is almost entirely due to extensive air
showers.

The footprint of one such event is shown in Fig. 1 and is used in the examples
that follow in later sections.

##  3 Reconstruction at the level of the station

For recorded events, the treatment of signals starts at the station level. In
this section, we describe the treatment of the FADC traces222 A few traces
(${<}3\%$) can in fact show anomalies due to hardware problems either in the
PMTs or in the electronics: these are identified prior to the analysis by
means of the monitoring data and excluded from the reconstruction. for each
PMT in stations passing the T1 or T2 trigger conditions. These procedures
provide a common standard to which all surface detector stations are
calibrated. Additionally, they provide a start time and a signal magnitude.
The uncertainties related to each of these are also discussed.

###  3.1 Calibration

Figure 2: The calibration of individual stations is performed by collecting 60
s of background signals. The histogram above is an example of the sizes of
integrated signals $Q$ in units of ADCT = ADC count $\times$ time bin. The
position of the peak (local maximum) at $Q={\sim}232$ ADCT, created by the
atmospheric background muons, is used as the reference value for conversion of
station signals into units of VEM.

The aim of the calibration is to enable a conversion of the signals recorded
by the three PMTs and electronics of each station from the unit of ADC-counts
to the reference unit of a vertical equivalent muon (VEM). The VEM is the
charge associated with a vertical muon passing through the center of the
station and also corresponds to ${\sim}240$ MeV of energy deposited due to
ionization losses of the passing muon. [18]. This conversion factor,
$Q^{\text{peak}}_{\text{VEM}}$, is evaluated online by constantly measuring
the properties of particles entering the detector volume at a rate of
${\sim}3/\text{s}$. The on-board SD electronics are used to record histograms
of the digitized pulse-heights and integrated charges from a minute’s worth of
background particles. An example of such a charge histogram is shown in Fig.
2. These histograms are stored alongside and exported with event data, thereby
enabling the off-line re-evaluation of the VEM unit with 3% accuracy.

The calibration factor $Q^{\text{peak}}_{\text{VEM}}$ is the position of the
local maximum corresponding to the most probable charge deposited by muons
coming from all directions, which is proportional to the VEM.

In each event, the FADC traces are independently calibrated for each PMT using
the background data recorded in the preceding minute. Dynamic station-by-
station conversion of all signals into units of VEM facilitates preservation
of a stable and uniform response to shower particles across the array.

Figure 3: _Top:_ Illustration of the baseline estimation (red), constructed
from the individual segments (highlighted in green) and interpolation, on a
high-gain FADC trace of a PMT. _Bottom:_ Identification of signal fragments
(magenta filled boxes) and the final merged signal window (large dotted cyan
box). Based on the number of peaks, this trace probably contains ${\sim}5$
muons with some additional electromagnetic component. Trace start and stop
times are denoted with the vertical dashed red and blue lines, respectively.

###  3.2 Determination of the baseline

The first step in the signal processing is the estimation of the baseline of
the high-gain and low-gain traces of each PMT separately. Since the baseline
of a trace can be affected by the signals themselves, it is not estimated as
one single value for the whole trace. Instead, the algorithm allows for a
changing baseline along the trace. Constant baseline segments are identified
as at least 40 consecutive bins (i.e. with a duration of at least 1 $\upmu$s)
that vary in amplitude by less than $\bar{y}\pm\sigma_{\text{b}}$, where
$\bar{y}$ is the mean of the baseline segment and $\sigma_{\text{b}}$ is the
tolerance needed to account for the possible noise in the baseline.
$\sigma_{\text{b}}$ is set to 2 ADC counts. In rare cases (${<}0.2$%), when
the PMT in question is noisy, the tolerance $\sigma_{\text{b}}$ is increased
in steps of 1 until a constant baseline segment is found. These baseline
segments, illustrated by the green boxes in Fig. 3-top are then combined to
form the total baseline. Each of the bins corresponding to a baseline segment
is inserted into the total baseline as a constant value given by $\bar{y}$ of
that baseline segment. Then, the parts between the baseline segments are
filled with values interpolated linearly in integrated signal, which results
in an estimation of the baseline covering all of the 768 bins constituting the
trace as shown with the red line in Fig. 3.

###  3.3 Identification of the time range of interest

The next step is the extraction of the relevant signal from the traces of
individual PMTs, i.e. identifying the start and stop time of the signal, which
is always performed on the high-gain channel due to its superior resolution.
As in the baseline case, the FADC traces are first scanned to identify
candidate signal _fragments_ , which consist of consecutive bins with
amplitudes of at least 3 ADC counts above the baseline. Examples of identified
signal fragments are indicated in Fig. 3-bottom with magenta boxes. After the
subtraction of the baseline in each PMT, traces are calibrated in VEM units
using the value of $Q^{\text{peak}}_{\text{VEM}}$.

To avoid inclusion of background signals (mostly single, unaccompanied muons),
the candidate fragments are additionally scrutinized before being merged into
the final signal _segment_. Two consecutive fragments A and B can only be
merged if they fulfill both of the following requirements:

  1. 1.

Fragment B must begin less than 500 ns (20 time bins) + $t_{\text{A}}$ after
the end of fragment A, where $t_{\text{A}}$ is the amount of time in fragment
A containing signal.

  2. 2.

The integral of fragment A must be greater than 30% of the integral of
fragment B, or fragment B must have an amplitude below
$5I^{\text{peak}}_{\text{VEM}}$, where $I^{\text{peak}}_{\text{VEM}}$ is the
mean amplitude of the background muons [18].

This procedure is repeated for all consecutive fragments until all fragments
are processed and signal segments are created. In Fig. 3-bottom, the final
signal segment resulting from the merging of the candidate fragments (magenta)
is indicated with a dotted cyan box.

The last step of the procedure consists in the averaging of the signal
segments of the three PMTs (or of all working PMTs333 Currently 12% (1%) of
stations have only two (one) functioning PMTs.) into the station-level signal
segment.

For the final station-level trace, only the segment with the largest signal is
selected and the corresponding start and stop times are adopted with a minimal
run of 250 ns to ensure containment of any potential sub-threshold signals.

With this start time, we obtain our best estimate of the beginning of the
passing shower front. With this stop time, we ensure that all particles
belonging to the shower in question are included while excluding as many
accidental signals as possible.

###  3.4 Estimation of the signal

The total signal of a station is obtained by integrating the final trace,
which consists of the bin-by-bin average of the high-gain (or low-gain, if the
high gain is saturated) traces of the working PMTs between the previously
determined start and stop times.

When the flux of particles entering a station is large (e.g. when a shower
lands close to a station), the signal can saturate the ADC range of the high-
gain channel. This normally occurs when the equivalent integrated signal is
greater than ${\sim}60$ VEM. When more than two FADC-trace bins reach the
maximum ADC count ($2^{10}-1=1023$), the final station-level trace is
constructed from the average of the low-gain PMT-traces instead, while the
start and stop times of the high-gain signal are kept.

When the core lands even closer to a station, the low-gain channels can also
saturate and the corresponding events are termed _saturated events_. For lower
energy showers, this occurs in 10% of events, for which the distance of this
station is closer than ${\sim}100$ m to the shower axis. For the highest
energy showers, saturation already occurs when a station is closer than
${\sim}500$ m, which occurs in more than 50% of events (see Section 5.2.2 for
more details). In the case of low-gain saturation, the equivalent integrated
signal is typically greater than ${\sim}1000$ VEM. Events with more than one
saturated station are extremely rare444 Out of the more than 400 000 events,
we have identified only three where two stations have saturated..

With the increasing signal, we first observe an overflow of the ADC range of
the low-gain channels and then, with even larger signals, we observe
saturation of the PMT itself (i.e. the departure from the linear behavior of
the photomultiplication). On such traces, a signal-recovery algorithm is
applied [19] to estimate the true magnitude of the low-gain signal from the
amplitude of the signal undershoot (the dip of the trace below baseline). In
addition, in a more detailed study it has been shown that the non-linear
response of the PMT, when measured precisely, can be used to recover the
saturated signals by fitting the non-saturated parts of the trace. For all
PMTs for which the saturation and non-linearity response has been measured
individually, the range of signal estimation can be increased up to $10^{6}$
VEM, with a resolution of 20% [19]. The latter recovery procedure is
implemented only in the _Observer_ framework.

When the high-gain channel is used to produce the calibrated station trace,
the high-gain signal is scaled with the value of
$Q^{\text{peak}}_{\text{VEM}}$ since this is the channel with which muon data
are continuously monitored by the online station calibration. In cases where
the high-gain channel saturates, the low-gain channel is used to create the
calibrated station trace. In this case, the low-gain signal is scaled with the
product of the $Q^{\text{peak}}_{\text{VEM}}$ factor and the high-to-low-gain
ratio, which is also constantly monitored by the online station calibration.

###  3.5 Uncertainties in the measurement of the signal

Figure 4: Comparison of the signal distributions of the two reconstructions.
_Left:_ Below 1000 VEM, above which the low-gain signal of some stations
begins to saturate, the distribution can be fitted with the following power
law
$S_{\text{Herald}}/\text{VEM}=1.02(S_{\text{Observer}}/\text{VEM})^{0.98}$.
The saturated stations are represented by the orange points. _Right:_
Distribution of relative signal difference $\Delta S/S_{\text{Observer}}$,
where $\Delta S=S_{\text{Herald}}-S_{\text{Observer}}$, for different ranges
of signal. For all ranges of signals, the average of the distribution is below
1%, and the standard deviation ranges from less than 2% to 5% for the full
data set.

The uncertainties on the integrated signal and on the start time of the signal
are obtained from studies using configurations with _twin_ stations, i.e.
pairs of stations placed only ${\sim}11$ m apart and thus sampling the same
part of a shower. In the studies laid out below we used all available data
recorded by twins. Two twins were already deployed in 2003 and remain
operational today. In 2006, 20 additional twins were added. Of these, 13
collected data for nearly 5 years, while a remaining hexagon of seven twins
remains operational.

The uncertainty in the start time strongly influences the angular resolution
of the reconstruction of the shower geometry. This uncertainty is composed of
the time resolution of the GPS time-tagging system (${\sim}10$ ns), the 40 MHz
FADC sampling (${\sim}7$ ns), and, more importantly, by the fluctuations of
the arrival times of the first particles. These depend on the number of
particles entering the station, i.e. a convolution between the thickness of
the shower front, the density of particles therein, and the cross-section of
the detector. At large radial distances only a few particles will enter a
detector and thus large start-time fluctuations can occur due to this sparse
sampling of the shower front. We can describe the effects above with a general
time variance model [20] for which the uncertainty of the start time
$t_{\text{start}}$ is estimated as

$\sigma_{t_{\text{start}}}^{2}=a^{2}\left(\frac{2t_{50}}{n}\right)^{2}\frac{n-1}{n+1}+b^{2},$
(3.1)

where $t_{50}$ is the length of the time interval that contains the first 50%
of the total signal, and the effective number of particles $n$ is obtained
from the signal $S$ as $n=S/\bar{\ell}(\theta)$, where $\bar{\ell}(\theta)$ is
the mean relative track-length of through-going particles crossing the
detector when arriving with a zenith angle $\theta$. Due to the cylindrical
shape of the water-Cherenkov tank, $1/\bar{\ell}(\theta)=\cos\theta+(2H/\pi
R)\sin\theta$, where $H$ and $R$ are the height and the radius of the
cylinder, respectively. The coefficients $a$ and $b$ (with typical magnitudes
of 0.7 and 16 ns, respectively) are parameterized as functions of the zenith
angle of the shower and were determined using twin stations.

For the signal uncertainty of non-saturated stations, a Poisson-like
parameterization,

$\sigma^{2}_{S}=f^{2}_{S}(\theta)\,S,\quad\text{where}\quad
f_{S}(\theta)=0.34+0.46\,\sec\theta,$ (3.2)

is a good description of the observed signal fluctuations. The values from
[21] have been updated as a result of studies with larger statistics.

However, for stations with saturated PMTs, the accuracy of the signal
measurement is degraded. In the _Observer_ framework, the signal uncertainties
are increased by the intrinsic uncertainties of the saturation-recovery
algorithm. The signal uncertainty for the recovered signals varies from 5% at
${\sim}1000$ VEM up to 60% at ${\sim}10\,000$ VEM. For very large signals,
i.e. above 10 000 VEM, the PMT response is highly non-linear, and only few
measurements of the saturation curve at large anode currents are available.
The algorithm uses average curves, which results in uncertainties greater than
50% on the recovered signal. In the _Herald_ framework, the algorithm for
signal recovery does not take into account the PMT non-linearity and therefore
only the timing information is used for signals larger than 2000 VEM.

For both the determination of the timing and the estimation of the signal and
its uncertainties, the two reconstruction frameworks exhibit some differences.
Nevertheless, this has only a limited impact on the final estimate of the
signal size. The differences between the signals are on average less than 1%
with a standard deviation varying from 2% to 4.9%. As can be seen in Fig.
4-right, they are mostly observed at the lowest signals, and are within the
estimated uncertainties. For the saturated signals, the large differences
observed in Fig. 4-left are explained by the additional recovery procedure
used in the _Observer_ framework. The impact of saturation on the event
reconstruction is discussed later in this paper.

##  4 Building events

Prior to reconstruction, each event is built from the data gathered by the T3
trigger. The sizes and start times of the signals of the triggered stations
are inspected to identify candidate showers. Details of the selection process
can be found in [14]. Here we summarize only the elements most relevant for
the reconstruction.

Since the highly sensitive SD stations may be triggered due to lightning
during thunderstorms, affected events must be identified and removed. This is
done using a Fourier-like algorithm to detect oscillations in traces, the
presence of which indicates that lightning has struck nearby. If even a single
PMT of any station participating in the event is found in this condition, the
whole event is discarded.

Next, a selection process is applied to the triggered stations of candidate
shower events to discard accidental stations, i.e. stations which triggered by
chance and are not, in fact, part of the event. To perform such a selection,
an estimate of the time the shower front was expected to have passed through
each station is required. This is obtained by means of a first order seed
reconstruction, which, in addition to facilitating the identification and
exclusion of accidental stations, also provides the initial estimates for the
geometric parameters of the proper event reconstruction detailed in Section 5.

The seed reconstruction is obtained as follows. The to-be-determined axis
$\hat{a}$ of the shower is anchored at the location of a barycenter
$\vec{x}_{\text{b}}$ (i.e. the signal-weighted555 Monte Carlo studies have
shown that weights proportional to the square-root of the signal in the
stations result in a barycenter on average closest to the impact point of the
shower core on the ground. center-of-mass of stations in an event), which
later666 Note that at this stage the axis $\hat{a}$ can be obtained without
any reference to the shower core. serves as a first estimate of the impact
position of the shower core at the ground. The assumed shower core travels in
the $-\hat{a}$ direction, impacting the ground at barycenter
$\vec{x}_{\text{b}}$ at the time $t_{\text{b}}$. Under the assumption that
particles in the shower front move in a plane perpendicular to the shower axis
with the same speed as the core of the shower (see Fig. 5-left), which is
assumed to be the speed of light $c$, the time $t_{\text{sh}}(\vec{x})$ when
the shower plane passes through some arbitrary point $\vec{x}$ (e.g. a station
on the ground) may be inferred through a simple projection onto the shower
axis as

$c\,t_{\text{sh}}(\vec{x})=c\,t_{\text{b}}-\hat{a}\cdot(\vec{x}-\vec{x}_{\text{b}}).$
(4.1)

The solution to this equation is obtained by identifying a reconstruction
_seed triangle_ consisting of three stations: a station with its two nearest-
neighbors in a non-aligned configuration, where all three stations must have
passed a station level trigger. The seed triangle with the highest sum of the
three station signals is used to analytically determine an exact solution,
thereby obtaining an estimate of the shower axis and an estimation of when the
shower front arrived at any position under the plane-front assumption.

Figure 5: _Left:_ Schematic description of the plane-front approximation for
air-shower development (see text). _Right:_ Schematic of the spherical shower-
front development. The apparent origin of the shower is at point
$\vec{x}_{\text{o}}$ at time $t_{\text{o}}$. The shower front reaches the
station $i$ at position $\vec{x}_{i}$ at time $t_{i}$.

We designate as accidental all stations where the onset of the signal is more
than $1\,\upmu$s prior to (early) or $2\,\upmu$s after the expected passage of
the plane shower front (late), which ensures that more than 99% of the
stations belonging to the shower are kept in the event. The criteria is more
strict for early stations since the seed triangle, which defines the shower
front estimate, is based on stations close to the shower axis with the highest
signals, which are more likely to attain signal already at the start of the
passage of the shower front. Stations more distant from the shower axis
measure smaller signals whose start times are biased to later times due to the
larger sampling fluctuations in the arrival time of the first particles as
well as due to the curvature of the shower front resulting in increasing delay
with respect to the plane shower front. Additionally, stations are removed
from an event when there are no triggered stations within 1.8 km (or 5 km) for
a single station (or two neighboring stations). Furthermore, active stations
which did not record enough signal to trigger and are located up to 5 km from
the shower core are considered as part of the event. They provide additional
constraints on the location of the impact point, as will be discussed in the
next section.

##  5 Reconstruction of events

The timing and the size of the signal measured in each station, as well as
accurate knowledge of the 3D positions of the stations, are the key inputs for
the reconstruction of the arrival directions and the sizes of the showers
selected according to the criteria described in the Section 4. To reconstruct
these quantities, we adopt a simplified model of air showers that allows us to
separate the process into two parts. First, from the timing information of the
stations, we determine the geometry of the shower, namely the direction of the
shower axis $\hat{a}$ and the position of the impact point of the shower core
on the ground $\vec{x}_{\text{c}}$. Using this geometry, the second step
consists of fitting the signal magnitudes as a function of distances from the
shower axis to an empirically-derived functional form describing the average
lateral distribution of particles. These two steps can be formally intertwined
and thus require iteration or a global fit, but can also be effectively
separated with a particular choice of the shower front.

###  5.1 Shower geometry

From the selection of the stations described in Section 4, the three stations
of the seed triangle are used to give the first rough estimate of the axis
$\hat{a}$ of the shower and the impact position $\vec{x}_{\text{c}}$ of the
core on the ground. The shower models used here describe the secondary
particles that “have traveled furthest in the forward direction”, also
referred to as the _extreme front_ [22], as moving with the speed of light in
a curved shower front. The radius of curvature is added as a free parameter
when five or more stations are participating in the event. Moreover, for the
reconstruction of the shower geometry, we consider the core to be moving in
the $-\hat{a}$ direction and intersecting with the ground at the impact point
of the shower core $\vec{x}_{\text{c}}$ at the time $t_{\text{c}}$. Thus four
parameters are fitted: the two directional cosines777 Here we try to avoid
fitting the angles $\theta$ and $\phi$ of the spherical coordinate system due
to the singularity at the zenith $\theta=0$. Nevertheless, they are used in
the fits of more inclined events with zenith angles $\theta>60^{\circ}$, c.f.
[11]. $u$ and $v$ of the unit vector $\hat{a}=(u,v,\sqrt{1-u^{2}-v^{2}})$, the
time $t_{\text{c}}$, and the radius of curvature of the shower front
$R_{\text{o}}$. Using this model for the arrival time of the shower front,
$t_{\text{sh}}(\vec{x})$, the shower geometry is fitted to the start-time of
the signals $t_{i}$ in each triggered station $i$ located at $\vec{x}_{i}$
using the Minuit framework [23]. The function which is minimized is the sum of
the squares of the differences in predicted and measured start times,

$\chi^{2}=\sum_{i}\frac{[t_{i}-t_{\text{sh}}(\vec{x}_{i})]^{2}}{\sigma^{2}_{t_{i}}},$
(5.1)

where $\sigma_{t_{i}}$ is the uncertainty of the start time $t_{i}$, as given
by Eq. 3.1, which derives, in part, from fluctuations in the arrival times of
the most forward particles in the shower front at a station $i$. The absolute
positions of the stations $\vec{x}_{i}$ were measured with the built-in GPS
receivers in dedicated campaigns to a precision of around 20 cm horizontally
and 50 cm vertically so that the positional uncertainty has only a negligible
effect on the event reconstruction.

As initial values for the fitting procedure, the signal-weighted barycenter
$\vec{x}_{\text{b}}$ from Section 4 is chosen as a suitable approximation for
the impact point $\vec{x}_{\text{c}}$, while the result of the plane fit from
Eq. 4.1 is used for the shower axis $\hat{a}$.

The two reconstruction frameworks include slightly different descriptions of
the shower front curvature. In the _Herald_ framework, the curvature is
considered constant, with the particles of the shower are distributed on a
spherical front moving along the direction of $-\hat{a}$. Compared to the
plane front in Eq. 4.1, the particles are thus delayed proportionally to
$R_{\text{o}}-\sqrt{R_{\text{o}}^{2}-[r_{\hat{a}}(\vec{x}-\vec{x}_{\text{c}})]^{2}}$,
where $R_{\text{o}}$ is the constant radius of the spherical front and where

$r_{\hat{a}}(\vec{x})=|\hat{a}\times\vec{x}|$ (5.2)

is the perpendicular distance of point $\vec{x}$ from the shower axis
$\hat{a}$. To keep $t_{\text{sh}}(\vec{x})$ linear in the fitted parameters,
only terms up to the second order in the $r/R_{\text{o}}$ expansion are used.
With a curvature parameter $k_{\text{o}}=1/2R_{\text{o}}$, we can then express
the shower timing as

$c\,t_{\text{sh}}(\vec{x})=c\,t_{\text{c}}-\vec{a}\cdot(\vec{x}-\vec{x}_{\text{c}})+k_{\text{o}}\,[r_{\hat{a}}(\vec{x}-\vec{x}_{\text{c}})]^{2},$
(5.3)

which is clearly a _paraboloidal_ extension888 The size of the quadratic term
is typically on the order of ${\sim}500$ ns at distances of ${\sim}2500$ m. of
Eq. 4.1.

The _Observer_ reconstruction approximates the shower development as starting
at time $t_{\text{o}}$ from a virtual point of origin $\vec{x}_{\text{o}}$
(see Fig. 5-right) and propagating towards the ground in the shape of a
_sphere_ , concentrically inflating with the speed of light999 The typical
difference in time between the paraboloidal and the spherical description of
the shower plane is of the order of 50 ns at distances of ${\sim}1500$ m.. The
arrival time of such a shower front at a point $\vec{x}$ is thus

$c\,t_{\text{sh}}(\vec{x})=c\,t_{\text{o}}+|\vec{x}-\vec{x}_{\text{o}}|.$
(5.4)

Note that contrary to the paraboloidal description in Eq. 5.3, this spherical
fit can be performed without any prior knowledge of the impact point
$\vec{x}_{\text{c}}$ or the shower axis. Once $\vec{x}_{\text{c}}$ is
determined at a later point, the shower axis may be obtained as a normalized
direction towards $\vec{x}_{\text{o}}$ as
$\hat{a}=(\vec{x}_{\text{o}}-\vec{x}_{\text{c}})/|\vec{x}_{\text{o}}-\vec{x}_{\text{c}}|$.
Due to the expansion of the sphere used in this model, the radius of curvature
of the shower front depends on time, $R(t)=c(t-t_{\text{o}})$. Nevertheless,
for consistency with Eq. 5.3, the radius of curvature is defined as the
distance between the virtual origin and the impact point,
$R_{\text{o}}=|\vec{x}_{\text{o}}-\vec{x}_{\text{c}}|$, i.e. as the radius of
the shower at the time of passage $t_{\text{c}}=t_{\text{o}}+R_{\text{o}}/c$
through the impact point $\vec{x}_{\text{c}}$.

For both of the models above, there are four free parameters to describe the
development of the shower front. Since low-energy events have a station
multiplicity of only three or four, they do not have enough degrees of freedom
in the timing data to solve for the shower-front curvature. For events with
less than five triggered stations, we therefore use the curved model with an
$R_{\text{o}}$ fixed to a parametrization optimized using events with a larger
number of stations.

As an example, a geometrical fit of the example event from Fig. 1 is shown in
Fig. 6-left for both reconstructions. The existence of curvature in the shower
front can be clearly observed in the increasing delays of signal start-times
with respect to the arrival of a plane front tangential at the shower core.

The angular differences between the axes of the two reconstructions are shown
in Fig. 6-right, from which we can conclude that 68%, 90%, and 95% of all the
events are reconstructed within an angular difference of less than
$0.40^{\circ}$, $0.83^{\circ}$, and $1.14^{\circ}$, respectively. The small,
non-zero difference of less than $0.1^{\circ}$ observed in $\Delta\theta$ will
be addressed in Section 5.2.1.

Figure 6: _Left:_ For the large example event from Fig. 1, the start times in
individual stations are shown as a function of the distance to the shower axis
for the _Herald_ (red triangles) and _Observer_ (blue circles)
reconstructions. The times are relative to the arrival time of the
corresponding plane front. The solid lines represent the predicted time for
the spherical model (Eq. 5.4, blue) and the parabolic model (Eq. 5.3, red).
_Right:_ Distribution of the angular difference $\eta$, where
$\sin\eta=|\hat{a}_{\text{p}}\times\hat{a}_{\text{s}}|$, between the axes of
the two curved shower-front models: $\hat{a}_{\text{p}}$ for the paraboloid
and $\hat{a}_{\text{s}}$ for the spherical model. The distribution has been
fitted with a weighted sum of two Rayleigh distributions from which we can
derive that 68% of reconstructed events have arrival directions which agree to
within $0.40^{\circ}$. _Inset:_ 2D distribution of the differences in the
zenith angle, $\Delta\theta=\theta_{\text{s}}-\theta_{\text{p}}$ and in the
azimuthal angle, $\Delta\varphi=\varphi_{\text{s}}-\varphi_{\text{p}}$.

###  5.2 Shower size

The shower size is estimated with a regression of the parameters of the
lateral distribution function (LDF), $S(r)$, to the signals in the triggered
stations and is additionally constrained by the absence of triggers in
stations which would have measured very little or no signal. The LDF is a
function of the perpendicular distance to the shower axis, $r$. However,
showers induced by identical primaries (i.e. with the same energy, mass, and
arrival direction) can be sampled at the ground at different stages of
development due to the shower-to-shower fluctuations arising from the
variability of the location and the nature of the leading interactions, and
the variability of the subsequent development of the cascade. This results in
a natural variability of the shower size, the mean value of which depends on
the mass of the primary particle.

For the majority of events, the multiplicity and spatial distribution of
triggered stations in the shower plane is insufficient to precisely estimate
the shower size and determine the shape of the LDF. Instead, station signals
of individual events are fitted with a scaled, data-derived _average_ LDF
shape $f_{\text{LDF}}(r)$ such that

$S(r)=S(r_{\text{opt}})\,f_{\text{LDF}}(r),$ (5.5)

where $S(r_{\text{opt}})$ is the shower-size estimator and $f_{\text{LDF}}(r)$
is normalized so that $f_{\text{LDF}}(r_{\text{opt}})\equiv 1$. The optimal
value for distance $r_{\text{opt}}$ has been chosen so that the variability of
the shower-size estimator $S(r_{\text{opt}})$ with respect to the
aforementioned shower-to-shower fluctuations is minimized. The value almost
exclusively101010 $r_{\text{opt}}$ depends only at ${\sim}5\%$ level on the
choice of the LDF form, the primary energy, and the zenith angle. depends on
the spacing and structure of the array [24]. In the case of the Auger SD,
where the array is an isometric triangular grid with a spacing of $1500$ m,
the optimal distance amounts to $r_{\text{opt}}\approx 1000$ m, so we have
decided to fix it at exactly 1000 m and denote the corresponding shorthand for
the shower-size estimator as $S(1000)$.

####  5.2.1 Lateral distribution function

Due to the lack of analytical solutions for the hadronic-cascade equations,
functional forms for $f_{\text{LDF}}(r)$ have traditionally been chosen
empirically (see [25] for an overview of functional forms chosen for previous
experiments). The _Herald_ reconstruction uses a log-log parabola,

$\ln f_{\text{LDF}}(r)=\beta\,\rho+\gamma\,\rho^{2},$ (5.6)

where $\rho=\ln(r/r_{\text{opt}})$ and $\beta$, $\gamma$ are the average slope
parameters of the LDF. For smaller axial distances, $r<r_{\text{c}}$ where
$r_{\text{c}}=300$ m, a tangential log-log linear function
$\beta\,\rho+\gamma(2\rho-\rho_{\text{c}})\rho_{\text{c}}$ is used where
$\rho_{\text{c}}=\ln(r_{\text{c}}/r_{\text{opt}})$.

The _Observer_ reconstruction instead uses a slightly modified NKG function
[26, 27, 28] of the form

$f_{\text{LDF}}(r)=\left(\frac{r}{r_{\text{opt}}}\right)^{\beta}\left(\frac{r+r_{\text{s}}}{r_{\text{opt}}+r_{\text{s}}}\right)^{\beta+\gamma},$
(5.7)

with a fixed $r_{\text{s}}=700$ m.

The average of the slope parameter $\beta$ is determined in both cases in a
data-driven way, by fitting $\beta$ for the subset of events with a
multiplicity and spatial distribution of stations providing a sufficient
_lever arm_. In these events, the slope $\beta$ is fitted given that there are
at least 2, 3, or more stations within a radius of $400<r/\text{m}<1600$ and
that at least two of the stations are separated by 900, 800, or 700 m,
respectively. This lever-arm criteria ensures that the lateral distribution is
sufficiently sampled around the distance of 1000 m to constrain the slope of
the LDF well.

Due to the strong correlation between $\beta$ and $\gamma$, it is more
difficult to determine the average of the slope parameter $\gamma$ in a data-
driven way in the case of both LDF functions. In Eqs. 5.6 and 5.7, the
parameter $\gamma$ is describing the deviation of the LDF from a simple power-
law function at large distances. In the _Observer_ reconstruction, similar
lever-arm criteria are used to fit $\gamma$ except the station search radius
is changed to $1000<r/\text{m}<2000$. In the _Herald_ reconstruction, the
parameter $\gamma$ has been determined from Monte-Carlo studies. Finally, both
parameters $\beta$ and $\gamma$ are parametrized as functions of the zenith
angle and $S(1000)$.

Both LDF models used within the Auger Collaboration assume that the deposited
signals in the stations are rotationally symmetric around the shower axis. In
truth, the showers are asymmetric due to a combination of the longitudinal
evolution and geometrical effects related to the angles of incidence of the
particles on the stations. The _Herald_ reconstruction corrects for these
asymmetries with an unbiasing model applied to the measured signals

$S_{\text{symm}}(r)=\frac{S(r)}{1+\alpha(r)\cos\zeta},$ (5.8)

where $\zeta$ is the azimuthal angle of the station (as measured in the shower
plane where $\zeta=0$ corresponds to the upstream direction) and where the
radial dependence of the amplitude $\alpha(r)$ has been obtained from
simulations. These corrections are only important in very specialized studies
as they only result in small shifts of ${\sim}40$ m in the estimate of the
position of the shower core, which is well within the typical uncertainties of
the reconstruction. These corrections ultimately have a ${<}0.2\%$ effect on
the value of $S(1000)$ and a ${<}0.2^{\circ}$ effect on the axial direction,
as estimated by performing the _Herald_ reconstruction with and without their
application. The exact magnitudes of these effects depend on the zenith angle
and state of development of the shower. Since no asymmetry correction is
applied in _Observer_ , the corresponding shift of the impact-point is also
responsible for the very small, yet systematic differences in arrival
directions reconstructed by the two reconstructions, where
$\Delta\theta<0.1^{\circ}$ as seen in Fig. 6-right.

####  5.2.2 Maximum-likelihood procedure

Using the results of the reconstruction of the shower geometry, the fit to the
LDF is, in the next step, maximized using the likelihood $\mathcal{L}$,
composed as the product of probabilities $P$ over the coordinates of the
shower impact point $\vec{x}_{\text{c}}$ and the size $S(1000)$, given the
observed signal sizes $S_{i}$ in the stations, located at $\vec{x}_{i}$. We
are thus maximizing the log-likelihood

$\ln\mathcal{L}=\sum_{i}\ln P(S(1000),\vec{x}_{c}\,|\,S_{i},\vec{x}_{i}).$
(5.9)

The two reconstructions use different models for calculation of the likelihood
$P$. In the _Observer_ framework, $P$ is a product of contributions from: (a)
stations with _small_ signals, (b) stations with _large_ signals, (c) stations
with _saturated_ signals, and (d) _non-triggered_ stations. The individual
terms are as follows:

  * (a)

_Small signals_. For a small signal $S_{i}<S_{\text{stat}}$ where
$S_{\text{stat}}\approx 20$ VEM, the likelihood is modeled as a Poissonian
probability $f_{\text{Poi}}(S_{i},S(r_{i}))$, with the signal scaled to a
particle count for which the corresponding Poisson variance matches the signal
variance from Eq. 3.2, and where $S(r_{i})$ is, according to the LDF of Eq.
5.5, the expected signal at the radial distance
$r_{i}=r_{\hat{a}}(\vec{x}_{i}-\vec{x}_{\text{c}})$ from Eq. 5.2.

  * (b)

_Large signals_. For a large number of particles the distribution of the
signal is modeled with a Gaussian approximation
$\mathcal{N}(S_{i};S(r_{i}),\sigma_{i})$ which is used for large signals
$S_{i}\geqslant S_{\text{stat}}$ and where $\sigma_{i}=\sigma_{S(r_{i})}$
according to the Eq. 3.2.

  * (c)

_Unrecoverable saturated signals_. For saturated signals that fail the signal
recovery procedure, mentioned in Section 3.4, $S_{i}$ is treated only as a
lower limit to the actual size of the signal. The Gaussian function
$\mathcal{N}$ for large signals is integrated over all possible values larger
than $S_{i}$ to get an estimate $F_{\text{sat}}(S_{i},S(r_{i}))$ of the
probability of detecting a signal larger than $S_{i}$, but only when no signal
recovery is applied.

  * (d)

_Non-triggered stations_. To trigger, a station has to acquire a certain
amount of signal, a process that can be modeled with a trigger probability
$p_{\text{trig}}(S(r_{i}))$, which is derived from data [14]. The
complementary probability for non-triggered stations is thus modeled as
$F_{\text{non}}(S(r_{i}))=1-p_{\text{trig}}(S(r_{i}))$.

This results in a likelihood function

$\mathcal{L}=\prod_{i}^{N_{\text{small}}}f_{\text{Poi}}(S_{i},S(r_{i}))\,\prod_{i}^{N_{\text{large}}}\mathcal{N}(S_{i};S(r_{i}),\sigma_{i})\,\prod_{i}^{N_{\text{sat}}}F_{\text{sat}}(S_{i},S(r_{i}))\,\prod_{i}^{N_{\text{non}}}F_{\text{non}}(S(r_{i})).$
(5.10)

For the _Herald_ reconstruction, the LDF likelihood function is split into two
parts based on whether stations have triggered or not. While for the triggered
stations, the signal $S_{i}$ is assumed to be a Gaussian fluctuation around an
average value $S(r_{i})$ given by the Eq. 5.6, for the non-triggered stations,
the relevant likelihood term is given, as in _Observer_ , by the complementary
probability. The _Herald_ reconstruction maximizes a likelihood function

$\mathcal{L}=\prod\limits_{i}^{N_{\text{sig}}}\mathcal{N}(S_{i};S(r_{i}),\sigma_{i})\,\prod\limits_{i}^{N_{\text{non}}}F_{\text{non}}(S(r_{i})).$
(5.11)

For showers with triggered stations satisfying the lever-arm conditions and
thus ensuring a strong constraint on the slope of the LDF in Eq. 5.6, the
parameter $\beta$ is always fitted in the _Herald_ reconstruction.

Figure 7: _Left_ : Relative difference of the shower size due to the treatment
of $\beta$ as a free or fixed parameter, obtained with the _Observer_
reconstruction where $\Delta
S(1000)=S(1000)_{\text{free}}-S(1000)_{\text{fix}}$. The relative difference
is shown for all events in addition to events where a station saturated.
_Right_ : Fraction of saturated events as a function of the shower size.

In the _Observer_ reconstruction, the parameter $\beta$ can be set as a free
parameter for the same set of showers. A comparison of the shower sizes
obtained with $\beta$ as a free or fixed parameter of the LDF in Eq. 5.7 is
shown in Fig. 7-left. The impact of freeing or fixing $\beta$ amounts to a
less than 3% (4%) bias in the reconstructed value of $S(1000)$ for non-
saturated (saturated) events and, therefore, fitting of $\beta$ is not
required to obtain an unbiased estimate of the shower size [6, 29]. The
fraction of all events which are saturated is ${\sim}$11% but strongly depends
on shower size as depicted in Fig. 7-right. The fraction is constant to within
2% across the full range of zenith angles.

Using the _Herald_ reconstruction, one final, global reconstruction step is
performed to account for the coupling between the geometry and the LDF fits.
Depending on the number of degrees of freedom, any parameters which were
previously fixed are at this stage allowed to vary. The global log-likelihood
function

$-2\ln\mathcal{L}_{\text{global}}=\chi^{2}-2\ln\mathcal{L}$ (5.12)

is minimized, where the two terms are given in Eqs. 5.1 and 5.11.

Figure 8: _Left:_ The lateral distribution of signals from the example event
in Fig. 1 is plotted as a function of the radial distance of stations from the
shower axis, shown for both _Herald_ (red triangles) and _Observer_ (blue
circles) reconstructions, along with their respective LDFs, Eqs. 5.6 and 5.7.
The reference signals at the radial distance of 1000 m are indicated by the
dashed line with values of $S(1000)=112$ VEM and 107 VEM for this particular
event. See Fig. 1 for the image of the corresponding footprint. _Right:_
Distribution of distances,
$|\Delta\vec{x}_{\text{c}}|=|\vec{x}_{\text{c}}^{\text{LLP}}-\vec{x}_{\text{c}}^{\text{NKG}}|$
between the impact points reconstructed in two different ways: with a Log-Log
Parabola (LLP) LDF in _Herald_ and a modified NKG in _Observer_. The
distribution has been fitted with a Rayleigh function (red line) indicating
that in 68% of the events, the two reconstructed impact points are at a
distance smaller than 92 m. _Inset:_ 2D distribution of $\Delta
x=x_{\text{LLP}}-x_{\text{NKG}}$ and $\Delta y=y_{\text{LLP}}-y_{\text{NKG}}$.
These coordinates have been computed in the reference frame of the _Observer_
shower and are oriented so that the positive/negative $\Delta y$ axis is in
the down/upstream direction.

The reconstructions described above successfully converge more than 99% of the
time. The convergence failures are mostly due to low-multiplicity events
formed by random coincidences which only by chance pass our selection
procedures. It is important to note that using the non-triggered stations adds
a significant geometrical constraint to the LDF fit, in particular for low-
multiplicity events. Without consideration of the non-triggered stations,
deviations of the shower size on the order of up to 8% are observed.

An example of the LDF fits obtained with the _Herald_ and _Observer_
reconstructions is shown in Fig. 8-left. In this example, the two LDF fits
give very similar results for the estimation of the shower size $S(1000)$. The
differences between the two LDFs (defined in Eqs. 5.6 and 5.7) are due to
uncertainty in the shape of the underlying lateral distribution close to and
far from the shower core. These differences lead to a distribution of
distances $|\Delta\vec{x}_{\text{c}}|$ between the shower impact points shown
in Fig. 8-right. For the two reconstructions, 68% of the events exhibit a
difference of less than 100 m. The systematic differences in the positions of
the impact points (as shown in the inset of Fig. 8-right) are largely due to
the asymmetry correction defined in Eq. 5.8, which is only employed in the
_Herald_ reconstruction.

Figure 9: _Left:_ Comparison of $S(1000)$ distributions for both
reconstructions. The distribution can be fitted with the following power law:
$S(1000)_{\text{Herald}}/\text{VEM}=1.154(S(1000)_{\text{Observer}}/\text{VEM})^{0.9904}$
(in black line). _Right:_ Relative difference $\Delta
S(1000)/S(1000)_{\text{Observer}}$, where $\Delta
S(1000)=S(1000)_{\text{Herald}}-S(1000)_{\text{Observer}}$, for all events and
also only those with $S(1000)$ greater than 10 or 100 VEM. The corresponding
means (standard deviations) of the distributions are 15.4% (10.8%), 11.2%
(9.2%), and 10.9% (5.3%), respectively.

In Fig. 9, a similar comparison is shown for the shower-size parameter
$S(1000)$. The skewness of the distributions seen in Fig. 9-right is a
consequence of a zenith-angle dependence of the systematic difference between
the two reconstructions of $S(1000)$. These systematic differences are
corrected for in the angular portion of the calibration procedure, which is
described in Section 8. The fit of $S(1000)_{\text{Herald}}$ vs.
$S(1000)_{\text{Observer}}$ shows that the correlation between the two
reconstructions follows a simple power law. Deviations from this power law do
not exceed 3% at any point across the range of shower sizes. Note that as long
as the correlation follows a power law, the energy calibration procedure
summarized in Section 8 will exactly reconcile any differences in the shower
size estimates of the two reconstructions.

##  6 Angular resolution

The angular accuracy of the reconstruction of the SD events is mostly driven
by the multiplicity of triggering stations and the precision of determination
of the time at which the shower front arrived at each station. The latter is
governed by an interplay between the number of particles, their time
distribution in the shower front, and the area of the detectors. The small
jitter of ${\sim}10$ ns in the GPS timing system [13] induces an angular
uncertainty111111 Increasing the jitter in the GPS timing to 25 ns (100 ns)
would induce an effect ${\sim}0.2^{\circ}$ (${\sim}1.5^{\circ}$) for the
smallest and ${<}0.1^{\circ}$ (${\sim}0.5^{\circ}$) for the largest shower
sizes. of less than ${\sim}0.1^{\circ}$. The uncertainty introduced by the 40
MHz sampling of the signals is of the same order of magnitude.

Deriving the angular resolution of cosmic-ray experiments usually relies on
the use of simulations where the reconstructed arrival direction is compared
with the injected one. We also present such a simulation-based derivation here
but additionally show results of measurement-driven studies of the angular
resolution. The angular resolution is obtained by comparing two arrival
directions $\hat{a}_{1}$ and $\hat{a}_{2}$, which are either (a) reconstructed
and true directions in the case of simulations, (b) two arrival directions
reconstructed by the same procedure and thus with the same resolution, or (c)
two arrival directions obtained by two different reconstruction methods, which
may have different resolutions. The angular difference $\eta$ is obtained from
$\sin\eta=|\hat{a}_{1}\times\hat{a}_{2}|$. For $\eta\ll 1$ and in a normal
approximation with one-dimensional variance $\sigma^{2}$, $\eta$ follows a
distribution $\mathrm{d}N/\mathrm{d}\eta=\mathcal{R}(\eta;\sigma)$ for the
case of (a) and $\mathrm{d}N/\mathrm{d}\eta=\mathcal{R}(\eta;\sqrt{2}\sigma)$
for the case of (b), where
$\mathcal{R}(r;\sigma)=(r/\sigma^{2})\exp(-r^{2}/2\sigma^{2})$ is the Rayleigh
distribution. We define the angular resolution (AR) as the angle at which the
value of the cumulative distribution of $\eta$ reaches 68.3%, this being in
turn equivalent to $\text{AR}\approx 1.52\sigma$ or $1.52\sigma/\sqrt{2}$ for
the two cases, respectively [30]. For the case (c) we assume that the two
measurements $\sigma_{1}$ and $\sigma_{2}$ are uncorrelated and thus can be
added in quadrature, i.e.
$\sigma_{\text{tot}}^{2}=\sigma_{1}^{2}+\sigma_{2}^{2}$.

Below, three different approaches are explored to derive the angular
resolution from measurements only, followed by a comparison of the results
with a derivation from simulations121212 Results from the _Observer_ framework
are shown for these analyses and are compatible with those from _Herald_..

###  6.1 Derivation from measurements

####  6.1.1 Hybrid data

In addition to the independent triggering systems of the FD and SD detectors,
a hybrid trigger [31] has been designed requiring a coincidence of at least
one SD station and an FD telescope. A sub-sample of these hybrid events can be
fully reconstructed by both the FD hybrid reconstruction [8] and,
independently, by the SD reconstruction. These events are used, among other
things, for the calibration of the SD energy estimator with the calorimetric
measurement performed by the FD [10]. From the direct comparison of the
arrival directions of the two reconstructions, $\hat{a}_{\text{SD}}$ and
$\hat{a}_{\text{hyb}}$, the angular resolution $\text{AR}_{\text{SD:hyb}}$ is
extracted from the angular difference $\eta_{\text{SD:hyb}}$ defined by
$\sin\eta_{\text{SD:hyb}}=|\hat{a}_{\text{SD}}\times\hat{a}_{\text{hyb}}|$.
Given the essential independence of the SD and hybrid reconstructions, the SD
angular resolution $\text{AR}_{\text{SD}}$ can be derived as
$\text{AR}_{\text{SD}}^{2}=\text{AR}_{\text{SD:hyb}}^{2}-\text{AR}_{\text{hyb}}^{2}$,
where $\text{AR}_{\text{hyb}}$ is the angular resolution of the hybrid events
only. With criteria to ensure the quality of the hybrid reconstruction, 29 344
events are selected. The details of the criteria applied can be found in [30,
32].

Moreover, 1049 events are seen by at least two telescopes belonging to two
different FD sites. For these _stereo_ hybrid events, a separate131313
Separate, but not completely independent since they share the same SD station
used for the timing information. hybrid reconstruction is performed for each
FD site. Computing the angular difference $\eta_{\text{hyb}}$ between these
arrival directions $\hat{a}^{\text{hyb}}_{1}$ and $\hat{a}^{\text{hyb}}_{2}$,
the resolution of the hybrid reconstruction of this particular selection of
stereo events is estimated to be $\text{AR}_{\text{hyb}}=(1.07\pm
0.05)^{\circ}$ assuming that resolutions for $\hat{a}^{\text{hyb}}_{1}$ and
$\hat{a}^{\text{hyb}}_{2}$ are the same141414 Note that stereo events have
different geometries than standard hybrid (mono) events (the latter being on
average closer to the telescopes) and therefore this number does not reflect
the average resolution of the hybrid reconstruction which was estimated to be
$0.5^{\circ}$ in [30].. With $\text{AR}_{\text{hyb}}$ at hand,
$\text{AR}_{\text{SD}}$ can be estimated as given above. The results are shown
in Fig. 11 (squares) as a function of the zenith angle and the shower size.

Figure 10: A map of the high-density sector of the SD array which is in the
North-West corner of the Observatory (see Fig. 1). Regular stations of the
1500 m grid are indicated with black markers. Stations of the 750 m grid are
indicated with blue markers, where the stations forming three additional 1500
m sub-grids are indicated with triangular, square, and pentagonal markers.
Twin and triplet stations in multiplets are indicated with red and green
circles, respectively.

####  6.1.2 Super-hexagon of multiplets

Located in the denser sector of the SD array, 12 twins and 7 triplets of
stations (two or three stations separated by ${\sim}11$ m) constitute the
_super-hexagon_ of the Observatory, as schematically depicted in Fig. 10 (red
and green circles). The twins (or triplets) of the super-hexagon can be
randomly divided to form two independent arrays and thereby enable two
independent reconstructions of events that land in this region and trigger at
least three twins. All other stations outside of the super-hexagon are not
included in this specialized event reconstruction. Note that for showers of
the same size, the multiplicity of stations used in this reconstruction can
be, relative to the normal reconstruction, greatly reduced given the limited
size of the super-hexagon.

Due to the small number of pairs available for this study, a selection on the
position of the reconstructed core is applied. Only 913 reconstructed
events151515 Requirement of having the same number of stations in each twin
set reduces the number of events to 301, nevertheless, the AR results remain
statistically compatible., for which the impact point lies within a seed
triangle constituted by twins only, are kept for the analysis. The estimate of
the angular resolution $\text{AR}_{\text{SD}}$ is obtained from the
distribution of the angular difference $\eta_{\text{tw}}$, defined as
$\sin\eta_{\text{tw}}=|\hat{a}^{\text{tw}}_{1}\times\hat{a}^{\text{tw}}_{2}|$.
The results of this procedure are shown in Fig. 11.

####  6.1.3 Sub-arrays

The stations comprising the 750 m array can be split into four small sub-
arrays with the usual 1500 m spacing, as illustrated in Fig. 10 (blue
symbols). Respectively, 841, 516, and 392 events are independently
reconstructed two, three, and four times, whereby requiring that the impact
point of each reconstruction lies inside the corresponding seed triangle.

For a given event with several successful reconstructions, the arrival
directions are compared in pairs, i.e. by computing the angular difference
$\eta_{\text{sub}_{ij}}$ from
$\sin\eta_{\text{sub}_{ij}}=|\hat{a}^{\text{sub}}_{i}\times\hat{a}^{\text{sub}}_{j}|$.
Here, we assume that all $\eta_{\text{sub}_{ij}}$ are independent, obtaining
4741 angles for the calculation of the estimate of angular resolution
$\text{AR}_{\text{SD}}$. The results are reported in Fig. 11 (circles).

###  6.2 Estimation using simulations

Using simulations of extensive air-showers, the angular resolution is obtained
as a function of zenith angle $\theta$ and the shower size $S(1000)$. Showers
initiated by a proton or an iron primary and developing according to the EPOS-
LHC [33] or QGSJet-II.04 [34] model of hadronic interactions are simulated
with energies between $10^{18.5}$ and $10^{20}$ eV and zenith angles between
$0^{\circ}$ and $60^{\circ}$. An SD simulation and event reconstruction is
performed 10 times for each shower simulation, where the impact position on
the SD array is randomly chosen according to a uniform probability
distribution across the complete array.

The reconstructed arrival directions $\hat{a}_{\text{rec}}$ are compared with
the true simulated axis $\hat{a}$ of the shower by computing the angular
difference $\eta$ from $\sin\eta=|\hat{a}\times\hat{a}_{\text{rec}}|$ and
estimating the AR from its distribution. For plotting purposes, the results
are empirically approximated to a good degree as a second-order polynomial in
$\sin^{2}\theta$ and as an exponential function in $\lg(S(1000)/\text{VEM})$.
The results are reported in Fig. 11 (lines and bands).

Figure 11: Angular resolution of the reconstructed arrival direction as
derived from a study of the hybrid events (filled black squares), a study of
events with at least three triggered twins (filled black circles), and a study
of three sub-arrays (open circles). The lines report the results from the
study with simulations produced with the hadronic model QGSJet-II.04 (full
lines) and EPOS-LHC (dashed lines) for iron (dark gray lines) and proton
(light gray lines) primaries. The angular resolution of the reconstructed
arrival direction is shown as a function of zenith angle (left) and the shower
size (right).

###  6.3 Comparison of results

The AR, as derived from both measured and simulated data, is shown to improve
with increasing zenith angle and shower size. The improvement for inclined
events is attributed to the increase in the number of stations participating
in an event due to the foreshortening of the array.

The slight discrepancies between the AR obtained from the measured and
simulated data can be explained by limitations in the reconstruction of
measurements in these specialized studies; considering that only a small
number of SD stations were available for the twin and sub-array studies, the
typical multiplicities of stations in events decreased and with them, the
angular accuracy of the performed reconstructions. Finally, computing the SD
angular resolution from events reconstructed by both the FD and the SD in a
data-driven way is limited by the knowledge of the angular resolution of the
hybrid reconstruction $\text{AR}_{\text{hyb}}$. There are fewer high-quality
events measured by the FD due to its limited duty cycle and the selection
criteria applied. Moreover, in order to derive $\text{AR}_{\text{hyb}}$ and
then $\text{AR}_{\text{SD}}$, we distinguish between the FD reconstructions
where the telescope was closer or further from the axis of the shower. The
geometrical reconstruction of the latter is thus expected to be worse, which
would lead to an overestimation of the accuracy of the SD.

It is known that there are some deficiencies in the simulations of air
showers, particularly with respect to the magnitude of the muonic component
[35]. Additionally, while simulations of the detector response have been
thoroughly validated [36], it is still possible that some experimental effects
are not reproduced. Performing comparisons of resolution estimates between
simulated and measured data in terms of shower size mitigates the impact of
the shortcomings of the simulations to the first-order. The actual angular
resolution of the SD is expected to lie somewhere between the measured data
points and the simulated lines. The angular resolution improves with both
increasing zenith angle and shower size but is always better than
$1.4^{\circ}$ and even approaches $0.7^{\circ}$ for the largest and most
inclined showers.

It is also important to note that the differences between the results obtained
from the _Herald_ and _Observer_ reconstructions are smaller than the angular
resolution of the SD. It is also worth mentioning that the angular differences
observed in Fig. 6-right are primarily due to the correction of the signal
asymmetry present in the _Herald_ reconstruction only. These asymmetry
corrections then induce a shift of the impact point of the shower core, as
observed in Fig. 8-right, which is, in turn, translated into the angular
differences observed. The resolution of the impact point depends on the
intrinsic properties of the shower and varies from ${\sim}100$ m at the lowest
energies to ${\sim}50$ m at the highest energies. The core resolution
(measured in the ground plane) worsens with increasing inclination, which can
be attributed to the projection effect.

##  7 Uncertainties in $S(1000)$

The uncertainty of the shower size estimator $S(1000)$ consists of statistical
and systematic contributions. While the statistical uncertainties of the
reconstructed shower size, $\sigma_{\text{stat}}(S(1000))$, are directly
related to the number of triggered stations and the uncertainties in their
signals, the systematic uncertainties, $\sigma_{\text{syst}}(S(1000))$, arise
from the lack of knowledge of the true shape of the LDF.

The uncertainties of $S(1000)$ are first estimated through propagation of
statistical and systematic errors in the fit of the LDF. These estimations are
complemented by two data-driven methods for the lower shower-size regime,
where sufficient data are available, and a comparison with uncertainty
estimates derived from simulations.

###  7.1 Fit uncertainties and model-dependent systematic errors

The statistical uncertainties in $S(1000)$ are directly estimated during the
fitting procedure. The derivation of the systematic uncertainties requires a
deeper study of the slope parameter $\beta$ from Eqs. 5.6 and 5.7. As a first
step, we use the same subset of events already used in the determination of
$\beta$ in Section 5.2.1, i.e. events for which $\beta$ can be reliably
fitted. The resulting standard deviation $\sigma_{\beta}$ of the difference
between the fixed and free $\beta$ is parametrized as a function of $S(1000)$.
In a second step, the systematic uncertainties are computed by propagating
$\sigma_{\beta}(S(1000))$ into the shower size estimator for all events.

This approach is applied to the full data set of events recorded by the SD,
while using the quality cuts described in the previous sections. The average
systematic uncertainties and statistical errors are reported in Fig. 12 as a
function of $\lg(S(1000)/\text{VEM})$ (left) and as a function of
$\sin^{2}\theta$ (right) for saturated and non-saturated events. In both
cases, the statistical uncertainty dominates the total uncertainty for smaller
shower sizes; however, for larger shower sizes, the systematic uncertainties
become significantly larger in the case of saturated events. This implies that
an accurate description of the shape of the lateral distribution is of
increasing relevance at the highest energies, where our estimates of the
shower size may benefit from an improved knowledge thereof. Note that the
uncertainties in Fig. 12-right are practically constant for the non-saturated
events when considering all energies together. On the other hand, the
resolution for the saturated events improves with increasing zenith angle.
This is a consequence of the fact that the saturated station for most of the
near-vertical saturated events is close to the shower core, resulting in a
clustering of the six nearest stations around a distance of ${\sim}1500$ m
from the shower axis, which is not the case for more inclined showers. At the
highest energies, these clustered configurations become responsible for the
systematic uncertainties being twice as high as the statistical uncertainties,
and the recovery procedure of the saturated signals (Section 3.4) has an
effect of only around 2% on the shower size estimator. A better knowledge of
the shape of the lateral distribution is thus needed to reduce the
uncertainties in $S(1000)$. Plans to improve this knowledge are discussed in
Section 9.

Figure 12: Relative systematic (syst.) and statistical (stat.) uncertainty
$\sigma(S(1000))/S(1000)$ of the shower-size reconstruction as a function of
shower size $\lg(S(1000)/\text{VEM})$ (left) and $\sin^{2}\theta$ (right),
obtained separately for saturated and non-saturated events.

###  7.2 Direct measurements

####  7.2.1 Super-hexagon of multiplets

As in the studies of the angular resolution discussed in Section 6.1.2, it is
also possible to perform independent reconstructions of $S(1000)$ for events
with at least three triggered twins located in the super-hexagon for which the
quantity $\delta_{\text{tw}}=\Delta S(1000)_{\text{tw}}/\sqrt{2}\langle
S(1000)\rangle_{\text{tw}}$ is computed, where $\Delta S(1000)_{\text{tw}}$
and $\langle S(1000)\rangle_{\text{tw}}$ are the difference and the average,
respectively, of the two reconstructed shower sizes $S(1000)_{1}^{\text{tw}}$
and $S(1000)_{2}^{\text{tw}}$. In each bin in $\lg(S(1000)/\text{VEM})$ or in
$\sin^{2}\theta$, the standard deviation $\sigma(\delta_{\text{tw}})$ is
calculated. The standard deviations of these distributions are plotted in Fig.
14 (as filled squares).

####  7.2.2 Sub-arrays

For events landing in the denser sector of the SD, a similar procedure as the
one used for twin study is applied. Using the events reconstructed two, three,
or four times with the stations from different sub-arrays, the quantity
$\delta_{\text{sub}}=\Delta S(1000)_{\text{sub}}/\sqrt{2}\langle
S(1000)\rangle_{\text{sub}}$ is derived, where $\Delta S(1000)_{\text{sub}}$
and $\langle S(1000)\rangle_{\text{sub}}$ are the difference and the average
of the reconstructed shower sizes $S(1000)^{\text{sub}}_{i}$ and
$S(1000)^{\text{sub}}_{j}$, respectively. For each pair of sub-arrays $(i,j)$,
the standard deviations $\sigma(\delta_{\text{sub}})_{ij}$ are calculated and
then averaged over all pairs $(i,j)$ in different bins of
$\lg(S(1000)/\text{VEM})$ and $\sin^{2}\theta$. The results are reported in
Fig. 14 (as open squares).

###  7.3 Estimation using simulations

A discrete set of simulations is used to determine the total uncertainty of
the energy estimator $S(1000)$. Showers initiated by a proton or an iron
primary and developing according to either the EPOS-LHC or the QGSJet-II.04
hadronic model are simulated with energies of $10^{18.5}$, $10^{19}$, and
$10^{19.5}$ eV and zenith angles of $0^{\circ}$, $12^{\circ}$, $22^{\circ}$,
$32^{\circ}$, $38^{\circ}$, $48^{\circ}$, and $56^{\circ}$. An SD simulation
and event reconstruction is performed 10 times for each simulated shower.

For simulations, the true value of $S(1000)$, denoted as
$S(1000)_{\text{true}}$, is determined by calculating the mean signal of 24
additional stations in a ring with a shower-plane radius of exactly 1000 m.
The shower sizes $S(1000)_{\text{true}}$ and $S(1000)_{\text{rec}}$, where the
latter is the value of the shower size reconstructed according to the
procedure detailed in Section 5.2.1, are compared for different primaries and
for different zenith angles by calculating the ratio
$\delta_{\text{sim}}=\Delta S(1000)_{\text{sim}}/S(1000)_{\text{true}}$ where
$\Delta S(1000)_{\text{sim}}=S(1000)_{\text{rec}}-S(1000)_{\text{true}}$. The
standard deviation $\sigma(\delta_{\text{sim}})$ of the relative differences
is reported in Fig. 13. The results show a total uncertainty in
$S(1000)_{\text{rec}}$, which decreases from 18% at the smallest shower sizes
to nearly 7% at the largest values of $S(1000)_{\text{true}}$. At larger
inclinations, the multiplicity of stations within a given distance from the
shower axis increases, whereas the magnitudes of the signals decrease due to
the larger attenuation of showers in the atmosphere. These competing effects,
which are shown in Fig. 13, result in the observed evolution of the standard
deviation $\sigma(\delta_{\text{sim}})$ with zenith angle.

Figure 13: Relative uncertainty of the reconstructed shower size
$\sigma(S(1000)_{\text{rec}})/S(1000)_{\text{true}}$ as obtained from QGSJet-
II.04 (left) and EPOS-LHC (right) simulations of proton (red) and iron (blue)
primaries for several initial Monte-Carlo energies and zenith angles of
$0^{\circ}$, $12^{\circ}$, $22^{\circ}$, $32^{\circ}$, $38^{\circ}$,
$48^{\circ}$, and $56^{\circ}$ (indicated with increasing transparency of the
markers).

###  7.4 Comparison of results

Figure 14: Relative total uncertainty of the shower-size reconstruction
$\sigma(S(1000))/S(1000)$ derived from a study of statistical and systematic
uncertainties separately (black circles), a study of events with at least
three triggered twins (black squares), and a study of three sub-grid arrays
(black open squares). The lines report the results from the study with
simulations produced with the hadronic models QGSJet-II.04 (full lines) and
EPOS-LHC (dashed lines) for iron (dark gray lines) and proton (light gray
lines) primaries. The relative total uncertainty is shown as a function of the
shower size (left) and zenith angle (right).

The statistical and systematic uncertainties of the reconstructed $S(1000)$
for all measured data are added in quadrature as
$\sigma^{2}(S(1000)_{\text{rec}})=\sigma^{2}_{\text{syst}}(S(1000)_{\text{rec}})+\sigma^{2}_{\text{stat}}(S(1000)_{\text{rec}})$,
where $\sigma_{\text{syst}}$ is from Fig. 7-left. The result is compared with
both data-driven and simulation-based estimates in Fig. 14, where the results
from simulations have been parameterized with simple polynomials.

While some discrepancies may be observed and can be attributed to the
systematic errors associated with each method of estimation, all methods
generally agree to within a few percentage points. The slight worsening of the
resolution from $\sim$6% at $\lg(S(1000)/\text{VEM})\approx 2.3$ to $\sim$8%
at the largest shower sizes is due to the fact that at a shower size of
$S(1000)=200$ VEM, 50% of the events have a saturated station and this
fraction increases with shower size. As seen in Fig. 14-right, the resolution
is relatively constant at 14% (since small shower sizes dominate calculation
of these estimates) over the zenith angle range with a slight decrease
observed for the most inclined shower sizes.

The impact of the resolution in shower size on the uncertainties in the
estimate of the primary energy is described in detail in [10].

##  8 Towards an efficient estimate of the energy

In the previous sections, we presented the reconstruction of the arrival
direction and the shower size from the timing and signal measurements of SD
stations. The estimation of the shower size is a first important step in the
determination of the energy of the primary cosmic ray. The second step
involves the correction of a number of biases to improve the energy
resolution. The shower size is further corrected for daily and yearly
atmospheric variations, which change the effective atmospheric overburden of
the array. This correction has been recently improved [37] using eleven years
of data from the Observatory. Next, a bias resulting from (oppositely) charged
shower particles experiencing deflection in the geomagnetic field is removed
[38]. This effect causes a dependence of the shower size on the azimuthal
orientation of the shower with respect to the orientation of the local
geomagnetic field. The shower size is corrected with a factor
$1+A_{\text{gm}}\,\cos^{k}\theta\,\sin^{2}\xi$, where $A_{\text{gm}}$ and $k$
are model parameters estimated from the Monte-Carlo data, and $\xi$ is the
angle between the geomagnetic field and the shower axis. The maximum value of
this correction stays below 4%.

The next steps are described in detail in [10], but are briefly summarized
here. Using the constant intensity cut (CIC) method [3, 39], the measured
shower size $S(1000)$ is converted to an angle-independent shower size
$S_{38}$. This is equivalent to the size of the shower had the primary
particle arrived with a zenith angle of $38^{\circ}$ (median of all events),
$S_{38}=S(1000)/f_{\text{CIC}}(\theta)$. This correction accounts for the
increasing slant depth of the array with increasing zenith angle. Since the
more inclined showers are sampled at a later shower age, the more-attenuated
electromagnetic cascade leads to a smaller observed shower size. With this
procedure, we obtain the minimally biased, zenith-independent energy estimator
$S_{38}$. This can be directly calibrated [40] by the nearly-calorimetric
energy measurement of the FD using hybrid events (see Section 6.1.1). The
power-law calibration curve $E=A\,(S_{38}/\text{VEM})^{B}$, obtained with
these _hybrid_ events, is then used to assign the energy to _all_ SD events.
All remaining lowest-order differences in reconstruction between the _Herald_
and _Observer_ frameworks are essentially removed by performing this energy
calibration step independently for each.

##  9 Conclusion and outlook

We have presented the methods used for reconstructing the properties of the
highest-energy cosmic rays from measurements of the surface-detector array of
the Pierre Auger Observatory. This description ranges from the station-level
calibration procedure to the event-level reconstruction of the arrival
direction and the shower size. At the core of these methods is an emphasis on
measurement-driven approaches.

At each step of the reconstruction procedure, whether it be modeling the
curvature of the impinging shower front or fitting the lateral distribution of
particles at the ground, a number of differing approaches may be developed and
applied, each with its own assumptions. Within the Pierre Auger Collaboration,
exploration and validation of different assumptions and methods has been
facilitated by the two reconstruction frameworks, _Herald_ and _Observer_. The
existence of these independent frameworks has provided a source of cross-
validation and has placed a spotlight onto the impact of different procedures
on the reconstructed observables of the primary cosmic-rays throughout the
development of our reconstruction procedures. It has also resulted in the
production of high-quality data sets of reconstructed events. Historically,
the _Herald_ data set has been used for anisotropy analyses and that of
_Observer_ for the energy spectrum and mass composition studies, where one or
the other was chosen out of practical reasons for publication. Nonetheless,
the consistency of physical results has been verified using both frameworks.
The angular resolution of the reconstructed arrival direction was shown to be
of the order of $1^{\circ}$ and approaches $0.5^{\circ}$ for the largest
shower sizes (i.e. above 200 VEM). The resolution in the reconstructed shower
size, from which the energy of the primary is inferred, was demonstrated to
improve from ${\sim}15\%$, for the smallest shower sizes down to ${\sim}6\%$
for the largest. These values correspond to primaries of approximately
$10^{18.5}$ eV and above $10^{19.5}$ eV, respectively (for more detailed
interpretation in terms of primary energy, see [10]). At the highest energies,
systematic uncertainties due to imperfect knowledge of the shape of the
lateral distribution within a few hundred meters of the shower core become of
the same order as statistical uncertainties, which will be addressed by the
extension of the dynamic range of the WCDs through the installation of a
smaller PMT [41].

Although the reconstruction methods presented here were developed and applied
to the 1500 m array of the Observatory, scaled adaptations thereof are also
successfully applied to the lower-energy extension of the surface detector
array with a denser detector spacing of 750 m [7]. In both cases, the results
are used for higher-level physics analyses, including studies of anisotropies
in arrival directions [9, 42], multimessenger [43] and mass-composition
studies [44, 45], and the cosmic-ray energy spectrum [10, 46, 47].

Looking to the future, the surface detector of the Pierre Auger Observatory is
currently undergoing a large-scale upgrade called AugerPrime [48, 49], with
the placement of a 3.8 m2 scintillator atop each water-Cherenkov detector as
one of the principal components. Exploitation of the differing responses of
the two detector types to the electromagnetic and muonic components of
extensive air showers will permit further improvements to the event
reconstruction algorithms described here and, most notably, add a primary-mass
estimate to the set of reconstructed observables.

## Acknowledgments

The successful installation, commissioning, and operation of the Pierre Auger
Observatory would not have been possible without the strong commitment and
effort from the technical and administrative staff in Malargüe. We are very
grateful to the following agencies and organizations for financial support:

Argentina – Comisión Nacional de Energía Atómica; Agencia Nacional de
Promoción Científica y Tecnológica (ANPCyT); Consejo Nacional de
Investigaciones Científicas y Técnicas (CONICET); Gobierno de la Provincia de
Mendoza; Municipalidad de Malargüe; NDM Holdings and Valle Las Leñas; in
gratitude for their continuing cooperation over land access; Australia – the
Australian Research Council; Brazil – Conselho Nacional de Desenvolvimento
Científico e Tecnológico (CNPq); Financiadora de Estudos e Projetos (FINEP);
Fundação de Amparo à Pesquisa do Estado de Rio de Janeiro (FAPERJ); São Paulo
Research Foundation (FAPESP) Grants No. 2019/10151-2, No. 2010/07359-6 and No.
1999/05404-3; Ministério da Ciência, Tecnologia, Inovações e Comunicações
(MCTIC); Czech Republic – Grant No. MSMT CR LTT18004, LM2015038, LM2018102,
CZ.02.1.01/0.0/0.0/16_013/0001402, CZ.02.1.01/0.0/0.0/18_046/0016010 and
CZ.02.1.01/0.0/0.0/17_049/0008422; France – Centre de Calcul IN2P3/CNRS;
Centre National de la Recherche Scientifique (CNRS); Conseil Régional Ile-de-
France; Département Physique Nucléaire et Corpusculaire (PNC-IN2P3/CNRS);
Département Sciences de l’Univers (SDU-INSU/CNRS); Institut Lagrange de Paris
(ILP) Grant No. LABEX ANR-10-LABX-63 within the Investissements d’Avenir
Programme Grant No. ANR-11-IDEX-0004-02; Germany – Bundesministerium für
Bildung und Forschung (BMBF); Deutsche Forschungsgemeinschaft (DFG);
Finanzministerium Baden-Württemberg; Helmholtz Alliance for Astroparticle
Physics (HAP); Helmholtz-Gemeinschaft Deutscher Forschungszentren (HGF);
Ministerium für Innovation, Wissenschaft und Forschung des Landes Nordrhein-
Westfalen; Ministerium für Wissenschaft, Forschung und Kunst des Landes Baden-
Württemberg; Italy – Istituto Nazionale di Fisica Nucleare (INFN); Istituto
Nazionale di Astrofisica (INAF); Ministero dell’Istruzione, dell’Universitá e
della Ricerca (MIUR); CETEMPS Center of Excellence; Ministero degli Affari
Esteri (MAE); México – Consejo Nacional de Ciencia y Tecnología (CONACYT) No.
167733; Universidad Nacional Autónoma de México (UNAM); PAPIIT DGAPA-UNAM; The
Netherlands – Ministry of Education, Culture and Science; Netherlands
Organisation for Scientific Research (NWO); Dutch national e-infrastructure
with the support of SURF Cooperative; Poland -Ministry of Science and Higher
Education, grant No. DIR/WK/2018/11; National Science Centre, Grants No.
2013/08/M/ST9/00322, No. 2016/23/B/ST9/01635 and No. HARMONIA
5–2013/10/M/ST9/00062, UMO-2016/22/M/ST9/00198; Portugal – Portuguese national
funds and FEDER funds within Programa Operacional Factores de Competitividade
through Fundação para a Ciência e a Tecnologia (COMPETE); Romania – Romanian
Ministry of Education and Research, the Program Nucleu within MCI
(PN19150201/16N/2019 and PN19060102) and project PN-
III-P1-1.2-PCCDI-2017-0839/19PCCDI/2018 within PNCDI III; Slovenia – Slovenian
Research Agency, grants P1-0031, P1-0385, I0-0033, N1-0111; Spain – Ministerio
de Economía, Industria y Competitividad (FPA2017-85114-P and FPA2017-85197-P),
Xunta de Galicia (ED431C 2017/07), Junta de Andalucía (SOMM17/6104/UGR), Feder
Funds, RENATA Red Nacional Temática de Astropartículas (FPA2015-68783-REDT)
and María de Maeztu Unit of Excellence (MDM-2016-0692); USA – Department of
Energy, Contracts No. DE-AC02-07CH11359, No. DE-FR02-04ER41300, No. DE-
FG02-99ER41107 and No. DE-SC0011689; National Science Foundation, Grant No.
0450696; The Grainger Foundation; Marie Curie-IRSES/EPLANET; European Particle
Physics Latin American Network; and UNESCO.

## References

  * [1] P. Bassi, G. Clark and B. Rossi, _Distribution of Arrival Times of Air Shower Particles_ , _Phys. Rev._ 92 (1953) 441.
  * [2] N. A. Porter, T. E. Cranshaw, J. F. D. Beer, A. G. Parham and A. C. Sherwood, _Observations on extensive air showers vii. the lateral distribution of energy in the electron-photon component_ , _The Philosophical Magazine: A Journal of Theoretical Experimental and Applied Physics_ 3 (1958) 826.
  * [3] J. Hersil, I. Escobar, D. Scott, G. Clark and S. Olbert, _Observations of Extensive Air Showers near the Maximum of Their Longitudinal Development_ , _Phys. Rev. Lett._ 6 (1961) 22.
  * [4] H. E. Bergeson, G. L. Cassiday, T.-W. Chiu, D. A. Cooper, J. W. Elbert, E. C. Loh et al., _Measurement of light emission from remote cosmic ray air showers_ , _Phys. Rev. Lett._ 39 (1977) 847.
  * [5] R. M. Baltrusaitis et al., _The Utah Fly’s Eye Detector_ , _Nucl. Instrum. Methods A_ 240 (1985) 410.
  * [6] A. M. Hillas, _Derivation of the EAS spectrum_ , _Conf. Proc. C_ 690825 (1969) 355.
  * [7] Pierre Auger collaboration, _The Pierre Auger Cosmic Ray Observatory_ , _Nucl. Instrum. Meth. A_ 798 (2015) 172 [1502.01323].
  * [8] Pierre Auger collaboration, _The Fluorescence Detector of the Pierre Auger Observatory_ , _Nucl. Instrum. Meth. A_ 620 (2010) 227 [0907.4282].
  * [9] Pierre Auger collaboration, _Observation of a Large-scale Anisotropy in the Arrival Directions of Cosmic Rays above $8{\times}10^{18}$ eV_, _Science_ 357 (2017) 1266 [1709.07321].
  * [10] Pierre Auger collaboration, _A measurement of the cosmic ray energy spectrum above $2.5{\times}10^{18}$ eV using the Pierre Auger Observatory_, _submitted to Phys. Rev. D_ (2020) [2008.06486].
  * [11] Pierre Auger collaboration, _Reconstruction of Inclined Air Showers Detected with the Pierre Auger Observatory_ , _JCAP_ 08 (2014) 019 [1407.3214].
  * [12] S. Argiró, S. L. C. Barroso, J. Gonzalez, L. Nellen, T. C. Paul, T. A. Porter et al., _The Offline Software Framework of the Pierre Auger Observatory_ , _Nucl. Instrum. Meth. A_ 580 (2007) 1485 [0707.1652].
  * [13] Pierre Auger collaboration, _Timing calibration and synchronization of surface and fluorescence detectors of the Pierre Auger Observatory_ , in _29th International Cosmic Ray Conference_ , 8, 2005\.
  * [14] Pierre Auger collaboration, _Trigger and Aperture of the Surface Detector Array of the Pierre Auger Observatory_ , _Nucl. Instrum. Meth. A_ 613 (2010) 29 [1111.6764].
  * [15] Pierre Auger collaboration, _Large-scale cosmic-ray anisotropies above 4 EeV measured by the Pierre Auger Observatory_ , _Astrophys. J._ 868 (2018) 4 [1808.03579].
  * [16] Pierre Auger collaboration, _Cosmic-ray anisotropies in right ascension measured by the Pierre Auger Observatory_ , _Astrophys. J._ 891 (2020) 142 [2002.06172].
  * [17] Pierre Auger collaboration, _The Pierre Auger Observatory Scaler Mode for the Study of Solar Activity Modulation of Galactic Cosmic Rays_ , _JINST_ 6 (2011) P01003.
  * [18] Pierre Auger collaboration, _Calibration of the surface array of the Pierre Auger Observatory_ , _Nucl. Instrum. Meth. A_ 568 (2006) 839.
  * [19] Pierre Auger collaboration, _Estimation of signal in saturated stations of Pierre Auger surface detector_ , in _Proc. of 33rd ICRC, Rio de Janeiro, Brazil_ , p. 0633, 2013, 1307.5059.
  * [20] Pierre Auger collaboration, _A model for the time uncertainty measurements in the Auger surface detector array_ , _Astropart. Phys._ 28 (2008) 523 [0705.1856].
  * [21] Pierre Auger collaboration, _The accuracy of signal measurement with the water Cherenkov detectors of the Pierre Auger Observatory_ , _Nucl. Instrum. Meth. A_ 578 (2007) 180.
  * [22] J. Linsley and L. Scarsi, _Arrival times of air shower particles at large distances from the axis_ , _Phys. Rev._ 128 (1962) 2384.
  * [23] F. James, _Minuit, function minimization and error analysis_ , long writeup D506, CERN, 1998.
  * [24] D. Newton, J. Knapp and A. A. Watson, _The Optimum Distance at which to Determine the Size of a Giant Air Shower_ , _Astropart. Phys._ 26 (2007) 414 [astro-ph/0608118].
  * [25] M. Nagano and A. A. Watson, _Observations and implications of the ultrahigh-energy cosmic rays_ , _Rev. Mod. Phys._ 72 (2000) 689.
  * [26] K. Kamata and J. Nishimura, _The lateral and the angular structure functions of electron showers_ , _Prog. Theor. Phys. Suppl._ 6 (1958) 93.
  * [27] K. Greisen, _The extensive air showers_ , _Prog. Cosmic Ray Phys._ 3 (1956) 1.
  * [28] K. Greisen, _Cosmic ray showers_ , _Ann. Rev. Nucl. Part. Sci._ 10 (1960) 63.
  * [29] A. M. Hillas, _Measurement of primary energy in air showers in the presence of fluctuations_ , in _Proc. of 12th ICRC, Hobart, Australia, 16–25 August, 1971_ , vol. 3, pp. 1001–1006, 1971.
  * [30] Pierre Auger collaboration, _Angular resolution of the Pierre Auger Observatory_ , in _29th International Cosmic Ray Conference_ , 8, 2005.
  * [31] Pierre Auger collaboration, _The Exposure of the Hybrid Detector of the Pierre Auger Observatory_ , _Astropart. Phys._ 34 (2011) 368 [1010.6162].
  * [32] Pierre Auger collaboration, _Hybrid Performance of the Pierre Auger Observatory_ , in _30th International Cosmic Ray Conference_ , vol. 4, pp. 425–428, 6, 2007, 0706.1105.
  * [33] T. Pierog, I. Karpenko, J. M. Katzy, E. Yatsenko and K. Werner, _EPOS LHC: Test of collective hadronization with data measured at the CERN Large Hadron Collider_ , _Phys. Rev. C_ 92 (2015) 034906 [1306.0121].
  * [34] S. Ostapchenko, _Monte Carlo treatment of hadronic interactions in enhanced Pomeron scheme: I. QGSJET-II model_ , _Phys. Rev. D_ 83 (2011) 014018 [1010.1869].
  * [35] Pierre Auger collaboration, _Testing Hadronic Interactions at Ultrahigh Energies with Air Showers Measured by the Pierre Auger Observatory_ , _Phys. Rev. Lett._ 117 (2016) 192001 [1610.08509].
  * [36] Pierre Auger collaboration, _Studies on the response of a water-Cherenkov detector of the Pierre Auger Observatory to atmospheric muons using an RPC hodoscope_ , 2007.04139.
  * [37] Pierre Auger collaboration, _Impact of Atmospheric Effects on the Energy Reconstruction of Air Showers Observed by the Surface Detectors of the Pierre Auger Observatory_ , _JINST_ 12 (2017) P02006 [1702.02835].
  * [38] Pierre Auger collaboration, _The effect of the geomagnetic field on cosmic ray energy estimates and large scale anisotropy searches on data from the Pierre Auger Observatory_ , _JCAP_ 11 (2011) 022 [1111.7122].
  * [39] T. K. Gaisser and A. M. Hillas, _Reliability of the method of constant intensity cuts for reconstructing the average development of vertical showers_ , in _Proc. of 15th ICRC, Plovdiv, Bulgaria, 13–26 August 1977_ , vol. 8, pp. 353–357, B’lgarska Akademiia na Naukite, Sofia, 1977.
  * [40] H. P. Dembinski, B. Kégl, I. C. Mariş, M. Roth and D. Veberič, _A likelihood method to cross-calibrate air-shower detectors_ , _Astropart. Phys._ 73 (2016) 44 [1503.09027].
  * [41] Pierre Auger collaboration, _The dynamic range of the AugerPrime Surface Detector: technical solution and physics reach_ , _PoS_ ICRC2017 (2018) 397.
  * [42] Pierre Auger collaboration, _An Indication of anisotropy in arrival directions of ultra-high-energy cosmic rays through comparison to the flux pattern of extragalactic gamma-ray sources_ , _Astrophys. J. Lett._ 853 (2018) L29 [1801.06160].
  * [43] Pierre Auger collaboration, _Multi-Messenger Physics with the Pierre Auger Observatory_ , _Front. Astron. Space Sci._ 6 (2019) 24 [1904.11918].
  * [44] Pierre Auger collaboration, _Evidence for a mixed mass composition at the ‘ankle’ in the cosmic-ray spectrum_ , _Phys. Lett. B_ 762 (2016) 288 [1609.08567].
  * [45] Pierre Auger collaboration, _Inferences on mass composition and tests of hadronic interactions from 0.3 to 100 EeV using the water-Cherenkov detectors of the Pierre Auger Observatory_ , _Phys. Rev. D_ 96 (2017) 122003 [1710.07249].
  * [46] Pierre Auger collaboration, _Features of the energy spectrum of cosmic rays above $2.5{\times}10^{18}$ eV using the Pierre Auger Observatory_, _submitted to Phys. Rev. Lett._ (2020) .
  * [47] Pierre Auger collaboration, _Measurement of the Cosmic Ray Flux near the Second Knee with the Pierre Auger Observatory_ , _PoS_ ICRC2019 (2019) 225.
  * [48] Pierre Auger collaboration, _The Pierre Auger Observatory Upgrade - Preliminary Design Report_ , 1604.03637.
  * [49] Pierre Auger collaboration, _Highlights from the Pierre Auger Observatory (ICRC2019)_ , _PoS_ ICRC2019 (2020) 004 [1909.10791].

$\bullet$

## The Pierre Auger Collaboration

A. Aab75, P. Abreu67, M. Aglietta50,49, J.M. Albury12, I. Allekotte1, A.
Almela8,11, J. Alvarez Castillo63, J. Alvarez-Muñiz74, R. Alves Batista75,
G.A. Anastasi58,49, L. Anchordoqui82, B. Andrada8, S. Andringa67, C. Aramo47,
P.R. Araújo Ferreira39, H. Asorey8, P. Assis67, G. Avila9,10, A.M. Badescu70,
A. Bakalova30, A. Balaceanu68, F. Barbato56,47, R.J. Barreira Luz67, K.H.
Becker35, J.A. Bellido12, C. Berat34, M.E. Bertaina58,49, X. Bertou1, P.L.
Biermannb, P. Billoir33, T. Bister39, J. Biteau32, A. Blanco67, J. Blazek30,
C. Bleve34, M. Boháčová30, D. Boncioli53,43, C. Bonifazi24, L. Bonneau
Arbeletche19, N. Borodai64, A.M. Botti8, J. Bracke, T. Bretz39, A.
Bridgeman37, F.L. Briechle39, P. Buchholz41, A. Bueno73, S. Buitink14, M.
Buscemi54,44, K.S. Caballero-Mora62, L. Caccianiga55,46, L. Calcagni4, A.
Cancio11,8, F. Canfora75,77, I. Caracas35, J.M. Carceller73, R. Caruso54,44,
A. Castellina50,49, F. Catalani17, G. Cataldi45, L. Cazon67, M. Cerda9, J.A.
Chinellato20, K. Choi74, J. Chudoba30, L. Chytka31, R.W. Clay12, A.C. Cobos
Cerutti7, R. Colalillo56,47, A. Coleman88, M.R. Coluccia52,45, R. Conceição67,
A. Condorelli42,43, G. Consolati46,51, F. Contreras9,10, F. Convenga52,45,
C.E. Covault80,h, S. Dasso5,3, K. Daumiller37, B.R. Dawson12, J.A. Day12, R.M.
de Almeida26, J. de Jesús8,37, S.J. de Jong75,77, G. De Mauro75,77, J.R.T. de
Mello Neto24,25, I. De Mitri42,43, J. de Oliveira26, D. de Oliveira Franco20,
V. de Souza18, E. De Vito52,45, J. Debatin36, M. del Río10, O. Deligny32, N.
Dhital64, A. Di Matteo49, M.L. Díaz Castro20, C. Dobrigkeit20, J.C. D’Olivo63,
Q. Dorosti41, R.C. dos Anjos23, M.T. Dova4, J. Ebr30, R. Engel36,37, I.
Epicoco52,45, M. Erdmann39, C.O. Escobarc, A. Etchegoyen8,11, H.
Falcke75,78,77, J. Farmer87, G. Farrar85, A.C. Fauth20, N. Fazzinic, F.
Feldbusch38, F. Fenu58,49, B. Fick84, J.M. Figueira8, A. Filipčič72,71, T.
Fodran75, M.M. Freire6, T. Fujii87,f, A. Fuster8,11, C. Galea75, C.
Galelli55,46, B. García7, A.L. Garcia Vegas39, H. Gemmeke38, F. Gesualdi8,37,
A. Gherghel-Lascu68, P.L. Ghia32, U. Giaccari75, M. Giammarchi46, M. Giller65,
J. Glombitza39, F. Gobbi9, F. Gollan8, G. Golup1, M. Gómez Berisso1, P.F.
Gómez Vitale9,10, J.P. Gongora9, N. González8, I. Goos1,37, D. Góra64, A.
Gorgi50,49, M. Gottowik35, T.D. Grubb12, F. Guarino56,47, G.P. Guedes21, E.
Guido49,58, S. Hahn37,8, R. Halliday80, M.R. Hampel8, P. Hansen4, D. Harari1,
V.M. Harvey12, A. Haungs37, T. Hebbeker39, D. Heck37, G.C. Hill12, C. Hojvatc,
J.R. Hörandel75,77, P. Horvath31, M. Hrabovský31, T. Huege37,14, J.
Hulsman8,37, A. Insolia54,44, P.G. Isar69, J.A. Johnsen81, J. Jurysek30, A.
Kääpä35, K.H. Kampert35, B. Keilhauer37, J. Kemp39, H.O. Klages37, M.
Kleifges38, J. Kleinfeller9, M. Köpke36, G. Kukec Mezek71, B.L. Lago16, D.
LaHurd80, R.G. Lang18, M.A. Leigui de Oliveira22, V. Lenok37, A. Letessier-
Selvon33, I. Lhenry-Yvon32, D. Lo Presti54,44, L. Lopes67, R. López59, R.
Lorek80, Q. Luce36, A. Lucero8, A. Machado Payeras20, M. Malacari87, G.
Mancarella52,45, D. Mandat30, B.C. Manning12, J. Manshanden40, P. Mantschc, S.
Marafico32, A.G. Mariazzi4, I.C. Mariş13, G. Marsella52,45, D. Martello52,45,
H. Martinez18, O. Martínez Bravo59, M. Mastrodicasa53,43, H.J. Mathes37, J.
Matthews83, G. Matthiae57,48, E. Mayotte35, P.O. Mazurc, G. Medina-Tanco63, D.
Melo8, A. Menshikov38, K.-D. Merenda81, S. Michal31, M.I. Micheletti6, L.
Miramonti55,46, D. Mockler13, S. Mollerach1, F. Montanet34, C. Morello50,49,
M. Mostafá86, A.L. Müller8,37, M.A. Muller20,d,24, K. Mulrey14, R. Mussa49, M.
Muzio85, W.M. Namasaka35, L. Nellen63, M. Niculescu-Oglinzanu68, M.
Niechciol41, D. Nitz84,g, D. Nosek29, V. Novotny29, L. Nožka31, A Nucita52,45,
L.A. Núñez28, M. Palatka30, J. Pallotta2, M.P. Panetta52,45, P. Papenbreer35,
G. Parente74, A. Parra59, M. Pech30, F. Pedreira74, J. Pȩkala64, R. Pelayo61,
J. Peña-Rodriguez28, J. Perez Armand19, M. Perlin8,37, L. Perrone52,45, C.
Peters39, S. Petrera42,43, T. Pierog37, M. Pimenta67, V. Pirronello54,44, M.
Platino8, B. Pont75, M. Pothast77,75, P. Privitera87, M. Prouza30, A.
Puyleart84, S. Querchfeld35, J. Rautenberg35, D. Ravignani8, M.
Reininghaus37,8, J. Ridky30, F. Riehn67, M. Risse41, P. Ristori2, V.
Rizi53,43, W. Rodrigues de Carvalho19, J. Rodriguez Rojo9, M.J. Roncoroni8, M.
Roth37, E. Roulet1, A.C. Rovero5, P. Ruehl41, S.J. Saffi12, A. Saftoiu68, F.
Salamida53,43, H. Salazar59, G. Salina48, J.D. Sanabria Gomez28, F. Sánchez8,
E.M. Santos19, E. Santos30, F. Sarazin81, R. Sarmento67, C. Sarmiento-Cano8,
R. Sato9, P. Savina52,45,32, C. Schäfer37, V. Scherini45, H. Schieler37, M.
Schimassek36,8, M. Schimp35, F. Schlüter37,8, D. Schmidt36, O. Scholten76,14,
P. Schovánek30, F.G. Schröder88,37, S. Schröder35, A. Schulz37, S.J. Sciutto4,
M. Scornavacche8,37, R.C. Shellard15, G. Sigl40, G. Silli8,37, O. Sima68,i, R.
Šmída87, P. Sommers86, J.F. Soriano82, J. Souchard34, R. Squartini9, M.
Stadelmaier37,8, D. Stanca68, S. Stanič71, J. Stasielak64, P. Stassi34, A.
Streich36,8, M. Suárez-Durán28, T. Sudholz12, T. Suomijärvi32, A.D.
Supanitsky8, J. Šupík31, Z. Szadkowski66, A. Taboada36, A. Tapia27, C.
Timmermans77,75, O. Tkachenko37, P. Tobiska30, C.J. Todero Peixoto17, B.
Tomé67, G. Torralba Elipe74, A. Travaini9, P. Travnicek30, C. Trimarelli53,43,
M. Trini71, M. Tueros4, R. Ulrich37, M. Unger37, M. Urban39, L. Vaclavek31, M.
Vacula31, J.F. Valdés Galicia63, I. Valiño42,43, L. Valore56,47, A. van
Vliet75, E. Varela59, B. Vargas Cárdenas63, A. Vásquez-Ramírez28, D.
Veberič37, C. Ventura25, I.D. Vergara Quispe4, V. Verzi48, J. Vicha30, L.
Villaseñor59, J. Vink79, S. Vorobiov71, H. Wahlberg4, A.A. Watsona, M.
Weber38, A. Weindl37, L. Wiencke81, H. Wilczyński64, T. Winchen14, M. Wirtz39,
D. Wittkowski35, B. Wundheiler8, A. Yushkov30, O. Zapparrata13, E. Zas74, D.
Zavrtanik71,72, M. Zavrtanik72,71, L. Zehrer71, A. Zepeda60, M. Ziolkowski41,
F. Zuccarello54,44

1

Centro Atómico Bariloche and Instituto Balseiro (CNEA-UNCuyo-CONICET), San
Carlos de Bariloche, Argentina

2

Centro de Investigaciones en Láseres y Aplicaciones, CITEDEF and CONICET,
Villa Martelli, Argentina

3

Departamento de Física and Departamento de Ciencias de la Atmósfera y los
Océanos, FCEyN, Universidad de Buenos Aires and CONICET, Buenos Aires,
Argentina

4

IFLP, Universidad Nacional de La Plata and CONICET, La Plata, Argentina

5

Instituto de Astronomía y Física del Espacio (IAFE, CONICET-UBA), Buenos
Aires, Argentina

6

Instituto de Física de Rosario (IFIR) – CONICET/U.N.R. and Facultad de
Ciencias Bioquímicas y Farmacéuticas U.N.R., Rosario, Argentina

7

Instituto de Tecnologías en Detección y Astropartículas (CNEA, CONICET,
UNSAM), and Universidad Tecnológica Nacional – Facultad Regional Mendoza
(CONICET/CNEA), Mendoza, Argentina

8

Instituto de Tecnologías en Detección y Astropartículas (CNEA, CONICET,
UNSAM), Buenos Aires, Argentina

9

Observatorio Pierre Auger, Malargüe, Argentina

10

Observatorio Pierre Auger and Comisión Nacional de Energía Atómica, Malargüe,
Argentina

11

Universidad Tecnológica Nacional – Facultad Regional Buenos Aires, Buenos
Aires, Argentina

12

University of Adelaide, Adelaide, S.A., Australia

13

Université Libre de Bruxelles (ULB), Brussels, Belgium

14

Vrije Universiteit Brussels, Brussels, Belgium

15

Centro Brasileiro de Pesquisas Fisicas, Rio de Janeiro, RJ, Brazil

16

Centro Federal de Educação Tecnológica Celso Suckow da Fonseca, Nova Friburgo,
Brazil

17

Universidade de São Paulo, Escola de Engenharia de Lorena, Lorena, SP, Brazil

18

Universidade de São Paulo, Instituto de Física de São Carlos, São Carlos, SP,
Brazil

19

Universidade de São Paulo, Instituto de Física, São Paulo, SP, Brazil

20

Universidade Estadual de Campinas, IFGW, Campinas, SP, Brazil

21

Universidade Estadual de Feira de Santana, Feira de Santana, Brazil

22

Universidade Federal do ABC, Santo André, SP, Brazil

23

Universidade Federal do Paraná, Setor Palotina, Palotina, Brazil

24

Universidade Federal do Rio de Janeiro, Instituto de Física, Rio de Janeiro,
RJ, Brazil

25

Universidade Federal do Rio de Janeiro (UFRJ), Observatório do Valongo, Rio de
Janeiro, RJ, Brazil

26

Universidade Federal Fluminense, EEIMVR, Volta Redonda, RJ, Brazil

27

Universidad de Medellín, Medellín, Colombia

28

Universidad Industrial de Santander, Bucaramanga, Colombia

29

Charles University, Faculty of Mathematics and Physics, Institute of Particle
and Nuclear Physics, Prague, Czech Republic

30

Institute of Physics of the Czech Academy of Sciences, Prague, Czech Republic

31

Palacky University, RCPTM, Olomouc, Czech Republic

32

Université Paris-Saclay, CNRS/IN2P3, IJCLab, Orsay, France

33

Laboratoire de Physique Nucléaire et de Hautes Energies (LPNHE), Universités
Paris 6 et Paris 7, CNRS-IN2P3, Paris, France

34

Univ. Grenoble Alpes, CNRS, Grenoble Institute of Engineering Univ. Grenoble
Alpes, LPSC-IN2P3, 38000 Grenoble, France

35

Bergische Universität Wuppertal, Department of Physics, Wuppertal, Germany

36

Karlsruhe Institute of Technology, Institute for Experimental Particle Physics
(ETP), Karlsruhe, Germany

37

Karlsruhe Institute of Technology, Institut für Kernphysik, Karlsruhe, Germany

38

Karlsruhe Institute of Technology, Institut für Prozessdatenverarbeitung und
Elektronik, Karlsruhe, Germany

39

RWTH Aachen University, III. Physikalisches Institut A, Aachen, Germany

40

Universität Hamburg, II. Institut für Theoretische Physik, Hamburg, Germany

41

Universität Siegen, Fachbereich 7 Physik – Experimentelle Teilchenphysik,
Siegen, Germany

42

Gran Sasso Science Institute, L’Aquila, Italy

43

INFN Laboratori Nazionali del Gran Sasso, Assergi (L’Aquila), Italy

44

INFN, Sezione di Catania, Catania, Italy

45

INFN, Sezione di Lecce, Lecce, Italy

46

INFN, Sezione di Milano, Milano, Italy

47

INFN, Sezione di Napoli, Napoli, Italy

48

INFN, Sezione di Roma “Tor Vergata”, Roma, Italy

49

INFN, Sezione di Torino, Torino, Italy

50

Osservatorio Astrofisico di Torino (INAF), Torino, Italy

51

Politecnico di Milano, Dipartimento di Scienze e Tecnologie Aerospaziali ,
Milano, Italy

52

Università del Salento, Dipartimento di Matematica e Fisica “E. De Giorgi”,
Lecce, Italy

53

Università dell’Aquila, Dipartimento di Scienze Fisiche e Chimiche, L’Aquila,
Italy

54

Università di Catania, Dipartimento di Fisica e Astronomia, Catania, Italy

55

Università di Milano, Dipartimento di Fisica, Milano, Italy

56

Università di Napoli “Federico II”, Dipartimento di Fisica “Ettore Pancini”,
Napoli, Italy

57

Università di Roma “Tor Vergata”, Dipartimento di Fisica, Roma, Italy

58

Università Torino, Dipartimento di Fisica, Torino, Italy

59

Benemérita Universidad Autónoma de Puebla, Puebla, México

60

Centro de Investigación y de Estudios Avanzados del IPN (CINVESTAV), México,
D.F., México

61

Unidad Profesional Interdisciplinaria en Ingeniería y Tecnologías Avanzadas
del Instituto Politécnico Nacional (UPIITA-IPN), México, D.F., México

62

Universidad Autónoma de Chiapas, Tuxtla Gutiérrez, Chiapas, México

63

Universidad Nacional Autónoma de México, México, D.F., México

64

Institute of Nuclear Physics PAN, Krakow, Poland

65

University of Łódź, Faculty of Astrophysics, Łódź, Poland

66

University of Łódź, Faculty of High-Energy Astrophysics,Łódź, Poland

67

Laboratório de Instrumentação e Física Experimental de Partículas – LIP and
Instituto Superior Técnico – IST, Universidade de Lisboa – UL, Lisboa,
Portugal

68

“Horia Hulubei” National Institute for Physics and Nuclear Engineering,
Bucharest-Magurele, Romania

69

Institute of Space Science, Bucharest-Magurele, Romania

70

University Politehnica of Bucharest, Bucharest, Romania

71

Center for Astrophysics and Cosmology (CAC), University of Nova Gorica, Nova
Gorica, Slovenia

72

Experimental Particle Physics Department, J. Stefan Institute, Ljubljana,
Slovenia

73

Universidad de Granada and C.A.F.P.E., Granada, Spain

74

Instituto Galego de Física de Altas Enerxías (IGFAE), Universidade de Santiago
de Compostela, Santiago de Compostela, Spain

75

IMAPP, Radboud University Nijmegen, Nijmegen, The Netherlands

76

KVI – Center for Advanced Radiation Technology, University of Groningen,
Groningen, The Netherlands

77

Nationaal Instituut voor Kernfysica en Hoge Energie Fysica (NIKHEF), Science
Park, Amsterdam, The Netherlands

78

Stichting Astronomisch Onderzoek in Nederland (ASTRON), Dwingeloo, The
Netherlands

79

Universiteit van Amsterdam, Faculty of Science, Amsterdam, The Netherlands

80

Case Western Reserve University, Cleveland, OH, USA

81

Colorado School of Mines, Golden, CO, USA

82

Department of Physics and Astronomy, Lehman College, City University of New
York, Bronx, NY, USA

83

Louisiana State University, Baton Rouge, LA, USA

84

Michigan Technological University, Houghton, MI, USA

85

New York University, New York, NY, USA

86

Pennsylvania State University, University Park, PA, USA

87

University of Chicago, Enrico Fermi Institute, Chicago, IL, USA

88

University of Delaware, Department of Physics and Astronomy, Bartol Research
Institute, Newark, DE, USA

—–

a

School of Physics and Astronomy, University of Leeds, Leeds, United Kingdom

b

Max-Planck-Institut für Radioastronomie, Bonn, Germany

c

Fermi National Accelerator Laboratory, USA

d

also at Universidade Federal de Alfenas, Poços de Caldas, Brazil

e

Colorado State University, Fort Collins, CO, USA

f

now at Hakubi Center for Advanced Research and Graduate School of Science,
Kyoto University, Kyoto, Japan

g

also at Karlsruhe Institute of Technology, Karlsruhe, Germany

h

also at Radboud Universtiy Nijmegen, Nijmegen, The Netherlands

i

also at University of Bucharest, Physics Department, Bucharest, Romania