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$ quarks 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}$.