We evaluate the production cross sections of $X(3872)$ at the LHC and Tevatron at NLO in $alpha_s$ in NRQCD by assuming that the short-distance production proceeds dominantly through its $chi_{c1}$ component in our $chi_{c1}mbox{-}D^0bar{D}^{*0}$ mixing model for $X(3872)$. The outcomes of the fits to the CMS $p_T$ distribution can well account for the recent ATLAS data in a much larger range of transverse momenta ($10~mbox{GeV}<p_T<70~mbox{GeV}$), and the CDF total cross section data, and are also consistent with the value of $k=Z_{cbar c}cdot Br(Xto J/psipi^+pi^-)$ constrained by the $B$-meson decay data. %It can also well describe the behavior of the CDF $psi(2S)$ data, which show a strong %resemblance to that of the X(3872). For LHCb the predicted X(3872) total cross section is larger than the data by a factor of 2, which is due to the problem of the fixed-order NRQCD calculation that may not be applicable for the region with small $p_T$ ($p_Tsim 5 ~mbox{GeV}$) and large forward rapidity $(2.5<y<4.5)$. In comparison, the prediction of molecule production mechanism for $X(3872)$ is inconsistent with both $p_T$ distributions and total cross sections of CMS and ATLAS, and the total cross section of CDF.