pyexec.f

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

C*********************************************************************
 
C...PYEXEC
C...Administrates the fragmentation and decay chain.
 
      SUBROUTINE pyexec
 
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)
      common/pydat3/mdcy(500,3),mdme(8000,2),brat(8000),kfdp(8000,5)
      common/pyint1/mint(400),vint(400)
      common/pyint4/mwid(500),wids(500,5)
      SAVE /pyjets/,/pydat1/,/pydat2/,/pydat3/,/pyint1/,/pyint4/
C...Local array.
      dimension ps(2,6),ijoin(100)
 
C...Initialize and reset.
      mstu(24)=0
      IF(mstu(12).NE.12345) CALL pylist(0)
      mstu(29)=0
      mstu(31)=mstu(31)+1
      mstu(1)=0
      mstu(2)=0
      mstu(3)=0
      IF(mstu(17).LE.0) mstu(90)=0
      mcons=1
 
C...Sum up momentum, energy and charge for starting entries.
      nsav=n
      DO 110 i=1,2
        DO 100 j=1,6
          ps(i,j)=0d0
  100   CONTINUE
  110 CONTINUE
      DO 130 i=1,n
        IF(k(i,1).LE.0.OR.k(i,1).GT.10) goto 130
        DO 120 j=1,4
          ps(1,j)=ps(1,j)+p(i,j)
  120   CONTINUE
        ps(1,6)=ps(1,6)+pychge(k(i,2))
  130 CONTINUE
      paru(21)=ps(1,4)
 
C...Start by all decays of coloured resonances involved in shower.
      norig=n
      DO 140 i=1,norig
        IF(k(i,1).EQ.3) THEN
          kc=pycomp(k(i,2))
          IF(mwid(kc).NE.0.AND.kchg(kc,2).NE.0) CALL pyresd(i)
        ENDIF
  140 CONTINUE
 
C...Prepare system for subsequent fragmentation/decay.
      CALL pyprep(0)
      IF(mint(51).NE.0) RETURN
 
C...Loop through jet fragmentation and particle decays.
      mbe=0
  150 mbe=mbe+1
      ip=0
  160 ip=ip+1
      kc=0
      IF(k(ip,1).GT.0.AND.k(ip,1).LE.10) kc=pycomp(k(ip,2))
      IF(kc.EQ.0) THEN
 
C...Deal with any remaining undecayed resonance
C...(normally the task of PYEVNT, so seldom used).
      ELSEIF(mwid(kc).NE.0) THEN
        ibeg=ip
        IF(kchg(kc,2).NE.0.AND.k(i,1).NE.3) THEN
          ibeg=ip+1
  170     ibeg=ibeg-1
          IF(ibeg.GE.2.AND.k(ibeg,1).EQ.2) goto 170
          IF(k(ibeg,1).NE.2) ibeg=ibeg+1
          iend=ip-1
  180     iend=iend+1
          IF(iend.LT.n.AND.k(iend,1).EQ.2) goto 180
          IF(iend.LT.n.AND.kchg(pycomp(k(iend,2)),2).EQ.0) goto 180
          njoin=0
          DO 190 i=ibeg,iend
            IF(kchg(pycomp(k(iend,2)),2).NE.0) THEN
              njoin=njoin+1
              ijoin(njoin)=i
            ENDIF
  190     CONTINUE
        ENDIF
        CALL pyresd(ip)
        CALL pyprep(ibeg)
        IF(mint(51).NE.0) RETURN
 
C...Particle decay if unstable and allowed. Save long-lived particle
C...decays until second pass after Bose-Einstein effects.
      ELSEIF(kchg(kc,2).EQ.0) THEN
        IF(mstj(21).GE.1.AND.mdcy(kc,1).GE.1.AND.(mstj(51).LE.0.OR.mbe
     &  .EQ.2.OR.pmas(kc,2).GE.parj(91).OR.iabs(k(ip,2)).EQ.311))
     &  CALL pydecy(ip)
 
C...Decay products may develop a shower.
        IF(mstj(92).GT.0) THEN
          ip1=mstj(92)
          qmax=sqrt(max(0d0,(p(ip1,4)+p(ip1+1,4))**2-(p(ip1,1)+p(ip1+1,
     &    1))**2-(p(ip1,2)+p(ip1+1,2))**2-(p(ip1,3)+p(ip1+1,3))**2))
          mint(33)=0
          CALL pyshow(ip1,ip1+1,qmax)
          CALL pyprep(ip1)
          IF(mint(51).NE.0) RETURN
          mstj(92)=0
        ELSEIF(mstj(92).LT.0) THEN
          ip1=-mstj(92)
          mint(33)=0
          CALL pyshow(ip1,-3,p(ip,5))
          CALL pyprep(ip1)
          IF(mint(51).NE.0) RETURN
          mstj(92)=0
        ENDIF
 
C...Jet fragmentation: string or independent fragmentation.
      ELSEIF(k(ip,1).EQ.1.OR.k(ip,1).EQ.2) THEN
        mfrag=mstj(1)
        IF(mfrag.GE.1.AND.k(ip,1).EQ.1) mfrag=2
        IF(mstj(21).GE.2.AND.k(ip,1).EQ.2.AND.n.GT.ip) THEN
          IF(k(ip+1,1).EQ.1.AND.k(ip+1,3).EQ.k(ip,3).AND.
     &    k(ip,3).GT.0.AND.k(ip,3).LT.ip) THEN
            IF(kchg(pycomp(k(k(ip,3),2)),2).EQ.0) mfrag=min(1,mfrag)
          ENDIF
        ENDIF
        IF(mfrag.EQ.1) CALL pystrf(ip)
        IF(mfrag.EQ.2) CALL pyindf(ip)
        IF(mfrag.EQ.2.AND.k(ip,1).EQ.1) mcons=0
        IF(mfrag.EQ.2.AND.(mstj(3).LE.0.OR.mod(mstj(3),5).EQ.0)) mcons=0
      ENDIF
 
C...Loop back if enough space left in PYJETS and no error abort.
      IF(mstu(24).NE.0.AND.mstu(21).GE.2) THEN
      ELSEIF(ip.LT.n.AND.n.LT.mstu(4)-20-mstu(32)) THEN
        goto 160
      ELSEIF(ip.LT.n) THEN
        CALL pyerrm(11,'(PYEXEC:) no more memory left in PYJETS')
      ENDIF
 
C...Include simple Bose-Einstein effect parametrization if desired.
      IF(mbe.EQ.1.AND.mstj(51).GE.1) THEN
        CALL pyboei(nsav)
        goto 150
      ENDIF
 
C...Check that momentum, energy and charge were conserved.
      DO 210 i=1,n
        IF(k(i,1).LE.0.OR.k(i,1).GT.10) goto 210
        DO 200 j=1,4
          ps(2,j)=ps(2,j)+p(i,j)
  200   CONTINUE
        ps(2,6)=ps(2,6)+pychge(k(i,2))
  210 CONTINUE
      pdev=(abs(ps(2,1)-ps(1,1))+abs(ps(2,2)-ps(1,2))+abs(ps(2,3)-
     &ps(1,3))+abs(ps(2,4)-ps(1,4)))/(1d0+abs(ps(2,4))+abs(ps(1,4)))
      IF(mcons.EQ.1.AND.pdev.GT.paru(11)) CALL pyerrm(15,
     &'(PYEXEC:) four-momentum was not conserved')
      IF(mcons.EQ.1.AND.abs(ps(2,6)-ps(1,6)).GT.0.1d0) CALL pyerrm(15,
     &'(PYEXEC:) charge was not conserved')
 
      RETURN
      END