  Long Gamma-Ray Bursts (LGRBs) have been shown to be powerful probes of the
Universe, in particular to study the star formation rate up to very high
redshift ($z \sim 9$). Since LGRBs are produced by only a small fraction of
massive stars, it is paramount to have a good understanding of their underlying
intrinsic population in order to use them as cosmological probes without
introducing any unwanted bias. The goal of this work is to constrain and
characterise this intrinsic population. We developed a Monte Carlo model where
each burst is described by its redshift and its properties at the peak of the
lightcurve. We derived the best fit parameters by comparing our synthetic
populations to carefully selected observational constraints based on the
CGRO/BATSE, Fermi/GBM and Swift/BAT samples with appropriate flux thresholds.
We explored different scenarios in terms of cosmic evolution of the luminosity
function and/or of the redshift distribution as well as including or not the
presence of intrinsic spectral-energetics (Ep-L) correlations. We find that the
existence of an intrinsic Ep-L correlation is preferred but with a shallower
slope than observed($\alpha_A \sim 0.3$) and a larger scatter ($\sim 0.4$ dex).
We find a strong degeneracy between the cosmic evolution of the luminosity and
of the LGRB rate, and show that a sample both larger and deeper than SHOALS by
a factor of three is needed to lift this degeneracy. The observed We conclude
that Ep-L correlation cannot be explained only by selection effects although
these do play a role in shaping the observed relation. The degeneracy between
cosmic evolution of the luminosity function and of the redshift distribution of
LGRBs should be included in the uncertainties of star formation rate estimates;
these amount to a factor of 10 at $z=6$ and up to a factor of 50 at $z=9$.
