  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. %that controls the proportion of isolated and
dynamical BBHs. 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.
