I present a deviation of the differential cross section, that (at least i hope so) makes things clearer. We rely on the following definition:
$$
\mathrm d P = \mathrm d \sigma \times F.
$$
Here $\mathrm d P$ is the transition probability per time and volume, with resulting particles $\{1,\dots,n\}$ carrying momenta $\{p_{f1},\dots, p_{fn}\} = p_f$. The incoming particles are labeled (as in P&S) with $A$ and $B$. F is the (relativistic) incoming flux density. We should express the flux density in terms of Lorentz covariants, in order to get a Lorentz covariant form of $\mathrm d \sigma$. We will come back to this later.
Lets start with the Transition probability. The incoming particles in the "infinite past" are considered as free particles (aswell as the final states in the "infinite future"). We define an initial state as the following wavepacket:
$$
\vert i \rangle = \vert \phi_A \phi_B \rangle_\text{in} = \int \dfrac{\mathrm d^3 k_A}{(2\pi)^3 2 k_A^0} \int \dfrac{\mathrm d^3 k_B}{(2\pi)^3 2 k_B^0} \tilde \phi_A(\vec k_A) \tilde \phi_B(\vec k_B) \vert \vec k_A \vec k_B \rangle_{\text{in}}.
$$
The typical restriction $k_i^0 = \sqrt{\vec k_i^2 + m_i^2}$ is made implicitly. Further, we use the abbreviation
$$
\mathrm{d} \tilde p = \dfrac{\mathrm{d}^3 p}{(2\pi)^3 2 p^0}
$$
from now on. The final state, can be a definite momentum eigenstate, as long as the detectors mainly measure momentum, and do not resolve positions at the level of de Broglie wavelenghts. E.g.
$$
\vert f \rangle = \vert \vec p_{f1},\dots,\vec p_{fn} \rangle_\text{out}.
$$
Before we continue, we make a few general statements on the initial wavepackets. The normalization of the momentum eigenstates is given as
$\langle \vec p \vert \vec p^\prime\rangle = (2\pi)^3 2 p^0 \delta(\vec p - \vec p^\prime)
$ and thus for a wavepacket
$$
\vert \phi \rangle = \int \mathrm d \tilde p \tilde \phi(\vec p) \vert \vec p \rangle,
$$
we obtain
$$
1 \stackrel{!}{=}\langle \phi \vert \phi \rangle = \int \mathrm d \tilde p \vert \tilde \phi(\vec p) \vert^2.
$$
We can define $\tilde \phi(\vec p)$ as Dirac sequence, e.g.
$$
\tilde \phi_\epsilon(\vec p) = \dfrac{1}{\sqrt{2\pi \epsilon}} \exp\left(-\dfrac{(\vec p - \vec p_A)^2}{2 \epsilon} \right) \stackrel{\epsilon \to 0}{\to} \delta(\vec p - \vec p_A).
$$
I.e. we can construct a wavepacket, with momenta tightly concentrated around a definite momentum $\vec p_A$:
$$
\vert \phi_\epsilon \rangle \stackrel{\epsilon \to 0}{\to} \vert \vec p_A \rangle.
$$
Let now $\vert \phi \rangle$ a wavepacket with concentrated momenta around $\vec p_A$. In QFT the position representation of such a wavepacket is defined as
$$
\phi(x) = \langle 0 \vert \Phi(x) \vert \phi \rangle = \int \mathrm d \tilde p e^{-ipx} \tilde \phi(\vec p),
$$
with the scalar field operator
$$
\Phi(x) = \int \mathrm d \tilde p \left(e^{-ipx} a(\vec p) + e^{ipx} b^\dagger(\vec p) \right).
$$
One could also use spinor or vector fields. For simplicity the further deviation is presented, using a scalar field operator. The final results are, however, valid for general fields. Now we define the following current (similar to the probability current in QM):
\begin{align}
j^\mu(x) &= i \left( \phi(x) \partial^\mu \phi^*(x) - \phi^*(x) \partial^\mu \phi(x) \right) \\
&= \int \mathrm d \tilde p \int \mathrm d \tilde p^\prime \tilde \phi^*(\vec p) \tilde \phi(\vec p^\prime) (p^\mu + p^{\prime \mu}) e^{i(p-p^\prime)x} \\
&\simeq 2 p_A^\mu \int \mathrm d \tilde p \int \mathrm d \tilde p^\prime \tilde \phi^*(\vec p) \tilde \phi(\vec p^\prime) e^{i(p-p^\prime)x}\\
&= 2 p^\mu_A \vert \phi(x) \vert^2.
\end{align}
We interpret $j^0(x) = \rho(x)$ as particle density and $\vec j(x) = (j^1(x),j^2(x),j^2(x))$ as particle flux density. The transition probability $P(A,B \to 1,\dots,n \text{ with momenta }p_f)$ is given by
\begin{align}
P(A,B \to 1,\dots,n \text{ with momenta }p_f) = \left(\prod_f \mathrm d \tilde p_f\right) \vert {}_{\text{out}}\langle \vec p_f \vert \phi_A \phi_B \rangle_{\text{in}}\vert^2.
\end{align}
It holds
\begin{align}
{}_{\text{out}}\langle \vec p_f \vert \vec k_A \vec k_B \rangle_\text{in} &= \langle \vec p_f \vert S \vert \vec k_A \vec k_B \rangle \\
&= \underbrace{\langle \vec p_f \vert \vec k_A \vec k_B \rangle}_{\text{unscattered part}} + \langle \vec p_f \vert \underbrace{S- \mathbb 1}_{T} \vert \vec k_A \vec k_B \rangle
\end{align}
(For a detailed definition of S, and the "in" and "out" states see P&S). The transition amplitude $\mathcal M$ is now defined as
\begin{align}
\langle \vec p_f \vert T \vert \vec k_A \vec k_B \rangle = i (2\pi)^4 \delta(k_A +k_B - {\textstyle\sum} p_f) \mathcal M(\{k_i\} \to p_f).
\end{align}
The physical interesting case is that the unscattered part vanishes (i.e. $\langle p_f \vert \phi_A \phi_B\rangle \simeq 0$). We furthermore assume, that the wavepackets $\vert \phi_A \phi_B \rangle$ are concentrated around $\vec p_A$ and $\vec p_B$. We obtain
\begin{align}
P(A,B \to 1,\dots,n \text{ with momenta }p_f) &=\left(\prod_f \mathrm d \tilde p_f\right) \left( \prod_{i=A,B} \int \mathrm d \tilde k_i \int \mathrm d \tilde k_i^\prime \tilde \phi_i^*(\vec k_i^\prime) \tilde \phi_i(\vec k_i) \right) \\
&\quad \times (2\pi)^4 \delta(k_A + k_B - {\textstyle\sum} p_f) (2\pi)^4 \delta(k_A^\prime + k_B^\prime - {\textstyle\sum} p_f) \\
&\quad \times \mathcal M^*(\{k_i^\prime\} \to p_f) \mathcal M(\{k_i\} \to p_f).
\end{align}
Since the momenta of the wavepackets are concentrated around $\vec p_A$ and $\vec p_B$ it holds
\begin{align}
\mathcal M^*(\{k_i^\prime\} \to p_f) \mathcal M(\{k_i\} \to p_f) \simeq \vert \mathcal M(\{\vec p_A ,\vec p_B\} \to p_f ) \vert^2,
\end{align}
and
\begin{align}
(2\pi)^4 \delta(k_A + k_B - {\textstyle\sum} p_f) (2\pi)^4 \delta(k_A^\prime + k_B^\prime - {\textstyle\sum} p_f) &= (2\pi)^4 \delta(k_A + k_B -(k_A^\prime + k_B^\prime)) (2\pi)^4 \delta(k_A^\prime + k_B^\prime-{\textstyle\sum} p_f) \\
&\simeq \int \mathrm d^4 x e^{ix(k_A^\prime + k_B^\prime - k_A -k_B)} (2\pi)^4 \delta(p_A + p_B - {\textstyle\sum} p_f) .
\end{align}
Back inserting yields
\begin{align}
P(A,B \to 1,\dots,n \text{ with momenta }p_f) &=\left(\prod_f \mathrm d \tilde p_f\right) (2\pi)^4 \delta(p_A + p_B -{\textstyle\sum} p_f) \\ &\quad \times \int \mathrm d^4 x \vert \phi_A(x) \vert^2 \vert \phi_B(x) \vert^2 \vert \mathcal M(\{\vec p_A ,\vec p_B\} \to p_f ) \vert^2.
\end{align}
Comparison with the definition of $\mathrm d P$ yields
\begin{align}
\mathrm d P = \mathrm d \Phi_f \vert \phi_A(x) \vert^2 \vert \phi_B(x) \vert^2 \vert \mathcal M(\{\vec p_A ,\vec p_B\} \to p_f ) \vert^2,
\end{align}
with
\begin{align}
\mathrm d \Phi_f = \left(\prod_f \mathrm d \tilde p_f\right) (2\pi)^4 \delta(p_A + p_B -{\textstyle\sum}p_f),
\end{align}
the LIPS (Lorentz invariant phase space) element. We define the incoming flux density
as
\begin{align}
F = \sqrt{(j_A \cdot j_B)^2 - j_A^2 j_B^2}.
\end{align}
This is manifest Lorentz invariant, since $j_i^\mu$ are Lorentz vectors by construction. With the result from above we obtain
\begin{align}
F &= \sqrt{(4 \vert \phi_A(x)\vert^2 \vert \phi_B(x)\vert^2 p_A \cdot p_B)^2 - 16 \vert \phi_A(x)\vert^4 \vert \phi_B(x)\vert^4 p_A^2 p_B^2 \vert} \\
&= 4 \vert \phi_A(x)\vert^2 \vert \phi_B(x)\vert^2 \sqrt{(p_A \cdot p_B)^2 - m_A^2 m_B^2}.
\end{align}
The final result for the differential cross section is
\begin{align}
\mathrm d \sigma = \mathrm d \Phi_f \dfrac{\vert \mathcal M(\{p_A, p_B\} \to p_f)\vert^2}{4 \sqrt{(p_A \cdot p_B)^2 - m_A^2 m_B^2}}.
\end{align}
This expression is manifest Lorentz invariant. Lets see, if we can restore the result of P&S, where the incoming particles have collinear velocities $\vec v_A, \vec v_B$, e.g. in the center of masse frame of the incoming particles. In this case $p_A$ and $p_B$ are given as
\begin{align}
p_A = (E_A,\vec p), \quad p_B = (E_B,-\vec p).
\end{align}
Inserting yields
\begin{align}
4 \sqrt{(p_A\cdot p_B)^2 - m_A^2 m_B^2} &= 4 \sqrt{(E_A E_B + \vec p^2)^2 - (E_A^2 - \vec p^2)(E_B^2 - \vec p^2)} \\
&= 4 E_A E_B \sqrt{\vec v_A^2 + \vec v_B^2 - 2 \vec v_A \vec v_B} \\
&= 2 E_A 2E_B \vert \vec v_A - \vec v_B \vert. \\
\Rightarrow \mathrm d \sigma &= \mathrm d \Phi_f\dfrac{\vert \mathcal M(\{p_A, p_B\} \to p_f)\vert^2}{2 E_A 2E_B \vert \vec v_A - \vec v_B \vert}.
\end{align}
This agrees with the formula for the differential cross section, given in (4.79) in P&S. This specific expression for the differential cross section, is only form invariant under boosts in the direction of $\vec p$. But the general result above is valid in all frames. The only problematic statement P&S made, in my opinion, is to call $\vert \vec v_A - \vec v_B \vert$ the relative velocity. In special relativity, $\vert \vec v_A - \vec v_B \vert$ cannot be interpreted as "relative velocity". E.g. with $\vec v_A = c \hat e$ and $\vec v_B = -c \hat e$, we obtain
\begin{align}
\vert \vec v_A - \vec v_B \vert = 2c.
\end{align}
A proper definition of the relativistic relative velocity is
\begin{align}
v_\text{rel} = \dfrac{\sqrt{(p_A \cdot p_B)^2 - m_A^2 m_B^2}}{p_A \cdot p_B} = \dfrac{\sqrt{(\vec v_A -\vec v_B)^2 - (\vec v_A \times \vec v_B)^2}}{1-\vec v_A \cdot \vec v_B}.
\end{align}
