The simple two-body processes

In the spirit of section , we want to select simple functions such that the true $\tau$, $y$ and $z$ dependence of the cross sections is approximately modelled. In particular, (almost) all conceivable kinematical peaks should be represented by separate terms in the approximate formulae. If this can be achieved, the ratio of the correct to the approximate cross sections will not fluctuate too much, but allow reasonable Monte Carlo efficiency.

Therefore the variables are generated according to the distributions $h_{\tau}(\tau)$, $h_y(y)$ and $h_z(z)$, where normally

$$h_{\tau}(\tau) = \frac{c_1}{{\cal I}_1} \, \frac{1}{\tau} + \frac{c_2}{{\cal I}_2}... ..._R)} + \frac{c_4}{{\cal I}_4} \, \frac{1}{(s\tau - m_R^2)^2 + m_R^2 \Gamma_R^2}$$

$$+ \frac{c_5}{{\cal I}_5} \, \frac{1}{\tau(\tau + \tau_{R'})} + \frac{c_6}{{\cal I}_6} \, \frac{1}{(s\tau - m_{R'}^2)^2 + m_{R'}^2 \Gamma_{R'}^2} ~,$$ (95)

$$h_y(y) = \frac{c_1}{{\cal I}_1} \, (y - y_{\mathrm{min}}) + \frac{c_2}{{\c... ...}_2} \, (y_{\mathrm{max}} - y) + \frac{c_3}{{\cal I}_3} \, \frac{1}{\cosh y} ~,$$ (96)

$$h_z(z) = \frac{c_1}{{\cal I}_1} + \frac{c_2}{{\cal I}_2} \, \frac{1}{a-z} ... ...cal I}_4} \, \frac{1}{(a-z)^2} + \frac{c_5}{{\cal I}_5} \, \frac{1}{(a+z)^2} ~.$$ (97)

Here each term is separately integrable, with an invertible primitive function, such that generation of $\tau$, $y$ and $z$ separately is a standard task, as described in section . In the following we describe the details of the scheme, including the meaning of the coefficients $c_i$ and ${\cal I}_i$, which are separate for $\tau$, $y$ and $z$.

The first variable to be selected is $\tau$. The range of allowed values, $\tau_{\mathrm{min}} \leq \tau \leq \tau_{\mathrm{max}}$, is generally constrained by a number of user-defined requirements. A cut on the allowed mass range is directly reflected in $\tau$, a cut on the $p_{\perp}$ range indirectly. The first two terms of $h_{\tau}$ are intended to represent a smooth $\tau$ dependence, as generally obtained in processes which do not receive contributions from $s$-channel resonances. Also $s$-channel exchange of essentially massless particles ($\gamma$, $\mathrm{g}$, light quarks and leptons) are accounted for, since these do not produce any separate peaks at non-vanishing $\tau$. The last four terms of $h_{\tau}$ are there to catch the peaks in the cross section from resonance production. These terms are only included when needed. Each resonance is represented by two pieces, a first to cover the interference with graphs which peak at $\tau =0$, plus the variation of parton distributions, and a second to approximate the Breit-Wigner shape of the resonance itself. The subscripts $R$ and $R'$ denote values pertaining to the two resonances, with $\tau_R = m_R^2/s$. Currently there is only one process where the full structure with two resonances is used, namely $\mathrm{f}\overline{\mathrm{f}}\to \gamma^*/\mathrm{Z}^0/\mathrm{Z}'^0$. Otherwise either one or no resonance peak is taken into account.

The kinematically allowed range of $y$ values is constrained by $\tau$, $\vert y\vert \leq - \frac{1}{2} \ln\tau$, and you may impose additional cuts. Therefore the allowed range $y_{\mathrm{min}} \leq y \leq y_{\mathrm{max}}$ is only constructed after $\tau$ has been selected. The first two terms of $h_y$ give a fairly flat $y$ dependence -- for processes which are symmetric in $y \leftrightarrow -y$, they will add to give a completely flat $y$ spectrum between the allowed limits. In principle, the natural subdivision would have been one term flat in $y$ and one forward-backward asymmetric, i.e. proportional to $y$. The latter is disallowed by the requirement of positivity, however. The $y - y_{\mathrm{min}}$ and $y_{\mathrm{max}} - y$ terms actually used give the same amount of freedom, but respect positivity. The third term is peaked at around $y = 0$, and represents the bias of parton distributions towards this region.

The allowed $z = \cos\hat{\theta}$ range is naïvely $-1 \leq z \leq 1$. However, most cross sections are divergent for $z \to \pm 1$, so some kind of regularization is necessary. Normally one requires $p_{\perp}\geq p_{\perp\mathrm{min}}$, which translates into $z^2 \leq 1 - 4 p_{\perp\mathrm{min}}^2/(\tau s)$ for massless outgoing particles. Since again the limits depend on $\tau$, the selection of $z$ is done after that of $\tau$. Additional requirements may constrain the range further. In particular, a $p_{\perp\mathrm{max}}$ constraint may split the allowed $z$ range into two, i.e. $z_{-\mathrm{min}} \leq z \leq z_{-\mathrm{max}}$ or $z_{+\mathrm{min}} \leq z \leq z_{+\mathrm{max}}$. An un-split range is represented by $z_{-\mathrm{max}} = z_{+\mathrm{min}} = 0$. For massless outgoing particles the parameter $a = 1$ in $h_z$, such that the five terms represent a piece flat in angle and pieces peaked as $1/\hat{t}$, $1/\hat{u}$, $1/\hat{t}^2$, and $1/\hat{u}^2$, respectively. For non-vanishing masses one has $a = 1 + 2 m_3^2 m_4^2/\hat{s}^2$. In this case, the full range $-1 \leq z \leq 1$ is therefore available -- physically, the standard $\hat{t}$ and $\hat{u}$ singularities are regularized by the masses $m_3$ and $m_4$.

For each of the terms, the ${\cal I}_i$ coefficients represent the integral over the quantity multiplying the coefficient $c_i$; thus, for instance:

$$h_{\tau}: {\cal I}_1 = \int \frac{\d\tau}{\tau} = \ln \left( \frac{\tau_{\mathrm{max}}}{\tau_{\mathrm{min}}} \right) ~,$$

$${\cal I}_2 = \int \frac{\d\tau}{\tau^2} = \frac{1}{\tau_{\mathrm{min}}} - \frac{1}{\tau_{\mathrm{max}}} ~;$$

$$h_y: {\cal I}_1 = \int (y - y_{\mathrm{min}}) \, \d y = \frac{1}{2} (y_{\mathrm{max}} - y_{\mathrm{min}})^2 ~;$$

$$h_z: {\cal I}_1 = \int \d z = (z_{-\mathrm{max}} - z_{-\mathrm{min}}) + (z_{+\mathrm{max}} - z_{+\mathrm{min}}) ,$$

$${\cal I}_2 = \int \frac{\d z}{a-z} = \ln \left( \frac{(a-z_{-\mat... ...)(a-z_{+\mathrm{min}})} {(a-z_{-\mathrm{max}})(a-z_{-\mathrm{min}})} \right) ~.$$ (98)

The $c_i$ coefficients are normalized to unit sum for $h_{\tau}$, $h_y$ and $h_z$ separately. They have a simple interpretation, as the probability for each of the terms to be used in the preliminary selection of $\tau$, $y$ and $z$, respectively. The variation of the cross section over the allowed phase space is explored in the initialization procedure of a PYTHIA run, and based on this knowledge the $c_i$ are optimized so as to give functions $h_{\tau}$, $h_y$ and $h_z$ that closely follow the general behaviour of the true cross section. For instance, the coefficient $c_4$ in $h_{\tau}$ is to be made larger the more the total cross section is dominated by the region around the resonance mass.

The phase-space points tested at initialization are put on a grid, with the number of points in each dimension given by the number of terms in the respective $h$ expression, and with the position of each point given by the median value of the distribution of one of the terms. For instance, the $\d\tau/\tau$ distribution gives a median point at $\sqrt{\tau_{\mathrm{min}}\tau_{\mathrm{max}}}$, and $\d\tau / \tau^2$ has the median $2 \tau_{\mathrm{min}} \tau_{\mathrm{max}} / (\tau_{\mathrm{min}} + \tau_{\mathrm{max}})$. Since the allowed $y$ and $z$ ranges depend on the $\tau$ value selected, then so do the median points defined for these two variables.

With only a limited set of phase-space points studied at the initialization, the `optimal' set of coefficients is not uniquely defined. To be on the safe side, 40% of the total weight is therefore assigned evenly between all allowed $c_i$, whereas the remaining 60% are assigned according to the relative importance surmised, under the constraint that no coefficient is allowed to receive a negative contribution from this second piece.

After a preliminary choice has been made of $\tau$, $y$ and $z$, it is necessary to find the weight of the event, which is to be used to determine whether to keep it or generate another one. Using the relation $\d\hat{t} = \hat{s} \, \beta_{34} \, \d z / 2$, eq. () may be rewritten as

$$\sigma = \int \int \int \frac{\d\tau}{\tau} \, \d y \, \frac{\hat{s} \beta... ...z \, x_1 f_1(x_1, Q^2) \, x_2 f_2(x_2, Q^2) \, \frac{\d\hat{\sigma}}{\d\hat{t}}$$

$$= \frac{\pi}{s} \int h_{\tau}(\tau) \, \d\tau \int h_y(y) \, \d y \... ..., h_y(y) \, 2 h_z(z)} \, \frac{\hat{s}^2}{\pi} \frac{\d\hat{\sigma}}{\d\hat{t}}$$

$$= \left\langle \frac{\pi}{s} \, \frac{\beta_{34}}{\tau^2 h_{\tau}(\... ...Q^2) \, \frac{\hat{s}^2}{\pi} \frac{\d\hat{\sigma}}{\d\hat{t}} \right\rangle ~.$$ (99)

In the middle line, a factor of $1 = h_{\tau}/h_{\tau}$ has been introduced to rewrite the $\tau$ integral in terms of a phase space of unit volume: $\int h_{\tau}(\tau) \, \d\tau = 1$ according to the relations above. Correspondingly for the $y$ and $z$ integrals. In addition, factors of $1 = \hat{s}/ (\tau s)$ and $1 = \pi / \pi$ are used to isolate the dimensionless cross section $(\hat{s}^2/\pi) \, \d\hat{\sigma} / \d\hat{t}$. The content of the last line is that, with $\tau$, $y$ and $z$ selected according to the expressions $h_{\tau}(\tau)$, $h_y(y)$ and $h_z(z)$, respectively, the cross section is obtained as the average of the final expression over all events. Since the $h$'s have been picked to give unit volume, there is no need to multiply by the total phase-space volume.

As can be seen, the cross section for a given Monte Carlo event is given as the product of four factors, as follows:

57.

The $\pi/s$ factor, which is common to all events, gives the overall dimensions of the cross section, in GeV$^{-2}$. Since the final cross section is given in units of mb, the conversion factor of 1 GeV $^{-2} = 0.3894$ mb is also included here.

58.

Next comes the Jacobian, which compensates for the change from the original to the final phase-space volume.

59.

The parton-distribution function weight is obtained by making use of the parton distribution libraries in PYTHIA or externally. The $x_1$ and $x_2$ values are obtained from $\tau$ and $y$ via the relations $x_{1,2} = \sqrt{\tau} \exp(\pm y)$.

60.

Finally, the dimensionless cross section $(\hat{s}^2/\pi) \, \d\hat{\sigma} / \d\hat{t}$ is the quantity that has to be coded for each process separately, and where the physics content is found.

Of course, the expression in the last line is not strictly necessary to obtain the cross section by Monte Carlo integration. One could also have used eq. () directly, selecting phase-space points evenly in $\tau$, $y$ and $\hat{t}$, and averaging over those Monte Carlo weights. Clearly this would be much simpler, but the price to be paid is that the weights of individual events could fluctuate wildly. For instance, if the cross section contains a narrow resonance, the few phase-space points that are generated in the resonance region obtain large weights, while the rest do not. With our procedure, a resonance would be included in the $h_{\tau}(\tau)$ factor, so that more events would be generated at around the appropriate $\tau_R$ value (owing to the $h_{\tau}$ numerator in the phase-space expression), but with these events assigned a lower, more normal weight (owing to the factor $1/h_{\tau}$ in the weight expression). Since the weights fluctuate less, fewer phase-space points need be selected to get a reasonable cross-section estimate.

In the program, the cross section is obtained as the average over all phase-space points generated. The events actually handed on to you should have unit weight, however (an option with weighted events exists, but does not represent the mainstream usage). At initialization, after the $c_i$ coefficients have been determined, a search inside the allowed phase-space volume is therefore made to find the maximum of the weight expression in the last line of eq. (). In the subsequent generation of events, a selected phase-space point is then retained with a probability equal to the weight in the point divided by the maximum weight. Only the retained phase-space points are considered further, and generated as complete events.

The search for the maximum is begun by evaluating the weight in the same grid of points as used to determine the $c_i$ coefficients. The point with highest weight is used as starting point for a search towards the maximum. In unfortunate cases, the convergence could be towards a local maximum which is not the global one. To somewhat reduce this risk, also the grid point with second-highest weight is used for another search. After initialization, when events are generated, a warning message will be given by default at any time a phase-space point is selected where the weight is larger than the maximum, and thereafter the maximum weight is adjusted to reflect the new knowledge. This means that events generated before this time have a somewhat erroneous distribution in phase space, but if the maximum violation is rather modest the effects should be negligible. The estimation of the cross section is not affected by any of these considerations, since the maximum weight does not enter into eq. ().

For $2 \to 2$ processes with identical final-state particles, the symmetrization factor of $1/2$ is explicitly included at the end of the $\d\hat{\sigma}/\d\hat{t}$ calculation. In the final cross section, a factor of 2 is retrieved because of integration over the full phase space (rather than only half of it). That way, no special provisions are needed in the phase-space integration machinery.