With the recent release of the Second Gravitational Wave Transient Catalog (GWTC-2), which introduced dozens of new detections, we are at a turning point of gravitational wave astronomy, as we are now able to directly infer constraints on the astrophysical population of compact objects. Here, we tackle the burning issue of understanding the origin of binary black hole (BBH) mergers. To this effect, we make use of state-of-the-art population synthesis and N-body simulations, to represent two distinct formation channels: BBHs formed in the field (isolated channel) and in young star clusters (dynamical channel). We then use a Bayesian hierarchical approach to infer the distribution of the mixing fraction f, with f = 0 (f = 1) in the pure dynamical (isolated) channel. We explore the effects of additional hyperparameters of the model, such as the spread in metallicity σZ and the parameter σsp, describing the distribution of spin magnitudes. We find that the dynamical model is slightly favoured with a median value of f = 0.26, when σsp = 0.1 and σZ = 0.4. Models with higher spin magnitudes tend to strongly favour dynamically formed BBHs (f ≤ 0.1 if σsp = 0.3). Furthermore, we show that hyperparameters controlling the rates of the model, such as σZ, have a large impact on the inference of the mixing fraction, which rises from 0.18 to 0.43 when we increase σZ from 0.2 to 0.6, for a fixed value of σsp = 0.1. Finally, our current set of observations is better described by a combination of both formation channels, as a pure dynamical scenario is excluded at the $99{{\ \rm per\ cent}}$ credible interval, except when the spin magnitude is high.
