next up previous contents
Next: PYTHIA as a generator Up: How to Include External Previous: Event information   Contents

An example

To exemplify the above discussion, consider the explicit case of $\mathrm{q}\overline{\mathrm{q}}$ or $\mathrm{g}\mathrm{g}\to \t\overline{\mathrm{t}}\to \b\mathrm{W}^+\overline{\mat...
...rline{\mathrm{q}}_2 \, \overline{\mathrm{b}}\mathrm{q}_3\overline{\mathrm{q}}_4$. 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

      END
There 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.000
Note 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 $k_{\perp}$; 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 $\mathrm{q}\overline{\mathrm{q}}$- and $\mathrm{g}\mathrm{g}$-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   

      END
Here 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
      END
However, 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 $\mathrm{g}\mathrm{g}\to \t\overline{\mathrm{t}}\to \b\mathrm{W}^+\overline{\mat...
...rline{\mathrm{q}}_2 \, \overline{\mathrm{b}}\mathrm{q}_3\overline{\mathrm{q}}_4$ event, only the $\mathrm{q}_1$, $\overline{\mathrm{q}}_2$, $\mathrm{q}_3$ and $\overline{\mathrm{q}}_4$ flavours need be given, assuming that the particles are always stored in the same order. For a $\mathrm{q}\overline{\mathrm{q}}$ initial state, the $\mathrm{q}$ 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 $\mathrm{g}\mathrm{g}$ 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 $\mathrm{W}$'s decayed leptonically to a $\tau$.) 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.


next up previous contents
Next: PYTHIA as a generator Up: How to Include External Previous: Event information   Contents
Stephen Mrenna 2007-10-30