py4jtw.f

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

 
C*********************************************************************
 
C...PY4JTW
C...Auxiliary to PY4JET, to evaluate weight of configuration.
 
      FUNCTION py4jtw(IA1,IA2,IA3,IA4)
 
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...First case: when both original partons radiate.
C...IA1 /= 0: N+1 -> IA1 + IA2, N+2 -> IA3 + IA4.
      IF(ia1.NE.0) THEN
        DO 100 j=1,4
          p(n+1,j)=p(ia1,j)+p(ia2,j)
          p(n+2,j)=p(ia3,j)+p(ia4,j)
  100   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))
        z1=p(ia1,4)/p(n+1,4)
        wt1=(4d0/3d0)*((1d0+z1**2)/(1d0-z1))/(p(n+1,5)**2-p(ia1,5)**2)
        z2=p(ia3,4)/p(n+2,4)
        wt2=(4d0/3d0)*((1d0+z2**2)/(1d0-z2))/(p(n+2,5)**2-p(ia3,5)**2)
 
C...Second case: when one original parton radiates to three.
C...IA1  = 0: N+1 -> IA2 + N+2, N+2 -> IA3 + IA4.
      ELSE
        DO 110 j=1,4
          p(n+2,j)=p(ia3,j)+p(ia4,j)
          p(n+1,j)=p(n+2,j)+p(ia2,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))
        IF(k(ia2,2).EQ.21) THEN
          z1=p(n+2,4)/p(n+1,4)
          wt1=(4d0/3d0)*((1d0+z1**2)/(1d0-z1))/(p(n+1,5)**2-
     &    p(ia3,5)**2)
        ELSE
          z1=p(ia2,4)/p(n+1,4)
          wt1=(4d0/3d0)*((1d0+z1**2)/(1d0-z1))/(p(n+1,5)**2-
     &    p(ia2,5)**2)
        ENDIF
        z2=p(ia3,4)/p(n+2,4)
        IF(k(ia2,2).EQ.21) THEN
          wt2=(4d0/3d0)*((1d0+z2**2)/(1d0-z2))/(p(n+2,5)**2-
     &    p(ia3,5)**2)
        ELSEIF(k(ia3,2).EQ.21) THEN
          wt2=3d0*((1d0-z2*(1d0-z2))**2/(z2*(1d0-z2)))/p(n+2,5)**2
        ELSE
          wt2=0.5d0*(z2**2+(1d0-z2)**2)
        ENDIF
      ENDIF
 
C...Total weight.
      py4jtw=wt1*wt2
 
      RETURN
      END