Getting Started with the Event Generation Machinery

A run with the full PYTHIA event generation machinery has to be more strictly organized than the simple examples above, in that it is necessary to initialize the generation before events can be generated, and in that it is not possible to change switches and parameters freely during the course of the run. A fairly precise recipe for how a run should be structured can therefore be given.

Thus, the nowadays normal usage of PYTHIA can be subdivided into three steps.

11.

The initialization step. It is here that all the basic characteristics of the coming generation are specified. The material in this section includes the following.

$\bullet$

Declarations for double precision and integer variables and integer functions:

 
    IMPLICIT DOUBLE PRECISION(A-H, O-Z)
    IMPLICIT INTEGER(I-N)
    INTEGER PYK,PYCHGE,PYCOMP

$\bullet$

Common blocks, at least the following, and maybe some more:

    COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
    COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
    COMMON/PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200)
    COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)

$\bullet$

Selection of required processes. Some fixed `menus' of subprocesses can be selected with different MSEL values, but with MSEL=0 it is possible to compose `à la carte', using the subprocess numbers. To generate processes 14, 18 and 29, for instance, one needs

    MSEL=0
    MSUB(14)=1
    MSUB(18)=1
    MSUB(29)=1

$\bullet$

Selection of kinematics cuts in the CKIN array. To generate hard scatterings with 5 GeV $\leq p_{\perp}\leq$ 10 GeV, for instance, use

    CKIN(3)=5D0
    CKIN(4)=10D0

Unfortunately, initial- and final-state radiation will shift around the kinematics of the hard scattering, making the effects of cuts less predictable. One therefore always has to be very careful that no desired event configurations are cut out.

$\bullet$

Definition of underlying physics scenario, e.g. Higgs mass.

$\bullet$

Selection of parton-distribution sets, $Q^2$ definitions, and all other details of the generation.

$\bullet$

Switching off of generator parts not needed for toy simulations, e.g. fragmentation for parton level studies.

$\bullet$

Initialization of the event generation procedure. Here kinematics is set up, maxima of differential cross sections are found for future Monte Carlo generation, and a number of other preparatory tasks carried out. Initialization is performed by PYINIT, which should be called only after the switches and parameters above have been set to their desired values. The frame, the beam particles and the energy have to be specified, e.g.

    CALL PYINIT('CMS','p','pbar',1800D0)

$\bullet$

Any other initial material required by you, e.g. histogram booking.

12.

The generation loop. It is here that events are generated and studied. It includes the following tasks:

$\bullet$

Generation of the next event, with

    CALL PYEVNT

$\bullet$

Printing of a few events, to check that everything is working as planned, with

    CALL PYLIST(1)

$\bullet$

An analysis of the event for properties of interest, either directly reading out information from the PYJETS common block or making use of the utility routines in PYTHIA.

$\bullet$

Saving of events on disk or tape, or interfacing to detector simulation.

13.

The finishing step. Here the tasks are:

$\bullet$

Printing a table of deduced cross sections, obtained as a by-product of the Monte Carlo generation activity, with the command

    CALL PYSTAT(1)

$\bullet$

Printing histograms and other user output.

To illustrate this structure, imagine a toy example, where one wants to simulate the production of a 300 GeV Higgs particle. In PYTHIA, a program for this might look something like the following.

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)
      COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
      COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
      COMMON/PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
      COMMON/PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200)
      COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
 
C...Number of events to generate. Switch on proper processes.
      NEV=1000
      MSEL=0
      MSUB(102)=1
      MSUB(123)=1
      MSUB(124)=1
 
C...Select Higgs mass and kinematics cuts in mass.
      PMAS(25,1)=300D0
      CKIN(1)=290D0
      CKIN(2)=310D0
 
C...For simulation of hard process only: cut out unnecessary tasks.
      MSTP(61)=0
      MSTP(71)=0
      MSTP(81)=0
      MSTP(111)=0
 
C...Initialize and list partial widths.
      CALL PYINIT('CMS','p','p',14000D0)
      CALL PYSTAT(2)
 
C...Book histogram.
      CALL PYBOOK(1,'Higgs mass',50,275D0,325D0)
 
C...Generate events. Look at first few.
      DO 200 IEV=1,NEV
        CALL PYEVNT
        IF(IEV.LE.3) CALL PYLIST(1)
 
C...Loop over particles to find Higgs and histogram its mass.
        DO 100 I=1,N
          IF(K(I,1).LT.20.AND.K(I,2).EQ.25) HMASS=P(I,5)
  100   CONTINUE
        CALL PYFILL(1,HMASS,1D0)
  200 CONTINUE
 
C...Print cross sections and histograms.
      CALL PYSTAT(1)
      CALL PYHIST
 
      END

Here 102, 123 and 124 are the three main Higgs production graphs $\mathrm{g}\mathrm{g}\rightarrow \mathrm{h}$, $\mathrm{Z}\mathrm{Z}\rightarrow \mathrm{h}$, and $\mathrm{W}\mathrm{W}\rightarrow \mathrm{h}$, and MSUB(ISUB)=1 is the command to switch on process ISUB. Full freedom to combine subprocesses `à la carte' is ensured by MSEL=0; ready-made `menus' can be ordered with other MSEL numbers. The PMAS command sets the mass of the Higgs, and the CKIN variables the desired mass range of the Higgs -- a Higgs with a 300 GeV nominal mass actually has a fairly broad Breit-Wigner type mass distribution. The MSTP switches that come next are there to modify the generation procedure, in this case to switch off initial- and final-state radiation, multiple interactions among beam jets, and fragmentation, to give only the `parton skeleton' of the hard process. The PYINIT call initializes PYTHIA, by finding maxima of cross sections, recalculating the Higgs decay properties (which depend on the Higgs mass), etc. The decay properties can be listed with PYSTAT(2).

Inside the event loop, PYEVNT is called to generate an event, and PYLIST(1) to list the event. The information used by PYLIST(1) is the event record, stored in the common block PYJETS. Here one finds all produced particles, both final and intermediate ones, with information on particle species and event history (K array), particle momenta (P array) and production vertices (V array). In the loop over all particles produced, 1 through N, the Higgs particle is found by its code, K(I,2)=25, and its mass is stored in P(I,5).

After all events have been generated, PYSTAT(1) gives a summary of the number of events generated in the various allowed channels, and the inferred cross sections.

In the run above, a typical event listing might look like the following.

                         Event listing (summary)
 
    I  particle/jet     KF    p_x      p_y      p_z       E        m
 
    1  !p+!           2212    0.000    0.000 8000.000 8000.000    0.938
    2  !p+!           2212    0.000    0.000-8000.000 8000.000    0.938
 ======================================================================
    3  !g!              21   -0.505   -0.229   28.553   28.558    0.000
    4  !g!              21    0.224    0.041 -788.073  788.073    0.000
    5  !g!              21   -0.505   -0.229   28.553   28.558    0.000
    6  !g!              21    0.224    0.041 -788.073  788.073    0.000
    7  !H0!             25   -0.281   -0.188 -759.520  816.631  300.027
    8  !W+!             24  120.648   35.239 -397.843  424.829   80.023
    9  !W-!            -24 -120.929  -35.426 -361.677  391.801   82.579
   10  !e+!            -11   12.922   -4.760 -160.940  161.528    0.001
   11  !nu_e!           12  107.726   39.999 -236.903  263.302    0.000
   12  !s!               3  -62.423    7.195 -256.713  264.292    0.199
   13  !cbar!           -4  -58.506  -42.621 -104.963  127.509    1.350
 ======================================================================
   14  (H0)             25   -0.281   -0.188 -759.520  816.631  300.027
   15  (W+)             24  120.648   35.239 -397.843  424.829   80.023
   16  (W-)            -24 -120.929  -35.426 -361.677  391.801   82.579
   17  e+              -11   12.922   -4.760 -160.940  161.528    0.001
   18  nu_e             12  107.726   39.999 -236.903  263.302    0.000
   19  s         A       3  -62.423    7.195 -256.713  264.292    0.199
   20  cbar      V      -4  -58.506  -42.621 -104.963  127.509    1.350
   21  ud_1      A    2103   -0.101    0.176 7971.328 7971.328    0.771
   22  d         V       1   -0.316    0.001  -87.390   87.390    0.010
   23  u         A       2    0.606    0.052   -0.751    0.967    0.006
   24  uu_1      V    2203    0.092   -0.042-7123.668 7123.668    0.771
 ======================================================================
                sum:  2.00     0.00     0.00     0.00 15999.98 15999.98

The above event listing is abnormally short, in part because some columns of information were removed to make it fit into this text, in part because all initial- and final-state QCD radiation, all non-trivial beam jet structure, and all fragmentation was inhibited in the generation. Therefore only the skeleton of the process is visible. In lines 1 and 2 one recognizes the two incoming protons. In lines 3 and 4 are incoming partons before initial-state radiation and in 5 and 6 after -- since there is no such radiation they coincide here. Line 7 shows the Higgs produced by $\mathrm{g}\mathrm{g}$ fusion, 8 and 9 its decay products and 10-13 the second-step decay products. Up to this point lines give a summary of the event history, indicated by the exclamation marks that surround particle names (and also reflected in the K(I,1) code, not shown). From line 14 onwards come the particles actually produced in the final states, first in lines 14-16 particles that subsequently decayed, which have their names surrounded by brackets, and finally the particles and partons left in the end, including beam remnants. Here this also includes a number of unfragmented partons, since fragmentation was inhibited. Ordinarily, the listing would have gone on for a few hundred more lines, with the particles produced in the fragmentation and their decay products. The final line gives total charge and momentum, as a convenient check that nothing unexpected happened. The first column of the listing is just a counter, the second gives the particle name and information on status and string drawing (the A and V), the third the particle-flavour code (which is used to give the name), and the subsequent columns give the momentum components.

One of the main problems is to select kinematics efficiently. Imagine for instance that one is interested in the production of a single $\mathrm{Z}$ with a transverse momentum in excess of 50 GeV. If one tries to generate the inclusive sample of $\mathrm{Z}$ events, by the basic production graphs $\mathrm{q}\overline{\mathrm{q}}\rightarrow \mathrm{Z}$, then most events will have low transverse momenta and will have to be discarded. That any of the desired events are produced at all is due to the initial-state generation machinery, which can build up transverse momenta for the incoming $\mathrm{q}$ and $\overline{\mathrm{q}}$. However, the amount of initial-state radiation cannot be constrained beforehand. To increase the efficiency, one may therefore turn to the higher-order processes $\mathrm{q}\mathrm{g}\rightarrow \mathrm{Z}\mathrm{q}$ and $\mathrm{q}\overline{\mathrm{q}}\rightarrow \mathrm{Z}\mathrm{g}$, where already the hard subprocess gives a transverse momentum to the $\mathrm{Z}$. This transverse momentum can be constrained as one wishes, but again initial- and final-state radiation will smear the picture. If one were to set a $p_{\perp}$ cut at 50 GeV for the hard-process generation, those events where the $\mathrm{Z}$ was given only 40 GeV in the hard process but got the rest from initial-state radiation would be missed. Not only therefore would cross sections come out wrong, but so might the typical event shapes. In the end, it is therefore necessary to find some reasonable compromise, by starting the generation at 30 GeV, say, if one knows that only rarely do events below this value fluctuate up to 50 GeV. Of course, most events will therefore not contain a $\mathrm{Z}$ above 50 GeV, and one will have to live with some inefficiency. It is not uncommon that only one event out of ten can be used, and occasionally it can be even worse.

If it is difficult to set kinematics, it is often easier to set the flavour content of a process. In a Higgs study, one might wish, for example, to consider the decay $\mathrm{h}^0 \rightarrow \mathrm{Z}^0 \mathrm{Z}^0$, with each $\mathrm{Z}^0 \rightarrow \mathrm{e}^+\mathrm{e}^-$ or $\mu^+ \mu^-$. It is therefore necessary to inhibit all other $\mathrm{h}^0$ and $\mathrm{Z}^0$ decay channels, and also to adjust cross sections to take into account this change, all of which is fairly straightforward. The same cannot be said for decays of ordinary hadrons, where the number produced in a process is not known beforehand, and therefore inconsistencies easily can arise if one tries to force specific decay channels.

In the examples given above, all run-specific parameters are set in the code (in the main program; alternatively it could be in a subroutine called by the main program). This approach is allowing maximum flexibility to change parameters during the course of the run. However, in many experimental collaborations one does not want to allow this freedom, but only one set of parameters, to be read in from an external file at the beginning of a run and thereafter never changed. This in particular applies when PYTHIA is to be linked with other libraries, such as GEANT [Bru89] and detector-specific software. While a linking of a normal-sized main program with PYTHIA is essentially instantaneous on current platforms (typically less than a second), this may not hold for other libraries. For this purpose one then needs a parser of PYTHIA parameters, the core of which can be provided by the PYGIVE routine.

As an example, consider a main program of the form

C...Double precision and integer declarations.
      IMPLICIT DOUBLE PRECISION(A-H, O-Z)
      IMPLICIT INTEGER(I-N)
      INTEGER PYK,PYCHGE,PYCOMP
C...Input and output strings.
      CHARACTER FRAME*12,BEAM*12,TARGET*12,PARAM*100

C...Read parameters for PYINIT call.
      READ(*,*) FRAME,BEAM,TARGET,ENERGY

C...Read number of events to generate, and to print.
      READ(*,*) NEV,NPRT

C...Loop over reading and setting parameters/switches.
  100 READ(*,'(A)',END=200) PARAM
      CALL PYGIVE(PARAM)
      GOTO 100

C...Initialize PYTHIA.
  200 CALL PYINIT(FRAME,BEAM,TARGET,ENERGY)

C...Event generation loop
      DO 300 IEV=1,NEV
        CALL PYEVNT
        IF(IEV.LE.NPRT) CALL PYLIST(1)
  300 CONTINUE   

C...Print cross sections.
      CALL PYSTAT(1)

      END

and a file indata with the contents

CMS,p,p,14000.
1000,3
! below follows commands sent to PYGIVE
MSEL=0           ! Mix processes freely
MSUB(102)=1      ! g + g -> h0 
MSUB(123)=1      ! Z0 + Z0 -> h0
MSUB(124)=1      ! W+ + W- -> h0 
PMAS(25,1)=300.  ! Higgs mass
CKIN(1)=290.     ! lower cutoff on mass
CKIN(2)=310.     ! upper cutoff on mass
MSTP(61)=0       ! no initial-state showers
MSTP(71)=0       ! no final-state showers    
MSTP(81)=0       ! no multiple interactions   
MSTP(111)=0      ! no hadronization

Here the text following the exclamation marks is interpreted as a comment by PYGIVE, and thus purely intended to allow better documentation of changes. The main program could then be linked to PYTHIA, to an executable a.out, and run e.g. with a Unix command line

  a.out < indata > output

to produce results on the file output. Here the indata could be changed without requiring a recompilation. Of course, the main program would have to be more realistic, e.g. with events saved to disk or tape, but the principle should be clear.