{'fraction_non_alphanumeric': 0.03057392862064214, 'fraction_numerical': 0.016139589361995315, 'mean_word_length': 4.634960691060856, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Cosmic Ray Composition and Energy Spectrum from 1-30 PeV Using the 40-String Configuration of IceTop and IceCube

The mass composition of high energy cosmic rays depends on their production, acceleration, and propagation. The study of cosmic ray composition can therefore reveal hints of the origin of these particles. At the South Pole, the IceCube Neutrino Observatory is capable of measuring two components of cosmic ray air showers in coincidence: the electromagnetic component at high altitude (2835 m) using the IceTop surface array, and the muonic component above ~1 TeV using the IceCube array. This unique detector arrangement provides an opportunity for precision measurements of the cosmic ray energy spectrum and composition in the region of the knee and beyond. We present the results of a neural network analysis technique to study the cosmic ray composition and the energy spectrum from 1 PeV to 30 PeV using data recorded using the 40-string/40-station configuration of the IceCube Neutrino Observatory.


Introduction
The flux of cosmic rays at Earth is known to follow a steep power-law spectrum over a large energy range. The index of this spectrum is approximately constant at energies lower than about 3 PeV, where the spectrum steepens in a feature known as the "knee". A further kink in the spectrum, known as the "ankle", occurs around 1 EeV where the spectrum becomes flatter again. The origins of these spectral changes are still uncertain. Currently, the most popular model predicts cosmic ray acceleration in shock fronts via the first order Fermi mechanism [1]. More specifically, at energies up to ∼ 10 17 eV, the source of this acceleration mechanism is often attributed to supernova remnants; a cut-off energy which depends upon nuclear charge (Z) of the particle accelerated at the source could be responsible for a mass-dependent knee; the ankle is then attributed to cosmic rays from extragalactic sources such as gamma-ray bursts or active galactic nuclei [2,3,4].
Various underlying source, acceleration, and propagation models, though tuned to predict similar energy spectra, differ considerably in energy-dependent composition in the region between the knee and the ankle [5]. This dependence on primary mass implies that a precise measurement of the composition of cosmic rays would provide important clues as to the origins of these particles. However, due to the decreasing flux of cosmic rays with increasing energy, measurements with high statistics are difficult to collect with satellite or balloon-based detectors above 100 TeV. Groundbased detectors can observe cosmic ray air showers above 100 TeV indirectly with high statistics. Indirect measurements of composition involve a close examination of the extensive air shower produced by a primary cosmic ray particle colliding with Earth's atmosphere. By using information from more than one component of the shower, such as the electromagnetic and muonic components, the energy and composition can be obtained for primary particles with much higher energies than those currently measurable through direct detection techniques.
At the South Pole, the IceCube Neutrino Observatory is sensitive to air showers with energies from below the knee to the ankle region of the energy spectrum. This region covers the predicted transition from galactic to extragalactic cosmic ray sources. The IceTop surface array measures a combination of the electromagnetic and the low energy muonic components of the cosmic ray shower, while the IceCube array measures the bundle of high-energy muons (>300 GeV) deep under the surface of the ice. In the following analysis, data from these two components are combined for a coincident composition measurement.

The IceCube Neutrino Observatory
The IceCube Neutrino Observatory consists of two parts: IceTop, a surface air shower array [6], and IceCube, a muon and neutrino telescope installed deep in the ice. These detectors are the successors to the SPASE [7] and AMANDA [8] experiments. Each array is comprised of light sensors called Digital Optical Modules (DOMs) [9], which detect Cherenkov photons emitted by relativistic charged particles passing through ice. Each DOM is a spherical, pressure-resistant glass shell containing a 25 cm diameter Hamamatsu photomultiplier tube (PMT), a mu-metal grid for magnetic shielding of the PMT, and electronics for operation and control of the PMT as well as amplification, digitization, filtering, and calibration. The colors represent the timing of the hits (red is earliest, blue is latest), and the size of the sphere around a DOM represents the amplitude of light seen by that DOM. In this large event, a big dust layer can clearly be distinguished as a "waist" in the amplitudes just over half-way down the IceCube array.
In IceCube, DOMs are frozen into the ice along strings which are placed in a 125 m triangular grid formation. The DOM's are vertically spaced 17 m apart, at depths from 1450 m to 2450 m below the surface, as shown in Fig. 1(a). The direction of muons (either from cosmic ray air showers above the surface, or neutrino interactions within the ice or bedrock) can be reconstructed from the pattern (position and timing) of hit DOMs. This feature is demonstrated in Fig. 1(b), which shows a coincident IceTop/IceCube event from the 40-string/40-station configuration of 2008. The minimum energy at the surface required for a muon to penetrate through the ice to the top of the IceCube array is around 300 GeV, while the threshold for the muon to pass through the whole detector volume is around 500 GeV.
Each IceCube string is associated with an IceTop "station", which consists of two cylindrical ice Cherenkov tanks separated by 10 m. Each tank has an inner diameter of 1.82 m, is 1.3 m high, and contains two down-facing DOMs, with center-to-center spacing of 58 cm, frozen in optically clear ice 90 cm in depth. The PMT in each DOM has an adjustable gain; to increase the dynamic range of the detector, one DOM in each tank is operated at low gain (LG) while its partner is operated at high gain (HG). The tanks are lined with a diffusely reflective coating of either Tyvek or, in most cases, zirconium fused polyethylene, and the surface of the ice is covered with perlite. IceTop measures mainly the electromagnetic component of incoming cosmic ray air showers. The signal response of a particle passing through an IceTop tank is measured in photoelectrons (PEs); however, each tank must be calibrated to obtain a uniform measurement. Therefore, the amount of charge deposited by mainly vertical 1 -10 GeV surface muons passing through a given tank is defined in terms of a common unit, called a "vertical equivalent muon", or VEM. A calibration process is then performed by comparing the charge spectra of each DOM -which shows a clear muon peak -in data and simulation. The altitude of the surface array is 2835 m, which corresponds to an atmospheric depth of ∼680 g/cm 2 . For showers with tens of PeV energies, this is just below the proton shower maximum; therefore, an excellent energy resolution is expected.
Construction of the IceCube Neutrino Observatory began at the geographic South Pole in the 2004-05 austral summer [10] and was completed during the 2010-11 austral summer. This work used the configuration of the detectors operational from April 2008 to May 2009, which consisted of 40 IceCube strings (2400 DOMs) and 40 IceTop stations (160 DOMs). Figure 1(a) shows the final detector configuration of 86 strings, while Fig. 1(b) is an example data event from the 40-string configuration studied in this work.

Data
This analysis uses data from August 2008, when the detector was in its 40-string/40-station configuration, for an overall detector livetime of 29.78 days. A study of the relationship between the electromagnetic and muonic air shower observables (K 70 and S 125 , as discussed in Section 4) throughout the full year of data from the 2008 configuration revealed significant seasonal variations in the relationship between the observables used for this analysis [11]. In particular, the effective temperature of the atmosphere (which is stable from day to day but varies dramatically between summer and winter) affects the production of muons in the upper atmosphere, changing the measured K 70 . Limiting the data set to August 2008, a mid-winter month which corresponds to our simulated data (described in Section 3.2), effectively mitigates these fluctuations.
Local coincidence (LC) and trigger requirements are the first steps in the data acquisition system of the detectors [9]. The LC requirement is satisfied in IceCube when two neighboring or next-tonearest-neighboring DOMs on a string pass a signal threshold of 0.25 photoelectrons (PE) within 1 µs. The LC requirement is satisfied in IceTop when either the HG DOMs in both tanks of a station, or a HG DOM and the LG DOM from the neighboring tank in the station have passed the signal threshold of about 20 PE within 1 µs. An event has triggered a detector when a certain number of DOMs record LC signals within a 5 µs sliding time window: eight DOMs for IceCube, and six DOMs for IceTop. Additionally, a "coincident event" is one where both detectors have triggered.

Simulation
Monte Carlo simulated events were produced for this analysis using the CORSIKA air shower generator [12] with the SIBYLL-2.1/FLUKA-2008 hadronic interaction models [13,14]. Five particle species (proton, helium, oxygen, silicon, and iron) were generated with an E −1 spectrum from 10 TeV to 50 PeV. The showers were generated uniformly over all azimuths and with a sin(θ)cos(θ) distribution in the zenith range from 0 • to 65 • , which is steeper than that reached by reconstructible coincident events. For each species, 3000 showers were simulated per third of a decade in energy. Each shower was oversampled 100 times over a circle with radius of 1200 m leading to a total of 16.5 million generated events. The atmospheric model chosen corresponds to the austral winter months at the South Pole. For initial comparison with experimental data, the simulation is reweighted, independent of primary mass, to an E −2.7 spectrum at energies below a knee at 3 PeV and an E −3.0 spectrum above the knee. The IceTop detector is simulated using parameterizations of the tank response obtained from a complete Geant4 [15,16] simulation. The muons are then propagated through the ice to the depth of the IceCube array using the muon propagator MuonMonteCarlo [17] which accounts for continuous and stochastic energy losses. The Cherenkov photons from the muon bundles passing through the volume of the array are simulated using the software package Photonics [18] which models the full structure of the ice properties according to the standard ice model used for the 40-string configuration of IceCube (known as the AHA ice model) [19]. This is followed by a simulation of each aspect of the DOM electronics and the trigger. The DOM signals are then processed in the same way as experimental data.
In IceCube, the parameterization for the light yield from muon bundles in the Photonics software has recently been improved. This improvement revealed an energy-dependent offset between our data and our standard simulation. Thus, a small set of simulation was generated using the new software and this discrepancy was parameterized using a comparison of the two simulation sets. This parameterization has been applied as a correction to the experimental data from this point forward.

Reconstruction with IceTop
Data collected by IceTop are analyzed to reconstruct the core position, direction and size of a cosmic ray air shower [20,21]. These parameters are determined by comparing the detected signal locations, charges, and times from hit stations (as well as the locations of unhit stations) to what is expected from a cosmic ray air shower, according to a likelihood function: where L s is the likelihood of signal charges from a lateral distribution function, L 0 is the likelihood of unhit stations, and L t is the likelihood of signal times from an expected timing profile. All three are described in detail in Ref. [21], but a brief overview is given here. The first term, L s , represents the likelihood of the array responding with measured signals in hit tanks S i , given a set of expected signals for each of those hit tanks S fit i . Signals are assumed to have normal distributions in log 10 (S i ), and are expressed in units of VEM (as discussed in Section 2).
The signal expectation value, S fit i , is a function of the perpendicular distance from the shower axis, r, and is given by a lateral distribution function which has been derived from simulations: where κ is a constant optimized from simulations, and β is fit as a free parameter to each event.
The shower size, S ref , is the signal expectation at a reference distance, R ref , from the shower axis.
The reference distance is chosen as the point in the fit which is the most stable or robust: multiple studies of the fit stability led to a choice of 125 m for this quantity [21]. Figure 2 shows an example of the distribution of charges for one shower from the 2008 IceTop configuration with 40 stations. The superimposed curve is the lateral distribution function of Eq. 2 with the best-fit S 125 for this particular event.
x( m) The amount of signal observed by each tank in a given station is denoted by the size of the hemisphere closest to it. Right: The lateral distribution of these hits, color-coded for timing in the same scheme as for the event display. The expected lateral distribution of signals is shown in black; S 125 is this function evaluated at 125 m from the core of the shower; in this case, 12.2 V EM .
The second likelihood term in Eq. (1) accounts for all stations j with the tanks (j, A) and (j, B) that did not trigger: where P hit j,A(B) is the probability that the "A" tank (or the "B" tank) within station j did not trigger, for a given signal expectation value S fit .
The third term of Eq. (1) describes the probability of measuring the observed signal arrival times, t i , given expected signal arrival times, t fit i . Each arrival time is assumed to be normallydistributed with standard deviations σ ti (r i ), depending on the radial distance, r i , of the tank to the shower core.
Studies of experimental data reveal that the shower front is not a flat plane. Rather, the expected arrival times of the signals can be described as those expected from a flat plane, modified by a "curvature term" which is symmetric around the shower axis and a function of radial distance. The expected signal arrival time at a tank with position x is thus parametrized as: Here, t 0 is the arrival time of the shower core at the surface of the ice (which is the surveyed level of the IceTop tanks), x c is the position of the shower core on the ground and n is the unit vector in the arrival direction of the shower. The last two terms in Eq. 4 (the parabola and the Gaussian) together are the "curvature term". The functional form of the curvature term and the values of the parameters a, b, and σ were obtained from a study of the time residuals in experimental data. Thus, in this likelihood function there are seven total possible free parameters: S 125 , β, and the track's core position (x c ,y c ) and direction and core time (θ, φ, t 0 ).

IceTop Fit Procedure.
Events in IceTop are initially reconstructed with two fast "first guess" algorithms. First, the core position is estimated using the center of gravity (COG) of the tanks which were hit, weighted with the square root of the charges: Then, the direction of the shower is estimated by fitting a plane (without curvature terms) to the arrival times of the signals. Both first guess core position and track direction together provide the seed track for the more advanced lateral distribution fit involving all three likelihood terms described in the previous section. Multiple lateral distribution fits are performed in succession for best results. At each step a full minimization is performed. First, the shower direction is fixed and only the core position, shower size, and β are allowed to vary as free parameters. Next, since Eq. (2) can exhibit unphysical behavior at small values of r, a second fit is made (again with shower direction fixed) in which all tanks closer than 11 m to the reconstructed shower core are excluded from the fit. If an 11 m circle around the new core location includes tanks which were not excluded before (which happens rarely), they are added to the list of excluded tanks and this second step is repeated until no new tanks are added to the list. In a third step, the track direction is allowed to vary, as well as the core location and shower variables with smaller stepsizes. Finally, the track direction is fixed again and a final fine-tuning of just the core location and shower variables (with small stepsizes) is performed.

Reconstruction with IceCube
Bundles of muons that reach the IceCube detector can have muon multiplicities ranging from 1 to more than 1000 at the highest energies, and an average lateral size of 20-30 meters. When these muons pass through the ice at a speed greater than the speed of light in the medium, they produce Cherenkov photons which can be detected by the DOMs of the IceCube array. Many techniques for reconstructing the direction and energy of single particle tracks exist, which vary track parameters until the likelihood of producing the measured amplitudes and arrival times of photons at the DOMs is maximized [22]. But when the source of the photons is a muon bundle rather than a single muon track, reconstruction can be improved with a likelihood function specific to muon bundles, which was developed using the prototype detectors SPASE-2 and AMANDA-B10 [23,24] and has been updated for use with IceCube. The expected signal for a given DOM can be described by an in-ice lateral distribution function (LDF), which depends upon its perpendicular distance, d, from the center of the muon bundle, as well as the slant depth, X, of the muon bundle at the point where it passes closest to the DOM, and the effective attenuation length, λ eff , of the ice surrounding the DOM. The overall expected signal, ξ, is where N µ,depth (X) is a "range-out" function which describes the number of muons as a function of slant depth: based on a model for continuous and stochastic energy loss of muons (parametrized by constants a and b respectively), combined with the muon energy spectrum index γ µ . As the bundle penetrates deeper, correspondingly less light should be expected at the deeper DOMs as muons "range out" (as shown conceptually in Fig. 3(a)): The effective attenuation length λ eff (z) is a combination of both scattering and absorption effects in the ice. Both of these properties vary with DOM location, due to horizontal dust layers. Therefore, λ eff is a function of each DOM's depth in the ice, z, and is parametrized as an average, or "bulk," λ 0 times a z-dependent correction factor: The correction factor c ice is based on scattering length data measured at many different depths. For DOMs which are far from the muon bundle (d > 80 m), light is likely to have traversed multiple dust layers before reaching the DOM. Thus, in Eq. (6) the average effective attenuation length of the bulk ice, λ 0 , (with no dust correction) is used for these DOMs, and the normalization, N , is modified by a factor chosen to ensure continuity of the function over all distances.
Therefore, the expected signal ξ for each DOM can be computed using Eq. (6), which already includes the effect of dust layers and muon range-out according to that DOM's position. This expected signal is then compared with the hits recorded by the DOMs using a likelihood method, as in Fig. 3(b). The overall normalization of the function, N , the average attenuation length, λ 0 , and the direction of the central track, are fitted as free parameters. Although λ 0 is known to be around 25 m [18], this length is fit as a free parameter because minor errors in the track reconstruction can change this fitted slope of the LDF of photons in IceCube, as shown in Fig. 4, for which the direction of a single simulated 12.3 PeV track was varied from the true direction by up to 1.1 • to show the effect of possible misreconstructions.
After the fitted track is found, a parameter is needed to quantify the overall amount of light (and number of muons) in the ice. Like the S 125 parameter (which is the IceTop lateral distribution function fit at a reference distance of 125 m), the K 70 parameter is the IceCube lateral distribution function of Eq. 6 evaluated at a reference slant depth of 1950 m, with a reference dust correction (c ice = 1), and a reference perpendicular distance D ref of 70 m: where A (the overall normalization) and λ 0 (the average attenuation length) are free parameters which were fit. The reference distance of 70 m was chosen to provide the most stable measurement under errors in the reconstructed track direction. A number of tests were performed to choose this reference point: one example is shown in Fig. 4, where the small spread around 70 m indicates the stability of this choice in reference distance.

Full Reconstruction: Multiple Iterations
The IceCube track reconstruction algorithm described in the previous section can clearly be used to estimate the multiplicity of the muon bundle, since this should scale with the overall normalization A of the LDF of Eq. 6. But since the track direction affects the expected ξ's at all the DOMs (by changing their perpendicular distances d and slant depths X), the same likelihood function can also be used to find the track direction. The core position and direction from the surface reconstruction are used as the "seed track" for the fit; the core position is anchored at the surface, while the direction of the track is allowed to vary in search of the best likelihood. Because the DOMs participating in the fit are far from where the track is fixed at the surface, a long "leverarm" gives improved angular resolution. Once this IceCube reconstruction has been performed, the resulting fitted track can be used as an improved seed for a second iteration of the surface reconstruction, leading to a more accurate core location. The improved surface reconstruction is then passed back as a seed to the IceCube reconstruction one last time. Thus, two iterations of both the lateral fit reconstruction and the muon bundle reconstruction are performed (IceTop → IceCube → IceTop → IceCube). Figure 5 shows the improvement to both the core location and angular reconstruction achieved through this iterative procedure. The two observables which will be used for composition (K 70 , which is related to the number of muons in IceCube, and S 125 , which is related to the number of electrons and photons in IceTop) are also drawn from the last iteration of fits, where the improved resolution allows for accuracy in measuring those parameters.

Snow on IceTop
During the course of this analysis, it was discovered that the presence of research buildings in the vicinity of the IceTop array have caused the snow to drift unevenly over the IceTop tanks. Some tanks are covered by more than 1 m of drifted snow, while others have a mere sprinkling: in particular, the part of the detector constructed first and nearest to structures (the "older half") tends to have a great deal more snow than the part constructed more recently (the "newer half"). Therefore, the signals in the tanks of the older half of the array are attenuated, and this effect is not present in the standard simulation used for this analysis. However, a small sample of simulated proton and iron showers with snow attenuation was used to study this effect and to develop a method of accounting for it in event reconstruction. For each tank, the signal expectation is lowered by a factor which depends upon the measured snow depth above each tank, as well as the snow density, which was measured to be about 0.4 g/cm 3 . Thus, a correction has been made to the data at this reconstruction stage so that it corresponds to the standard simulation described above. The error introduced by this correction will be addressed in Section 9.

Event Selection
Experimentally observed and simulated events have been selected in multiple steps in order to ensure a reliable reconstruction of direction and core location, to retain as many events as possible, and to avoid introduction of any bias by primary mass and energy: : The angular resolution (left) and distance between true and reconstructed core position (right) for simulated data. The "ideal limit" (black) is the maximum resolution that could be achieved by the IceCube track reconstruction (left) and by the IceTop core reconstruction (right), obtained by using the fixed true core location and true direction respectively. After two iterations of both IceTop and IceCube reconstructions, the optimal core resolution is already found, so the direction reconstruction is also limited to two iterations. The second iteration through the two reconstructions clearly provides an improvement to both the angular reconstruction of the track and the shower core position at the surface [25].  : The angular resolution (left) and core position resolution (right) between the true track and the final reconstructed track in proton (red dashed) and iron (blue solid) simulation after the event selection described in Section 5. An angular resolution of 0.5 • and a core resolution of 9 m are achieved (where resolution is defined such that 68%, or one gaussian sigma, of all events are contained within that range).
• Step 1: Initial filter and reconstruction requirements. For this analysis, a simple filter was applied requiring 5 IceTop stations and 5 IceCube DOMs for an event to pass on to the next stage. (Note that although 8 IceCube DOMs were required to trigger IceCube, some of those hits will have been removed in a cleaning phase where hits are clustered into events. This process can leave less than 5 DOMs in the final event.) Furthermore, for some events either the IceTop or IceCube reconstruction (as discussed in Section 4.3) fails to find a successful minimum. At this stage, an event is excluded if it does not have successful reconstructions in both IceTop and IceCube. • Step 2: The shower and the track are reconstructed as passing within the area and volume of the two arrays. For an event to be considered "contained" within IceTop, the shower core must be reconstructed within the two-dimensional outer boundary of the array. To satisfy IceCube containment, the track must be reconstructed within the three-dimensional outer boundary of IceCube. This selection is applied using three different reconstruction algorithms: the final muon bundle reconstruction (which is used for the analysis), the first surface reconstruction (which has not been biased by any IceCube track reconstructions), and an independent IceCube track reconstruction (which has not been biased by taking as input a surface reconstruction as a first guess) must all satisfy IceTop and IceCube containment. As shown in Table 1, this basic selection eliminates the largest fraction of events, but is the most important one for improving core reconstruction and angular resolution. showers with small track lengths are likely to be events with a poorly reconstructed direction and/or core position. Thus, a requirement that a track length -estimated by calculating the length of the track segment between the first and last direct hit if within a -15 ns to +75 ns residual time window -is greater than 850 m for showers with S 125 >5.625, linearly decreasing to 400 m at S 125 = 0, is applied. • Step 5: The fitted slope of the LDF function is within a reasonable range. The reconstructed effective propagation length, λ 0 , should resemble the independent measurements for Antarctic ice, though a range of values can still provide a robust K 70 , as described in Section 4.2 and shown in Fig. 4. As this parameter is a free parameter in the fit, it does have a slight dependence on the mass of the primary particle; thus, to avoid introducing a bias into the event sample very loose requirements are made: •   Section 5). In this table, simulation includes proton and iron only, of which 6.6 million total events were generated. The event selection refers to the steps in Section 5. "Contained" shows the number of simulated events at each step for which the true simulated track was contained within the area and volume of the two arrays. ∆ψ is the opening angle between the true and the reconstructed direction, and ∆R is the difference between the true and the reconstructed core location. In each case, at each cut level, the number of events is given as well as the percentage they comprise of the total number of events at that cut level. After all six selection levels, at least 94% of the events used for this analysis are well-contained and have a core resolution within 25 m and an angular resolution within 1 • . Simulation events listed here are not re-weighted from their E −1 generation spectrum.
parameter for the arrival time difference of the track hypothesis at the surface of the ice: where z IT is the altitude of the reconstructed core position of the shower in IceTop, z IC is the point on the IceCube track closest to the center of gravity of the hits, t is the reconstructed arrival time of the shower at depth z and θ is the zenith angle, as reconstructed by independent IceCube (IC) and IceTop (IT) reconstructions. It has been found that these uncorrelated events can be distinguished from the truly coincident sample by requiring ∆t less than 3 µs.
With this event selection, the resolution in core position is better than 9 m and in track direction is less than 0.5 • (where resolution is defined such that 68.3% of all events are contained within that range), as shown in Fig. 6. The final event sample contains 239797 events from the August 2008 experimental data and 20289 for all five simulated primaries from the Monte Carlo simulation.

Neural Network Mapping Technique
After the event selection criteria from Section 5 have been applied, a strong relationship between muon and electron observables (the vertical axis K 70 , and the horizontal axis S 125 ), and cosmic ray primary mass and energy (A and E 0 ) can be seen in Fig. 7. Black lines approximating contours of energy (in log 10 (E 0 /GeV)) guide the eye. Mass is indicated by color; each colored bin in the K 70 -S 125 space is shaded according to the percentage of proton found in that bin for simulated proton and iron primaries only for clarity. In the analysis intermediate species were included. The separation between proton and iron is clearly visible in Fig. 7, especially at the highest energies. However, the energy contour lines are not parallel and the proton and iron distributions are also not parallel, which implies that the relation between the K 70 -S 125 space and the mass-energy space (S125/VEM) is non-linear, as the proton/iron distribution becomes distinctly more orthogonal to S 125 at higher energies.
In light of the non-linear correlation between the first two-dimensional space (K 70 -S 125 ) and the second (A-E 0 ), a mapping technique is required to step between them. For flexibility in testing combinations of multiple new input parameters under development in IceCube and IceTop, a neural network technique was developed. It is important to emphasize that both outputs are continuous in nature (as opposed to quantized); therefore, this style of neural network is referred to as either an interpolation or a regression problem [26].
A neural network consists of a set of input parameters which are connected to a set of output parameters through a series of nodes which are arranged in layers. Each node in a layer is connected to all nodes of the previous and subsequent layer via a series of weights, and is assigned an activation function which, for an interpolation problem, is typically a sigmoid function for the internal layers and a linear function for the output layer. Each node then represents the activation function acting on the summed, weighted input parameters. The structure used for this analysis can be found in Fig. 8.
The neural network used here, shown in Fig. 8, is a feed-forward multilayer-perceptron as found in the ROOT TMultilayerPerceptron package [27]. The network takes log 10 (S 125 ) and log 10 (K 70 ) as input parameters and is trained on Monte Carlo simulation data to find output parameters of log 10 (E 0 ) and ln(A) [11]. Both the two inputs and the two outputs are renormalized so that they are numbers between 0 and 1.
Over the course of a number of training cycles, the network "learns" to find the best fit to the expected outputs. During each cycle, the error between the output parameters (obtained from the given network structure) and the true parameters is calculated. The weights in the network are updated, as prescribed by a learning algorithm, to reduce this error. The error is minimized after a sufficient number of learning cycles. This iterative training procedure is performed using 1/4 of the final Monte Carlo sample (∼1000 events of each species). This "training" sample provides the network with K 70 and S 125 as well as the true mass and energy outputs.
A "testing" sample, consisting of an independent set of 1/4 of the final Monte Carlo samples, is used to check the progress of the neural network training. After each learning cycle, the testing sample is also passed through the network, the output parameters are calculated and then compared with the true values, to verify that the network is not becoming too specific to the training sample (in the case of "overtraining" [26]). The testing set is also used to optimize a number of other network parameters, such as learning method, architecture, and activation functions. For this analysis a broad range of combinations were evaluated and the network which was chosen provided the testing sample with the best energy resolution and the best mass reconstruction for the energy range accessible (1 PeV-30 PeV), within the fewest learning cycles.
The final optimized network had the following characteristics: • Architecture: One method for choosing network architecture is to incorporate the fewest hidden layers and nodes possible to provide results of the desired accuracy [26]. Following this prescription, the network chosen had two hidden layers of five nodes each (for the 2:5:5:2 setup as seen in Fig. 8).
• Activation Functions: The architecture is defined by nodes with sigmoid activation functions for the hidden layers and linear activation functions for the output layer.
• Learning Algorithm: The quasi-Newtonian learning algorithm of Broyden, Fletcher, Goldfarb and Shanno (known as BFGS) was chosen based on its efficiency in returning a desirable network [27].
• Training Cycles: Using 947 training cycles provided the smallest errors without overtraining for this particular structure, learning algorithm, and simulated data sample.
To avoid biases, the training dataset was not included in the simulation data sample used for the analysis. Furthermore, while the testing sample was not used to train the network, it was used to optimize the choice of the network structure, learning algorithm and the number of learning cycles, and it was used to optimize the minimization process (as will be discussed in Section 8), all of which could introduce a bias. Thus, the testing sample was also removed from the final sample. Half of the simulated data was independent of the training and testing samples (∼2000 events of each species) and is called the "analysis" sample. This half was used for the final analysis steps described below.
The normalized energy output from the neural network can be translated directly back into units of log 10 (E/GeV). The energy resolution of this measurement ranges from 18 − 20% in the threshold region (∼1.5 PeV) to 6 − 8% at 30 PeV, for an average resolution better than 14% over the full range of energies, with a bias which is well within the resolution, as seen in Fig. 9. In the analysis of energy spectrum and composition described in the following sections, events are placed in energy bins which are wider than this energy resolution. The normalized mass output from the neural network, by contrast, exhibits much more spread in mass-space and is not well-resolved, as seen in Fig. 11. The most likely average composition of the population of cosmic rays can still be measured, with a further minimization step to be described in Section 8.

Energy Spectrum
At this point it is important to note that the correction described in Section 4.4 accounts for the signal attenuation due to the snow, but cannot account for the change in the composition-dependent energy threshold -which is higher in energy for the older half of the array due to the deep snow. Therefore, from this point forward, in both simulated and experimental data events with a shower core located on the older side of the array (x-position more than 200 m in Fig. 2(a)) have been removed to minimize this effect. The number of events surviving as a function of reconstructed energy can be found in Table 4.
For a determination of the energy spectrum, the data were divided into energy bins. Nine bins for each decade in energy are chosen so that the bin width is sufficiently larger than the resolution and bias in the reconstructed energy.
The flux, Ψ, as a function of energy, E 0 , is given by: where e is Euler's constant (the base of the natural logarithm); η(E 0 ) is the efficiency (Fig. 10(a)), defined as the ratio of simulated events left after the final event selection to the number generated; A = π(1200 m) 2 is the geometrical area over which the CORSIKA air showers were generated; τ = 29.78 days is the exposure time, or livetime, of the detector; and Ω is the solid angle weighted by the generated zenith distribution: where θ is the zenith angle. In other words, the flux can be found by dividing the energy distribution by the detector aperture, A(E 0 ) = η(E 0 )AΩ, and the exposure time, τ . The efficiency can be found in Fig. 10(a) for proton and iron separately. To account for the fluctuations caused by the comparatively small statistics found in the simulated aperture distribution the aperture, A, is fitted to a sigmoid function modified by a line, of the form: as can be found in Fig. 10(b), which shows the aperture for all generated particle types. The slight loss of efficiency at the highest energies is due to events for which the likelihood minimization algorithm does not converge. The resulting measured spectrum is shown in Fig. 10(c), where the minimum energy shown corresponds to the point of maximum efficiency from the fit Eq. (14), and the maximum shown is chosen to be two energy bins below the maximum energy simulated to avoid edge effects. A phenomenological evaluation was developed [28,29], where the differential cosmic ray flux was expressed as: where Z is the charge of a particle, Φ 0 Z is the absolute flux at 1 TeV, γ Z and γ c are the indices below and above the cut-off energyÊ Z , respectively, and c is a parameter defining the curvature of the knee (where the knee becomes sharper as c increases). The measured spectrum in Fig. 10(c) has been fit to the flux model in Eq. (15), with the resulting power law indices of 2.62 ± 0.05 at energies below and 3.26 ± 0.06 above a slowly turning knee around 4.98 ± 0.37 PeV, where the errors given are statistical. The data points are given in Table 2, and the fit parameters are shown in Table 3. The systematic error bars will be discussed in Section 9.

Composition
As mentioned in Section 6, mass composition requires a further minimization step. Within each slice in energy and for each simulated species, the neural network produces a distribution of mass outputs which are referred to as "template histograms". Examples of template histograms for all five species are shown in Fig. 11(a)(right), for a particular slice in energy selected from Fig. 11(a)(left). Experimental data, when passed through the same neural network, also has a histogram of mass outputs which can be decomposed into a linear combination of the template histograms of the individual nuclear species (proton, oxygen, etc.). A minimizer finds the optimal mixture of simulated species to match the experimental data.
Note that the distributions of the templates in Fig. 11(a)(right) very rarely approach either 0 or 1 on the output axis: this is due to the overlap in the distributions, as shown in violet in Fig. 7. The type output from the network can be seen as a probability of being a certain species of particle: as the separation between species is never complete the probability that a particle is any given type  (c) Figure 10: Upper left: The efficiency (as defined in the text) with respect to the true energy, where (red) circles mark proton and (blue) squares denote iron. The efficiency is much smaller than 1: this is a combined effect of the requirement of coincidence between the two detectors, and the large radius within which the showers are generated.
Upper right: The aperture with respect to the energy reconstructed by the neural network for all simulated data (all particle types), fitted as described in Eq. (14). Lower: The reconstructed energy spectrum for the August 2008 data, which uses the aperture shown above. The additional symbols show alternate results from simulations using alternate models, as discussed in Section 9.    Table 4: Number of events of experimental data and simulation, for the six energy ranges analyzed for composition, after application of all cuts (including the selection for the newer half of the array).
will never be 1. Furthermore, the higher mass particles have much more overlap between types, thereby causing probability of reconstructing an iron particle to be slightly lower.
To find the best mixture of species, the contribution of each species is quantified as a fraction, with the constraints that each fraction must be between 0 and 1, and the fractions of all nuclei must sum to unity. The template histograms C for each species are weighted according to these fractions and summed to produce a histogram of their mixture C mix . Two species "A" and "B" (i.e. proton and iron only, as shown in Fig. 11(b)), are mixed using: where f A is the fraction of the mixture which is species A. Expanding to three species "A", "B", and "C", the mixture histogram is formed by: Here, f A is the fraction which is species A, f B is the fraction of what remains (which is not A) which is of species B. Mixtures of four or more species can be constructed with similar logic. A minimizer can vary f A , f B , etc. as free parameters. Constructing the mixture in this way makes it easy to enforce in the minimizer that all fractions are physical (that is, between 0 and 1, and properly normalized). The mixture of species that best matches the data is found using a χ 2 test and the software package MINUIT [30], where: Here, C data is the measured data histogram, as reconstructed by the neural network output, and σ 2 data is the variance in the data. C mix is a mixture of simulated species as described above. An example of a best-fit mixture of template histograms for a slice in energy using only two species is Energy Range Reconstructed (log 10 (E reco /GeV)) ln A ± stat. ± sys.  shown in Fig. 11(b), where the "testing" sample from Section 6 (including all species) was treated as though it were data, and compared to template histograms of proton and iron only (for easier viewing). Simulated nuclei from the testing sample (described in Section 6) were "hand-mixed" in different proportions at different energies to test the effectiveness of various minimization strategies, ranging from using only the two most extreme template histograms (protons and iron) with only one free parameter, to using all five template histograms with four free parameters. With too few nuclei, the minimizer is confident about the best mixture, but the mixture is not a very good fit. With too many, the minimizer can find many options for achieving a good fit (many of them unphysical mixtures, such as iron and oxygen present with no silicon in between), so it is less confident about the best one. Using three template histograms: protons, "intermediate nuclei" (a 50-50 mixture of helium and oxygen), and iron, is a balanced approach: A list of the number of events for each species as a function of energy is given in Table 4. The mass output distributions with respect to energy are shown for the final simulated "analysis" sample in Fig. 12(a). The results of a check of this method using a mixed composition sub-set of the "testing" sample of simulated data are shown in Fig. 12(b), where the true mass value (gold bars) is compared with the reconstructed mass value (red stars). This minimization technique reconstructs the true mean logarithmic mass within the error bars for each slice in energy. Furthermore, Fig. 12(c) and Fig. 12(d) show this technique applied to pure proton and iron subsets of the "testing" sample of simulated data, respectively, to demonstrate the capability of the method to reconstruct both extremes. As discussed above, as the probability of finding an iron is never 1, the iron reconstruction never quite reaches the expected value but remains within statistical error bars.
The composition results of this combined neural network and minimization technique applied to the August 2008 data are shown in Fig. 13, with data points and statistical error bars in red. The systematic errors (shown as alternate symbols) are discussed below and are summarized in Table 5.  Figure 11: (a), left: the mass output from the neural network is shown as RMS profiles in slices of the reconstructed energy with colors corresponding to the different species, as labeled. This plot is analogous to Fig. 7 after the neural network transformation to the new axes; therefore, the large overlap at the lowest energies corresponds to the similarly large overlap in Fig. 7. (a), right: The underlying distributions, or template histograms, for a single slice in energy for all 5 particle species. (b): An example of template histogram fitting: the testing simulation, or "test data", sample (black with statistical error bars), which includes all five species, was fit to the proton and iron template histograms only -for easier viewing -using the minimization technique described in the text. The grey "mixing ratio", or combination of proton and iron comprising the test data, matches the data quite well even though only two template histograms were used to fit a sample containing five species.    Figure 13: The mean logarithmic mass vs primary energy for August 2008 data as calculated using the analysis technique described in this work. The error bars on the data points are statistical, while the additional symbols indicate results after systematic shifts estimated using alternate models.

Systematic Uncertainties
This measurement of cosmic ray composition and energy spectrum depends crucially on the two independent observables K 70 and S 125 . Systematic errors, particularly in the scaling of one observable (independent of the other) can make the composition seem heavier or lighter. To estimate the magnitude of potential systematic errors and their effect on the measurement of composition, the K 70 -S 125 parameter space (as depicted in Fig. 7) was studied using samples of simulated events with alternate models, listed below.
• Hadronic Interaction Model: As the forward-direction physics of the first interaction of highenergy cosmic rays has not yet been measured by accelerator experiments on Earth, cosmic ray physicists depend upon extrapolations from the known physics at lower energies. There are, therefore, a number of different high-energy hadronic interaction models. Model used: SIBYLL-2.1 [13]. Alternate models: EPOS-1.99 [31] and QGSJET-II-03 [32]. • Optical Properties of the Ice: For the 40-string configuration of IceCube there were two models available for measuring the depth and wavelength-dependence of the scattering and absorption in the ice. In general, the absorption lengths vary from 100 m to 200 m and effective scattering lengths range from 20 m to 70 m due to the concentration of dust in the glacial ice, which has a strong depth-dependence [19,33,34]. Model used: AHA [19], the standard ice model for the 40-string configuration of IceCube and was developed using data from the predecessor of IceCube (AMANDA). Alternate model: SPICE-1, which measured the ice properties over the full depth range of IceCube using the in-situ LEDs present on every DOM main board [33]. • DOM Efficiency: The efficiency of the DOMs in the ice depends upon a number of factors including the PMT quantum efficiency, the transmission of the glass and the gel, and the refrozen ice from the deployment hole (which is less clear than the bulk ice). Measurements of the DOM efficiency in a controlled setting have revealed an estimated uncertainty of approximately 8% [35], which is equivalent to a shift in log 10 (K 70 ) of 0.035 in log 10 (K 70 ). This was treated as a constant systematic uncertainty across all values of S 125 which could either raise or lower the measured K 70 in IceCube. IceTop is not affected because its signals are calibrated to Vertical Equivalent Muons, as described in Section 2. • Snow correction in reconstruction: As discussed in Section 4.4, snow accumulation at the surface can affect signals in IceTop, and reconstruction methods can correct for this effect in data. However, the correction applied during reconstruction will not work perfectly for all events at all energies. Potential systematic errors in S 125 due to imprecision in the technique itself were studied by comparing simulations in which snow was "removed" during reconstruction to simulations in which there was no snow at all. • Atmosphere: effect on high-energy muons: Atmospheric conditions high in the atmosphere, and the effective temperature of the atmosphere, affect the number of high-energy muons and K 70 . With the selection of one month of data for this analysis, the effective temperature was stable; K70 exhibits negligible daily or weekly variation within that month. In addition, the atmospheric profile high in the atmosphere was compared to that simulated by CORSIKA, and the systematic error due to this mismatch is also negligible compared to other systematic effects on K 70 discussed in this section. The actual average atmospheric conditions during August 2008 differs from that simulated in CORSIKA. The systematic effect on S 125 due to this mismatch was estimated by comparing CORSIKA showers (10 PeV and 10 degree zenith angle) made using the standard South Pole winter atmosphere profile to those made using a customized profile built from measurements taken during August 2008. Distributions of log 10 (S125) differ by less than 0.01, and this potential shift was treated as constant across all energies. Understanding these effects in IceTop and IceCube is the subject of ongoing study.
Some of these effects (such as ice model and DOM efficiency) affect only the in-ice measurement of K 70 , while others (such as snow removal and atmospheric pressure) affect only the surface measurement of S 125 . The hadronic interaction model potentially affects both. A summary of how K 70 and S 125 are observed to shift when alternate models are explored is shown in Fig. 14 on the left and right, respectively. Potential shifts in log 10 (K 70 ) are shown as a function of log 10 (S 125 ), while potential shifts in log 10 (S 125 ) are shown as a function of log 10 (K 70 ).
To explore how the final results (composition and spectrum) would change under an alternate model, values of either log 10 (K 70 ) or log 10 (S 125 ) in the data were artificially shifted according to the estimated curves in Fig. 14. (For alternate hadronic interaction models, events were shifted half in K 70 and half in S 125 .) For each alternate model, these shifted data sets were then run through the full analysis machinery (neural network reconstruction of A and E 0 and minimization with template histograms) to yield "alternate" composition and energy spectrum results. Although primitive, this procedure allowed for exploration of many alternate models quickly, without having to generate large sets of simulated data for each one. Figure 13 shows the composition result (with standard simulations) together with these alternate results, after applying a shift either to K 70 or to S 125 (or both). These results, including errors, are also summarized numerically in Table 5. The "systematic errors" are drawn from the maximum and minimum ln A yielded by any of the alternate models at each energy. Figure 10(c) shows a similar collection of alternate results for the energy spectrum. The numeric values are given in Table 2.  Figure 14: Left: Estimates of systematic uncertainties affecting the in-ice measurement (as shifts in log 10 K 70 for a given log 10 S 125 ). Right: Estimates of systematic uncertainties affecting the surface measurement (as shifts in log 10 S 125 for a given log 10 K 70 ).

Results and Discussion
Data from the combined IceTop and IceCube detectors in their 2008 40-string/40-station configuration have been used to develop a technique to measure the cosmic ray energy spectrum and composition at energies between 1 PeV and 30 PeV. Using measurements of the electromagnetic component of the air shower at the surface (S 125 ) and the muonic component of the air shower in the ice (K 70 ), a neural network in conjunction with a χ 2 minimization algorithm have been employed to extract the relevant physical parameters: ln A and E 0 .
This technique provided an energy resolution better than 14% with small offset due to primary mass, which led to an energy spectrum that compares favorably with other recent measurements shown in Fig. 15. In these figures, the shaded (red) region indicate the systematic errors for this analysis. Comparisons are made to other experiments (with statistical error bars only) in Fig. 15(a). Figure 15(b) highlights a comparison (including systematic errors) to two spectra measured by the 26-station configuration of IceTop only (IT-26) [21], one assuming protons and one assuming a two-component model. This work (using a coincidence detector and a neural network based on both surface and deep observables) and the IT-26 work (using the IceTop surface detector only and an unfolding technique) yield slightly different spectra, but are consistent within systematic errors. When fit to the flux model in Eq. (15), the spectrum presented here indicates a power law of indices 2.61 ± 0.07 below and 3.23 ± 0.09 above a slowly turning knee around 4.75 ± 0.59 PeV, where the errors given are statistical.
This method also provided a reconstructed mass parameter and, in turn, a measurement of cosmic ray composition: the mean logarithmic mass, ln A . The mass composition is compared with other experimental measurements in Fig. 16 where, additionally, systematic error bars are E / GeV  Figure 15: The all-particle cosmic ray spectrum found using data from 40-station IceTop/IceCube as described in the text. For this analysis the results are shown as (red) stars, with statistical error bars as solid (red) lines and systematic errors as the shaded (red) region. Upper: this result compared to other experiments with their statistical errors only. Lower: this result compared to two spectra from IceTop-26, with systematic errors also shown. Fluxes are multiplied by E 2.7 , and the scales of the two plots are the same. The overall flux and shape found in this work are very similar to those reported by other experiments. Data points from [21].   15), but dissimilar to others (in particular, optical measurements). Data points compiled from [28,36].
shown for IceCube/IceTop. In this analysis, using the 40-string/40-station configuration of IceTop and IceCube from 2008, a mass is observed which similar to some other measured results (especially those using a similar electons vs. muons stretegy as in this work), and in disagreement with others (especially the optical measurements, which emply a very different strategy). The slope of the strong increase in mass through the knee region is nearly identical to the model in Eq. (15). Additionally, this new analysis technique with IceCube-40/IceTop-40 shows results consistent within systematic errors to those from a different technique applied to data from its predecessor, SPASE-2/AMANDA-B10 [24].
In the future this technique can be expanded to include new composition-sensitive input parameters, as well as employ the larger IceCube detector and a greater quantity of data. From looking at several simulated alternate models, it is clear that systematic uncertainties (both in the in-ice measurement and surface measurement) can greatly affect the measured composition and spectrum, regardless of the detector's size or livetime. Although the mean logarithmic mass itself is difficult to measure absolutely because of systematics, this work shows an unmistakable trend of increasing mass, regardless of the differing absolute scale of ln A measured with respect to different models. The systematic studies performed for this analysis have led to many improvements which will allow for more precise measurements with new data. Therefore, in the coming years the complete IceCube Neutrino Observatory will be the only detector of its kind able to provide both composition and energy spectrum measurements from energies overlapping with direct measurements below the knee to energies nearing the ankle.