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 10 parameters in total: 3 cosmological and 7 bias and counterterms parameters. Amongst them, the 4 parameters of interest are: the cold dark matter density, ωc, the Hubble parameter, h, the primordial amplitude of the power spectrum, As, and the linear HI bias, b1. 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, with all parameters being >5σ away from the true values. We find that scale cuts kmin∼0.03h/Mpc are required to return accurate estimates for ωc and h, at the price of a decrease in the precision, while As and b1 remain biased. We comment on the implications of these results for current and forthcoming real data analyses.