py2frm.f

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

 
C*********************************************************************
 
C...PY2FRM
C...An interface from a two-fermion generator to include
C...parton showers and hadronization.
 
      SUBROUTINE py2frm(IRAD,ITAU,ICOM)
 
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)
      SAVE /pyjets/,/pydat1/
C...Local arrays.
      dimension ijoin(2),intau(2)
 
C...Call PYHEPC to convert input from HEPEVT to PYJETS common.
      IF(icom.EQ.0) THEN
        mstu(28)=0
        CALL pyhepc(2)
      ENDIF
 
C...Loop through entries and pick up all final fermions/antifermions.
      i1=0
      i2=0
      DO 100 i=1,n
      IF(k(i,1).LE.0.OR.k(i,1).GT.10) goto 100
      kfa=iabs(k(i,2))
      IF((kfa.GE.1.AND.kfa.LE.6).OR.(kfa.GE.11.AND.kfa.LE.16)) THEN
        IF(k(i,2).GT.0) THEN
          IF(i1.EQ.0) THEN
            i1=i
          ELSE
            CALL pyerrm(16,'(PY2FRM:) more than one fermion')
          ENDIF
        ELSE
          IF(i2.EQ.0) THEN
            i2=i
          ELSE
            CALL pyerrm(16,'(PY2FRM:) more than one antifermion')
          ENDIF
        ENDIF
      ENDIF
  100 CONTINUE
 
C...Check that event is arranged according to conventions.
      IF(i1.EQ.0.OR.i2.EQ.0) THEN
        CALL pyerrm(16,'(PY2FRM:) event contains too few fermions')
      ENDIF
      IF(i2.LT.i1) THEN
        CALL pyerrm(6,'(PY2FRM:) fermions arranged in wrong order')
      ENDIF
 
C...Check whether fermion pair is quarks or leptons.
      IF(iabs(k(i1,2)).LT.10.AND.iabs(k(i2,2)).LT.10) THEN
        iql12=1
      ELSEIF(iabs(k(i1,2)).GT.10.AND.iabs(k(i2,2)).GT.10) THEN
        iql12=2
      ELSE
        CALL pyerrm(16,'(PY2FRM:) fermion pair inconsistent')
      ENDIF
 
C...Decide whether to allow or not photon radiation in showers.
      mstj(41)=2
      IF(irad.EQ.0) mstj(41)=1
 
C...Do colour joining and parton showers.
      ip1=i1
      ip2=i2
      IF(iql12.EQ.1) THEN
        ijoin(1)=ip1
        ijoin(2)=ip2
        CALL pyjoin(2,ijoin)
      ENDIF
      IF(iql12.EQ.1.OR.irad.EQ.1) THEN
        pm12s=(p(ip1,4)+p(ip2,4))**2-(p(ip1,1)+p(ip2,1))**2-
     &  (p(ip1,2)+p(ip2,2))**2-(p(ip1,3)+p(ip2,3))**2
        CALL pyshow(ip1,ip2,sqrt(max(0d0,pm12s)))
      ENDIF
 
C...Do fragmentation and decays. Possibly except tau decay.
      IF(itau.EQ.0) THEN
        ntau=0
        DO 110 i=1,n
        IF(iabs(k(i,2)).EQ.15.AND.k(i,1).EQ.1) THEN
          ntau=ntau+1
          intau(ntau)=i
          k(i,1)=11
        ENDIF
  110   CONTINUE
      ENDIF
      CALL pyexec
      IF(itau.EQ.0) THEN
        DO 120 i=1,ntau
        k(intau(i),1)=1
  120   CONTINUE
      ENDIF
 
C...Call PYHEPC to convert output from PYJETS to HEPEVT common.
      IF(icom.EQ.0) THEN
        mstu(28)=0
        CALL pyhepc(1)
      ENDIF
 
      END