next up previous contents
Next: Interfaces to Other Generators Up: How to Include External Previous: PYTHIA as a generator   Contents


Further comments

This section contains additional information on a few different topics: cross section interpretation, negative-weight events, relations with other PYTHIA switches and routines, and error conditions.

In several IDWTUP options, the XWGTUP variable is supposed to give the differential cross section of the current event, times the phase-space volume within which events are generated, expressed in picobarns. (Converted to millibarns inside PYTHIA.) This means that, in the limit that many events are generated, the average value of XWGTUP gives the total cross section of the simulated process. Of course, the tricky part is that the differential cross section usually is strongly peaked in a few regions of the phase space, such that the average probability to accept an event, $\langle$XWGTUP$\rangle /$XMAXUP(i) is small. It may then be necessary to find a suitable set of transformed phase-space coordinates, for which the correspondingly transformed differential cross section is better behaved.

To avoid confusion, here is a more formal version of the above paragraph. Call $\d X$ the differential phase space, e.g. for a $2 \to 2$ process $\d X = \d x_1 \, \d x_2 \, \d\hat{t}$, where $x_1$ and $x_2$ are the momentum fractions carried by the two incoming partons and $\hat{t}$ the Mandelstam variable of the scattering (see section [*]). Call $\d\sigma / \d X$ the differential cross section of the process, e.g. for $2 \to 2$: $\d\sigma / \d X = \sum_{ij}
f_i(x_1,Q^2) \, f_j(x_2,Q^2) \, \d\hat{\sigma}_{ij} / \d\hat{t}$, i.e. the product of parton distributions and hard-scattering matrix elements, summed over all allowed incoming flavours $i$ and $j$. The physical cross section that one then wants to generate is $\sigma = \int (\d\sigma / \d X) \, \d X$, where the integral is over the allowed phase-space volume. The event generation procedure consists of selecting an $X$ uniformly in $\d X$ and then evaluating the weight $\d\sigma / \d X$ at this point. XWGTUP is now simply XWGTUP $ = \d\sigma / \d X \, \int \d X$, i.e. the differential cross section times the considered volume of phase space. Clearly, when averaged over many events, XWGTUP will correctly estimate the desired cross section. If XWGTUP fluctuates too much, one may try to transform to new variables $X'$, where events are now picked accordingly to $\d X'$ and XWGTUP $ = \d\sigma / \d X' \, \int \d X'$.

A warning. It is important that $X$ is indeed uniformly picked within the allowed phase space, alternatively that any Jacobians are properly taken into account. For instance, in the case above, one approach would be to pick $x_1$, $x_2$ and $\hat{t}$ uniformly in the ranges $0 < x_1 < 1$, $0 < x_2 < 1$, and $-s < \hat{t} < 0$, with full phase space volume $\int \d X = s$. The cross section would only be non-vanishing inside the physical region given by $-s x_1 x_2 < \hat{t}$ (in the massless case), i.e. Monte Carlo efficiency is likely to be low. However, if one were to choose $\hat{t}$ values only in the range $-\hat{s} < \hat{t} < 0$, small $\hat{s}$ values would be favoured, since the density of selected $\hat{t}$ values would be larger there. Without the use of a compensating Jacobian $\hat{s}/s$, an incorrect answer would be obtained. Alternatively, one could start out with a phase space like $\d X = \d x_1 \, \d x_2 \, \d (\cos\hat{\theta})$, where the limits decouple. Of course, the $\cos\hat{\theta}$ variable can be translated back into a $\hat{t}$, which will then always be in the desired range $-\hat{s} < \hat{t} < 0$. The transformation itself here gives the necessary Jacobian.

At times, it is convenient to split a process into a discrete set of subprocesses for the parton-level generation, without retaining these in the IDPRUP classification. For instance, the cross section above contains a summation over incoming partons. An alternative would then have been to let each subprocess correspond to one unique combination of incoming flavours. When an event of process type $i$ is to be generated, first a specific subprocess $ik$ is selected with probability $f^{ik}$, where $\sum_k f^{ik} = 1$. For this subprocess an XWGTUP$^k$ is generated as above, except that there is no longer a summation over incoming flavours. Since only a fraction $f^{ik}$ of all events now contain this part of the cross section, a compensating factor $1/f^{ik}$ is introduced, i.e. XWGTUP$ = $XWGTUP$^k/f^{ik}$. Further, one has to define XMAXUP(i) $ = \mathrm{max}_k$ XMAXUP$^{ik}/f^{ik}$ and XSECUP(i)$ = \sum_k$ XSECUP$^{ik}$. The generation efficiency will be maximized for the $f^{ik}$ coefficients selected proportional to XMAXUP$^{ik}$, but this is no requirement.

The standard allows external parton-level events to come with negative weights, unlike the case for internal PYTHIA processes. In order to avoid indiscriminate use of this option, some words of caution are in place. In next-to-leading-order calculations, events with negative weights may occur as part of the virtual corrections. In any physical observable quantity, the effect of such events should cancel against the effect of real events with one more parton in the final state. For instance, the next-to-leading order calculation of gluon scattering contains the real process $\mathrm{g}\mathrm{g}\to \mathrm{g}\mathrm{g}\mathrm{g}$, with a positive divergence in the soft and collinear regions, while the virtual corrections to $\mathrm{g}\mathrm{g}\to \mathrm{g}\mathrm{g}$ are negatively divergent. Neglecting the problems of outgoing gluons collinear with the beams, and those of soft gluons, two nearby outgoing gluons in the $\mathrm{g}\mathrm{g}\to \mathrm{g}\mathrm{g}\mathrm{g}$ process can be combined into one effective one, such that the divergences can be cancelled.

If rather widely separated gluons can be combined, the remaining negative contributions are not particularly large. Different separation criteria could be considered; one example would be $\Delta R = \sqrt{(\Delta \eta)^2 + (\Delta \varphi)^2} \approx 1$. The recombination of well separated partons is at the price of an arbitrariness in the choice of clustering algorithm, when two gluons of nonvanishing invariant mass are to be combined into one massless one, as required to be able to compare with the kinematics of the massless $\mathrm{g}\mathrm{g}\to \mathrm{g}\mathrm{g}$ process when performing the divergence cancellations. Alternatively, if a smaller $\Delta R$ cut is used, where the combining procedure is less critical, there will be more events with large positive and negative weights that are to cancel.

Without proper care, this cancellation could easily be destroyed by the subsequent showering description, as follows. The standard for external processes does not provide any way to pass information on the clustering algorithm used, so the showering routine will have to make its own choice what region of phase space to populate with radiation. One choice could be to allow a cone defined by the nearest colour-connected parton (see section [*] for a discussion). There could then arise a significant mismatch in shower description between events where two gluons are just below or just above the $\Delta R$ cut for being recombined, equivalently between $\mathrm{g}\mathrm{g}\to \mathrm{g}\mathrm{g}$ and $\mathrm{g}\mathrm{g}\to \mathrm{g}\mathrm{g}\mathrm{g}$ events. Most of the phase space may be open for the former, while only the region below $\Delta R$ may be it for the latter. Thus the average `two-parton' events may end up containing significantly more jet activity than the corresponding `three-parton' ones. The smaller the $\Delta R$ cut, the more severe the mismatch, both on an event-by-event basis and in terms of the event rates involved.

One solution would be to extend the standard also to specify which clustering algorithm has been used in the matrix-element calculation, and with what parameter values. Any shower emission that would give rise to an extra jet, according to this algorithm, would be vetoed. If a small $\Delta R$ cut is used, this is almost equivalent to allowing no shower activity at all. (That would still leave us with potential mismatch problems in the hadronization description. Fortunately the string fragmentation approach is very powerful in this respect, with a smooth transition between two almost parallel gluons and a single one with the full energy [Sjö84].) But we know that the unassisted matrix-element description cannot do a good job of the internal structure of jets on an event-by-event basis, since multiple-gluon emission is the norm in this region. Therefore a $\Delta R \sim 1$ will be required, to let the matrix elements describe the wide-angle emission and the showers the small-angle one. This again suggests a picture with only a small contribution from negative-weight events. In summary, the appearance of a large fraction of negative-weight events should be a sure warning sign that physics is likely to be incorrectly described.

The above example illustrates that it may, at times, be desirable to sidestep the standard and provide further information directly in the PYTHIA common blocks. (Currently there is no exact match to the clustering task mentioned above, although the already-described UPVETO routine, section [*], could be constructed to make use of such information. Here we concentrate on a few simpler ways to intervene, possibly to be used in conjunction with UPVETO.) Then it is useful to note that, apart from the hard-process generation machinery itself, the external processes are handled almost exactly as the internal ones. Thus essentially all switches and parameter values related to showers, underlying events and hadronization can be modified at will. This even applies to alternative listing modes of events and history pointers, as set by MSTP(128). Also some of the information on the hard scattering is available, such as MSTI(3), MSTI(21) - MSTI(26), and PARI(33) - PARI(38). Before using them, however, it is prudent to check that your variables of interest do work as intended for the particular process you study. Several differences do remain between internal and external processes, in particular related to final-state showers and resonance decays. For internal processes, the PYRESD routine will perform a shower (if relevant) directly after each decay. A typical example would be that a $\t\to \b\mathrm{W}$ decay is immediately followed by a shower, which could change the momentum of the $\mathrm{W}$ before it decays in its turn. For an external process, this decay chain would presumably already have been carried out. When the equivalent shower to the above is performed, it is therefore now necessary also to boost the decay products of the $\mathrm{W}$. The special sequence of showers and boosts for external processes is administrated by the PYADSH routine. Should the decay chain not have been carried out, e.g if HEPEUP event record contains an undecayed $\mathrm{Z}^0$, then PYRESD will be called to let it decay. The decay products will be visible also in the documentation section, as for internal processes.

You are free to make use of whatever tools you want in your UPINIT and UPEVNT routines, and normally there would be little or no contact with the rest of PYTHIA, except as described above. However, several PYTHIA tools can be used, if you so wish. One attractive possibility is to use PYPDFU for parton-distribution-function evaluation. Other possible tools could be PYR for random-number generation, PYALPS for $\alpha_{\mathrm{s}}$ evaluation, PYALEM for evaluation of a running $\alpha_{\mathrm{em}}$, and maybe a few more.

We end with a few comments on anomalous situations. As already described, you may put NUP = 0 inside UPEVNT, e.g. to signal the end of the file from which events are read. If the program encounters this value at a return from UPEVNT, then it will also exit from PYEVNT, without incrementing the counters for the number of events generated. It is then up to you to have a check on this condition in your main event-generation loop. This you do either by looking at NUP or at MSTI(51); the latter is set to 1 if no event was generated.

It may also happen that a parton-level configuration fails elsewhere in the PYEVNT call. For instance, the beam-remnant treatment occasionally encounters situations it cannot handle, wherefore the parton-level event is rejected and a new one generated. This happens also with ordinary (not user-defined) events, and usually comes about as a consequence of the initial-state radiation description leaving too little energy for the remnant. If the same hard scattering were to be used as input for a new initial-state radiation and beam-remnant attempt, it could then work fine. There is a possibility to give events that chance, as follows. MSTI(52) counts the number of times a hard-scattering configuration has failed to date. If you come in to UPEVNT with MSTI(52) non-vanishing, this means that the latest configuration failed. So long as the contents of the HEPEUP common block are not changed, such an event may be given another try. For instance, a line

      IF(MSTI(52).GE.1.AND.MSTI(52).LE.4) RETURN
at the beginning of UPEVNT will give each event up to five tries; thereafter a new one would be generated as usual. Note that the counter for the number of events is updated at each new try. The fraction of failed configurations is given in the bottom line of the PYSTAT(1) table.

The above comment only refers to very rare occurrences (less than one in a hundred), which are not errors in a strict sense; for instance, they do not produce any error messages on output. If you get warnings and error messages that the program does not understand the flavour codes or cannot reconstruct the colour flows, it is due to faults of yours, and giving such events more tries is not going to help.


next up previous contents
Next: Interfaces to Other Generators Up: How to Include External Previous: PYTHIA as a generator   Contents
Stephen Mrenna 2007-10-30