py2ent.f

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

 
C*********************************************************************
 
C...PY2ENT
C...Stores two partons/particles in their CM frame,
C...with the first along the +z axis.
 
      SUBROUTINE py2ent(IP,KF1,KF2,PECM)
 
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)-1) CALL pyerrm(21,
     &'(PY2ENT:) writing outside PYJETS memory')
      kc1=pycomp(kf1)
      kc2=pycomp(kf2)
      IF(kc1.EQ.0.OR.kc2.EQ.0) CALL pyerrm(12,
     &'(PY2ENT:) 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)
      DO 110 i=ipa,ipa+1
        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)
      IF(mstu(19).EQ.1) THEN
        mstu(19)=0
      ELSE
        IF(kq1+kq2.NE.0.AND.kq1+kq2.NE.4) CALL pyerrm(2,
     &  '(PY2ENT:) unphysical flavour combination')
      ENDIF
      k(ipa,2)=kf1
      k(ipa+1,2)=kf2
 
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) k(ipa,1)=2
        k(ipa+1,1)=1
 
C...Store partons in K vectors for parton shower evolution.
      ELSE
        k(ipa,1)=3
        k(ipa+1,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)
      ENDIF
 
C...Check kinematics and store partons/particles in P vectors.
      IF(pecm.LE.pm1+pm2) CALL pyerrm(13,
     &'(PY2ENT:) energy smaller than sum of masses')
      pa=sqrt(max(0d0,(pecm**2-pm1**2-pm2**2)**2-(2d0*pm1*pm2)**2))/
     &(2d0*pecm)
      p(ipa,3)=pa
      p(ipa,4)=sqrt(pm1**2+pa**2)
      p(ipa,5)=pm1
      p(ipa+1,3)=-pa
      p(ipa+1,4)=sqrt(pm2**2+pa**2)
      p(ipa+1,5)=pm2
 
C...Set N. Optionally fragment/decay.
      n=ipa+1
      IF(ip.EQ.0) CALL pyexec
 
      RETURN
      END