In relativistic nuclear collisions the production of hadrons with light (u,d,s) quarks is quantitatively described in the framework of the Statistical Hadronization Model (SHM). Charm quarks are dominantly produced in initial hard collisions but interact strongly in the hot fireball and thermalize. Therefore charmed hadrons can be incorporated into the SHM by treating charm quarks as 'impurities' with thermal distributions, while the total charm content of the fireball is fixed by the measured open charm cross section. We call this model SHMc and demonstrate that with SHMc the measured multiplicities of single charm hadrons in lead-lead collisions at LHC energies can be well described with the same thermal parameters as for (u,d,s) hadrons. Furthermore, transverse momentum distributions are computed in a blast-wave model, which includes the resonance decay kinematics. SHMc is extended to lighter collision systems down to oxygen-oxygen and includes doubly- and triply-charmed hadrons. We show predictions for production probabilities of such states exhibiting a characteristic and quite spectacular enhancement hierarchy.