bin/HARDCOL/hardcol.f

0001 C...Extract parton-level events according to Les Houches Accord.
0002 C...The Les Houches Event File produced here can then be used
0003 C...as input for hadron-level event simulation, also in Pythia8.
0004
0005 C...Note: you need to create two temporary files for MSTP(161) and MSTP(162).
0006 C...The final call to PYLHEF will pack them into a standard-compliant
0007 C...Les Houches Event File on your unit MSTP(163), and erase the two
0008 C...temporary files (unless you set MSTP(164)=1).
0009
0010 C...IMPORTANT: the PYLHEF routine attached below is necessary if you run
0011 C...with PYTHIA versions up to and including 6.403. Starting with 6.404
0012 C...this routine is already part of the standard distribution, so you
0013 C...should remove the copy in this file.
0014
0015 C...Double precision and integer declarations.
0016       SUBROUTINE HARDCOL
0017       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0018       IMPLICIT INTEGER(I-N)
0019
0020 C...Three PYTHIA functios retunr integers, so need declaring.
0021       INTEGER PYK,PYCHGE,PYCOMP
0022
0023 C...EXTERNAL statement links PYDATA on most machines.
0024       EXTERNAL PYDATA
0025
0026 C...Commonblocks.
0027       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
0028
0029 c...pythia common block.
0030       PARAMETER (MAXNUP=500)
0031       COMMON/HEPEUP/NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,IDUP(MAXNUP),
0032      &ISTUP(MAXNUP),MOTHUP(2,MAXNUP),ICOLUP(2,MAXNUP),PUP(5,MAXNUP),
0033      &VTIMUP(MAXNUP),SPINUP(MAXNUP)
0034       SAVE /HEPEUP/
0035
0036       PARAMETER (MAXPUP=100)
0037       INTEGER PDFGUP,PDFSUP,LPRUP
0038       COMMON/HEPRUP/IDBMUP(2),EBMUP(2),PDFGUP(2),PDFSUP(2),
0039      &IDWTUP,NPRUP,XSECUP(MAXPUP),XERRUP(MAXPUP),XMAXUP(MAXPUP),
0040      &LPRUP(MAXPUP)
0041       SAVE /HEPRUP/
0042
0043 c...user process event common block.
0044       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0045       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0046       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0047       COMMON/PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200)
0048       COMMON/PYDATR/MRPY(6),RRPY(100)
0049       COMMON/HCLPAR/ECM,NEV
0050
0051
0052 C...Temporary files for initialization/event output.
0053       MSTP(161)=21
0054       OPEN(21,FILE='HARDCOL.init',STATUS='unknown')
0055       MSTP(162)=22
0056       OPEN(22,FILE='HARDCOL.evnt',STATUS='unknown')
0057
0058 C...Final Les Houches Event File, obtained by combining above two.
0059       MSTP(163)=23
0060       OPEN(23,FILE='HARDCOL.lhe',STATUS='unknown')
0061
0062 c*********************************************
0063 c... HARDCOL
0064 c*********************************************
0065       CALL SETPARAMETERS
0066       CALL read_hcs_file
0067
0068 C...Initialize.
0069       CALL EVNTINIT
0070
0071 C...Fills the HEPRUP commonblock with info on incoming beams and allowed
0072 C...processes, and optionally stores that information on file.
0073 c      CALL HARDCOL_PYUPIN
0074 c      CALL PYUPIN
0075
0076 C...Event loop. List first few events.
0077       DO 200 IEV=1,NEV
0078         CALL PYUPEV
0079         IF(IEV.LE.2) THEN
0080            CALL PYLIST(2)
0081            CALL PYLIST(7)
0082         ENDIF
0083  200  CONTINUE
0084
0085 C...Final statistics.
0086       CALL PYSTAT(1)
0087       CALL PYUPIN
0088
0089 C...Produce final Les Houches Event File.
0090       CALL PYLHEF
0091
0092       END
0093
0094 C*********************************************************************
0095
0096 C...Combine the two old-style Pythia initialization and event files
0097 C...into a single Les Houches Event File.
0098
0099       SUBROUTINE PYLHEF
0100
0101 C...Double precision and integer declarations.
0102       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0103       IMPLICIT INTEGER(I-N)
0104
0105 C...PYTHIA commonblock: only used to provide read/write units and version.
0106       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0107       SAVE /PYPARS/
0108
0109 C...User process initialization commonblock.
0110       INTEGER MAXPUP
0111       PARAMETER (MAXPUP=100)
0112       INTEGER IDBMUP,PDFGUP,PDFSUP,IDWTUP,NPRUP,LPRUP
0113       DOUBLE PRECISION EBMUP,XSECUP,XERRUP,XMAXUP
0114       COMMON/HEPRUP/IDBMUP(2),EBMUP(2),PDFGUP(2),PDFSUP(2),
0115      &IDWTUP,NPRUP,XSECUP(MAXPUP),XERRUP(MAXPUP),XMAXUP(MAXPUP),
0116      &LPRUP(MAXPUP)
0117       SAVE /HEPRUP/
0118
0119 C...User process event common block.
0120       INTEGER MAXNUP
0121       PARAMETER (MAXNUP=500)
0122       INTEGER NUP,IDPRUP,IDUP,ISTUP,MOTHUP,ICOLUP
0123       DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,PUP,VTIMUP,SPINUP
0124       COMMON/HEPEUP/NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,IDUP(MAXNUP),
0125      &ISTUP(MAXNUP),MOTHUP(2,MAXNUP),ICOLUP(2,MAXNUP),PUP(5,MAXNUP),
0126      &VTIMUP(MAXNUP),SPINUP(MAXNUP)
0127       SAVE /HEPEUP/
0128
0129 C...Lines to read in assumed never longer than 200 characters.
0130       PARAMETER (MAXLEN=200)
0131       CHARACTER*(MAXLEN) STRING
0132
0133 C...Format for reading lines.
0134       CHARACTER*6 STRFMT
0135       STRFMT='(A000)'
0136       WRITE(STRFMT(3:5),'(I3)') MAXLEN
0137
0138 C...Rewind initialization and event files.
0139       REWIND MSTP(161)
0140       REWIND MSTP(162)
0141
0142 C...Write header info.
0143       WRITE(MSTP(163),'(A)') '<LesHouchesEvents version="1.0">'
0144       WRITE(MSTP(163),'(A)') '<!--'
0145       WRITE(MSTP(163),'(A,I1,A1,I3)') 'File generated with PYTHIA ',
0146      &MSTP(181),'.',MSTP(182)
0147       WRITE(MSTP(163),'(A)') '-->'
0148
0149 C...Read first line of initialization info and get number of processes.
0150       READ(MSTP(161),'(A)',END=300,ERR=300) STRING
0151       READ(STRING,*,ERR=300) IDBMUP(1),IDBMUP(2),EBMUP(1),
0152      &EBMUP(2),PDFGUP(1),PDFGUP(2),PDFSUP(1),PDFSUP(2),IDWTUP,NPRUP
0153
0154 C...Copy initialization lines, omitting trailing blanks.
0155 C...Embed in <init> ... </init> block.
0156       WRITE(MSTP(163),'(A)') '<init>'
0157       DO 120 IPR=0,NPRUP
0158         IF(IPR.GT.0) READ(MSTP(161),'(A)',END=300,ERR=300) STRING
0159         LEN=MAXLEN+1
0160   110   LEN=LEN-1
0161         IF(LEN.GT.1.AND.STRING(LEN:LEN).EQ.' ') GOTO 110
0162         WRITE(MSTP(163),'(A)',ERR=300) STRING(1:LEN)
0163   120 CONTINUE
0164       WRITE(MSTP(163),'(A)') '</init>'
0165
0166 C...Begin event loop.
0167   200 CONTINUE
0168
0169 C...Read first line of event info and get number of particles.
0170       READ(MSTP(162),'(A)',END=280,ERR=300) STRING
0171       READ(STRING,*,ERR=300) NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP
0172
0173 C...Copy event lines, omitting trailing blanks.
0174 C...Embed in <event> ... </event> block.
0175       WRITE(MSTP(163),'(A)') '<event>'
0176       DO 220 I=0,NUP
0177         IF(I.GT.0) READ(MSTP(162),'(A)',END=300,ERR=300) STRING
0178         LEN=MAXLEN+1
0179   210   LEN=LEN-1
0180         IF(LEN.GT.1.AND.STRING(LEN:LEN).EQ.' ') GOTO 210
0181         WRITE(MSTP(163),'(A)',ERR=300) STRING(1:LEN)
0182   220 CONTINUE
0183       WRITE(MSTP(163),'(A)') '</event>'
0184
0185 C...Loop back to look for next event.
0186       GOTO 200
0187
0188 C...Successfully reached end of event loop: write closing tag
0189 C...and remove temporary intermediate files (unless asked not to).
0190   280 WRITE(MSTP(163),'(A)') '</LesHouchesEvents>'
0191       IF(MSTP(164).EQ.1) RETURN
0192       CLOSE(MSTP(161),ERR=300,STATUS='DELETE')
0193       CLOSE(MSTP(162),ERR=300,STATUS='DELETE')
0194       RETURN
0195
0196 C...Error exit.
0197   300 WRITE(*,*) ' PYLHEF file joining failed!'
0198
0199       RETURN
0200       END