  Angular momentum and spin precession are expected to be generic features of a
significant fraction of binary black hole systems. As such, it is essential to
have waveform models that faithfully incorporate the effects of precession.
Here, we assess how well the current state of the art models achieve this for
waveform strains constructed only from the $\ell=2$ multipoles.Specifically, we
conduct a survey on the faithfulness of the waveform models %(approximants)
\texttt{SEOBNRv5PHM}, \texttt{TEOBResumS}, \texttt{IMRPhenomTPHM},
\texttt{IMRPhenomXPHM} to the numerical relativity (NR) surrogate
\texttt{NRSur7dq4} and to NR waveforms from the \texttt{SXS} catalog. The
former assessment involves systems with mass ratios up to six and dimensionless
spins up to 0.8. The latter employs $317$ short and $23$ long \texttt{SXS}
waveforms. For all cases, we use reference inclinations of zero and $90^\circ$.
We find that all four models become more faithful as the mass ratio approaches
unity and when the merger-ringdown portion of the waveforms are excluded. We
also uncover a correlation between the co-precessing $(2,\pm2)$ multipole
mismatches and the overall strain mismatch. We additionally find that for high
inclinations, precessing $(2,\pm 1)$ multipoles that are more faithful than
their $(2,\pm2)$ counterparts, and comparable in magnitude, improve waveform
faithfulness. As a side note, we show that use of uniformly-filled parameter
spaces may lead to an overestimation of precessing model faithfulness. We
conclude our survey with a parameter estimation study in which we inject two
precessing \texttt{SXS} waveforms (at low and high masses) and recover the
signal with \texttt{SEOBNRv5PHM}, \texttt{IMRPhenomTPHM} and
\texttt{IMRPhenomXPHM}. As a bonus, we present preliminary multidimensional
fits to model unfaithfulness for Bayesian model selection in parameter
estimation studies.
