We re-analyse photometric near-infrared data in order to investigate why it is so hard to get a consensus for the shape and density law of the bulge, as seen from the literature. To solve the problem we use the Besancon Galaxy Model to provide a scheme for parameter fitting of the structural characteristics of the bulge region. The fitting process allows the determination of the global shape of the bulge main structure. We explore various parameters and shape for the bulge/bar structure based on Ferrers ellipsoids and fit the shape of the inner disc in the same process. The results show that the main structure is a quite standard triaxial boxy bar/bulge with an orientation of about 13 degree with respect to the Sun-centre direction. But the fit is greatly improved when we add a second structure, which is a longer and thicker ellipsoid. We emphasize that our first ellipsoid represent the main boxy bar of the Galaxy, and that the thick bulge could be either a classical bulge slightly flattened by the effect of the bar potential, or a inner thick disc counterpart. We show that the double clump seen at intermediate latitudes can be reproduced by adding a slight flare to the bar. In order to better characterize the populations, we further simulate several fields which have been surveyed in spectroscopy and for which metallicity distribution function (MDF) are available. The model is in good agreement with these MDF along the minor axis if we assume that the main bar has a mean solar metallicity and the second thicker population has a lower metallicity. It then creates naturally a vertical metallicity gradient by the mixing of the two poulations. (abridged)