  We present a novel approach to derive constraints on neutrino masses from
cosmological data, while taking into account our ignorance of the neutrino mass
ordering. We derive constraints from a combination of current and future
cosmological datasets on the total neutrino mass $M_\nu$ and on the mass
fractions carried by each of the mass eigenstates, after marginalizing over the
(unknown) neutrino mass ordering, either normal (NH) or inverted (IH). The
bounds take therefore into account the uncertainty related to our ignorance of
the mass hierarchy. This novel approach is carried out in the framework of
Bayesian analysis of a typical hierarchical problem. In this context, the
choice of the neutrino mass ordering is modeled via the discrete hyperparameter
$h_{type}$. The preference for either the NH or the IH scenarios is then
encoded in the posterior distribution of $h_{type}$ itself. Current CMB
measurements assign equal odds to the two hierarchies, and are thus unable to
distinguish between them. However, after the addition of BAO measurements, a
weak preference for NH appears, with odds of 4:3 from Planck temperature and
large-scale polarization in combination with BAO (3:2 if small-scale
polarization is also included). Forecasts suggest that the combination of
upcoming CMB (COrE) and BAO surveys (DESI) may determine the neutrino mass
hierarchy at a high statistical significance if the mass is very close to the
minimal value allowed by oscillations, as for NH and $M_\nu=0.06$ eV there is a
9:1 preference of NH vs IH. On the contrary, if $M_\nu$ is of the order of 0.1
eV or larger, even future cosmological observations will be inconclusive. The
unbiased limit on $M_\nu$ we obtain with this innovative statistical strategy
is crucial for ongoing and planned neutrinoless double beta decay searches.
