py4jet.f

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

 
C*********************************************************************
 
C...PY4JET
C...An interface from a four-parton generator to include
C...parton showers and hadronization.
 
      SUBROUTINE py4jet(PMAX,IRAD,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),ptot(4),beta(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 partons.
      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.EQ.21) THEN
        IF(k(i,2).GT.0.AND.k(i,2).LE.6) THEN
          IF(i1.EQ.0) THEN
            i1=i
          ELSEIF(i3.EQ.0) THEN
            i3=i
          ELSE
            CALL pyerrm(16,'(PY4JET:) more than two quarks')
          ENDIF
        ELSEIF(k(i,2).LT.0) THEN
          IF(i2.EQ.0) THEN
            i2=i
          ELSEIF(i4.EQ.0) THEN
            i4=i
          ELSE
            CALL pyerrm(16,'(PY4JET:) more than two antiquarks')
          ENDIF
        ELSE
          IF(i3.EQ.0) THEN
            i3=i
          ELSEIF(i4.EQ.0) THEN
            i4=i
          ELSE
            CALL pyerrm(16,'(PY4JET:) more than two gluons')
          ENDIF
        ENDIF
      ENDIF
  100 CONTINUE
 
C...Check that event is arranged according to conventions.
      IF(i1.EQ.0.OR.i2.EQ.0.OR.i3.EQ.0.OR.i4.EQ.0) THEN
        CALL pyerrm(16,'(PY4JET:) event contains too few partons')
      ENDIF
      IF(i2.LT.i1.OR.i3.LT.i2.OR.i4.LT.i3) THEN
        CALL pyerrm(6,'(PY4JET:) partons arranged in wrong order')
      ENDIF
 
C...Check whether second pair are quarks or gluons.
      IF(iabs(k(i3,2)).LT.10.AND.iabs(k(i4,2)).LT.10) THEN
        iqg34=1
      ELSEIF(k(i3,2).EQ.21.AND.k(i4,2).EQ.21) THEN
        iqg34=2
      ELSE
        CALL pyerrm(16,'(PY4JET:) second parton pair inconsistent')
      ENDIF
 
C...Boost partons to their cm frame.
      DO 110 j=1,4
        ptot(j)=p(i1,j)+p(i2,j)+p(i3,j)+p(i4,j)
  110 CONTINUE
      ecm=sqrt(max(0d0,ptot(4)**2-ptot(1)**2-ptot(2)**2-ptot(3)**2))
      DO 120 j=1,3
        beta(j)=ptot(j)/ptot(4)
  120 CONTINUE
      CALL pyrobo(i1,i1,0d0,0d0,-beta(1),-beta(2),-beta(3))
      CALL pyrobo(i2,i2,0d0,0d0,-beta(1),-beta(2),-beta(3))
      CALL pyrobo(i3,i3,0d0,0d0,-beta(1),-beta(2),-beta(3))
      CALL pyrobo(i4,i4,0d0,0d0,-beta(1),-beta(2),-beta(3))
      nsav=n
 
C...Decide and set up shower history for q qbar q' qbar' events.
      IF(iqg34.EQ.1) THEN
        w1=py4jtw(0,i1,i3,i4)
        w2=py4jtw(0,i2,i3,i4)
        IF(w1.GT.pyr(0)*(w1+w2)) THEN
          CALL py4jts(0,i1,i3,i4,i2,qmax)
        ELSE
          CALL py4jts(0,i2,i3,i4,i1,qmax)
        ENDIF
 
C...Decide and set up shower history for q qbar g g events.
      ELSE
        w1=py4jtw(i1,i3,i2,i4)
        w2=py4jtw(i1,i4,i2,i3)
        w3=py4jtw(0,i3,i1,i4)
        w4=py4jtw(0,i4,i1,i3)
        w5=py4jtw(0,i3,i2,i4)
        w6=py4jtw(0,i4,i2,i3)
        w7=py4jtw(0,i1,i3,i4)
        w8=py4jtw(0,i2,i3,i4)
        wr=(w1+w2+w3+w4+w5+w6+w7+w8)*pyr(0)
        IF(w1.GT.wr) THEN
          CALL py4jts(i1,i3,i2,i4,0,qmax)
        ELSEIF(w1+w2.GT.wr) THEN
          CALL py4jts(i1,i4,i2,i3,0,qmax)
        ELSEIF(w1+w2+w3.GT.wr) THEN
          CALL py4jts(0,i3,i1,i4,i2,qmax)
        ELSEIF(w1+w2+w3+w4.GT.wr) THEN
          CALL py4jts(0,i4,i1,i3,i2,qmax)
        ELSEIF(w1+w2+w3+w4+w5.GT.wr) THEN
          CALL py4jts(0,i3,i2,i4,i1,qmax)
        ELSEIF(w1+w2+w3+w4+w5+w6.GT.wr) THEN
          CALL py4jts(0,i4,i2,i3,i1,qmax)
        ELSEIF(w1+w2+w3+w4+w5+w6+w7.GT.wr) THEN
          CALL py4jts(0,i1,i3,i4,i2,qmax)
        ELSE
          CALL py4jts(0,i2,i3,i4,i1,qmax)
        ENDIF
      ENDIF
 
C...Boost back original partons and mark them as deleted.
      CALL pyrobo(i1,i1,0d0,0d0,beta(1),beta(2),beta(3))
      CALL pyrobo(i2,i2,0d0,0d0,beta(1),beta(2),beta(3))
      CALL pyrobo(i3,i3,0d0,0d0,beta(1),beta(2),beta(3))
      CALL pyrobo(i4,i4,0d0,0d0,beta(1),beta(2),beta(3))
      k(i1,1)=k(i1,1)+10
      k(i2,1)=k(i2,1)+10
      k(i3,1)=k(i3,1)+10
      k(i4,1)=k(i4,1)+10
 
C...Rotate shower initiating partons to be along z axis.
      phi=pyangl(p(nsav+1,1),p(nsav+1,2))
      CALL pyrobo(nsav+1,nsav+6,0d0,-phi,0d0,0d0,0d0)
      the=pyangl(p(nsav+1,3),p(nsav+1,1))
      CALL pyrobo(nsav+1,nsav+6,-the,0d0,0d0,0d0,0d0)
 
C...Set up copy of shower initiating partons as on mass shell.
      DO 140 i=n+1,n+2
        DO 130 j=1,5
          k(i,j)=0
          p(i,j)=0d0
          v(i,j)=v(i1,j)
  130   CONTINUE
        k(i,1)=1
        k(i,2)=k(i-6,2)
  140 CONTINUE
      IF(k(nsav+1,2).EQ.k(i1,2)) THEN
        k(n+1,3)=i1
        p(n+1,5)=p(i1,5)
        k(n+2,3)=i2
        p(n+2,5)=p(i2,5)
      ELSE
        k(n+1,3)=i2
        p(n+1,5)=p(i2,5)
        k(n+2,3)=i1
        p(n+2,5)=p(i1,5)
      ENDIF
      pabs=sqrt(max(0d0,(ecm**2-p(n+1,5)**2-p(n+2,5)**2)**2-
     &(2d0*p(n+1,5)*p(n+2,5))**2))/(2d0*ecm)
      p(n+1,3)=pabs
      p(n+1,4)=sqrt(pabs**2+p(n+1,5)**2)
      p(n+2,3)=-pabs
      p(n+2,4)=sqrt(pabs**2+p(n+2,5)**2)
      n=n+2
 
C...Decide whether to allow or not photon radiation in showers.
C...Connect up colours.
      mstj(41)=2
      IF(irad.EQ.0) mstj(41)=1
      ijoin(1)=n-1
      ijoin(2)=n
      CALL pyjoin(2,ijoin)
 
C...Decide on maximum virtuality and do parton shower.
      IF(pmax.LT.parj(82)) THEN
        pqmax=qmax
      ELSE
        pqmax=pmax
      ENDIF
      CALL pyshow(nsav+1,-100,pqmax)
 
C...Rotate and boost back system.
      CALL pyrobo(nsav+1,n,the,phi,beta(1),beta(2),beta(3))
 
C...Do fragmentation and decays.
      CALL pyexec
 
C...Call PYHEPC to convert output from PYJETS to HEPEVT common.
      IF(icom.EQ.0) THEN
        mstu(28)=0
        CALL pyhepc(1)
      ENDIF
 
      RETURN
      END