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.
IMPLICIT DOUBLE PRECISION(A-H, O-Z) IMPLICIT INTEGER(I-N) INTEGER PYK,PYCHGE,PYCOMP
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)
MSEL=0 MSUB(14)=1 MSUB(18)=1 MSUB(29)=1
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 , , and , 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.98The 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 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 with a transverse momentum in excess of 50 GeV. If one tries to generate the inclusive sample of events, by the basic production graphs , 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 and . However, the amount of initial-state radiation cannot be constrained beforehand. To increase the efficiency, one may therefore turn to the higher-order processes and , where already the hard subprocess gives a transverse momentum to the . 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 cut at 50 GeV for the hard-process generation, those events where the 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 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 , with each or . It is therefore necessary to inhibit all other and 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) ENDand 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 hadronizationHere 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 > outputto 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.