To exemplify the above discussion, consider the explicit case of or . These two processes are already available in PYTHIA, but without full spin correlations. One might therefore wish to include them from some external generator. A physics analysis would then most likely involve angular correlations intended to set limits on (or reveal traces of) anomalous couplings. However, so as to give a simple self-contained example, instead consider the analysis of the charged multiplicity distribution. This actually offers a simple first cross-check between the internal and external implementations of the same process. The main program might then look something like
IMPLICIT DOUBLE PRECISION(A-H, O-Z) IMPLICIT INTEGER(I-N) C...User process event common block. INTEGER MAXNUP PARAMETER (MAXNUP=500) INTEGER NUP,IDPRUP,IDUP,ISTUP,MOTHUP,ICOLUP DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,PUP,VTIMUP,SPINUP COMMON/HEPEUP/NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,IDUP(MAXNUP), &ISTUP(MAXNUP),MOTHUP(2,MAXNUP),ICOLUP(2,MAXNUP),PUP(5,MAXNUP), &VTIMUP(MAXNUP),SPINUP(MAXNUP) SAVE /HEPEUP/ C...PYTHIA common block. COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5) SAVE /PYJETS/ C...Initialize. CALL PYINIT('USER',' ',' ',0D0) C...Book histogram. Reset event counter. CALL PYBOOK(1,'Charged multiplicity',100,-1D0,199D0) NACC=0 C...Event loop; check that not at end of run; list first events. DO 100 IEV=1,1000 CALL PYEVNT IF(NUP.EQ.0) GOTO 110 NACC=NACC+1 IF(IEV.LE.3) CALL PYLIST(7) IF(IEV.LE.3) CALL PYLIST(2) C...Analyse event; end event loop. CALL PYEDIT(3) CALL PYFILL(1,DBLE(N),1D0) 100 CONTINUE C...Statistics and histograms. 110 CALL PYSTAT(1) CALL PYFACT(1,1D0/DBLE(NACC)) CALL PYHIST ENDThere PYINIT is called with 'USER' as first argument, implying that the rest is dummy. The event loop itself looks fairly familiar, but with two additions. One is that NUP is checked after each event, since NUP = 0 would signal a premature end of the run, with the external generator unable to return more events. This would be the case e.g. if events have been stored on file, and the end of this file is reached. The other is that CALL PYLIST(7) can be used to list the particle content of the HEPEUP common block (with some information omitted, on vertices and spin), so that one can easily compare this input with the output after PYTHIA processing, CALL PYLIST(2). An example of a PYLIST(7) listing would be
Event listing of user process at input (simplified) I IST ID Mothers Colours p_x p_y p_z E m 1 -1 21 0 0 101 109 0.000 0.000 269.223 269.223 0.000 2 -1 21 0 0 109 102 0.000 0.000 -225.566 225.566 0.000 3 2 6 1 2 101 0 72.569 153.924 -10.554 244.347 175.030 4 2 -6 1 2 0 102 -72.569 -153.924 54.211 250.441 175.565 5 1 5 3 0 101 0 56.519 33.343 53.910 85.045 4.500 6 2 24 3 0 0 0 16.050 120.581 -64.464 159.302 80.150 7 1 -5 4 0 0 102 44.127 -60.882 25.507 79.527 4.500 8 2 -24 4 0 0 0 -116.696 -93.042 28.705 170.914 78.184 9 1 2 6 0 103 0 -8.667 11.859 16.063 21.766 0.000 10 1 -1 6 0 0 103 24.717 108.722 -80.527 137.536 0.000 11 1 -2 8 0 0 104 -33.709 -22.471 -26.877 48.617 0.000 12 1 1 8 0 104 0 -82.988 -70.571 55.582 122.297 0.000Note the reverse listing of ID(UP) and IST(UP) relative to the HEPEUP order, to have better agreement with the PYJETS one. (The ID column is wider in real life, to allow for longer codes, but has here been reduced to fit the listing onto the page.)
The corresponding PYLIST(2) listing of course would be considerably longer, containing a complete event as it does. Also the particles above would there appear boosted by the effects of initial-state radiation and primordial ; copies of them further down in the event record would also include the effects of final-state radiation. The full story is available with MSTP(125) = 2, while the default listing omits some of the intermediate steps.
The PYINIT call will generate a call to the user-supplied routine UPINIT. It is here that we need to specify the details of the generation model. Assume, for instance, that - and -initiated events have been generated in two separate runs for Tevatron Run II, with weighted events stored in two separate files. By the end of each run, cross section and maximum weight information has also been obtained, and stored on separate files. Then UPINIT could look like
SUBROUTINE UPINIT C...Double precision and integer declarations. IMPLICIT DOUBLE PRECISION(A-H, O-Z) IMPLICIT INTEGER(I-N) C...User process initialization common block. INTEGER MAXPUP PARAMETER (MAXPUP=100) INTEGER IDBMUP,PDFGUP,PDFSUP,IDWTUP,NPRUP,LPRUP DOUBLE PRECISION EBMUP,XSECUP,XERRUP,XMAXUP COMMON/HEPRUP/IDBMUP(2),EBMUP(2),PDFGUP(2),PDFSUP(2), &IDWTUP,NPRUP,XSECUP(MAXPUP),XERRUP(MAXPUP),XMAXUP(MAXPUP), &LPRUP(MAXPUP) SAVE /HEPRUP/ C....Pythia common block - needed for setting PDF's; see below. COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200) SAVE /PYPARS/ C...Set incoming beams: Tevatron Run II. IDBMUP(1)=2212 IDBMUP(2)=-2212 EBMUP(1)=1000D0 EBMUP(2)=1000D0 C...Set PDF's of incoming beams: CTEQ 5L. C...Note that Pythia will not look at PDFGUP and PDFSUP. PDFGUP(1)=4 PDFSUP(1)=46 PDFGUP(2)=PDFGUP(1) PDFSUP(2)=PDFSUP(1) C...Set use of CTEQ 5L in internal Pythia code. MSTP(51)=7 C...If you want Pythia to use PDFLIB, you have to set it by hand. C...(You also have to ensure that the dummy routines C...PDFSET, STRUCTM and STRUCTP in Pythia are not linked.) C MSTP(52)=2 C MSTP(51)=1000*PDFGUP(1)+PDFSUP(1) C...Decide on weighting strategy: weighted on input, cross section known. IDWTUP=2 C...Number of external processes. NPRUP=2 C...Set up q qbar -> t tbar. OPEN(21,FILE='qqtt.info',FORM='unformatted',ERR=100) READ(21,ERR=100) XSECUP(1),XERRUP(1),XMAXUP(1) LPRUP(1)=661 OPEN(22,FILE='qqtt.events',FORM='unformatted',ERR=100) C...Set up g g -> t tbar. OPEN(23,FILE='ggtt.info',FORM='unformatted',ERR=100) READ(23,ERR=100) XSECUP(2),XERRUP(2),XMAXUP(2) LPRUP(2)=662 OPEN(24,FILE='ggtt.events',FORM='unformatted',ERR=100) RETURN C...Stop run if file operations fail. 100 WRITE(*,*) 'Error! File open or read failed. Program stopped.' STOP ENDHere unformatted read/write is used to reduce the size of the event files, but at the price of a platform dependence. Formatted files are preferred if they are to be shipped elsewhere. The rest should be self-explanatory.
Inside the event loop of the main program, PYEVNT will call UPEVNT to obtain the next parton-level event. In its simplest form, only a single READ statement would be necessary to read information on the next event, e.g. what is shown in the event listing earlier in this section, with a few additions. Then the routine could look like
SUBROUTINE UPEVNT C...Double precision and integer declarations. IMPLICIT DOUBLE PRECISION(A-H, O-Z) IMPLICIT INTEGER(I-N) C...User process event common block. INTEGER MAXNUP PARAMETER (MAXNUP=500) INTEGER NUP,IDPRUP,IDUP,ISTUP,MOTHUP,ICOLUP DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,PUP,VTIMUP,SPINUP COMMON/HEPEUP/NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,IDUP(MAXNUP), &ISTUP(MAXNUP),MOTHUP(2,MAXNUP),ICOLUP(2,MAXNUP),PUP(5,MAXNUP), &VTIMUP(MAXNUP),SPINUP(MAXNUP) SAVE /HEPEUP/ C...Pick file to read from, based on requested event type. IUNIT=22 IF(IDPRUP.EQ.662) IUNIT=24 C...Read event from this file. (Except that NUP and IDPRUP are known.) NUP=12 READ(IUNIT,ERR=100,END=100) XWGTUP,SCALUP,AQEDUP,AQCDUP, &(IDUP(I),ISTUP(I),MOTHUP(1,I),MOTHUP(2,I),ICOLUP(1,I), &ICOLUP(2,I),(PUP(J,I),J=1,5),VTIMUP(I),SPINUP(I),I=1,NUP) C...Return, with NUP=0 if read failed. RETURN 100 NUP=0 RETURN ENDHowever, in reality one might wish to save disk space by not storing redundant information. The XWGTUP and SCALUP numbers are vital, while AQEDUP and AQCDUP are purely informational and can be omitted. In a event, only the , , and flavours need be given, assuming that the particles are always stored in the same order. For a initial state, the flavour should be added to the list. The ISTUP, MOTHUP and ICOLUP information is the same in all events of a given process, except for a twofold ambiguity in the colour flow for initial states. All VTIMUP vanish and the SPINUP are uninteresting since spins have already been taken into account by the kinematics of the final fermions. (It would be different if one of the 's decayed leptonically to a .) Only the PUP values of the six final fermions need be given, since the other momenta and masses can be reconstructed from this, remembering that the two initial partons are massless and along the beam pipe. The final fermions are on the mass shell, so their masses need not be stored event by event, but can be read from a small table. The energy of a particle can be reconstructed from its momentum and mass. Overall transverse momentum conservation removes two further numbers. What remains is thus 5 integers and 18 real numbers, where the reals could well be stored in single precision. Of course, the code needed to unpack information stored this way would be lengthy but trivial. Even more compact storage strategies could be envisaged, e.g. only to save the weight and the seed of a dedicated random-number generator, to be used to generate the next parton-level event. It is up to you to find the optimal balance between disk space and coding effort.