Why a factor of 1/2 in cosmological Boltzmann equation for the phase space density?

In "The Early Universe" book by Kolb and Turner the Boltzmann equation is written for FRWL metric. The equation states:

begin{align}
Efrac{partial f}{partial t}- frac{dot{R}}{R}|vec{p}|^2frac{partial f}{partial E}=mathbf{C}[f]
end{align}
Now to get the result for the number density $n=g_iint frac{d^3p}{(2pi)^3} f(p)$ we divide by $E$ (and not by $2E$) and integrate the equation with the measure $g_ifrac{d^3p}{(2pi)^3}=2E_i dPi_i$.
We get the familiar result for processes ${a_i}rightarrow {a_f}$ (assuming the species $x$ considered is in the inital state):
$$dot{n}+3Hn=-int prod_{iin {inital}}d Pi_i f_i(p_i)prod_{fin {final}}d Pi_f (1pm f_f(p_f))(2pi)^4delta^4(sum_ip_i-sum_f p_f)|mathcal{M}_{{inital}rightarrow{final}}|^2+int prod_{iin {inital}}d Pi_i (1pm f_i(p_i))prod_{fin {final}}d Pi_f f_f(p_f)(2pi)^4delta^4(sum_ip_i-sum_f p_f)|mathcal{M}_{{inital}leftarrow{final}}|^2$$
But what is the expression for the phase space density $f(p)$. My naive expectation would say the same expression but without integrating over the species under consideration $dPi_x$ but this would mean that the collision term is, because we need to cancel the factor from the Lorentz invariant measure:
$$mathbf{C}[f]=-frac{1}{2}f_x(p)int prod_{iin {inital}, ineq x}d Pi_i f_i(p_i)prod_{fin {final}}d Pi_f (1pm f_f(p_f))(2pi)^4delta^4(sum_ip_i-sum_f p_f)|mathcal{M}_{{inital}rightarrow{final}}|^2+frac{1}{2}(1pm f_x(p))int prod_{iin {inital}}d Pi_i (1pm f_i(p_i))prod_{fin {final}}d Pi_f f_f(p_f)(2pi)^4delta^4(sum_ip_i-sum_f p_f)|mathcal{M}_{{inital}leftarrow{final}}|^2$$
I would have not expected this factor of $1/2$ in the above (but also according to this paper link) the factor appears. How can we understand this intuitively?
I always thought of collision term of the Boltzmann equation as the distribution version of the cross section.