  We have performed the first $n_f = 2+1+1$ lattice QCD computations of the
properties (masses and decay constants) of ground-state charmonium mesons. Our
calculation uses the HISQ action to generate quark-line connected two-point
correlation functions on MILC gluon field configurations that include $u/d$
quark masses going down to the physical point, tuning the $c$ quark mass from
$M_{J/\psi}$ and including the effect of the $c$ quark's electric charge
through quenched QED. We obtain $M_{J/\psi}-M_{\eta_c}$ (connected) =
120.3(1.1) MeV and interpret the difference with experiment as the impact on
$M_{\eta_c}$ of its decay to gluons, missing from the lattice calculation. This
allows us to determine $\Delta M_{\eta_c}^{\mathrm{annihiln}}$ =+7.3(1.2) MeV,
giving its value for the first time. Our result of $f_{J/\psi}=$ 0.4104(17)
GeV, gives $\Gamma(J/\psi \rightarrow e^+e^-)$=5.637(49) keV, in agreement
with, but now more accurate than experiment. At the same time we have improved
the determination of the $c$ quark mass, including the impact of quenched QED
to give $\overline{m}_c(3\,\mathrm{GeV})$ = 0.9841(51) GeV. We have also used
the time-moments of the vector charmonium current-current correlators to
improve the lattice QCD result for the $c$ quark HVP contribution to the
anomalous magnetic moment of the muon. We obtain $a_{\mu}^c = 14.638(47) \times
10^{-10}$, which is 2.5$\sigma$ higher than the value derived using moments
extracted from some sets of experimental data on $R(e^+e^- \rightarrow
\mathrm{hadrons})$. This value for $a_{\mu}^c$ includes our determination of
the effect of QED on this quantity, $\delta a_{\mu}^c = 0.0313(28) \times
10^{-10}$.
