Solution

First, let us go to the rest frame of the $\pi^0$ mesons. In this frame, denote the angles to the velocity $\beta c$ as $\theta$. Obviously, for the conservation of momentum, the photons of each pair fly apart strictly in opposite directions, but overall isotropically – all $\theta$ for the line along which the photons fly apart are equally probable.

Now let us switch to the laboratory frame. We need to find how the angles change. This can be done using the ready‑made formula for light aberration, but I will show the derivation using the Lorentz transformations for the 4‑momentum. Taking into account that for a photon $E=pc$, it is not necessary to write the full transformations; it is sufficient to write
\begin{pmatrix} p' \\ p_x' \end{pmatrix}= \begin{pmatrix} p \\ p_x \end{pmatrix}\cdot \begin{pmatrix} \gamma & \beta\gamma \\ \beta\gamma & \gamma \end{pmatrix}

$$p' \begin{pmatrix} 1\\ \cos\theta' \end{pmatrix}=\gamma p \begin{pmatrix} 1 \\ \cos\theta \end{pmatrix}\cdot \begin{pmatrix} 1 & \beta \\ \beta& 1 \end{pmatrix}=\gamma p \begin{pmatrix} 1+\beta\cos\theta \\ \beta+\cos\theta \end{pmatrix}$$
$$\cos\theta'=\frac{p_x'}{p'}=\frac{\beta+\cos\theta}{1+\beta\cos\theta }$$

Working with the full matrix gives the same result.

The condition $\theta'<\frac{\pi}{2}$ is equivalent to $\cos\theta'>0$.
$$0<\beta<1;\ \cos\theta\ge -1\quad \to \quad 1+\beta\cos\theta>0 \text{ – always.}$$
Then the condition can be rewritten as
$$\cos\theta>-\beta.$$
The probability, taking into account the isotropy of $\theta$ , is then found as follows ($k$ is a normalisation factor chosen so that the total probability equals 1):
$$N=\frac{\int_{-\beta}^1 k\, d(\cos\theta)}{\int_{-1}^1 k\, d(\cos\theta)}=\frac{1+\beta}{2}.$$