This work presents the differential charmonia production cross sections in high energy p+p collisions calculated using NRQCD formalism. The NRQCD formalism, factorizes the quarkonia production cross sections in terms of short distance QCD cross sections and long distance matrix elements (LDMEs). The short distance crosssections are calculated in terms of perturbative QCD and LDMEs are obtained by fitting the experimental data. Measured transverse momentum distributions of \(\chi_{\rm c}\), \(\psi\)(2S) and J/\(\psi\) in p +{\(\bar {\rm p}\)} collisions at \(\sqrt{s}=\) 1.8, 1.96 TeV and in p+p collisions at \(\sqrt{s}=\) 7, 8 and 13 TeV are used to constrain LDMEs. The feed-down contribution to each state from the higher states are taken into account. The formalism provides a very good description of the data in a wide energy range. The values of LDMEs are used to predict the charmonia cross sections in p+p collisions at 13 TeV in kinematic bins relevant for the LHC detectors.