  With the recent release of the second gravitational-wave transient catalogue
(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
hyper-parameters of the model, such as the spread in metallicity
$\sigma_{\text{Z}}$ and the parameter $\sigma_{\text{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 $\sigma_{\text{sp}}=0.1$ and
$\sigma_{\text{Z}}=0.4$. Models with higher spin magnitudes tend to strongly
favour dynamically formed BBHs ($f\le{}0.1$ if $\sigma_{\text{sp}}=0.3$).
Furthermore, we show that hyper-parameters controlling the rates of the model,
such as $\sigma_{\rm Z}$, have a large impact on the inference of the mixing
fraction, which rises from $0.18$ to $0.43$ when we increase
$\sigma_{\text{Z}}$ from 0.2 to 0.6, for a fixed value of
$\sigma_{\text{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\%$ credible interval, except when the spin
magnitude is high.
