These identities should work:
\begin{align*}
s&=4E^2
\\
t&=-2E^2(1-\cos\theta)=-4E^2\sin^2(\theta/2)
\\
u&=-2E^2(1+\cos\theta)=-4E^2\cos^2(\theta/2)
\end{align*}
Derivation of probability density for Bhabha scattering:
In a typical collider experiment the momentum vectors are
\begin{equation*}
p_1=\begin{pmatrix}E\\0\\0\\p\end{pmatrix}\qquad
p_2=\begin{pmatrix}E\\0\\0\\-p\end{pmatrix}\qquad
p_3=\begin{pmatrix}
E\\
p\sin\theta\cos\phi\\
p\sin\theta\sin\phi\\
p\cos\theta
\end{pmatrix}
\qquad
p_4=\begin{pmatrix}
E\\
-p\sin\theta\cos\phi\\
-p\sin\theta\sin\phi\\
-p\cos\theta
\end{pmatrix}
\end{equation*}
where $p=\sqrt{E^2-m^2}$.
The spinors are
\begin{gather*}
v_{11}=\begin{pmatrix}p\\0\\E+m\\0\end{pmatrix}\quad
u_{21}=\begin{pmatrix}E+m\\0\\-p\\0\end{pmatrix}\quad
v_{31}=\begin{pmatrix}p_3^z\\p_3^x+ip_3^y\\E+m\\0\end{pmatrix}\quad
u_{41}=\begin{pmatrix}E+m\\0\\p_4^z\\p_4^x+ip_4^y\end{pmatrix}
\\
v_{12}=\begin{pmatrix}0\\-p\\0\\E+m\end{pmatrix}\quad
u_{22}=\begin{pmatrix}0\\E+m\\0\\p\end{pmatrix}\quad
v_{32}=\begin{pmatrix}p_3^x-ip_3^y\\-p_3^z\\0\\E+m\end{pmatrix}\quad
u_{42}=\begin{pmatrix}0\\E+m\\p_4^x-ip_4^y\\-p_4^z\end{pmatrix}
\end{gather*}
The last digit in a spinor subscript is 1 for spin up and 2 for spin down.
Note that the spinors are not individually normalized.
Instead, a combined spinor normalization constant $N=(E+m)^4$
will be used where needed.
This is the probability density for Bhabha scattering.
The formula is from Feynman diagrams.
$$
|\mathcal{M}(s_1,s_2,s_3,s_4)|^2=\frac{e^4}{N}
\left|
-\frac{1}{t}(\bar{v}_1\gamma^\mu v_3)(\bar{u}_4\gamma_\mu u_2)
+\frac{1}{s}(\bar{v}_1\gamma^\nu u_2)(\bar{u}_4\gamma_\nu v_3)
\right|^2
$$
Symbol $s_j$ selects the spin (up or down) of spinor $j$.
Symbol $e$ is electron charge.
Symbols $s$ and $t$ are Mandelstam variables $s=(p_1+p_2)^2$ and $t=(p_1-p_3)^2$.
Let
\begin{equation*}
a_1=(\bar{v}_1\gamma^\mu v_3)(\bar{u}_4\gamma_\mu u_2)
\qquad
a_2=(\bar{v}_1\gamma^\nu u_2)(\bar{u}_4\gamma_\nu v_3)
\end{equation*}
Then
\begin{align*}
|\mathcal{M}(s_1,s_2,s_3,s_4)|^2
&=
\frac{e^4}{N}\left|{-\frac{a_1}{t}} + \frac{a_2}{s}\right|^2\\
&=
\frac{e^4}{N}\left(-\frac{a_1}{t} + \frac{a_2}{s}\right)\left(-\frac{a_1}{t} + \frac{a_2}{s}\right)^*\\
&=
\frac{e^4}{N}
\left(
\frac{a_1a_1^*}{t^2} - \frac{a_1a_2^*}{st} -
\frac{a_1^*a_2}{st} + \frac{a_2a_2^*}{s^2}
\right)
\end{align*}
The expected probability density $\langle|\mathcal{M}|^2\rangle$ is computed
by summing $|\mathcal{M}|^2$ over all spin states and dividing by the number of inbound states.
There are four inbound states.
\begin{align*}
\langle|\mathcal{M}|^2\rangle
&=
\frac{1}{4}\sum_{s_1=1}^2\sum_{s_2=1}^2\sum_{s_3=1}^2\sum_{s_4=1}^2
|\mathcal{M}(s_1,s_2,s_3,s_4)|^2\\
&=
\frac{e^4}{4}\sum_{s_1=1}^2\sum_{s_2=1}^2\sum_{s_3=1}^2\sum_{s_4=1}^2
\frac{1}{N}
\left(
\frac{a_1a_1^*}{t^2} - \frac{a_1a_2^*}{st} -
\frac{a_1^*a_2}{st} + \frac{a_2a_2^*}{s^2}
\right)
\end{align*}
Use the Casimir trick to replace sums over spins with matrix products.
\begin{align*}
f_{11}&=\frac{1}{N}\sum_\text{spins}a_1a_1^*=
\mathop{\rm Tr}\left(
(\not p_1-m)\gamma^\mu(\not p_3-m)\gamma^\nu
\right)
\mathop{\rm Tr}\left(
(\not p_4+m)\gamma_\mu(\not p_2+m)\gamma_\nu
\right)
\\
f_{12}&=\frac{1}{N}\sum_{\rm spins}a_1a_2^*=
\mathop{\rm Tr}\left(
(\not p_1-m)\gamma^\mu(\not p_2+m)\gamma^\nu
(\not p_4+m)\gamma_\mu(\not p_3-m)\gamma_\nu
\right)
\\
f_{22}&=\frac{1}{N}\sum_\text{spins}a_2a_2^*=
\mathop{\rm Tr}\left(
(\not p_1-m)\gamma^\mu(\not p_2+m)\gamma^\nu
\right)
\mathop{\rm Tr}\left(
(\not p_4+m)\gamma_\mu(\not p_3-m)\gamma_\nu
\right)
\end{align*}
Hence
\begin{equation*}
\langle|\mathcal{M}|^2\rangle
=
\frac{e^4}{4}
\left(
\frac{f_{11}}{t^2} - \frac{f_{12}}{st} -
\frac{f_{12}^*}{st} + \frac{f_{22}}{s^2}
\right)
\end{equation*}
These formulas compute probability densities directly from momentum vectors.
\begin{align*}
f_{11}&=
32(p_1\cdot p_2)(p_3\cdot p_4)+32(p_1\cdot p_4)(p_2\cdot p_3)
-32 m^2(p_1\cdot p_3)-32 m^2(p_2\cdot p_4)+64 m^4
\\
f_{12}&=
-32 (p_1\cdot p_4) (p_2\cdot p_3)
-16 m^2 (p_1\cdot p_2) + 16 m^2 (p_1\cdot p_3) - 16 m^2 (p_1\cdot p_4)\\
&\phantom{=}\qquad
{}- 16 m^2 (p_2\cdot p_3) + 16 m^2 (p_2\cdot p_4) - 16 m^2 (p_3\cdot p_4) - 32 m^4
\\
f_{22}&=
32(p_1\cdot p_3)(p_2\cdot p_4)+32(p_1\cdot p_4)(p_2\cdot p_3)
+32 m^2(p_1\cdot p_2)+32 m^2(p_3\cdot p_4)+64 m^4
\end{align*}
In Mandelstam variables $s=(p_1+p_2)^2$, $t=(p_1-p_3)^2$, $u=(p_1-p_4)^2$ the formulas are
\begin{align*}
f_{11} &= 8 s^2 + 8 u^2 - 64 s m^2 - 64 u m^2 + 192 m^4
\\
f_{12} &= -8 u^2 + 64 u m^2 - 96 m^4
\\
f_{22} &= 8 t^2 + 8 u^2 - 64 t m^2 - 64 u m^2 + 192 m^4
\end{align*}
When $E\gg m$ a useful approximation is to set $m=0$ and obtain
\begin{align*}
f_{11}&= 8 s^2 + 8 u^2\\
f_{12}&= -8 u^2\\
f_{22}&= 8 t^2 + 8 u^2
\end{align*}
For $m=0$ the Mandelstam variables are
\begin{align*}
s&=4E^2
\\
t&=-2E^2(1-\cos\theta)=-4E^2\sin^2(\theta/2)
\\
u&=-2E^2(1+\cos\theta)=-4E^2\cos^2(\theta/2)
\end{align*}
The corresponding expected probability density is
\begin{align*}
\langle|\mathcal{M}|^2\rangle
&=
\frac{e^4}{4}
\left(
\frac{8s^2+8u^2}{t^2}+\frac{16u^2}{st}+\frac{8t^2+8u^2}{s^2}
\right)
\\
&=2e^4\left(\frac{s^2+u^2}{t^2}+\frac{2u^2}{st}+\frac{t^2+u^2}{s^2}\right)
\\
&=2e^4
\left(
\frac{1+\cos^4(\theta/2)}{\sin^4(\theta/2)}
-\frac{2\cos^4(\theta/2)}{\sin^2(\theta/2)}
+\frac{1+\cos^2\theta}{2}
\right)
\end{align*}
