py4jts.f

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

 
C*********************************************************************
 
C...PY4JTS
C...Auxiliary to PY4JET, to set up chosen configuration.
 
      SUBROUTINE py4jts(IA1,IA2,IA3,IA4,IA5,QMAX)
 
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)
      SAVE /pyjets/
 
C...Reset info.
      DO 110 i=n+1,n+6
        DO 100 j=1,5
          k(i,j)=0
          v(i,j)=v(ia2,j)
  100   CONTINUE
        k(i,1)=16
  110 CONTINUE
 
C...First case: when both original partons radiate.
C...N+1 -> (IA1=N+3) + (IA2=N+4), N+2 -> (IA3=N+5) + (IA4=N+6).
      IF(ia1.NE.0) THEN
 
C...Set up flavour and history pointers for new partons.
        k(n+1,2)=k(ia1,2)
        k(n+2,2)=k(ia3,2)
        k(n+3,2)=k(ia1,2)
        k(n+4,2)=k(ia2,2)
        k(n+5,2)=k(ia3,2)
        k(n+6,2)=k(ia4,2)
        k(n+1,3)=ia1
        k(n+1,4)=n+3
        k(n+1,5)=n+4
        k(n+2,3)=ia3
        k(n+2,4)=n+5
        k(n+2,5)=n+6
        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...Set up momenta for new partons.
        DO 120 j=1,5
          p(n+1,j)=p(ia1,j)+p(ia2,j)
          p(n+2,j)=p(ia3,j)+p(ia4,j)
          p(n+3,j)=p(ia1,j)
          p(n+4,j)=p(ia2,j)
          p(n+5,j)=p(ia3,j)
          p(n+6,j)=p(ia4,j)
  120   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))
        qmax=min(p(n+1,5),p(n+2,5))
 
C...Second case: q radiates twice.
C...N+1 -> (IA2=N+4) + N+3, N+3 -> (IA3=N+5) + (IA4=N+6),
C...IA5=N+2 does not radiate.
      ELSEIF(k(ia2,2).EQ.21) THEN
 
C...Set up flavour and history pointers for new partons.
        k(n+1,2)=k(ia3,2)
        k(n+2,2)=k(ia5,2)
        k(n+3,2)=k(ia3,2)
        k(n+4,2)=k(ia2,2)
        k(n+5,2)=k(ia3,2)
        k(n+6,2)=k(ia4,2)
        k(n+1,3)=ia3
        k(n+1,4)=n+3
        k(n+1,5)=n+4
        k(n+2,3)=ia5
        k(n+3,3)=n+1
        k(n+3,4)=n+5
        k(n+3,5)=n+6
        k(n+4,3)=n+1
        k(n+5,3)=n+3
        k(n+6,3)=n+3
 
C...Set up momenta for new partons.
        DO 130 j=1,5
          p(n+1,j)=p(ia2,j)+p(ia3,j)+p(ia4,j)
          p(n+2,j)=p(ia5,j)
          p(n+3,j)=p(ia3,j)+p(ia4,j)
          p(n+4,j)=p(ia2,j)
          p(n+5,j)=p(ia3,j)
          p(n+6,j)=p(ia4,j)
  130   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+3,5)=sqrt(max(0d0,p(n+3,4)**2-p(n+3,1)**2-p(n+3,2)**2-
     &  p(n+3,3)**2))
        qmax=p(n+3,5)
 
C...Third case: q radiates g, g branches.
C...N+1 -> (IA2=N+3) + N+4, N+4 -> (IA3=N+5) + (IA4=N+6),
C...IA5=N+2 does not radiate.
      ELSE
 
C...Set up flavour and history pointers for new partons.
        k(n+1,2)=k(ia2,2)
        k(n+2,2)=k(ia5,2)
        k(n+3,2)=k(ia2,2)
        k(n+4,2)=21
        k(n+5,2)=k(ia3,2)
        k(n+6,2)=k(ia4,2)
        k(n+1,3)=ia2
        k(n+1,4)=n+3
        k(n+1,5)=n+4
        k(n+2,3)=ia5
        k(n+3,3)=n+1
        k(n+4,3)=n+1
        k(n+4,4)=n+5
        k(n+4,5)=n+6
        k(n+5,3)=n+4
        k(n+6,3)=n+4
 
C...Set up momenta for new partons.
        DO 140 j=1,5
          p(n+1,j)=p(ia2,j)+p(ia3,j)+p(ia4,j)
          p(n+2,j)=p(ia5,j)
          p(n+3,j)=p(ia2,j)
          p(n+4,j)=p(ia3,j)+p(ia4,j)
          p(n+5,j)=p(ia3,j)
          p(n+6,j)=p(ia4,j)
  140   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+4,5)=sqrt(max(0d0,p(n+4,4)**2-p(n+4,1)**2-p(n+4,2)**2-
     &  p(n+4,3)**2))
        qmax=p(n+4,5)
 
      ENDIF
      n=n+6
 
      RETURN
      END