We analyse the high-quality Hertzsprung-Russell diagram (HRD) derived from Gaia data release 2 for the Solar Neighbourhood. We start building an almost complete sample within 200 pc and for |b| > 25(circle), so as to limit the impact of known errors and artefacts in the Gaia catalogue. Particular effort is then put into improving the modelling of population of binaries, which produce two marked features in the HRD: the sequence of near-equal mass binaries along the lower main sequence, and the isolated group of hot subdwarfs. We describe a new tool, BinaPSE, to follow the evolution of interacting binaries in a way that improves the consistency with PARSEC evolutionary tracks for single stars. BinaPSE is implemented into the TRILEGAL code for the generation of 'partial models' for both single and binary stellar populations, taking into account the presence of resolved and unresolved binaries. We then fit the Gaia HRD via Markov chain Monte Carlo methods that search for the star formation history and an initial binary fraction (by mass) that maximize the likelihood. The main results are (i) the binary fraction derived from the lower MS is close to 0.4, while twice larger values are favoured when the upper part of the HRD is fitted; (ii) present models predict the observed numbers of hot subdwarfs to within a factor of 2; and (iii) irrespective of the prescription for the binaries, the star formation rate peaks at values similar to 1.5 x 10(-4) M(circle dot)yr(-1) at ages slightly above 2 Gyr, and then decreases to similar to 0.8 x 10(-4) M(circle dot)yr(-1) at very old ages.
