We present a self-consistent Bayesian formalism to sample the primordial density fields compatible with a set of dark matter density tracers after a cosmic evolution observed in redshift space. Previous works on density reconstruction did not self-consistently consider redshift space distortions or included an additional iterative distortion correction step. We present here the analytic solution of coherent flows within a Hamiltonian Monte Carlo posterior sampling of the primordial density field. We test our method within the Zel'dovich approximation, presenting also an analytic solution including tidal fields and spherical collapse on small scales. Our resulting reconstructed fields are isotropic and their power spectra are unbiased compared to the true field defined by our mock observations. Novel algorithmic implementations are introduced regarding the mass assignment kernels when defining the dark matter density field and optimization of the time-step in the Hamiltonian equations of motions. Our algorithm, dubbed BARCODE, promises to be specially suited for analysis of the dark matter cosmic web down to scales of a few megaparsecs. This large-scale structure is implied by the observed spatial distribution of galaxy clusters - such as obtained from X-ray, Sunyaev-Zel'dovich, or weak lensing surveys - as well as that of the intergalactic medium sampled by the Ly α forest or perhaps even by deep hydrogen intensity mapping. In these cases, virialized motions are negligible, and the tracers cannot be modelled as point-like objects. It could be used in all of these contexts as a baryon acoustic oscillation reconstruction algorithm.