Neutrino Modelling in Friedmann Equation
The energy density of a Fermi gas is
$$\rho_{\nu}= \int \rho(p)\ dp = \int E(p)F(p)g(p)\ dp $$
$$\rho_{\nu}  = \int \left(\sqrt{p^2 c^2 + m^2 c^4}\right)\left(\exp (E/k_BT) + 1\right)^{-1} \left(g_s 4\pi p^2/h^3\right)\ dp$$
in units of energy per unit volume.

Before neutrino decoupling at $k_B T \sim 1$ MeV, the neutrinos are ultrarelativistic with $pc \gg m_{\nu}c^2$. After decoupling, the shape of the occupation index function $F(p)$ does not change - so $F(p) = \left(pc/k_BT_{\nu} +1\right)^{-1}$ in subsequent evolution.

Thus
$$\rho_{\nu} = \frac{g_s c}{h^3} \int \frac{ \sqrt{p^2 + m_{\nu}^2c^2}}
{\exp(pc/k_BT_{\nu}) +1}\ 4\pi p^2\ dp$$

I don't understand where your $(2\pi)^3$ comes from, other than to suggest that the unit system is actually $\hbar = 1$.