Context. The ability to detect and characterise an increasing variety of exoplanets has been made possible by the continuous development of stable, high-resolution spectrographs and the Doppler radial velocity (RV) method. The cross-correlation function (CCF) method is one of the traditional approaches used to derive RVs. More recently, template matching has been introduced as an advantageous alternative for M-dwarf stars. Aims. We describe a new implementation of the template matching technique for stellar RV estimation within a semi-Bayesian framework, providing a more statistically principled characterisation of the RV measurements and associated uncertainties. This methodology, named the Semi-Bayesian Approach for RVs with Template matching, S-BART, can currently be applied to HARPS and ESPRESSO data. We first validate its performance with respect to other template matching pipelines using HARPS data. We then apply S-BART to ESPRESSO observations, comparing the scatter and uncertainty of the derived RV time series with those obtained using the CCF method. We leave a full analysis of the planetary and activity signals present in the considered datasets for future work. Methods. In the context of a semi-Bayesian framework, a common RV shift is assumed to describe the difference between each spectral order of a given stellar spectrum and a template built from the available observations. Posterior probability distributions are obtained for the relative RV associated with each spectrum using the Laplace approximation, after marginalization with respect to the continuum. We also implemented, for validation purposes, a traditional template matching approach, where a RV shift is estimated individually for each spectral order and the final RV estimate is calculated as a weighted average of the RVs of the individual orders. Results. The application of our template-based methods to HARPS archival observations of Barnard’s star allowed us to validate our implementation against other template matching methods. Although we find similar results, the standard deviation of the RVs derived with S-BART is smaller than that obtained with the HARPS-TERRA and SERVAL pipelines. We believe this is due to differences in the construction of the stellar template and the handling of telluric features. After validating S-BART, we applied it to 33 ESPRESSO GTO targets, evaluating its performance and comparing it to the CCF method as implemented in ESO’s official pipeline. We find a decrease in the median RV scatter of ~10 and ~4% for M- and K-type stars, respectively. Our semi-Bayesian framework yields more precise RV estimates than the CCF method, in particular in the case of M-type stars where S-BART achieves a median uncertainty of ~15 cm s⁻¹ over 309 observations of 16 targets. Further, with the same data we estimated the nightly zero point (NZP) of the instrument, finding a weighted NZP scatter of below ~0.7 m s⁻¹. Given that this includes stellar variability, photon noise, and potential planetary signals, it should be taken as an upper limit on the RV precision attainable with ESPRESSO data.