I wonder whether there is a systematic approach to find (or at least whether there are criteria for the existence of) vectors $P_0, P_1, \dots, P_n$, say in $n$-dimensional Minkowski space of signature $(1, n-1)$, such that:
$$P_0{}^2 = m_0, \quad P_i{}^2 = p_i{}^2 - 2 m_i - m_0, \\ P_0 \cdot P_i = m_i, \quad P_i \cdot P_j = p_i \cdot p_j - (m_i + m_j) - m_0,$$
where $m_0$, $m_1$, ..., $m_n$ are given constants and $p_1$, ..., $p_n$ are given (independent) vectors. My gut tells me that for $n \geq 2$ there ought to be such $\{ P_0, P_1, \dots, P_n \}$, since they have $(n+1)n$ degrees of freedom, while there are only $\frac{(n+1)(n+2)}{2}$ constraints, however, I struggle to find a way to express them in terms of $\{ m_0, m_1, \dots, m_n \}$ and $\{ p_1, \dots, p_n \}$?
Since these are all scalar products, they're invariant under Lorentz transformations of the $P_\alpha$. You can use this to choose a coordinate system in which $P_\alpha$ only has components up to $\alpha$, e.g. $P_0$ would be $(\sqrt{m_0},0,0,\ldots)$, $P_1$ would be $(m_1/\sqrt{m_0},\sqrt{m_1^2/m_0+2m_1+m_0-p_1^2},0,0,\ldots)$, and so on. Of course whether the radicands are non-negative depends on the parameters; if they aren't, there's no solution (since any solution could be transformed to this form).
Since you're looking for $n+1$ vectors, not $n$, there won't be a solution in general, unless the parameters are such that $P_n$ has the right length without requiring an $(n+1)$-th component.