  We describe a model of the short gamma-ray burst (SGRB) population under a
`quasi-universal jet' scenario in which jets can differ in their on-axis peak
prompt emission luminosity $L_c$, but share a universal angular luminosity
profile $\ell(\theta_v)=L(\theta_v)/L_c$ as a function of the viewing angle
$\theta_v$. The model is fitted, through a Bayesian hierarchical approach
inspired by gravitational wave (GW) population analyses, to 3 observed SGRB
samples simultaneously: the Fermi/GBM sample of SGRBs with spectral information
in the catalogue (367 events); a flux-complete sample of 16 Swift/BAT SGRBs
also detected by GBM, with a measured redshift; and a sample of SGRBs with a
binary neutron star (BNS) merger counterpart, which only includes GRB~170817A
at present. The results favour a narrow jet core with half-opening angle
$\theta_c=2.1_{-1.4}^{+2.4}$ deg (90\% credible intervals from our fiducial
`full sample' analysis) whose on-axis peak luminosity is distributed as $p(L_c)
\propto L_c^{-A}$ with $A=3.2_{-0.4}^{+0.7}$ above a minimum luminosity
$L_c^\star = 5_{-2}^{+11}\times 10^{51}$ erg s$^{-1}$. For $\theta_v>\theta_c$,
the luminosity scales as a power law $\ell\propto \theta_v^{-\alpha_L}$ with
$\alpha_L=4.7_{-1.4}^{+1.2}$, with no evidence for a break. While the model
implies an intrinsic `Yonetoku' correlation between $L$ and the peak photon
energy $E_p$, its slope is somewhat shallower $E_p\propto L^{0.4\pm 0.2}$ than
the apparent one, and the normalization is offset towards larger $E_p$, due to
selection effects. The implied local rate density of SGRBs is between about 100
up to several thousands of events per Gpc$^{3}$ yr, in line with the BNS merger
rate density inferred from GW observations. Based on the model, we predict 0.2
to 1.3 joint GW+SGRB detections per year by the Advanced GW detector network
and Fermi/GBM during the O4 observing run.
