The Physics Routines

Once the initial parton/particle configuration has been specified and default parameter values changed, if so desired, only a PYEXEC call is necessary to simulate the whole fragmentation and decay chain. Therefore a normal user will not directly see any of the other routines in this section. Some of them could be called directly, but the danger of faulty usage is then non-negligible.

The PYTAUD routine provides an optional interface to an external $\tau$ decay library, where polarization effects could be included. It is up to you to write the appropriate calls, as explained at the end of this section.

Purpose:

to administrate the fragmentation and decay chain. PYEXEC may be called several times, but only entries which have not yet been treated (more precisely, which have $1 \leq$K(I,1)$\leq 10$) can be affected by further calls. This may apply if more partons/particles have been added by you, or if particles previously considered stable are now allowed to decay. The actions that will be taken during a PYEXEC call can be tailored extensively via the PYDAT1-PYDAT3 common blocks, in particular by setting the MSTJ values suitably.

SUBROUTINE PYPREP(IP) :

to rearrange parton shower end products (marked with K(I,1)=3) sequentially along strings; also to (optionally) allow small parton systems to collapse into two particles or one only, in the latter case with energy and momentum to be shuffled elsewhere in the event; also to perform checks that e.g. flavours of colour-singlet systems make sense.

SUBROUTINE PYSTRF(IP) :

to generate the fragmentation of an arbitrary colour-singlet parton system according to the Lund string fragmentation model. One of the absolutely central routines of PYTHIA.

SUBROUTINE PYINDF(IP) :

to handle the fragmentation of a parton system according to independent fragmentation models, and implement energy, momentum and flavour conservation, if so desired. Also the fragmentation of a single parton, not belonging to a parton system, is considered here (this is of course physical nonsense, but may sometimes be convenient for specific tasks).

SUBROUTINE PYDECY(IP) :

to perform a particle decay, according to known branching ratios or different kinds of models, depending on our level of knowledge. Various matrix elements are included for specific processes.

SUBROUTINE PYKFDI(KFL1,KFL2,KFL3,KF) :

to generate a new quark or diquark flavour and to combine it with an existing flavour to give a hadron.

KFL1:

incoming flavour.

KFL2:

extra incoming flavour, e.g. for formation of final particle, where the flavours are completely specified. Is normally 0.

KFL3:

newly created flavour; is 0 if KFL2 is non-zero.

KF:

produced hadron. Is 0 if something went wrong (e.g. inconsistent combination of incoming flavours).

SUBROUTINE PYPTDI(KFL,PX,PY) :

to give transverse momentum, e.g. for a $\mathrm{q}\overline{\mathrm{q}}$ pair created in the colour field, according to independent Gaussian distributions in $p_x$ and $p_y$.

SUBROUTINE PYZDIS(KFL1,KFL3,PR,Z) :

to generate the longitudinal scaling variable $z$ in jet fragmentation, either according to the Lund symmetric fragmentation function, or according to a choice of other shapes.

SUBROUTINE PYBOEI :

to include Bose-Einstein effects according to a simple parameterization. By default, this routine is not called. If called from PYEXEC, this is done after the decay of short-lived resonances, but before the decay of long-lived ones. This means the routine should never be called directly by you, nor would effects be correctly simulated if decays are switched off. See MSTJ(51) - MSTJ(57) for switches of the routine.

FUNCTION PYMASS(KF) :

to give the mass for a parton/particle.

SUBROUTINE PYNAME(KF,CHAU) :

to give the parton/particle name (as a string of type CHARACTER CHAU*16). The name is read out from the CHAF array.

FUNCTION PYCHGE(KF) :

to give three times the charge for a parton/particle. The value is read out from the KCHG(KC,1) array.

FUNCTION PYCOMP(KF) :

to give the compressed parton/particle code KC for a given KF code, as required to find entry into mass and decay data tables. Also checks whether the given KF code is actually an allowed one (i.e. known by the program), and returns 0 if not. Note that KF may be positive or negative, while the resulting KC code is never negative.
Internally PYCOMP uses a binary search in a table, with KF codes arranged in increasing order, based on the KCHG(KC,4) array. This table is constructed the first time PYCOMP is called, at which time MSTU(20) is set to 1. In case of a user change of the KCHG(KC,4) array one should reset MSTU(20)=0 to force a re-initialization at the next PYCOMP call (this is automatically done in PYUPDA calls). To speed up execution, the latest (KF,KC) pair is kept in memory and checked before the standard binary search.

SUBROUTINE PYERRM(MERR,MESSAG) :

to keep track of the number of errors and warnings encountered, write out information on them, and abort the program in case of too many errors.

FUNCTION PYANGL(X,Y) :

to calculate the angle from the $x$ and $y$ coordinates.

SUBROUTINE PYLOGO :

to write a title page for the PYTHIA programs. Called by PYLIST(0).

SUBROUTINE PYTIME(IDATI) :

to give the date and time, for use in PYLOGO and elsewhere. Since Fortran 77 does not contain a standard way of obtaining this information, the routine is dummy, to be replaced by you. Some commented-out examples are given, e.g. for Fortran 90 or the GNU Linux libU77. The output is given in an integer array ITIME(6), with components year, month, day, hour, minute and second. If there should be no such information available on a system, it is acceptable to put all the numbers above to 0.

Purpose:

to act as an interface between the standard decay routine PYDECY and a user-supplied $\tau$ lepton decay library such as TAUOLA [Jad91]. The latter library would normally know how to handle polarized $\tau$'s, given the $\tau$ helicity as input, so one task of the interface routine is to construct the $\tau$ polarization/helicity from the information available. Input to the routine (from PYDECY) is provided in the first three arguments, while the last argument and some event record information have to be set before return. To use this facility you have to set the switch MSTJ(28), include your own interface routine PYTAUD and see to it that the dummy routine PYTAUD is not linked. The dummy routine is there only to avoid unresolved external references when no user-supplied interface is linked.

ITAU :

line number in the event record where the $\tau$ is stored. The four-momentum of this $\tau$ has first been boosted back to the rest frame of the decaying mother and thereafter rotated to move out along the $+z$ axis. It would have been possible to also perform a final boost to the rest frame of the $\tau$ itself, but this has been avoided so as not to suppress the kinematics aspect of close-to-threshold production (e.g. in $\mathrm{B}$ decays) vs. high-energy production (e.g. in real $\mathrm{W}$ decays). The choice of frame should help the calculation of the helicity configuration. After the PYTAUD call the $\tau$ and its decay products will automatically be rotated and boosted back. However, seemingly, the event record does not conserve momentum at this intermediate stage.

IORIG :

line number where the mother particle to the $\tau$ is stored. Is 0 if the mother is not stored. This does not have to mean the mother is unknown. For instance, in semileptonic $\mathrm{B}$ decays the mother is a $\mathrm{W}^{\pm}$ with known four-momentum $p_{\mathrm{W}} = p_{\tau} + p_{\nu_{\tau}}$, but there is no $\mathrm{W}$ line in the event record. When several copies of the mother is stored (e.g. one in the documentation section of the event record and one in the main section), IORIG points to the last. If a branchings like $\tau \to \tau\gamma$ occurs, the `grandmother' is given, i.e. the mother of the direct $\tau$ before branching.

KFORIG :

flavour code for the mother particle. Is 0 if the mother is unknown. The mother would typically be a resonance such as $\gamma^* / \mathrm{Z}^0$ (23), $\mathrm{W}^{\pm}$ ($\pm 24$), $\mathrm{h}^0$ (25), or $\H ^{\pm}$ ($\pm 37$). Often the helicity choice would be clear just by the knowledge of this mother species, e.g., $\mathrm{W}^{\pm}$ vs. $\H ^{\pm}$. However, sometimes further complications may exist. For instance, the KF code 23 represents a mixture of $\gamma^*$ and $\mathrm{Z}^0$; a knowledge of the mother mass (in P(IORIG,5)) would here be required to make the choice of helicities. Further, a $\mathrm{W}^{\pm}$ or $\mathrm{Z}^0$ may either be (predominantly) transverse or longitudinal, depending on the production process under study.

NDECAY :

the number of decay products of the $\tau$; to be given by the user routine. You must also store the KF flavour codes of those decay products in the positions K(I,2), N+1$\leq$I$\leq$N+NDECAY, of the event record. The corresponding five-momentum (momentum, energy and mass) should be stored in the associated P(I,J) positions, 1$\leq$J$\leq$5. The four-momenta are expected to add up to the four-momentum of the $\tau$ in position ITAU. You should not change the N value or any of the other K or V values (neither for the $\tau$ nor for its decay products) since this is automatically done in PYDECY.