py4ent.f

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

 
C*********************************************************************
 
C...PY4ENT
C...Stores four partons or particles in their CM frame, with
C...the first along the +z axis, the last in the xz plane with x > 0
C...and the second having y < 0 and y > 0 with equal probability.
 
      SUBROUTINE py4ent(IP,KF1,KF2,KF3,KF4,PECM,X1,X2,X4,X12,X14)
 
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)-3) CALL pyerrm(21,
     &'(PY4ENT:) writing outside PYJETS momory')
      kc1=pycomp(kf1)
      kc2=pycomp(kf2)
      kc3=pycomp(kf3)
      kc4=pycomp(kf4)
      IF(kc1.EQ.0.OR.kc2.EQ.0.OR.kc3.EQ.0.OR.kc4.EQ.0) CALL pyerrm(12,
     &'(PY4ENT:) 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)
      pm4=0d0
      IF(mstu(10).EQ.1) pm4=p(ipa+3,5)
      IF(mstu(10).GE.2) pm4=pymass(kf4)
      DO 110 i=ipa,ipa+3
        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)
      kq4=kchg(kc4,2)*isign(1,kf4)
      IF(mstu(19).EQ.1) THEN
        mstu(19)=0
      ELSEIF(kq1.EQ.0.AND.kq2.EQ.0.AND.kq3.EQ.0.AND.kq4.EQ.0) THEN
      ELSEIF(kq1.NE.0.AND.kq2.EQ.2.AND.kq3.EQ.2.AND.(kq1+kq4.EQ.0.OR.
     &  kq1+kq4.EQ.4)) THEN
      ELSEIF(kq1.NE.0.AND.kq1+kq2.EQ.0.AND.kq3.NE.0.AND.kq3+kq4.EQ.0d0)
     &  THEN
      ELSE
        CALL pyerrm(2,'(PY4ENT:) unphysical flavour combination')
      ENDIF
      k(ipa,2)=kf1
      k(ipa+1,2)=kf2
      k(ipa+2,2)=kf3
      k(ipa+3,2)=kf4
 
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.OR.kq4.NE.0)) k(ipa,1)=2
        k(ipa+1,1)=1
        IF(kq2.NE.0.AND.kq1+kq2.NE.0.AND.(kq3.NE.0.OR.kq4.NE.0))
     &  k(ipa+1,1)=2
        k(ipa+2,1)=1
        IF(kq3.NE.0.AND.kq4.NE.0) k(ipa+2,1)=2
        k(ipa+3,1)=1
 
C...Store partons for parton shower evolution from q-g-g-qbar or
C...g-g-g-g event.
      ELSEIF(kq1+kq2.NE.0) THEN
        k(ipa,1)=3
        k(ipa+1,1)=3
        k(ipa+2,1)=3
        k(ipa+3,1)=3
        kcs=4
        IF(kq1.EQ.-1) kcs=5
        k(ipa,kcs)=mstu(5)*(ipa+1)
        k(ipa,9-kcs)=mstu(5)*(ipa+3)
        k(ipa+1,kcs)=mstu(5)*(ipa+2)
        k(ipa+1,9-kcs)=mstu(5)*ipa
        k(ipa+2,kcs)=mstu(5)*(ipa+3)
        k(ipa+2,9-kcs)=mstu(5)*(ipa+1)
        k(ipa+3,kcs)=mstu(5)*ipa
        k(ipa+3,9-kcs)=mstu(5)*(ipa+2)
 
C...Store partons for parton shower evolution from q-qbar-q-qbar event.
      ELSE
        k(ipa,1)=3
        k(ipa+1,1)=3
        k(ipa+2,1)=3
        k(ipa+3,1)=3
        k(ipa,4)=mstu(5)*(ipa+1)
        k(ipa,5)=k(ipa,4)
        k(ipa+1,4)=mstu(5)*ipa
        k(ipa+1,5)=k(ipa+1,4)
        k(ipa+2,4)=mstu(5)*(ipa+3)
        k(ipa+2,5)=k(ipa+2,4)
        k(ipa+3,4)=mstu(5)*(ipa+2)
        k(ipa+3,5)=k(ipa+3,4)
      ENDIF
 
C...Check kinematics.
      mkerr=0
      IF(0.5d0*x1*pecm.LE.pm1.OR.0.5d0*x2*pecm.LE.pm2.OR.
     &0.5d0*(2d0-x1-x2-x4)*pecm.LE.pm3.OR.0.5d0*x4*pecm.LE.pm4)
     &mkerr=1
      pa1=sqrt(max(1d-10,(0.5d0*x1*pecm)**2-pm1**2))
      pa2=sqrt(max(1d-10,(0.5d0*x2*pecm)**2-pm2**2))
      pa4=sqrt(max(1d-10,(0.5d0*x4*pecm)**2-pm4**2))
      x24=x1+x2+x4-1d0-x12-x14+(pm3**2-pm1**2-pm2**2-pm4**2)/pecm**2
      cthe4=(x1*x4-2d0*x14)*pecm**2/(4d0*pa1*pa4)
      IF(abs(cthe4).GE.1.002d0) mkerr=1
      cthe4=max(-1d0,min(1d0,cthe4))
      sthe4=sqrt(1d0-cthe4**2)
      cthe2=(x1*x2-2d0*x12)*pecm**2/(4d0*pa1*pa2)
      IF(abs(cthe2).GE.1.002d0) mkerr=1
      cthe2=max(-1d0,min(1d0,cthe2))
      sthe2=sqrt(1d0-cthe2**2)
      cphi2=((x2*x4-2d0*x24)*pecm**2-4d0*pa2*cthe2*pa4*cthe4)/
     &max(1d-8*pecm**2,4d0*pa2*sthe2*pa4*sthe4)
      IF(abs(cphi2).GE.1.05d0) mkerr=1
      cphi2=max(-1d0,min(1d0,cphi2))
      IF(mkerr.EQ.1) CALL pyerrm(13,
     &'(PY4ENT:) 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+3,1)=pa4*sthe4
      p(ipa+3,3)=pa4*cthe4
      p(ipa+3,4)=sqrt(pa4**2+pm4**2)
      p(ipa+3,5)=pm4
      p(ipa+1,1)=pa2*sthe2*cphi2
      p(ipa+1,2)=pa2*sthe2*sqrt(1d0-cphi2**2)*(-1d0)**int(pyr(0)+0.5d0)
      p(ipa+1,3)=pa2*cthe2
      p(ipa+1,4)=sqrt(pa2**2+pm2**2)
      p(ipa+1,5)=pm2
      p(ipa+2,1)=-p(ipa+1,1)-p(ipa+3,1)
      p(ipa+2,2)=-p(ipa+1,2)
      p(ipa+2,3)=-p(ipa,3)-p(ipa+1,3)-p(ipa+3,3)
      p(ipa+2,4)=sqrt(p(ipa+2,1)**2+p(ipa+2,2)**2+p(ipa+2,3)**2+pm3**2)
      p(ipa+2,5)=pm3
 
C...Set N. Optionally fragment/decay.
      n=ipa+3
      IF(ip.EQ.0) CALL pyexec
 
      RETURN
      END