  $J/\psi$ observables of p+p collisions at $\sqrt{s}=$ 200 GeV and 2.76 TeV
are calculated using the PYTHIA event generator for the production of $c$ and
$\bar c$ and the Wigner density formalism. The charm quark momentum from PYTHIA
is tuned in such a way that the transverse momentum and rapidity distribution
correspond to the Fixed-Order Next-to-Leading Logarithm (FONLL) calculations
and hence to the experimental data for open charm mesons. Using the Wigner
density of the charmonium wave function we calculate for each charm quark pair
the probability to form a charmonium. As a result, the experimental data on the
total $J/\psi$ yield, its rapidity and transverse momentum distribution as well
as the fractions of feed-down from excited states of charmonium are reproduced.
Applying the same approach to relativistic heavy-ion collisions, we find, if
quark gluon plasma (QGP) effects are ignored, a considerable enhancement of the
primordial $J/\psi$'s.
