The question how well nonrelativistic QCD (NRQCD) factorization can describe quarkonium production has been subject to debate since its invention. We review our recent reanalysis of a classic next-to-leading order color octet (CO) long distance matrix element (LDME) fit to large transverse momentum \(p_T\) \(J/\psi\) and \(\eta_c\) LHC production data. Our analysis differs from previous analysis of this kind not only by implementing for the first time a systematic treatment of scale uncertainties, but also by scrutinizing a much broader range of observables. Surprisingly, \(J/\psi\) hadroproduction is well described up to the highest measured values of \(p_T\). Potential NRQCD based relations nontrivially lead to a perfect description of \(\Upsilon(nS)\) production data. Furthermore, \(J/\psi\) production in \(\gamma p\) and \(\gamma \gamma\) collisions is, contrary to prevailing conceptions, reproduced down to \(p_T=1\) GeV, as long as the region of large inelasticity \(z\) is excluded. The overall picture is much rosier than usually perceived, the more so as the remaining discrepancies appear in phase space regions where solutions via varying kinds of resummations have been proposed.