Cosmic ray composition study using machine learning at the IceCube Neutrino Observatory

Matthias Plum (for the IceCube Collaboration)

REVIEW 3 major objections 3 minor 15 references

Cosmic ray composition study using machine learning at the IceCube Neutrino Observatory

T0 review · 3 major / 3 minor · reviewed 2026-08-14 · deepseek-v4-flash

Challenge this review Re-run · record.json Download PDF Read on arXiv ↗ Review JSON

Pith's one-line read Adding shower age $\beta$ improves cosmic-ray composition resolution.

desk verdict A modest, internally consistent simulation study showing that adding IceTop's beta parameter slightly improves IceCube's ML-based mass-template resolution, with the caveat that the evidence is a closure test under a single hadronic model. read the letter →

arxiv 1908.06433 v1 pith:RHLHB67L submitted 2019-08-18 astro-ph.HE

Matthias Plum (for the IceCube Collaboration) This is my paper · ORCID

classification astro-ph.HE

keywords cosmic-raycompositionkneeregionmachinelearningrandomforestIceCubeNeutrinoObservatoryTopshowerageparameterbetatemplatefit

verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

This paper argues that cosmic-ray mass composition in the knee region around 3 PeV can be measured more precisely by feeding the IceTop shower-age parameter $\beta$, a measure of the shower's development stage from the Double Logarithmic Parabola fit, into a machine-learning mass regressor alongside the existing IceCube/IceTop observables. Using full Monte Carlo air-shower simulations, the author trains a random-forest regressor for protons, helium, oxygen, and iron, converts its output into kernel-density templates per energy bin, and fits weighted template mixtures to mock data. The claimed result is that adding $\beta$ improves the mass resolution over the whole studied energy range relative to the baseline analysis, while preserving accurate reconstruction of the input composition in both a maximum-mixing scenario and a realistic H4a scenario. A sympathetic reader would care because the knee is where the cosmic-ray source population is thought to transition from galactic to extragalactic, and better composition resolution directly tightens the statistical mapping of that transition. The demonstration is entirely simulation-based, so the claim is about what the method can do, not yet what real data show.

What carries the argument

The load-bearing machinery is a random-forest regressor whose inputs are the IceCube/IceTop coincidence observables; the new ingredient is the IceTop shower-age parameter $\beta$, derived from the Double Logarithmic Parabola fit to the lateral charge distribution, which is sensitive to how early or late in its development a shower is observed and therefore to the primary mass. The regressor's output for each simulated primary is converted into kernel-density-estimated probability densities that serve as templates, and an extended likelihood fit over the four-element weighted mixture $P_{\mathrm{mass}}(X) = \sum_i w_i P_i(X)$ extracts the fraction of each element in a data set. The $\beta$ parameter is what distinguishes the improved analysis from the baseline; the templates, fits, and mock-data studies are the machinery that demonstrates the improvement.

What would settle it

Retrain the entire pipeline with an alternative high-energy hadronic model, such as EPOS-LHC, and compare the $\beta$-based resolution improvement; if the gain shrinks or reverses, the claimed improvement is tied to SIBYLL-2.1 rather than to actual air showers. Alternatively, apply the same templates to real IceCube/IceTop coincidence data and check whether the reconstructed fractions agree with the published composition from the same detector.

Watch

Extended reading notes

Core claim

The paper's central claim is that the shower-age parameter $\beta$, obtained from the Double Logarithmic Parabola fit to the IceTop lateral charge distribution, carries composition information that the baseline observables—the logarithmic IceTop signal at 125 m, the reconstructed zenith angle, and the muon-bundle energy deposit at 1500 m slant depth—do not fully capture. Adding $\beta$ to a random-forest regressor produces per-element mass-output templates that are visibly more separated for hydrogen, helium, oxygen, and iron, and the resulting template fit reconstructs input fractions within statistical uncertainties in both the 25/25/25/25 maximum-mixing scenario and the H4a scenario. Across the energy range $\log_{10}(E/\mathrm{GeV})$ from roughly 6.8 to 8.0, the improved analysis shows a slight but consistent improvement in mass resolution relative to the baseline, with the largest gains in the intermediate helium and oxygen groups, which are hardest to separate because their template distributions overlap.

Load-bearing premise

The load-bearing premise is that the full Monte Carlo chain—CORSIKA with FLUKA and SIBYLL-2.1 plus the IceCube/IceTop detector simulation including snow-depth corrections—reproduces the real air-shower observables and their composition dependence, because every template, mock data set, and resolution number is generated from that simulation.

Editorial extensions

If this is right

Applied to real IceCube/IceTop coincidence data, the improved templates should return mass fractions in the knee region with smaller statistical uncertainties than the baseline analysis used in [6,7].
The template-fitting chain can be rerun for any proposed composition scenario; the H4a test shows that realistic non-equal fractions are recoverable within statistical errors.
Because the improvement appears across the whole energy range, the method sharpens the measurement of where the galactic-to-extragalactic transition shows up in the composition.
The resolution gain is largest for helium and oxygen, the intermediate groups whose template distributions overlap most, so future work that separates those two elements will yield the largest improvement.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

A testable extension the paper does not report: retrain with an alternative high-energy hadronic model such as EPOS-LHC; if the beta gain persists, the improvement is physical rather than a SIBYLL-2.1 artifact.
The same regressor-plus-template structure can be pointed at other shower-classification tasks, such as separating muon-rich from electromagnetic showers, since beta already encodes longitudinal development.
Before real-data application, one should check how strongly beta is correlated with the IceTop signal strength S125; if the two are nearly redundant, the improvement may shrink once detector noise is included.
Since the demonstration is entirely in simulation, the immediate next step is to apply the templates to actual IceCube/IceTop data and compare the reconstructed fractions with the published composition analysis.

Share X Bluesky LinkedIn Reddit HN

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Desk Editor's Note

read the letter

Short version: this is an honest, incremental method paper from the IceCube collaboration. The genuinely new element is the shower-age parameter beta from IceTop's double-log-parabola fit, added to the random-forest mass templates the group already uses in [6,7]. The paper reports a small resolution improvement over the whole knee energy range from that extra variable. The simulation study is internally consistent: in the maximum-mixing and H4a scenarios, the template fit recovers the injected fractions, and Figure 9 does show the beta-augmented analysis with narrower fitted distributions. The text is appropriately modest, everything is labeled preliminary, and the analysis is framed as under development. As far as I can tell, the citation pattern is normal for an IceCube proceeding; the listed prior work is the relevant baseline, not padding.

The soft spot is the one you would expect from a short proceedings: all the quantitative support comes from a closure test. The fast Monte Carlo data in Section 4 are sampled from the same KDE templates the fit then uses, so the recovered fractions being close to truth is mostly the fitter inverting its own input. That validates the machinery but does not test whether the templates match real events. The templates themselves come from one simulation chain, CORSIKA with FLUKA and SIBYLL-2.1, so the claimed beta sensitivity in Figure 4 and the improvement in Figure 9 depend on that one hadronic model being right in exactly the new observable. There is no real-data comparison, no alternative hadronic model, and no error bars on the resolution curves. Those are genuine limitations, but they are limitations of scope rather than signs of carelessness; the paper does not oversell.

If I were working on composition template methods, I would read this as a useful data point: the beta variable is worth trying, and the result suggests a modest gain if it survives real data. But I would not cite the resolution improvement as established, and I would not treat this as a measurement. The right audience is people tracking IceCube's composition program or ML-based air-shower analyses.

For peer review: I would not desk reject it. If it came to a journal as a full paper, I would send it to a referee and expect the referee to ask for either real-data validation or a second hadronic model before the improvement claim is accepted. For an ICRC proceedings, it is a fair, clearly-scoped contribution.

Referee Report

3 major / 3 minor

Summary. The paper describes a machine-learning-based template fit for cosmic-ray mass composition using IceCube/IceTop coincidence events. A random forest regressor is trained on full Monte Carlo simulations with the baseline observables S125, cos(theta), and log10(dE/dX1500m), and an 'improved analysis' additionally includes the IceTop shower age parameter beta. KDE templates per primary element (H, He, O, Fe) and per energy bin are combined into a weighted response PDF, and an extended likelihood fit extracts composition fractions. The method is tested on fast Monte Carlo pseudo-experiments for a maximum-mixing scenario (25%:25%:25%:25%) and an H4a scenario. The paper reports that the improved analysis reconstructs the input composition accurately and slightly improves the mass resolution over the whole energy range compared with the baseline.

Significance. The paper is a well-scoped method contribution: it applies a standard random-forest/KDE template technique and tests a new observable (shower age beta) in a controlled simulation environment. The comparison between baseline and improved analyses is clearly presented, and the use of full Monte Carlo training plus fast pseudo-experiments is methodologically clean as a closure test. If the claimed resolution improvement survived a cross-check with an independent hadronic model or with real data, it would be a useful step toward improving IceCube/IceTop composition measurements in the knee region. At present, however, the results are strictly simulation-only and self-consistent by construction, so the quantitative significance is limited. The paper does not provide data validation, an alternative hadronic interaction model, or uncertainties on the quoted resolutions.

major comments (3)

[Section 4 (Figs. 6-10) and Section 5] The validation is a closure test: the fast mock data sets are generated by sampling from the same KDE template PDFs (the weighted sum in Section 3, Pmass(X) = sum_i w_i P_i(X)) that the extended likelihood fit subsequently uses. Recovering the input fractions and observing an improved resolution in Fig. 9 therefore verifies only that the fitter can invert its own generative model; it does not test whether the templates, including the composition dependence of the shower age parameter beta from Fig. 4, match real IceCube/IceTop events. Because the entire analysis chain rests on CORSIKA with FLUKA and SIBYLL-2.1 (Section 2), the central claim in Section 5 that the improved analysis 'shows a slight improvement over the whole energy range' is conditional on that hadronic model being accurate in exactly the observable added. Please either add a cross-check with an independent hadronic interaction model or a data-based validation, or explicitly reframe the result as a simulation-only closure test.
[Section 4, Figs. 8 and 9] The resolution curves in Fig. 9 and the reconstructed fractions in Fig. 8 are shown without statistical uncertainties. The resolution is obtained from Gaussian fits to distributions such as Fig. 7, so each resolution value carries a sampling uncertainty determined by the number of pseudo-experiments and the fitted event counts; without error bars or a significance test, the reported 'slight improvement' cannot be distinguished from statistical fluctuation, especially for helium and oxygen, which the text notes have large overlap with neighboring distributions. Similarly, the statement that both analyses reconstruct the true composition 'inside the statistical uncertainties' (Section 4) is not quantitatively supported in the figures as presented. Please add uncertainties to Figs. 8 and 9, or provide the underlying numerical values.
[Section 4, H4a scenario (Fig. 10)] The H4a test inherits the same closure-test limitation as the maximum-mixing test, and it is presented only as average reconstructed fractions in Fig. 10, with no uncertainties, no resolution comparison, and no goodness-of-fit measure. The claim that 'both analyses are capable of reconstructing the primary mass composition in a realistic source scenario' is therefore not quantitatively supported as stated. Please add per-bin uncertainties or pull distributions for the H4a reconstruction, or soften the claim accordingly.

minor comments (3)

[Section 3] The abbreviation 'RFT' is defined as 'random forest tree'; since a random forest is already an ensemble of trees, consider using 'random forest regressor' or 'RF' for clarity.
[Section 2] The simulation description would benefit from stating the CORSIKA version and the versions of FLUKA and SIBYLL-2.1 used, because composition-sensitive results are known to depend on hadronic interaction model versions.
[References] Reference [7] is cited as 'arXiv:1906.04317' without a title or author list; if this is an IceCube conference paper, please provide the full citation, such as the corresponding PoS(ICRC2019) number.

Circularity Check

1 steps flagged · score 6.0 of 10

Composition reconstruction claim reduces to a closure test because the fast mock data are generated from the same KDE templates that the fit uses.

self definitional [Section 4, 'Reconstruction Capabilities of Different Composition Scenarios'; see Eq. (1) in Section 3 and Figure 6 caption.]
"The mass composition response PDF for each energy bin can be used to generate fast mock scenario data sets where the number of events and the fractions of each primary group are artificially selected. In this way the reconstruction capabilities of the method are tested for several Monte Carlo scenarios."

The mock data are drawn from the same Pmass(X) = sum_i w_i P_i(X) model (Eq. 1) that is then used as the fit template. This makes the likelihood correctly specified by construction, so the fitted fractions returning the input fractions within statistical error is a property of the estimator, not an independent test of composition sensitivity. It cannot detect a mismatch between the KDE templates and real air showers; the only external input is the CORSIKA/FLUKA/SIBYLL-2.1 simulation chain of Section 2. The paper's accompanying statement that 'the input PDF is reconstructed within the statistical uncertainties' (Section 4, Figure 6) is therefore a closure test rather than a prediction.

full rationale

The paper's central demonstration that the template fit 'reconstructs the true composition' rests on fast Monte Carlo data sets generated from the same KDE response PDF that the extended-likelihood fit uses as its template. Because the generating model and the fitting model are identical by construction, recovering the input fractions in Figures 8 and 10 verifies internal consistency of the fitter, not the power of the observables to separate masses in real data. This is the main circular step. The beta-based improvement in Figure 9 is less formally circular, since beta is an additional simulated observable and both analyses are treated symmetrically; however, the resolution numbers are also produced by generating and fitting fast Monte Carlo from the same templates, so they quantify only how well the simulated templates can be inverted under one hadronic model. No real-data comparison or alternative hadronic model check is presented, so the beta improvement remains an in-sample result. Because the composition-reconstruction claim reduces by construction to the fitting model, a partial-circularity score of 6 is appropriate; the beta comparison retains some independent content within the simulation and therefore does not warrant a higher score.

Assumptions & free parameters 3 free parameters · 4 assumptions · 0 invented entities

The central result rests on a single simulation chain and on template PDFs that are reused to generate and fit mock data, which is the main source of circularity. No free parameter is reported numerically, and no systematic variations such as hadronic model, snow depth, or selection effects are shown.

free parameters (3)

Template composition fractions w_i = Not reported; constrained so that the sum of w_i equals 1
The mass-response PDF Pmass(X) = sum_i w_i P_i(X) is fit to data or mock data, and the fitted fractions are the extracted composition. The fractions are not independent of the template shapes and are highly correlated with each other, as noted in Section 3.
KDE bandwidth = Not stated
Kernel density templates such as those shown in Figure 5 require a bandwidth choice that affects template smoothness and resolution; the paper does not specify the bandwidth or any optimization procedure in Section 3.
Random forest hyperparameters = Not stated
The random forest regressor is mentioned in Section 3, but tree count, depth, feature set, and other hyperparameters are not given. These choices affect the mass-output separation and hence the measured resolution.

assumptions (4)

domain assumption CORSIKA with FLUKA and SIBYLL-2.1 accurately simulates air showers and detector response, including snow coverage.
All training, templates, and mock data come from this simulation chain; no data or alternative hadronic models are used for validation, so the whole analysis inherits the simulation's accuracy. The chain is introduced in Section 2 and used throughout Sections 3 and 4.
domain assumption The four primaries H, He, O, and Fe, equally spaced in ln A, form a sufficient training basis for arbitrary compositions.
The template PDFs are built only for these four elements, and Pmass linearly combines them via Equation (1). Any real composition is assumed to be representable as a weighted sum of these four templates, which is an assumption about the smoothness of composition sensitivity in the observable space.
domain assumption The selected high-quality full Monte Carlo events and the vertical zenith cut of theta <= 30 degrees are representative of the final data sample.
Only vertical coincident events are used, and only roughly 2000 events per energy bin are selected for training and testing. Selection efficiency and bias are not quantified, and the paper does not show how the cut affects the composition resolution. This is stated in Sections 2 and 3.
standard math KDE template PDFs and the extended-likelihood fit yield unbiased estimates in the statistical-only limit.
The method's statistical consistency is assumed; the paper tests only one realization with fast Monte Carlo generated from the same templates, so fitting biases, template binning effects, and KDE smoothing biases are not separately characterized. This appears in Sections 3 and 4.

how reviews work

0 comments

Cite this review

Pith. "Pith review of Cosmic ray composition study using machine learning at the IceCube Neutrino Observatory." pith.science (2026). https://pith.science/paper/RHLHB67L

@misc{pith2026190806433,
  author       = {Pith},
  title        = {Pith review of: Cosmic ray composition study using machine learning at the IceCube Neutrino Observatory},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/RHLHB67L}},
  note         = {Machine review of arXiv:1908.06433}
}

abstract

The evaluation of mass composition of cosmic rays in the knee region ($\sim 3$ PeV) is critical to understanding the transition in the origin of cosmic rays from galactic to extragalactic sources. The IceCube Neutrino Observatory at the South Pole is a multi-component detector consisting of the surface IceTop array and the deep in-ice IceCube detector. By applying modern machine-learning techniques to cosmic-ray air showers reconstructed coincidentally in both detector components of IceCube observatory, the energy and the mass of primary cosmic rays in this transition region can be measured. In this contribution, we will discuss the reconstruction performance and composition sensitivity of IceCube observables presently under development.

Figures

Figures reproduced from arXiv: 1908.06433 by the authors.

**Figure 2.** Example coincidence event in IceCube and IceTop. The color represents the timing information from early (red) to late (blue). The size of the circles corresponds to the size of the measured signal in the PMTs. The red line is the reconstructed shower axis. In total, IceCube covers an effective detector volume of 1 km3 . On top of most IceCube strings, two IceTop [2] tanks are placed, separated by roughly 10 m. Due t… view at source ↗

**Figure 3.** Mean and standard deviation of reconstructed log10(dE/dX1500m)in IceCube as a function of energy for different primaries. 6.4 6.6 6.8 7.0 7.2 7.4 7.6 7.8 8.0 log10(ETrue/GeV) 2.8 3.0 3.2 3.4 3.6 rec H He O Fe ICECUBE PRELIMINARY [PITH_FULL_IMAGE:figures/full_fig_p003_3.png] view at source ↗

**Figure 5.** Example bin with the KDE [9] mass template PDFs generated with a RFT regressor of the baseline (solid) and improved (dashed) analysis in the energy bin log10(E/GeV)=6.8-6.9. tively different from the templates of the baseline (solid) analysis. These template PDF’s for each energy bin of the four elementary groups are combined to a joined mass composition response PDF by introducing a weight factor for each primary g… view at source ↗

Figures from the paper (5 more)

**Figure 6.** Figure 6: Example fast Monte Carlo data set from the improved analysis from the maximum mixing scenario of the four elementary groups (Fraction: 25%:25%:25%:25%) in the energy bin log10(E/GeV)=6.8- 6.9 with in total 3000 MC events. The generator PDF is shown as a solid black li…

**Figure 7.** Figure 7: Example bin of the reconstructed number of events from the Monte Carlo reconstruction study for the maximum mixing scenario of the four elementary groups (Fractions: 25%:25%:25%:25%) in the energy bin log10(E/GeV)=6.8-6.9 with in total 3000 MC events per bin. The basel…

**Figure 8.** Figure 8: Comparison of average reconstructed mass fractions by the baseline and improved analysis with each bin containing 3000 events. Both analyses reconstruct on average the Monte Carlo truth, represented by the solid line. Note the zoom on the y-axis around 25% to emphasize…

**Figure 9.** Figure 9: Comparison of mass composition resolution derived from a Monte Carlo study with 3000 MC events per bin. The improved analysis shows a slight improvement over the baseline analysis. 0.00 0.25 0.50 0.75 1.00 H Baseline Analysis H Improved Analysis 0.00 0.25 0.50 0.75 1.0…

**Figure 10.** Figure 10: Comparison of mass reconstruction based on the H4a [11] model of the baseline and the improved analysis with 3000 MC events per bin. 7 [PITH_FULL_IMAGE:figures/full_fig_p007_10.png]

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

15 extracted references · 5 canonical work pages

[1]

IceCube Collaboration, M. G. Aartsen et al., JINST 12 (2017) P03012

2017
[2]

Abbasi et al., Nucl

IceCube Collaboration, R. Abbasi et al., Nucl. Instr . and Meth. A 700 (2013) 188–220

2013
[3]

Heck et al., CORSIKA: A Monte Carlo code to simulate extensive air showers, Report FZKA 6019, F orschungszentrum Karlsruhe, 1998

D. Heck et al., CORSIKA: A Monte Carlo code to simulate extensive air showers, Report FZKA 6019, F orschungszentrum Karlsruhe, 1998

1998
[4]

Battistoni et al., AIP Conference Proceedings 896 (2007) 31–49

G. Battistoni et al., AIP Conference Proceedings 896 (2007) 31–49

2007
[5]

E. Ahn, R. Engel, T. Gaisser, P. Lipari, and T. Stanev, Physical Review D 80 (2009) 94003

2009
[6]

Andeen and M

IceCube Collaboration, K. Andeen and M. Plum, PoS(ICRC2019)172 (these proceedings)

work page
[7]

IceCube Collaboration, arXiv:1906.04317

arXiv 1906
[8]

Pedregosa et al., Journal of Machine Learning Research 12 (2011) 2825–2830

F. Pedregosa et al., Journal of Machine Learning Research 12 (2011) 2825–2830

2011

Show all 15 references

[9]

Cranmer, Computer Physics Communications 136 (2001) 198 – 207

K. Cranmer, Computer Physics Communications 136 (2001) 198 – 207

2001
[10]

Verkerke and D

W. Verkerke and D. Kirkby, arXiv e-prints (2003) physics/0306116

2003 arXiv
[11]

T. K. Gaisser, Astroparticle Physics 35 (2012) 801

2012
[12]

Xinhua and E

IceCube Collaboration, B. Xinhua and E. Dvorak, PoS(ICRC2019)244 (these proceedings)
[13]

Kauer, PoS(ICRC2019)309 (these proceedings)

IceCube Collaboration, M. Kauer, PoS(ICRC2019)309 (these proceedings)
[14]

Leszczy´nska and M

IceCube Collaboration, A. Leszczy´nska and M. Plum, PoS(ICRC2019)332 (these proceedings)
[15]

Schaufel and K

IceCube Collaboration, M. Schaufel and K. Andeen, PoS(ICRC2019)179 (these proceedings). 8

Pith tools

Reviewed August 14, 2026 · model on record in the stance chip above.

[1] [1]

IceCube Collaboration, M. G. Aartsen et al., JINST 12 (2017) P03012

2017

[2] [2]

Abbasi et al., Nucl

IceCube Collaboration, R. Abbasi et al., Nucl. Instr . and Meth. A 700 (2013) 188–220

2013

[3] [3]

Heck et al., CORSIKA: A Monte Carlo code to simulate extensive air showers, Report FZKA 6019, F orschungszentrum Karlsruhe, 1998

D. Heck et al., CORSIKA: A Monte Carlo code to simulate extensive air showers, Report FZKA 6019, F orschungszentrum Karlsruhe, 1998

1998

[4] [4]

Battistoni et al., AIP Conference Proceedings 896 (2007) 31–49

G. Battistoni et al., AIP Conference Proceedings 896 (2007) 31–49

2007

[5] [5]

E. Ahn, R. Engel, T. Gaisser, P. Lipari, and T. Stanev, Physical Review D 80 (2009) 94003

2009

[6] [6]

Andeen and M

IceCube Collaboration, K. Andeen and M. Plum, PoS(ICRC2019)172 (these proceedings)

work page

[7] [7]

IceCube Collaboration, arXiv:1906.04317

arXiv 1906

[8] [8]

Pedregosa et al., Journal of Machine Learning Research 12 (2011) 2825–2830

F. Pedregosa et al., Journal of Machine Learning Research 12 (2011) 2825–2830

2011