py6frm.f

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

 
C*********************************************************************
 
C...PY6FRM
C...An interface from a six-fermion generator to include
C...parton showers and hadronization.
 
      SUBROUTINE py6frm(P12,P13,P21,P23,P31,P32,PTOP,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)
      SAVE /pyjets/,/pydat1/
C...Local arrays.
      dimension ijoin(2),intau(6),beta(3),betao(3),betan(3)
 
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
      i5=0
      i6=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
          ELSEIF(i5.EQ.0) THEN
            i5=i
          ELSE
            CALL pyerrm(16,'(PY6FRM:) more than three fermions')
          ENDIF
        ELSE
          IF(i2.EQ.0) THEN
            i2=i
          ELSEIF(i4.EQ.0) THEN
            i4=i
          ELSEIF(i6.EQ.0) THEN
            i6=i
          ELSE
            CALL pyerrm(16,'(PY6FRM:) more than three antifermions')
          ENDIF
        ENDIF
      ENDIF
  100 CONTINUE
 
C...Check that event is arranged according to conventions.
      IF(i5.EQ.0.OR.i6.EQ.0) THEN
        CALL pyerrm(16,'(PY6FRM:) event contains too few fermions')
      ENDIF
      IF(i2.LT.i1.OR.i3.LT.i2.OR.i4.LT.i3.OR.i5.LT.i4.OR.i6.LT.i5) THEN
        CALL pyerrm(6,'(PY6FRM:) 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,'(PY6FRM:) 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,'(PY6FRM:) second fermion pair inconsistent')
      ENDIF
      IF(iabs(k(i5,2)).LT.10.AND.iabs(k(i6,2)).LT.10) THEN
        iql56=1
      ELSEIF(iabs(k(i5,2)).GT.10.AND.iabs(k(i6,2)).GT.10) THEN
        iql56=2
      ELSE
        CALL pyerrm(16,'(PY6FRM:) third 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...Allow dipole pairings only among leptons and quarks separately.
      p12d=p12
      p13d=0d0
      IF(iql34.EQ.iql56) p13d=p13
      p21d=0d0
      IF(iql12.EQ.iql34) p21d=p21
      p23d=0d0
      IF(iql12.EQ.iql34.AND.iql12.EQ.iql56) p23d=p23
      p31d=0d0
      IF(iql12.EQ.iql34.AND.iql12.EQ.iql56) p31d=p31
      p32d=0d0
      IF(iql12.EQ.iql56) p32d=p32
 
C...Decide whether t+tbar.
      itop=0
      IF(pyr(0).LT.ptop) THEN
        itop=1
 
C...If t+tbar: reconstruct t's.
        it=n+1
        itb=n+2
        DO 110 j=1,5
          k(it,j)=0
          k(itb,j)=0
          p(it,j)=p(i1,j)+p(i3,j)+p(i4,j)
          p(itb,j)=p(i2,j)+p(i5,j)+p(i6,j)
          v(it,j)=0d0
          v(itb,j)=0d0
  110   CONTINUE
        k(it,1)=1
        k(itb,1)=1
        k(it,2)=6
        k(itb,2)=-6
        p(it,5)=sqrt(max(0d0,p(it,4)**2-p(it,1)**2-p(it,2)**2-
     &  p(it,3)**2))
        p(itb,5)=sqrt(max(0d0,p(itb,4)**2-p(itb,1)**2-p(itb,2)**2-
     &  p(itb,3)**2))
        n=n+2
 
C...If t+tbar: colour join t's and let them shower.
        ijoin(1)=it
        ijoin(2)=itb
        CALL pyjoin(2,ijoin)
        pmtts=(p(it,4)+p(itb,4))**2-(p(it,1)+p(itb,1))**2-
     &  (p(it,2)+p(itb,2))**2-(p(it,3)+p(itb,3))**2
        CALL pyshow(it,itb,sqrt(max(0d0,pmtts)))
 
C...If t+tbar: pick up the t's after shower.
        itnew=it
        itbnew=itb
        DO 120 i=itb+1,n
          IF(k(i,2).EQ.6) itnew=i
          IF(k(i,2).EQ.-6) itbnew=i
  120   CONTINUE
 
C...If t+tbar: loop over two top systems.
        DO 200 it1=1,2
          IF(it1.EQ.1) THEN
            ito=it
            itn=itnew
            ibo=i1
            iw1=i3
            iw2=i4
          ELSE
            ito=itb
            itn=itbnew
            ibo=i2
            iw1=i5
            iw2=i6
          ENDIF
          IF(iabs(k(ibo,2)).NE.5) CALL pyerrm(6,
     &    '(PY6FRM:) not b in t decay')
 
C...If t+tbar: find boost from original to new top frame.
          DO 130 j=1,3
            betao(j)=p(ito,j)/p(ito,4)
            betan(j)=p(itn,j)/p(itn,4)
  130     CONTINUE
 
C...If t+tbar: boost copy of b by t shower and connect it in colour.
          n=n+1
          ib=n
          k(ib,1)=3
          k(ib,2)=k(ibo,2)
          k(ib,3)=itn
          DO 140 j=1,5
            p(ib,j)=p(ibo,j)
            v(ib,j)=0d0
  140     CONTINUE
          CALL pyrobo(ib,ib,0d0,0d0,-betao(1),-betao(2),-betao(3))
          CALL pyrobo(ib,ib,0d0,0d0,betan(1),betan(2),betan(3))
          k(ib,4)=mstu(5)*itn
          k(ib,5)=mstu(5)*itn
          k(itn,4)=k(itn,4)+ib
          k(itn,5)=k(itn,5)+ib
          k(itn,1)=k(itn,1)+10
          k(ibo,1)=k(ibo,1)+10
 
C...If t+tbar: construct W recoiling against b.
          n=n+1
          iw=n
          DO 150 j=1,5
            k(iw,j)=0
            v(iw,j)=0d0
  150     CONTINUE
          k(iw,1)=1
          kchw=pychge(k(iw1,2))+pychge(k(iw2,2))
          IF(iabs(kchw).EQ.3) THEN
            k(iw,2)=isign(24,kchw)
          ELSE
            CALL pyerrm(16,'(PY6FRM:) fermion pair inconsistent with W')
          ENDIF
          k(iw,3)=iw1
 
C...If t+tbar: construct W momentum, including boost by t shower.
          DO 160 j=1,4
            p(iw,j)=p(iw1,j)+p(iw2,j)
  160     CONTINUE
          p(iw,5)=sqrt(max(0d0,p(iw,4)**2-p(iw,1)**2-p(iw,2)**2-
     &    p(iw,3)**2))
          CALL pyrobo(iw,iw,0d0,0d0,-betao(1),-betao(2),-betao(3))
          CALL pyrobo(iw,iw,0d0,0d0,betan(1),betan(2),betan(3))
 
C...If t+tbar: boost b and W to top rest frame.
          DO 170 j=1,3
            beta(j)=(p(ib,j)+p(iw,j))/(p(ib,4)+p(iw,4))
  170     CONTINUE
          CALL pyrobo(ib,ib,0d0,0d0,-beta(1),-beta(2),-beta(3))
          CALL pyrobo(iw,iw,0d0,0d0,-beta(1),-beta(2),-beta(3))
 
C...If t+tbar: let b shower and pick up modified W.
          pmts=(p(ib,4)+p(iw,4))**2-(p(ib,1)+p(iw,1))**2-
     &    (p(ib,2)+p(iw,2))**2-(p(ib,3)+p(iw,3))**2
          CALL pyshow(ib,iw,sqrt(max(0d0,pmts)))
          DO 180 i=iw,n
            IF(iabs(k(i,2)).EQ.24) iwm=i
  180     CONTINUE
 
C...If t+tbar: take copy of W decay products.
          DO 190 j=1,5
            k(n+1,j)=k(iw1,j)
            p(n+1,j)=p(iw1,j)
            v(n+1,j)=v(iw1,j)
            k(n+2,j)=k(iw2,j)
            p(n+2,j)=p(iw2,j)
            v(n+2,j)=v(iw2,j)
  190     CONTINUE
          k(iw1,1)=k(iw1,1)+10
          k(iw2,1)=k(iw2,1)+10
          k(iwm,1)=k(iwm,1)+10
          k(iwm,4)=n+1
          k(iwm,5)=n+2
          k(n+1,3)=iwm
          k(n+2,3)=iwm
          IF(it1.EQ.1) THEN
            i3=n+1
            i4=n+2
          ELSE
            i5=n+1
            i6=n+2
          ENDIF
          n=n+2
 
C...If t+tbar: boost W decay products, first by effects of t shower,
C...then by those of b shower. b and its shower simple boost back.
          CALL pyrobo(n-1,n,0d0,0d0,-betao(1),-betao(2),-betao(3))
          CALL pyrobo(n-1,n,0d0,0d0,betan(1),betan(2),betan(3))
          CALL pyrobo(n-1,n,0d0,0d0,-beta(1),-beta(2),-beta(3))
          CALL pyrobo(n-1,n,0d0,0d0,-p(iw,1)/p(iw,4),
     &    -p(iw,2)/p(iw,4),-p(iw,3)/p(iw,4))
          CALL pyrobo(n-1,n,0d0,0d0,p(iwm,1)/p(iwm,4),
     &    p(iwm,2)/p(iwm,4),p(iwm,3)/p(iwm,4))
          CALL pyrobo(ib,ib,0d0,0d0,beta(1),beta(2),beta(3))
          CALL pyrobo(iw,n,0d0,0d0,beta(1),beta(2),beta(3))
  200   CONTINUE
      ENDIF
 
C...Decide on dipole pairing.
      ip1=i1
      ip3=i3
      ip5=i5
      prn=pyr(0)*(p12d+p13d+p21d+p23d+p31d+p32d)
      IF(itop.EQ.1.OR.prn.LT.p12d) THEN
        ip2=i2
        ip4=i4
        ip6=i6
      ELSEIF(prn.LT.p12d+p13d) THEN
        ip2=i2
        ip4=i6
        ip6=i4
      ELSEIF(prn.LT.p12d+p13d+p21d) THEN
        ip2=i4
        ip4=i2
        ip6=i6
      ELSEIF(prn.LT.p12d+p13d+p21d+p23d) THEN
        ip2=i4
        ip4=i6
        ip6=i2
      ELSEIF(prn.LT.p12d+p13d+p21d+p23d+p31d) THEN
        ip2=i6
        ip4=i2
        ip6=i4
      ELSE
        ip2=i6
        ip4=i4
        ip6=i2
      ENDIF
 
C...Do colour joinings and parton showers
C...(except ones already made for t+tbar).
      IF(itop.EQ.0) THEN
        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
      ENDIF
      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
      IF(iql56.EQ.1) THEN
        ijoin(1)=ip5
        ijoin(2)=ip6
        CALL pyjoin(2,ijoin)
      ENDIF
      IF(iql56.EQ.1.OR.irad.EQ.1) THEN
        pm56s=(p(ip5,4)+p(ip6,4))**2-(p(ip5,1)+p(ip6,1))**2-
     &  (p(ip5,2)+p(ip6,2))**2-(p(ip5,3)+p(ip6,3))**2
        CALL pyshow(ip5,ip6,sqrt(max(0d0,pm56s)))
      ENDIF
 
C...Do fragmentation and decays. Possibly except tau decay.
      IF(itau.EQ.0) THEN
        ntau=0
        DO 210 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
  210   CONTINUE
      ENDIF
      CALL pyexec
      IF(itau.EQ.0) THEN
        DO 220 i=1,ntau
        k(intau(i),1)=1
  220   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