  We present the first effective-one-body (EOB) model for generic-spins
quasi-circular black-hole--neutron-star (BHNS) inspiral-merger-ringdown
gravitational waveforms (GWs). Our model is based on a new numerical-relativity
(NR) informed expression of the BH remnant and its ringdown. It reproduces the
NR $(\ell,m)=(2,2)$ waveform with typical phase agreement of
${\lesssim0.5}\,$rad (${\lesssim1}\,$rad) to merger (ringdown). The maximum
(minimum) mismatch between the $(2,2)$ and the NR data is 4% (0.6%). Higher
modes (HMs) $(2,1)$, $(3,2)$, $(3,3)$, $(4,4)$ and $(5,5)$ are included and
their mismatch with the available NR waveforms are up to (down to) a 60% (1%)
depending on the inclination. Phase comparison with a 16 orbit precessing
simulation shows differences within the NR uncertainties. We demonstrate the
applicability of the model in GW parameter estimation by perfoming the first
BHNS Bayesian analysis with HMs (and non-precessing spins) of the event
GW190814, together with new $(2,2)$-mode analysis of GW200105 and GW200115. For
the GW190814 study, the inclusion of HMs gives tighter parameter posteriors.
The Bayes factors of our analyses on this event show decisive evidence for the
presence of HMs, but no clear preference for a BHNS or a binary black hole
(BBH) source. Similarly, we confirm GW200105 and GW200115 show no evidence for
tidal effects.
