We present a novel approach to derive constraints on neutrino masses , as well as on other cosmological parameters , from cosmological data , while taking into account our ignorance of the neutrino mass ordering . We derive constraints from a combination of current as well as future cosmological datasets on the total neutrino mass M _ { \nu } and on the mass fractions f _ { \nu,i } = m _ { i } / M _ { \nu } ( where the index i = 1 , 2 , 3 indicates the three mass eigenstates ) carried by each of the mass eigenstates m _ { i } , after marginalizing over the ( unknown ) neutrino mass ordering , either normal ordering ( NH ) or inverted ordering ( IH ) . The bounds on all the cosmological parameters , including those on the total neutrino mass , take therefore into account the uncertainty related to our ignorance of the mass hierarchy that is actually realized in nature . This novel approach is carried out in the framework of Bayesian analysis of a typical hierarchical problem , where the distribution of the parameters of the model depends on further parameters , the hyperparameters . In this context , the choice of the neutrino mass ordering is modeled via the discrete hyperparameter h _ { \mathrm { type } } , which we introduce in the usual Markov chain analysis . The preference from cosmological data for either the NH or the IH scenarios is then simply encoded in the posterior distribution of the hyperparameter itself . Current cosmic microwave background ( CMB ) measurements assign equal odds to the two hierarchies , and are thus unable to distinguish between them . However , after the addition of baryon acoustic oscillation ( BAO ) measurements , a weak preference for the normal hierarchical scenario 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 ) . Concerning next-generation cosmological experiments , 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 oscillation experiments , as for NH and a fiducial value of M _ { \nu } = 0.06 \mathrm { eV } there is a 9:1 preference of normal versus inverted hierarchy . On the contrary , if the sum of the masses is of the order of 0.1 \mathrm { eV } or larger , even future cosmological observations will be inconclusive . The innovative statistical strategy exploited here represents a very simple , efficient and robust tool to study the sensitivity of present and future cosmological data to the neutrino mass hierarchy , and a sound competitor to the standard Bayesian model comparison . The unbiased limit on M _ { \nu } we obtain is crucial for ongoing and planned neutrinoless double beta decay searches .