Getting Started with the Simple Routines

Normally PYTHIA is expected to take care of the full event generation process. At times, however, one may want to access the more simple underlying routines, which allow a large flexibility to `do it yourself'. We therefore start with a few cases of this kind, at the same time introducing some of the more frequently used utility routines.

As a first example, assume that you want to study the production of $\u\overline{\mathrm{u}}$ 2-jet systems at 20 GeV energy. To do this, write a main program

      IMPLICIT DOUBLE PRECISION(A-H, O-Z)
      CALL PY2ENT(0,2,-2,20D0)
      CALL PYLIST(1)
      END

and run this program, linked together with PYTHIA. The routine PY2ENT is specifically intended for storing two entries (partons or particles). The first argument (0) is a command to perform fragmentation and decay directly after the entries have been stored, the second and third that the two entries are $\u$ (2) and $\overline{\mathrm{u}}$ ($-2$), and the last that the c.m. energy of the pair is 20 GeV, in double precision. When this is run, the resulting event is stored in the PYJETS common block. This information can then be read out by you. No output is produced by PY2ENT itself, except for a title page which appears once for every PYTHIA run.

Instead the second command, to PYLIST, provides a simple visible summary of the information stored in PYJETS. The argument (1) indicates that the short version should be used, which is suitable for viewing the listing directly on an 80-column terminal screen. It might look as shown here.

                        Event listing (summary)
 
  I  particle/jet KS     KF orig   p_x     p_y     p_z      E       m
 
  1  (u)       A  12      2    0   0.000   0.000  10.000  10.000   0.006
  2  (ubar)    V  11     -2    0   0.000   0.000 -10.000  10.000   0.006
  3  (string)     11     92    1   0.000   0.000   0.000  20.000  20.000
  4  (rho+)       11    213    3   0.098  -0.154   2.710   2.856   0.885
  5  (rho-)       11   -213    3  -0.227   0.145   6.538   6.590   0.781
  6  pi+           1    211    3   0.125  -0.266   0.097   0.339   0.140
  7  (Sigma0)     11   3212    3  -0.254   0.034  -1.397   1.855   1.193
  8  (K*+)        11    323    3  -0.124   0.709  -2.753   2.968   0.846
  9  p~-           1  -2212    3   0.395  -0.614  -3.806   3.988   0.938
 10  pi-           1   -211    3  -0.013   0.146  -1.389   1.403   0.140
 11  pi+           1    211    4   0.109  -0.456   2.164   2.218   0.140
 12  (pi0)        11    111    4  -0.011   0.301   0.546   0.638   0.135
 13  pi-           1   -211    5   0.089   0.343   2.089   2.124   0.140
 14  (pi0)        11    111    5  -0.316  -0.197   4.449   4.467   0.135
 15  (Lambda0)    11   3122    7  -0.208   0.014  -1.403   1.804   1.116
 16  gamma         1     22    7  -0.046   0.020   0.006   0.050   0.000
 17  K+            1    321    8  -0.084   0.299  -2.139   2.217   0.494
 18  (pi0)        11    111    8  -0.040   0.410  -0.614   0.751   0.135
 19  gamma         1     22   12   0.059   0.146   0.224   0.274   0.000
 20  gamma         1     22   12  -0.070   0.155   0.322   0.364   0.000
 21  gamma         1     22   14  -0.322  -0.162   4.027   4.043   0.000
 22  gamma         1     22   14   0.006  -0.035   0.422   0.423   0.000
 23  p+            1   2212   15  -0.178   0.033  -1.343   1.649   0.938
 24  pi-           1   -211   15  -0.030  -0.018  -0.059   0.156   0.140
 25  gamma         1     22   18  -0.006   0.384  -0.585   0.699   0.000
 26  gamma         1     22   18  -0.034   0.026  -0.029   0.052   0.000
                 sum:  0.00        0.000   0.000   0.000  20.000  20.000

(A few blanks have been removed between the columns to make it fit into the format of this text.) Look in the particle/jet column and note that the first two lines are the original $\u$ and $\overline{\mathrm{u}}$. The parentheses enclosing the names, `(u)' and `(ubar)', are there as a reminder that these partons actually have been allowed to fragment. The partons are still retained so that event histories can be studied. Also note that the KF (flavour code) column contains 2 in the first line and $-2$ in the second. These are the codes actually stored to denote the presence of a $\u$ and a $\overline{\mathrm{u}}$, cf. the PY2ENT call, while the names written are just conveniences used when producing visible output. The A and V near the end of the particle/jet column indicate the beginning and end of a string (or cluster, or independent fragmentation) parton system; any intermediate entries belonging to the same system would have had an I in that column. (This gives a poor man's representation of an up-down arrow, $\updownarrow$.)

In the orig (origin) column, the zeros indicate that $\u$ and $\overline{\mathrm{u}}$ are two initial entries. The subsequent line, number 3, denotes the fragmenting $\u\overline{\mathrm{u}}$ string system as a whole, and has origin 1, since the first parton of this string system is entry number 1. The particles in lines 4-10 have origin 3 to denote that they come directly from the fragmentation of this string. In string fragmentation it is not meaningful to say that a particle comes from only the $\u$ quark or only the $\overline{\mathrm{u}}$ one. It is the string system as a whole that gives a $\rho^+$, a $\rho^-$, a $\pi^+$, a $\Sigma^0$, a $\mathrm{K}^{*+}$, a $\overline{\mathrm{p}}^-$, and a $\pi^-$. Note that some of the particle names are again enclosed in parentheses, indicating that these particles are not present in the final state either, but have decayed further. Thus the $\pi^-$ in line 13 and the $\pi^0$ in line 14 have origin 5, as an indication that they come from the decay of the $\rho^-$ in line 5. Only the names not enclosed in parentheses remain at the end of the fragmentation/decay chain, and are thus experimentally observable. The actual status code used to distinguish between different classes of entries is given in the KS column; codes in the range 1-10 correspond to remaining entries, and those above 10 to those that have fragmented or decayed.

The columns with p_x, p_y, p_z, E and m are quite self-explanatory. All momenta, energies and masses are given in units of GeV, since the speed of light is taken to be $c = 1$. Note that energy and momentum are conserved at each step of the fragmentation/decay process (although there exist options where this is not true). Also note that the $z$ axis plays the rôle of preferred direction, along which the original partons are placed. The final line is intended as a quick check that nothing funny happened. It contains the summed charge, summed momentum, summed energy and invariant mass of the final entries at the end of the fragmentation/decay chain, and the values should agree with the input implied by the PY2ENT arguments. (In fact, warnings would normally appear on the output if anything untoward happened, but that is another story.)

The above example has illustrated roughly what information is to be had in the event record, but not so much about how it is stored. This is better seen by using a 132-column format for listing events. Try e.g. the following program

      IMPLICIT DOUBLE PRECISION(A-H, O-Z)
      CALL PY3ENT(0,1,21,-1,30D0,0.9D0,0.7D0)
      CALL PYLIST(2)
      CALL PYEDIT(3)
      CALL PYLIST(2)
      END

where a 3-jet $\d\mathrm{g}\overline{\mathrm{d}}$ event is generated in the first executable line and listed in the second. This listing will contain the numbers as directly stored in the common block PYJETS

      COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)

For particle I, K(I,1) thus gives information on whether or not a parton or particle has fragmented or decayed, K(I,2) gives the particle code, K(I,3) its origin, K(I,4) and K(I,5) the position of fragmentation/decay products, and P(I,1)-P(I,5) momentum, energy and mass. The number of lines in current use is given by N, i.e. 1 $\leq$ I $\leq$ N. The V matrix contains decay vertices; to view those PYLIST(3) has to be used. NPAD is a dummy, needed to avoid some compiler troubles. It is important to learn the rules for how information is stored in PYJETS.

The third executable line in the program illustrates another important point about PYTHIA: a number of routines are available for manipulating and analyzing the event record after the event has been generated. Thus PYEDIT(3) will remove everything except stable charged particles, as shown by the result of the second PYLIST call. More advanced possibilities include things like sphericity or clustering routines. PYTHIA also contains some simple routines for histogramming, used to give self-contained examples of analysis procedures.

Apart from the input arguments of subroutine calls, control on the doings of PYTHIA may be imposed via many common blocks. Here sensible default values are always provided. A user might want to switch off all particle decays by putting MSTJ(21)=0 or increase the $\mathrm{s}/ \u$ ratio in fragmentation by putting PARJ(2)=0.40D0, to give but two examples. It is by exploring the possibilities offered here that PYTHIA can be turned into an extremely versatile tool, even if all the nice physics is already present in the default values.

As a final, semi-realistic example, assume that the $p_{\perp}$ spectrum of $\pi^+$ particles is to be studied in 91.2 GeV $\mathrm{e}^+\mathrm{e}^-$ annihilation events, where $p_{\perp}$ is to be defined with respect to the sphericity axis. Using the internal routines for simple histogramming, a complete program might look like

C...Double precision and integer declarations.
      IMPLICIT DOUBLE PRECISION(A-H, O-Z)
      IMPLICIT INTEGER(I-N)
      INTEGER PYK,PYCHGE,PYCOMP

C...Common blocks.
      COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)

C...Book histograms.
      CALL PYBOOK(1,'pT spectrum of pi+',100,0D0,5D0)

C...Number of events to generate. Loop over events.
      NEVT=100
      DO 110 IEVT=1,NEVT

C...Generate event. List first one.
        CALL PYEEVT(0,91.2D0)
        IF(IEVT.EQ.1) CALL PYLIST(1)

C...Find sphericity axis.
        CALL PYSPHE(SPH,APL)

C...Rotate event so that sphericity axis is along z axis.
        CALL PYEDIT(31)

C...Loop over all particles, but skip if not pi+.
        DO 100 I=1,N
          IF(K(I,2).NE.211) GOTO 100

C...Calculate pT and fill in histogram.
          PT=SQRT(P(I,1)**2+P(I,2)**2)
          CALL PYFILL(1,PT,1D0)

C...End of particle and event loops.
  100   CONTINUE
  110 CONTINUE

C...Normalize histogram properly and list it.
      CALL PYFACT(1,20D0/NEVT)
      CALL PYHIST

      END

Study this program, try to understand what happens at each step, and run it to check that it works. You should then be ready to look at the relevant sections of this report and start writing your own programs.