  Knowledge of the shape of the mass spectrum of compact objects can be used to
help break the degeneracy between the mass and redshift of the gravitational
wave (GW) sources, and thus can be used to infer cosmological parameters in the
absence of redshift measurements obtained from electromagnetic observations. In
this paper, we study extensively different aspects of this approach, including
its computational limits and achievable accuracy. We focus on ground-based
detectors with current and future sensitivities, we first perform the analysis
of an extensive set of simulated data with a hierarchical Bayesian scheme
inferring population and cosmological parameters. We consider a population
model (power-law plus Gaussian) which exhibits characteristic scales (extremes
of the mass spectrum, presence of an accumulation point) that allows an
indirect estimate of the source redshift. Our analysis of this catalog
highlights and quantifies the tight interplay between source population and
cosmological parameters, as well as the influence of initial assumptions
(whether formulated on the source or cosmological parameters). We then validate
our results by an "end-to-end" analysis using simulated GW data and posterior
samples generated from Bayesian samplers used for GW parameter estimation, thus
mirroring the analysis chain used for observational data for the first time in
literature. Our results then lead us to re-examine the estimation of $H_0$
obtained with GWTC-1, and we show explicitly how population assumptions impact
the final $H_0$ result. Our results underline the importance of inferring
population and cosmological parameters jointly (and not separately as is often
assumed). The only exception, as we discuss, is if an electromagnetic
counterpart was to be observed for all the BBH events: then the population
assumptions have less impact on the estimation of cosmological parameters.
