pyhepc.f

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pyhepc.f

 
C*********************************************************************
 
C...PYHEPC
C...Converts PYTHIA event record contents to or from
C...the standard event record commonblock.
 
      SUBROUTINE pyhepc(MCONV)
 
C...Double precision and integer declarations.
      IMPLICIT DOUBLE PRECISION(a-h, o-z)
      IMPLICIT INTEGER(i-n)
      INTEGER pyk,pychge,pycomp
C...Commonblocks.
      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)
      SAVE /pyjets/,/pydat1/,/pydat2/
C...HEPEVT commonblock.
      parameter(nmxhep=4000)
      common/hepevt/nevhep,nhep,isthep(nmxhep),idhep(nmxhep),
     &jmohep(2,nmxhep),jdahep(2,nmxhep),phep(5,nmxhep),vhep(4,nmxhep)
      DOUBLE PRECISION phep,vhep
      SAVE /hepevt/

C...Store HEPEVT commonblock size (for interfacing issues).
      mstu(8)=nmxhep
 
C...Conversion from PYTHIA to standard, the easy part.
      IF(mconv.EQ.1) THEN
        nevhep=0
        IF(n.GT.nmxhep) CALL pyerrm(8,
     &  '(PYHEPC:) no more space in /HEPEVT/')
        nhep=min(n,nmxhep)
        DO 150 i=1,nhep
          isthep(i)=0
          IF(k(i,1).GE.1.AND.k(i,1).LE.10) isthep(i)=1
          IF(k(i,1).GE.11.AND.k(i,1).LE.20) isthep(i)=2
          IF(k(i,1).GE.21.AND.k(i,1).LE.30) isthep(i)=3
          IF(k(i,1).GE.31.AND.k(i,1).LE.100) isthep(i)=k(i,1)
          idhep(i)=k(i,2)
          jmohep(1,i)=k(i,3)
          jmohep(2,i)=0
          IF(k(i,1).NE.3.AND.k(i,1).NE.13.AND.k(i,1).NE.14) THEN
            jdahep(1,i)=k(i,4)
            jdahep(2,i)=k(i,5)
          ELSE
            jdahep(1,i)=0
            jdahep(2,i)=0
          ENDIF
          DO 100 j=1,5
            phep(j,i)=p(i,j)
  100     CONTINUE
          DO 110 j=1,4
            vhep(j,i)=v(i,j)
  110     CONTINUE
 
C...Check if new event (from pileup).
          IF(i.EQ.1) THEN
            inew=1
          ELSE
            IF(k(i,1).EQ.21.AND.k(i-1,1).NE.21) inew=i
          ENDIF
 
C...Fill in missing mother information.
          IF(i.GE.inew+2.AND.k(i,1).EQ.21.AND.k(i,3).EQ.0) THEN
            imo1=i-2
  120       IF(imo1.GT.inew.AND.k(imo1+1,1).EQ.21.AND.k(imo1+1,3).EQ.0)
     &      THEN
              imo1=imo1-1
              goto 120
            ENDIF
            jmohep(1,i)=imo1
            jmohep(2,i)=imo1+1
          ELSEIF(k(i,2).GE.91.AND.k(i,2).LE.93) THEN
            i1=k(i,3)-1
  130       i1=i1+1
            IF(i1.GE.i) CALL pyerrm(8,
     &      '(PYHEPC:) translation of inconsistent event history')
            IF(i1.LT.i.AND.k(i1,1).NE.1.AND.k(i1,1).NE.11) goto 130
            kc=pycomp(k(i1,2))
            IF(i1.LT.i.AND.kc.EQ.0) goto 130
            IF(i1.LT.i.AND.kchg(kc,2).EQ.0) goto 130
            jmohep(2,i)=i1
          ELSEIF(k(i,2).EQ.94) THEN
            njet=2
            IF(nhep.GE.i+3.AND.k(i+3,3).LE.i) njet=3
            IF(nhep.GE.i+4.AND.k(i+4,3).LE.i) njet=4
            jmohep(2,i)=mod(k(i+njet,4)/mstu(5),mstu(5))
            IF(jmohep(2,i).EQ.jmohep(1,i)) jmohep(2,i)=
     &      mod(k(i+1,4)/mstu(5),mstu(5))
          ENDIF
 
C...Fill in missing daughter information.
          IF(k(i,2).EQ.94.AND.mstu(16).NE.2) THEN
            DO 140 i1=jdahep(1,i),jdahep(2,i)
              i2=mod(k(i1,4)/mstu(5),mstu(5))
              jdahep(1,i2)=i
  140       CONTINUE
          ENDIF
          IF(k(i,2).GE.91.AND.k(i,2).LE.94) goto 150
          i1=jmohep(1,i)
          IF(i1.LE.0.OR.i1.GT.nhep) goto 150
          IF(k(i1,1).NE.13.AND.k(i1,1).NE.14) goto 150
          IF(jdahep(1,i1).EQ.0) THEN
            jdahep(1,i1)=i
          ELSE
            jdahep(2,i1)=i
          ENDIF
  150   CONTINUE
        DO 160 i=1,nhep
          IF(k(i,1).NE.13.AND.k(i,1).NE.14) goto 160
          IF(jdahep(2,i).EQ.0) jdahep(2,i)=jdahep(1,i)
  160   CONTINUE
 
C...Conversion from standard to PYTHIA, the easy part.
      ELSE
        IF(nhep.GT.mstu(4)) CALL pyerrm(8,
     &  '(PYHEPC:) no more space in /PYJETS/')
        n=min(nhep,mstu(4))
        nkq=0
        kqsum=0
        DO 190 i=1,n
          k(i,1)=0
          IF(isthep(i).EQ.1) k(i,1)=1
          IF(isthep(i).EQ.2) k(i,1)=11
          IF(isthep(i).EQ.3) k(i,1)=21
          k(i,2)=idhep(i)
          k(i,3)=jmohep(1,i)
          k(i,4)=jdahep(1,i)
          k(i,5)=jdahep(2,i)
          DO 170 j=1,5
            p(i,j)=phep(j,i)
  170     CONTINUE
          DO 180 j=1,4
            v(i,j)=vhep(j,i)
  180     CONTINUE
          v(i,5)=0d0
          IF(isthep(i).EQ.2.AND.phep(4,i).GT.phep(5,i)) THEN
            i1=jdahep(1,i)
            IF(i1.GT.0.AND.i1.LE.nhep) v(i,5)=(vhep(4,i1)-vhep(4,i))*
     &      phep(5,i)/phep(4,i)
          ENDIF
 
C...Fill in missing information on colour connection in jet systems.
          IF(isthep(i).EQ.1) THEN
            kc=pycomp(k(i,2))
            kq=0
            IF(kc.NE.0) kq=kchg(kc,2)*isign(1,k(i,2))
            IF(kq.NE.0) nkq=nkq+1
            IF(kq.NE.2) kqsum=kqsum+kq
            IF(kq.NE.0.AND.kqsum.NE.0) THEN
              k(i,1)=2
            ELSEIF(kq.EQ.2.AND.i.LT.n) THEN
              IF(k(i+1,2).EQ.21) k(i,1)=2
            ENDIF
          ENDIF
  190   CONTINUE
        IF(nkq.EQ.1.OR.kqsum.NE.0) CALL pyerrm(8,
     &  '(PYHEPC:) input parton configuration not colour singlet')
      ENDIF
 
      END