Upper bounds on phase space momenta

Suppose I wish to calculate the phase space volume for the process $overline{X}X to A_1 A_2 A_3 A_4 A_5$ in the CM frame of $overline{X}, X$ so that $sqrt{s} = 2m_X$. The volume is given by

$$
V propto
int d^3 p_1
int d^3 p_2
int d^3 p_3
int d^3 p_4
int d^3 p_5
frac{1}{E_1} frac{1}{E_2} frac{1}{E_3} frac{1}{E_4} frac{1}{E_5}
delta[2m_X – (sum_{f=1}^5 E_f)] delta^{(3)}(sum_{f=1}^5 mathbf{p}_f).
$$

Define $mathbf{P} = sum_{f=1}^5 mathbf{p}_i$. Integrating out the 3D delta function we obtain

$$
V= (2pi)
int d^3 p_1
int d^3 p_2
int d^3 p_3
int d p_4 d(costheta_{34})
frac{1}{E_1} frac{1}{E_2} frac{1}{E_3} frac{p_4^2}{E_4} frac{1}{E_5(mathbf{P} = 0)}
frac{delta(costheta_{34} -z_0)}{|partial_{costheta_{34}} E_5(mathbf{P}=0)|_{z_0} }
$$

where $costheta_{34} := mathbf{p_3cdot p_4}/p_3p_4$ and $z_0$ is the value of $costheta_{34}$ such that

$$2m_X – E_1 – E_2 – E_3 -E_4 – E_5(mathbf{P}=0) = 0$$

and the factor of $1/{|partial_{costheta_{34}} E_5(mathbf{P}=0)|_{z_0} }$ is just the Jacobian we pick up from the identity

$$
delta(g(x)) = sum_{x_0}frac{delta(x-x_0)}{|g'(x_0)|}
$$

---

Now, after integrating out the delta function we must enforce the fact that the integral is to be $0$ whenever $|z_0|>1$. This amounts to forcing conservation on energy on the integrand.

Define the region is phase space $Xi$ as

$$
Xi(p_1, p_2, p_3, p_4, costheta_{14}, costheta_{12}, costheta_{23}) = begin{cases} 1 , -1 leq z_0 leq 1
0 , text{otherwise}
end{cases}
$$

Therefore, an equivalent expression for the phase space volume $V$ is given by

begin{align}
V &= (2pi)
{int_{mathbb{R}^3} d^3 p_1
int_{mathbb{R}^3} d^3 p_2
int_{mathbb{R}^3} d^3 p_3
int_0^infty d p_4 }
& Xi cdot frac{1}{E_1} frac{1}{E_2} frac{1}{E_3} frac{p_4^2}{E_4} frac{1}{E_5(mathbf{P} = 0)_{z_0}}
frac{1}{|partial_{costheta_{34}} E_5(mathbf{P}=0)|_{z_0} }.
end{align}

---

Doing this integral numerically proves to be very difficult, since I need to integrate over all of space. Hence every numerical integration scheme will be fooled into thinking this is zero since it doesn’t know that all the action is happening in the regions where $Xi =1$.

Question: Is there a finite upper bound on the limit of integration with respect to $dp_i$, such that I may entirely cover the region $Xi = 1$?

My Conjecture: My reasoning is that, yes, after enforcing conservation of energy via the $1D$-energy delta function the momenta $p_i$ are bound above by $p_i^{max} leq sqrt{(2 m_X )^2 – m_i^2}$. Which leads me to purpose the following conjecture that

begin{align}
frac{dV}{dOmega_1 dOmega_2 dOmega_3} &= (2pi)
{int_0^{p_1^{max}} d p_1
int_0^{p_2^{max}} d p_2
int_0^{p_3^{max}} d p_3
int_0^{p_4^{max}} d p_4 }
& Xi cdot frac{p_1^2}{E_1} frac{p_2^2}{E_2} frac{p_3^2}{E_3} frac{p_4^2}{E_4} frac{1}{E_5(mathbf{P} = 0)_{z_0}}
frac{ 1}{|partial_{costheta_{34}} E_5(mathbf{P}=0)|_{z_0} }.
end{align}