py3ent.f

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

 
C*********************************************************************
 
C...PY3ENT
C...Stores three partons or particles in their CM frame,
C...with the first along the +z axis and the third in the (x,z)
C...plane with x > 0.
 
      SUBROUTINE py3ent(IP,KF1,KF2,KF3,PECM,X1,X3)
 
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/pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
      SAVE /pyjets/,/pydat1/,/pydat2/
 
C...Standard checks.
      mstu(28)=0
      IF(mstu(12).NE.12345) CALL pylist(0)
      ipa=max(1,iabs(ip))
      IF(ipa.GT.mstu(4)-2) CALL pyerrm(21,
     &'(PY3ENT:) writing outside PYJETS memory')
      kc1=pycomp(kf1)
      kc2=pycomp(kf2)
      kc3=pycomp(kf3)
      IF(kc1.EQ.0.OR.kc2.EQ.0.OR.kc3.EQ.0) CALL pyerrm(12,
     &'(PY3ENT:) unknown flavour code')
 
C...Find masses. Reset K, P and V vectors.
      pm1=0d0
      IF(mstu(10).EQ.1) pm1=p(ipa,5)
      IF(mstu(10).GE.2) pm1=pymass(kf1)
      pm2=0d0
      IF(mstu(10).EQ.1) pm2=p(ipa+1,5)
      IF(mstu(10).GE.2) pm2=pymass(kf2)
      pm3=0d0
      IF(mstu(10).EQ.1) pm3=p(ipa+2,5)
      IF(mstu(10).GE.2) pm3=pymass(kf3)
      DO 110 i=ipa,ipa+2
        DO 100 j=1,5
          k(i,j)=0
          p(i,j)=0d0
          v(i,j)=0d0
  100   CONTINUE
  110 CONTINUE
 
C...Check flavours.
      kq1=kchg(kc1,2)*isign(1,kf1)
      kq2=kchg(kc2,2)*isign(1,kf2)
      kq3=kchg(kc3,2)*isign(1,kf3)
      IF(mstu(19).EQ.1) THEN
        mstu(19)=0
      ELSEIF(kq1.EQ.0.AND.kq2.EQ.0.AND.kq3.EQ.0) THEN
      ELSEIF(kq1.NE.0.AND.kq2.EQ.2.AND.(kq1+kq3.EQ.0.OR.
     &  kq1+kq3.EQ.4)) THEN
      ELSE
        CALL pyerrm(2,'(PY3ENT:) unphysical flavour combination')
      ENDIF
      k(ipa,2)=kf1
      k(ipa+1,2)=kf2
      k(ipa+2,2)=kf3
 
C...Store partons/particles in K vectors for normal case.
      IF(ip.GE.0) THEN
        k(ipa,1)=1
        IF(kq1.NE.0.AND.(kq2.NE.0.OR.kq3.NE.0)) k(ipa,1)=2
        k(ipa+1,1)=1
        IF(kq2.NE.0.AND.kq3.NE.0) k(ipa+1,1)=2
        k(ipa+2,1)=1
 
C...Store partons in K vectors for parton shower evolution.
      ELSE
        k(ipa,1)=3
        k(ipa+1,1)=3
        k(ipa+2,1)=3
        kcs=4
        IF(kq1.EQ.-1) kcs=5
        k(ipa,kcs)=mstu(5)*(ipa+1)
        k(ipa,9-kcs)=mstu(5)*(ipa+2)
        k(ipa+1,kcs)=mstu(5)*(ipa+2)
        k(ipa+1,9-kcs)=mstu(5)*ipa
        k(ipa+2,kcs)=mstu(5)*ipa
        k(ipa+2,9-kcs)=mstu(5)*(ipa+1)
      ENDIF
 
C...Check kinematics.
      mkerr=0
      IF(0.5d0*x1*pecm.LE.pm1.OR.0.5d0*(2d0-x1-x3)*pecm.LE.pm2.OR.
     &0.5d0*x3*pecm.LE.pm3) mkerr=1
      pa1=sqrt(max(1d-10,(0.5d0*x1*pecm)**2-pm1**2))
      pa2=sqrt(max(1d-10,(0.5d0*(2d0-x1-x3)*pecm)**2-pm2**2))
      pa3=sqrt(max(1d-10,(0.5d0*x3*pecm)**2-pm3**2))
      cthe2=(pa3**2-pa1**2-pa2**2)/(2d0*pa1*pa2)
      cthe3=(pa2**2-pa1**2-pa3**2)/(2d0*pa1*pa3)
      IF(abs(cthe2).GE.1.001d0.OR.abs(cthe3).GE.1.001d0) mkerr=1
      cthe3=max(-1d0,min(1d0,cthe3))
      IF(mkerr.NE.0) CALL pyerrm(13,
     &'(PY3ENT:) unphysical kinematical variable setup')
 
C...Store partons/particles in P vectors.
      p(ipa,3)=pa1
      p(ipa,4)=sqrt(pa1**2+pm1**2)
      p(ipa,5)=pm1
      p(ipa+2,1)=pa3*sqrt(1d0-cthe3**2)
      p(ipa+2,3)=pa3*cthe3
      p(ipa+2,4)=sqrt(pa3**2+pm3**2)
      p(ipa+2,5)=pm3
      p(ipa+1,1)=-p(ipa+2,1)
      p(ipa+1,3)=-p(ipa,3)-p(ipa+2,3)
      p(ipa+1,4)=sqrt(p(ipa+1,1)**2+p(ipa+1,3)**2+pm2**2)
      p(ipa+1,5)=pm2
 
C...Set N. Optionally fragment/decay.
      n=ipa+2
      IF(ip.EQ.0) CALL pyexec
 
      RETURN
      END