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

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.


Introduction
The intensity of cosmic rays with energies above ∼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 ∼10 km 2 , 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 ∼1 per km 2 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 θ < 60 • , 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.

The surface detector
The SD consists of more than 1600 water-Cherenkov stations arranged on a isometric triangular grid covering 3000 km 2 . 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 ∼8, stations. Within the boundary of the 1500 m array, there are two denser regions: one extending over 28 km 2 within which the stations are spaced 750 m apart from one another, and the other covering 2 km 2 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 m 2 . 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 ∼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.
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 NT5, 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 events1 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/s, due mainly to single uncorrelated atmospheric muons [17], down to about 3×10 −5 /s for high quality physics Example of a shower event detected by the SD, which is reconstructed with θ = 55.2 • 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.
events. Consequently, a rate of events across the complete array of about 0.05/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.

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 traces2 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.

Calibration
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 ∼240 MeV of energy deposited due to ionization losses of the passing muon. [18]. This conversion factor, Q peak VEM , is evaluated online by constantly measuring the properties of particles entering the detector volume at a rate of ∼3/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 peak 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.

Determination of the baseline
The first step in the signal processing is the estimation of the baseline of the high-gain and lowgain 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 µs) that vary in amplitude by less thanȳ ± σ b , whereȳ is the mean of the baseline segment and σ b is the tolerance needed to account for the possible noise in the baseline. σ b is set to 2 ADC counts. In rare cases (<0.2%), when the PMT in question is noisy, the tolerance σ 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ȳ of that baseline segment. Then, the parts between the

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 peak 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. Fragment B must begin less than 500 ns (20 time bins) + t A after the end of fragment A, where t A is the amount of time in fragment A containing signal.
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 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 PMTs3) 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.

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 ∼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. For lower energy showers, this occurs most of the time when the distance of this station is closer than 100 m to the shower axis. For the highest energy showers, this typically already occurs when a station is closer than 500 m. In the case of low-gain saturation, the equivalent integrated signal is typically greater than ∼1000 VEM. The fraction of such events with a saturated station4, declared saturated events, is less than 20% for the lowest energy showers but increases to ∼80% for the highest energies. 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 peak 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 peak VEM factor and the high-to-low-gain ratio, which is also constantly monitored by the online station calibration. 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.

Uncertainties in the measurement of the signal
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 ∼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 (∼10 ns), the 40 MHz FADC sampling (∼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 start is estimated as 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/¯ (θ), where¯ (θ) is the mean relative track-length of through-going particles crossing the detector when arriving with a zenith angle θ. Due to the cylindrical shape of the water-Cherenkov tank, 1/¯ (θ) = cos θ + (2H/πR) sin θ, 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, 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 ∼1000 VEM up to 60% at ∼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.

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â of the shower is anchored at the location of a barycenter ì x b (i.e. the signal-weighted5 center-of-mass of stations in an event), which later6 serves as a first estimate of the impact position of the shower core at the ground. The assumed shower core travels in the −â direction, impacting the ground at barycenter ì x b at the time t 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 sh ( ì x) when the shower plane passes through some arbitrary point ì x (e.g. a station on the ground) may be inferred through a simple projection onto the shower axis as 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.
We designate as accidental all stations where the onset of the signal is more than 1 µs prior to (early) or 2 µ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.

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â and the position of the impact point of the shower core on the ground ì x 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.

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â of the shower and the impact position ì x 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 −â direction and intersecting with the ground at the impact point of the shower core ì x c at the time t c . Thus four parameters are fitted: the two directional cosines7 u and v of the unit vectorâ = (u, v, , the time t c , and the radius of curvature of the shower front R o . Using this model for the arrival time of the shower front, t sh ( ì x), the shower geometry is fitted to the start-time of the signals t i in each triggered station i located at ì x i using the M framework [23]. The function which is minimized is the sum of the squares of the differences in predicted and measured start times, where σ 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 ì 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 ì x b from Section 4 is chosen as a suitable approximation for the impact point ì x c , while the result of the plane fit from Eq. (4.1) is used for the shower axisâ.
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 −â. Compared to the plane front in Eq. (4.1), the particles are thus delayed proportionally to where R o is the constant radius of the spherical front and where is the perpendicular distance of point ì x from the shower axisâ. To keep t sh ( ì x) linear in the fitted parameters, only terms up to the second order in the r/R o expansion are used. With a curvature parameter k o = 1/2R o , we can then express the shower timing as which is clearly a paraboloidal extension8 of Eq. (4.1).
The Observer reconstruction approximates the shower development as starting at time t o from a virtual point of origin ì x o (see Fig. 5-right) and propagating towards the ground in the shape of a sphere, concentrically inflating with the speed of light9. The arrival time of such a shower front at a point ì x is thus 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 ì x c or the shower axis. Once ì x c is determined at a later point, the shower axis may be obtained as a normalized direction towards ì 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 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, as the radius of the shower at the time of passage t c = t o + R o /c through the impact point ì x 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 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 • , 0.83 • , and 1.14 • , respectively. The small, non-zero difference of less than 0.1 • observed in ∆θ will be addressed in Section 5.2.1.

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 LDF (r) such that S(r) = S(r opt ) f LDF (r), where S(r opt ) is the shower-size estimator and f LDF (r) is normalized so that f LDF (r opt ) ≡ 1. The optimal value for distance r opt has been chosen so that the variability of the shower-size estimator S(r opt ) with respect to the aforementioned shower-to-shower fluctuations is minimized. The value almost exclusively10 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 opt ≈ 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).

Lateral distribution function
Due to the lack of analytical solutions for the hadronic-cascade equations, functional forms for f 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, where ρ = ln(r/r opt ) and β, γ are the average slope parameters of the LDF. For smaller axial distances, r < r c where r c = 300 m, a tangential log-log linear function β ρ + γ(2ρ − ρ c )ρ c is used where ρ c = ln(r c /r opt ).
The Observer reconstruction instead uses a slightly modified NKG function [26][27][28] of the form with a fixed r s = 700 m. The average of the slope parameter β is determined in both cases in a data-driven way, by fitting β for the subset of events with a multiplicity and spatial distribution of stations providing a sufficient lever arm. In these events, the slope β is fitted given that there are at least 2, 3, or more stations within a radius of 400 < r/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 β and γ, it is more difficult to determine the average of the slope parameter γ in a data-driven way in the case of both LDF functions. In Eqs. (5.6) and (5.7), the parameter γ 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 γ except the station search radius is changed to 1000 < r/m < 2000. In the Herald reconstruction, the parameter γ has been determined from Monte-Carlo studies. Finally, both parameters β and γ 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 magnitude of these effects depends on the zenith angle and state of development of the shower. Using a rotationally-symmetric LDF approximation results in a systematic shift of the core by ∼40 m in the direction opposite to the shower development (upstream) and a corresponding shift in the axis of less than <0.1 • , although the exact values for both depend on the zenith angle. In the Herald reconstruction, the signal asymmetry is corrected with an unbiasing model for the measured signals S symm (r) = S(r) 1 + α(r) cos ζ , (5.8) where ζ is the azimuthal angle of the station (as measured in the shower plane and not in the ground plane), and where ζ = 0 corresponds to the upstream direction, and an amplitude α(r), the radial dependence of which has been obtained from simulations. These corrections are only important in very specialized studies as such small shifts of the shower axis are well within the typical uncertainties of the reconstruction and ultimately have a <0.2% effect on the value of S(1000) and a <0.2 • effect on the axial direction. 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 ∆θ < 0.1 • as seen in Fig. 6-right.

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 L, composed as the product of probabilities P over the coordinates of the shower impact point ì x c and the size S(1000), given the observed signal sizes S i in the stations, located at ì x i . We are thus maximizing the log-likelihood The two reconstructions use different models for calculation of the likelihood P. In the Observer framework, P is a product of contributions from: (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 N for large signals is integrated over all possible values larger than S i to get an estimate F 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 trig (S(r i )), which is derived from data [14]. The complementary probability for non-triggered stations is thus modeled as F non (S(r i )) = 1 − p trig (S(r i )).

This results in a likelihood function
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 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 β is always fitted in the Herald reconstruction.
In the Observer reconstruction, the parameter β can be set as a free parameter for the same set of showers. A comparison of the shower sizes obtained with β as a free or fixed parameter of the LDF in Eq. (5.7) is shown in Fig. 7. The impact of freeing or fixing β amounts to a less than 3% (4%) bias in the reconstructed value of S(1000) for non-saturated (saturated) events and, therefore, fitting of β is not required to obtain an unbiased estimate of the shower size [6,29].
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 | 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 ∆x = x LLP − x NKG and ∆y = y LLP − y NKG . These coordinates have been computed in the reference frame of the Observer shower and are oriented so that the positive/negative ∆y axis is in the down/upstream direction. freedom, any parameters which were previously fixed are at this stage allowed to vary. The global log-likelihood function − 2 ln L global = χ 2 − 2 ln L (5.12) is minimized, where the two terms are given in Eqs. (5.1) and (5.11). 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 |∆ ì x 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.
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) Herald vs. S(1000) 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.

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 ∼10 ns in the GPS timing system [13] induces an angular uncertainty11 of less than ∼0.1 • . 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 directionŝ a 1 andâ 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 η is obtained from sin η = |â 1 ×â 2 |. For η 1 and in a normal approximation with one-dimensional variance σ 2 , η follows a distribution dN/dη = R(η; σ) for the case of (a) and dN/dη = R(η; √ 2σ) for the case of (b), where R(r; σ) = (r/σ 2 ) exp(−r 2 /2σ 2 ) is the Rayleigh distribution. We define the angular resolution (AR) as the angle at which the value of the cumulative distribution of η reaches 68.3%, this being in turn equivalent to AR ≈ 1.52σ or 1.52σ/ √ 2 for the two cases, respectively [30]. For the case (c) we assume that the two measurements σ 1 and σ 2 are uncorrelated and thus can be added in quadrature, i.e. σ 2 tot = σ 2 1 + σ 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 simulations12.

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 subsample 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,â SD and a hyb , the angular resolution AR SD:hyb is extracted from the angular difference η SD:hyb defined by sin η SD:hyb = |â SD ×â hyb |. Given the essential independence of the SD and hybrid reconstructions, the SD angular resolution AR SD can be derived as AR 2 SD = AR 2 SD:hyb − AR 2 hyb , where AR 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 separate13 hybrid reconstruction is performed for each FD site. Computing the angular difference η hyb between these arrival directionsâ hyb 1 andâ hyb 2 , the resolution of the hybrid reconstruction of this particular selection of stereo events is estimated to be AR hyb = (1.07 ± 0.05) • assuming that resolutions forâ hyb 1 andâ hyb 2 are the same14. With AR hyb at hand, AR 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.

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 ∼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, 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. 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 events15, 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 AR SD is obtained from the distribution of the angular difference η tw , defined as sin η tw = |â tw 1 ×â tw 2 |. The results of this procedure are shown in Fig. 11.

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 η sub i j from sin η sub i j = |â sub i ×â sub j |. Here, we assume that all η sub i j are independent, obtaining 4741 angles for the calculation of the estimate of angular resolution AR SD . The results are reported in Fig. 11 (circles).   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).

Estimation using simulations
Using simulations of extensive air-showers, the angular resolution is obtained as a function of zenith angle θ 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 • and 60 • . 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â rec are compared with the true simulated axisâ of the shower by computing the angular difference η from sin η = |â ×â 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 θ and as an exponential function in lg(S(1000)/VEM). The results are reported in Fig. 11 (lines and bands).

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 AR 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 AR hyb and then AR 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. Acknowledging that simulations do not emulate all experimental aspects [35], the estimate of the angular resolution derived therefrom may be optimistic. Therefore, 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 • and even approaches 0.7 • 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 ∼100 m at the lowest energies to ∼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.

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, σ stat (S(1000)), are directly related to the number of triggered stations and the uncertainties in their signals, the systematic uncertainties, σ 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.

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 β from Eqs. (5.6) and (5.7). As a first step, we use the same subset of events already used in the determination of β in Section 5.2.1, i.e. events for which β can be reliably fitted. The resulting standard deviation σ β of the difference between the fixed and free β is parametrized as a function of S(1000). In a second step, the systematic uncertainties are computed by propagating σ β (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  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 ∼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.

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 δ tw = ∆S(1000) tw / √ 2 S(1000) tw is computed, where ∆S(1000) tw and S(1000) tw are the difference and the average, respectively, of the two reconstructed shower sizes S(1000) tw 1 and S(1000) tw 2 . In each bin in lg(S(1000)/VEM) or in sin 2 θ, the standard deviation σ(δ tw ) is calculated. The standard deviations of these distributions are plotted in Fig. 14 (as filled squares).

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 δ sub = ∆S(1000) sub / √ 2 S(1000) sub is derived, where ∆S(1000) sub and S(1000) sub are the difference and the average of the reconstructed shower sizes S(1000) sub i and S(1000) sub j , respectively. For each pair of sub-arrays (i, j), the standard deviations σ(δ sub ) i j are calculated and then averaged over all pairs (i, j) in different bins of lg(S(1000)/VEM) and sin 2 θ. The results are reported in Fig. 14 (as open squares).

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 • , 12 • , 22 • , 32 • , 38 • , 48 • , and 56 • . 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) 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) true and S(1000) 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 δ sim = ∆S(1000) sim /S(1000) true where ∆S(1000) sim = S(1000) rec − S(1000) true . The standard deviation σ(δ sim ) of the relative differences is reported in Fig. 13. The results show a total uncertainty in S(1000) rec , which decreases from 18% at the smallest shower sizes to nearly 7% at the largest values of S(1000) 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 σ(δ sim ) with zenith angle.

Comparison of results
The statistical and systematic uncertainties of the reconstructed S(1000) for all measured data are added in quadrature as σ 2 (S(1000) rec ) = σ 2 syst (S(1000) rec ) + σ 2 stat (S(1000) rec ), where σ syst is from Fig. 7. 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 ∼6% at lg(S(1000)/VEM) ≈ 2.3 to ∼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 Figure 14. Relative total uncertainty of the shower-size reconstruction σ(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 impact of the resolution in shower size on the uncertainties in the estimate of the primary energy is described in detail in [10].

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 [36] 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 [37]. 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 gm cos k θ sin 2 ξ, where A gm and k are model parameters estimated from the Monte-Carlo data, and ξ 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,38], 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 • (median of all events), S 38 = S(1000)/ f CIC (θ). 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 [39] by the nearlycalorimetric energy measurement of the FD using hybrid events (see Section 6.1.1). The power-law calibration curve E = A (S 38 /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.

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 • and approaches 0.5 • 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 ∼15%, for the smallest shower sizes down to ∼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 [40].
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, 41], multimessenger [42] and mass-composition studies [43,44], and the cosmic-ray energy spectrum [10,45,46].
Looking to the future, the surface detector of the Pierre Auger Observatory is currently undergoing a large-scale upgrade called AugerPrime [47,48], with the placement of a 3.8 m 2 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. [10] P A collaboration, A measurement of the cosmic ray energy spectrum above 2.5×10 18 eV using the Pierre Auger Observatory, submitted to Phys. Rev. D (2020) .