  Orbital eccentricity is a crucial physical effect to unveil the origin of
compact-object binaries detected by ground- and spaced-based gravitational-wave
(GW) observatories. Here, we perform for the first time a Bayesian inference
study of inspiral-merger-ringdown eccentric waveforms for binary black holes
with non-precessing spins using two (instead of one) eccentric parameters:
eccentricity and relativistic anomaly. We employ for our study the multipolar
effective-one-body (EOB) waveform model SEOBNRv4EHM, and use initial conditions
such that the eccentric parameters are specified at an orbit-averaged
frequency. We show that this new parametrization of the initial conditions
leads to a more efficient sampling of the parameter space. We also assess the
impact of the relativistic-anomaly parameter by performing mock-signal
injections, and we show that neglecting such a parameter can lead to
significant biases in several binary parameters. We validate our model with
mock-signal injections based on numerical-relativity waveforms, and we
demonstrate the ability of the model to accurately recover the injected
parameters. Finally, using standard stochastic samplers employed by the
LIGO-Virgo-KAGRA Collaboration, we analyze a set of real GW signals observed by
the LIGO-Virgo detectors during the first and third runs. We do not find clear
evidence of eccentricity in the signals analyzed, more specifically we measure
$e^{\text{GW150914}}_{\text{gw, 10Hz}}= 0.08^{+0.09}_{-0.06}$,
$e^{\text{GW151226}}_{\text{gw, 20Hz}}= {0.04}^{+0.05}_{-0.04} $, and
$e^{\text{GW190521}}_{\text{gw, 5.5Hz}}= 0.15^{+0.12}_{-0.12}$.
