production of $n$ pions as $$N_n = \int_{s_{\rm min}}^{4E_0^2} ds\,P_n (s) \,
\frac{dN_{\rm in}^{\rm hard}}{ds}\,.
\label{number}$$ The lower limit of integration is fixed by kinematics and the upper limit is determined by the beam energy.

We shall approximate the pion number distribution function $P_n(s)$ with a binomial [@comment] and choose the parameters of this binomial such that we have some rough agreement with experiment [@topo]. $$P_n (s) = \left( \begin{array}{c} n_{\rm max}\\n \end{array} \right)
\xi^n
(1-\xi)^{n_{\rm max}-n}
\label{bino}$$ The maximum number of pions produced in a nucleon-nucleon collision is determined by kinematics. $$n_{\rm max} (s)= {\rm Integer} \left(\frac{\sqrt{s}-2m_N)}{m_{\pi}}\right)
\label{nc}$$ The parameter $\xi$ is related to the mean multiplicity by $$\xi(s) = \frac{\langle n \rangle}{n_{\rm max}}
= \frac{3}{n_{\rm max}}\,\left(\frac{1}{4} \langle n_{pp}^- \rangle
+\frac{1}{2} \langle n_{pn}^- \rangle
+\frac{1}{4} \langle n_{nn}^- \rangle\right)\, .
\label{q}$$ Here $\langle n \rangle$ is the average pion multiplicity averaged over $pp$, $pn$ and $nn$ collisions while $\langle n_{pp}^- \rangle$, $\langle n_{pn}^- \rangle$ and $\langle n_{nn}^- \rangle$ represent the average negative pion multiplicity in those collisions. All average multiplicities are functions of $s$, of course. The factor of 3 is due to isospin averaging.

Experimental data were compiled and parametrized in [@multi] as $$\begin{aligned}
\langle n_{pp}^- \rangle &=& -0.41 