  We calculate the relativistic corrections to prompt $J/\psi$ pair and
$J/\psi+\psi(2S)$ hadroproduction through the color-singlet channel within the
framework of nonrelativistic QCD (NRQCD) factorization. The short-distance
coefficients are obtained by matching full-QCD and NRQCD calculations at the
partonic level, in which both squared amplitude and phase space are expanded in
$v^2$. We find that such an expansion of the phase space spoils the convergence
of NRQCD factorization near the production threshold. To fix this problem, we
propose to modify the matching between full QCD and NRQCD by adopting the
physical phase space. In this modified approach, the theoretical uncertainties
due to the choice of charm quark mass $m_c$ are largely reduced and the overall
agreement of our predictions with LHC data is significantly improved, both as
for total and differential cross sections.
