Learning the Composition of Ultra High Energy Cosmic Rays

We apply statistical inference on the Pierre Auger Open Data to discern for the first time the full mass composition of cosmic rays at different energies. Working with longitudinal electromagnetic profiles of cosmic ray showers, in particular their peaking depths $X_{\rm max}$, we employ central moments of the $X_{\rm max}$ distributions as features to discriminate between different shower compositions. We find that already the first few moments entail the most relevant information to infer the primary cosmic ray mass spectrum. Our approach, based on an unbinned likelihood, allows us to consistently account for sources of statistical uncertainties due to finite datasets, both measured and simulated, as well as systematic effects. Finally, we provide a quantitative comparison of different high energy hadronic interaction models available in the atmospheric shower simulation codes.