  We provide perturbation theory predictions for the HI intensity mapping power
spectrum multipoles using the Effective Field Theory of Large Scale Structure
(EFTofLSS), which should allow us to constrain cosmological parameters
exploiting mildly nonlinear scales. Assuming survey specifications typical of
proposed interferometric HI intensity mapping experiments like CHORD and PUMA,
and realistic ranges of validity for the perturbation theory modelling, we run
mock full shape MCMC analyses at a redshift bin centred at $z=0.5$, and compare
with Stage-IV optical galaxy surveys. We include the impact of 21cm foreground
removal using simulations-based prescriptions, and quantify the effects on the
precision and accuracy of the parameter estimation. We vary 11 parameters in
total: 3 cosmological parameters, 7 bias and counterterms parameters, and the
HI brightness temperature. Amongst them, the 4 parameters of interest are: the
cold dark matter density, $\omega_{\rm c}$, the Hubble parameter, $h$, the
primordial amplitude of the power spectrum, $A_{\rm s}$, and the linear HI
bias, $b_1$. For the best case scenario, we obtain unbiased constraints on all
parameters with $<3\%$ errors at $68\%$ confidence level. When we include the
foreground removal effects, the parameter estimation becomes strongly biased
for $\omega_{\rm c}, h$, and $b_1$, while $A_{\rm s}$ is less biased ($<
2\sigma$). We find that scale cuts $k_{\rm min} \geq 0.03 \, h/{\mathrm{Mpc}}$
are required to return accurate estimates for $\omega_{\rm c}$ and $h$, at the
price of a decrease in the precision, while $b_1$ remains strongly biased. We
comment on the implications of these results for real data analyses.
