Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:29:54

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