py4frm.f

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

 
C*********************************************************************
 
C...PY4FRM
C...An interface from a four-fermion generator to include
C...parton showers and hadronization.
 
      SUBROUTINE py4frm(ATOTSQ,A1SQ,A2SQ,ISTRAT,IRAD,ITAU,ICOM)
 
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/pypars/mstp(200),parp(200),msti(200),pari(200)
      common/pyint1/mint(400),vint(400)
      SAVE /pyjets/,/pydat1/,/pypars/,/pyint1/
C...Local arrays.
      dimension ijoin(2),intau(4)
 
C...Call PYHEPC to convert input from HEPEVT to PYJETS common.
      IF(icom.EQ.0) THEN
        mstu(28)=0
        CALL pyhepc(2)
      ENDIF
 
C...Loop through entries and pick up all final fermions/antifermions.
      i1=0
      i2=0
      i3=0
      i4=0
      DO 100 i=1,n
      IF(k(i,1).LE.0.OR.k(i,1).GT.10) goto 100
      kfa=iabs(k(i,2))
      IF((kfa.GE.1.AND.kfa.LE.6).OR.(kfa.GE.11.AND.kfa.LE.16)) THEN
        IF(k(i,2).GT.0) THEN
          IF(i1.EQ.0) THEN
            i1=i
          ELSEIF(i3.EQ.0) THEN
            i3=i
          ELSE
            CALL pyerrm(16,'(PY4FRM:) more than two fermions')
          ENDIF
        ELSE
          IF(i2.EQ.0) THEN
            i2=i
          ELSEIF(i4.EQ.0) THEN
            i4=i
          ELSE
            CALL pyerrm(16,'(PY4FRM:) more than two antifermions')
          ENDIF
        ENDIF
      ENDIF
  100 CONTINUE
 
C...Check that event is arranged according to conventions.
      IF(i3.EQ.0.OR.i4.EQ.0) THEN
        CALL pyerrm(16,'(PY4FRM:) event contains too few fermions')
      ENDIF
      IF(i2.LT.i1.OR.i3.LT.i2.OR.i4.LT.i3) THEN
        CALL pyerrm(6,'(PY4FRM:) fermions arranged in wrong order')
      ENDIF
 
C...Check which fermion pairs are quarks and which leptons.
      IF(iabs(k(i1,2)).LT.10.AND.iabs(k(i2,2)).LT.10) THEN
        iql12=1
      ELSEIF(iabs(k(i1,2)).GT.10.AND.iabs(k(i2,2)).GT.10) THEN
        iql12=2
      ELSE
        CALL pyerrm(16,'(PY4FRM:) first fermion pair inconsistent')
      ENDIF
      IF(iabs(k(i3,2)).LT.10.AND.iabs(k(i4,2)).LT.10) THEN
        iql34=1
      ELSEIF(iabs(k(i3,2)).GT.10.AND.iabs(k(i4,2)).GT.10) THEN
        iql34=2
      ELSE
        CALL pyerrm(16,'(PY4FRM:) second fermion pair inconsistent')
      ENDIF
 
C...Decide whether to allow or not photon radiation in showers.
      mstj(41)=2
      IF(irad.EQ.0) mstj(41)=1
 
C...Decide on dipole pairing.
      ip1=i1
      ip2=i2
      ip3=i3
      ip4=i4
      IF(iql12.EQ.iql34) THEN
        r1sq=a1sq
        r2sq=a2sq
        delta=atotsq-a1sq-a2sq
        IF(istrat.EQ.1) THEN
          IF(delta.GT.0d0) r1sq=r1sq+delta
          IF(delta.LT.0d0) r2sq=max(0d0,r2sq+delta)
        ELSEIF(istrat.EQ.2) THEN
          IF(delta.GT.0d0) r2sq=r2sq+delta
          IF(delta.LT.0d0) r1sq=max(0d0,r1sq+delta)
        ENDIF
        IF(r2sq.GT.pyr(0)*(r1sq+r2sq)) THEN
          ip2=i4
          ip4=i2
        ENDIF
      ENDIF
 
C...If colour reconnection then bookkeep W+W- or Z0Z0
C...and copy q qbar q qbar consecutively.
      IF(mstp(115).GE.1.AND.iql12.EQ.1.AND.iql34.EQ.1) THEN
        k(n+1,1)=11
        k(n+1,3)=ip1
        k(n+1,4)=n+3
        k(n+1,5)=n+4
        k(n+2,1)=11
        k(n+2,3)=ip3
        k(n+2,4)=n+5
        k(n+2,5)=n+6
        IF(k(ip1,2)+k(ip2,2).EQ.0) THEN
          k(n+1,2)=23
          k(n+2,2)=23
          mint(1)=22
        ELSEIF(pychge(k(ip1,2)).GT.0) THEN
          k(n+1,2)=24
          k(n+2,2)=-24
          mint(1)=25
        ELSE
          k(n+1,2)=-24
          k(n+2,2)=24
          mint(1)=25
        ENDIF
        DO 110 j=1,5
          k(n+3,j)=k(ip1,j)
          k(n+4,j)=k(ip2,j)
          k(n+5,j)=k(ip3,j)
          k(n+6,j)=k(ip4,j)
          p(n+1,j)=p(ip1,j)+p(ip2,j)
          p(n+2,j)=p(ip3,j)+p(ip4,j)
          p(n+3,j)=p(ip1,j)
          p(n+4,j)=p(ip2,j)
          p(n+5,j)=p(ip3,j)
          p(n+6,j)=p(ip4,j)
          v(n+1,j)=v(ip1,j)
          v(n+2,j)=v(ip3,j)
          v(n+3,j)=v(ip1,j)
          v(n+4,j)=v(ip2,j)
          v(n+5,j)=v(ip3,j)
          v(n+6,j)=v(ip4,j)
  110   CONTINUE
        p(n+1,5)=sqrt(max(0d0,p(n+1,4)**2-p(n+1,1)**2-p(n+1,2)**2-
     &  p(n+1,3)**2))
        p(n+2,5)=sqrt(max(0d0,p(n+2,4)**2-p(n+2,1)**2-p(n+2,2)**2-
     &  p(n+2,3)**2))
        k(n+3,3)=n+1
        k(n+4,3)=n+1
        k(n+5,3)=n+2
        k(n+6,3)=n+2
C...Remove original q qbar q qbar and update counters.
        k(ip1,1)=k(ip1,1)+10
        k(ip2,1)=k(ip2,1)+10
        k(ip3,1)=k(ip3,1)+10
        k(ip4,1)=k(ip4,1)+10
        iw1=n+1
        iw2=n+2
        nsd1=n+2
        ip1=n+3
        ip2=n+4
        ip3=n+5
        ip4=n+6
        n=n+6
      ENDIF
 
C...Do colour joinings and parton showers.
      IF(iql12.EQ.1) THEN
        ijoin(1)=ip1
        ijoin(2)=ip2
        CALL pyjoin(2,ijoin)
      ENDIF
      IF(iql12.EQ.1.OR.irad.EQ.1) THEN
        pm12s=(p(ip1,4)+p(ip2,4))**2-(p(ip1,1)+p(ip2,1))**2-
     &  (p(ip1,2)+p(ip2,2))**2-(p(ip1,3)+p(ip2,3))**2
        CALL pyshow(ip1,ip2,sqrt(max(0d0,pm12s)))
      ENDIF
      naft1=n
      IF(iql34.EQ.1) THEN
        ijoin(1)=ip3
        ijoin(2)=ip4
        CALL pyjoin(2,ijoin)
      ENDIF
      IF(iql34.EQ.1.OR.irad.EQ.1) THEN
        pm34s=(p(ip3,4)+p(ip4,4))**2-(p(ip3,1)+p(ip4,1))**2-
     &  (p(ip3,2)+p(ip4,2))**2-(p(ip3,3)+p(ip4,3))**2
        CALL pyshow(ip3,ip4,sqrt(max(0d0,pm34s)))
      ENDIF
 
C...Optionally do colour reconnection.
      mint(32)=0
      msti(32)=0
      IF(mstp(115).GE.1.AND.iql12.EQ.1.AND.iql34.EQ.1) THEN
        CALL pyreco(iw1,iw2,nsd1,naft1)
        msti(32)=mint(32)
      ENDIF
 
C...Do fragmentation and decays. Possibly except tau decay.
      IF(itau.EQ.0) THEN
        ntau=0
        DO 120 i=1,n
        IF(iabs(k(i,2)).EQ.15.AND.k(i,1).EQ.1) THEN
          ntau=ntau+1
          intau(ntau)=i
          k(i,1)=11
        ENDIF
  120   CONTINUE
      ENDIF
      CALL pyexec
      IF(itau.EQ.0) THEN
        DO 130 i=1,ntau
        k(intau(i),1)=1
  130   CONTINUE
      ENDIF
 
C...Call PYHEPC to convert output from PYJETS to HEPEVT common.
      IF(icom.EQ.0) THEN
        mstu(28)=0
        CALL pyhepc(1)
      ENDIF
 
      END