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 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) ENDand 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 (2) and (), 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 and . 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 in the second. These are the codes actually stored to denote the presence of a and a , 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, .)
In the orig (origin) column, the zeros indicate that and are two initial entries. The subsequent line, number 3, denotes the fragmenting 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 quark or only the one. It is the string system as a whole that gives a , a , a , a , a , a , and a . 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 in line 13 and the in line 14 have origin 5, as an indication that they come from the decay of the 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 . 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 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) ENDwhere a 3-jet 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 I 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 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 spectrum of particles is to be studied in 91.2 GeV annihilation events, where 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 ENDStudy 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.