In this work, we extended our statistical model with charmed and bottomed hadrons, and fit the quark creational probabilities for the heavy quarks, using low energy inclusive charmonium and bottomonium data. With the finalized fit for all the relevant types of quarks (up, down, strange, charm, bottom) at the energy range from a few GeV up to a few tens of GeV’s, the model is now considered complete. Some examples are also given for proton–proton, pion–proton, and proton–antiproton collisions with charmonium, bottomonium, and open charm hadrons in the final state.