pyprep.f

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

 
C*********************************************************************
 
C...PYPREP
C...Rearranges partons along strings.
C...Special considerations for systems with junctions, with
C...possibility of junction-antijunction annihilation.
C...Allows small systems to collapse into one or two particles.
C...Checks flavours and colour singlet invariant masses.
 
      SUBROUTINE pyprep(IP)
 
C...Double precision and integer declarations.
      IMPLICIT DOUBLE PRECISION(a-h, o-z)
      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/pypars/mstp(200),parp(200),msti(200),pari(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)
C...The common block of colour tags.
      common/pyctag/nct,mct(4000,2)
      SAVE /pyjets/,/pydat1/,/pydat2/,/pydat3/,/pyint1/,/pyctag/,
     &/pypars/
      DATA nerrpr/0/
      SAVE nerrpr
C...Local arrays.
      dimension dps(5),dpc(5),ue(3),pg(5),e1(3),e2(3),e3(3),e4(3),
     &ecl(3),ijunc(10,0:4),ipiece(30,0:4),kfend(4),kfq(4),
     &ijur(4),pju(4,6),irng(4,2),tjj(2,5),t(5),pul(3,5),
     &ijcp(0:6),tjuold(5)
      CHARACTER chtmp*6
 
C...Function to give four-product.
      four(i,j)=p(i,4)*p(j,4)-p(i,1)*p(j,1)-p(i,2)*p(j,2)-p(i,3)*p(j,3)
 
C...Rearrange parton shower product listing along strings: begin loop.
      mstu(24)=0
      nold=n
      i1=n
      njunc=0
      npiece=0
      njjstr=0
      mstu32=mstu(32)+1
      DO 100 i=max(1,ip),n
C...First store junction positions.
        IF(k(i,1).EQ.42) THEN
          njunc=njunc+1
          ijunc(njunc,0)=i
          ijunc(njunc,4)=0
        ENDIF
  100 CONTINUE
 
      DO 250 mqgst=1,3
        DO 240 i=max(1,ip),n
C...Special treatment for junctions
          IF (k(i,1).LE.0) goto 240
          IF(k(i,1).EQ.42) THEN
C...MQGST=2: Look for junction-junction strings (not detected in the
C...main search below).
            IF (mqgst.EQ.2.AND.npiece.NE.3*njunc) THEN
              IF (njjstr.EQ.0) THEN
                njjstr = (3*njunc-npiece)/2
              ENDIF
C...Check how many already identified strings end on this junction
              ilc=0
              DO 110 j=1,npiece
                IF (ipiece(j,4).EQ.i) ilc=ilc+1
  110         CONTINUE
C...If less than 3, remaining must be to another junction
              IF (ilc.LT.3) THEN
                IF (ilc.NE.2) THEN
C...Multiple j-j connections not handled yet.
                  CALL pyerrm(2,
     &            '(PYPREP:) Too many junction-junction strings.')
                  mint(51)=1
                  RETURN
                ENDIF
C...The colour information in the junction is unreadable for the
C...colour space search further down in this routine, so we must
C...start on the colour mother of this junction and then "artificially"
C...prevent the colour mother from connecting here again.
                itjunc=mod(k(i,4)/mstu(5),mstu(5))
                kcs=4
                IF (mod(itjunc,2).EQ.0) kcs=5
C...Switch colour if the junction-junction leg is presumably a
C...junction mother leg rather than a junction daughter leg.
                IF (itjunc.GE.3) kcs=9-kcs
                IF (mint(33).EQ.0) THEN
C...Find the unconnected leg and reorder junction daughter pointers so
C...MOD(K(I,4),MSTU(5)) always points to the junction-junction string
C...piece.
                  ia=mod(k(i,4),mstu(5))
                  IF (k(ia,kcs)/mstu(5)**2.GE.2) THEN
                    itmp=mod(k(i,5),mstu(5))
                    IF (k(itmp,kcs)/mstu(5)**2.GE.2) THEN
                      itmp=mod(k(i,5)/mstu(5),mstu(5))
                      k(i,5)=k(i,5)+(ia-itmp)*mstu(5)
                    ELSE
                      k(i,5)=k(i,5)+(ia-itmp)
                    ENDIF
                    k(i,4)=k(i,4)+(itmp-ia)
                    ia=itmp
                  ENDIF
                  IF (itjunc.LE.2) THEN
C...Beam baryon junction
                    k(ia,kcs)   = k(ia,kcs) + 2*mstu(5)**2
                    k(i,kcs)    = k(i,kcs) + 1*mstu(5)**2
C...Else 1 -> 2 decay junction
                  ELSE
                    k(ia,kcs)   = k(ia,kcs) + mstu(5)**2
                    k(i,kcs)    = k(i,kcs) + 2*mstu(5)**2
                  ENDIF
                  i1beg = i1
                  nstp = 0
                  goto 170
C...Alternatively use colour tag information.
                ELSE
C...Find a final state parton with appropriate dangling colour tag.
                  jct=0
                  ia=0
                  ijumo=k(i,3)
                  DO 140 j1=max(1,ip),n
                    IF (k(j1,1).NE.3) goto 140
C...Check for matching final-state colour tag
                    imatch=0
                    DO 120 j2=max(1,ip),n
                      IF (k(j2,1).NE.3) goto 120
                      IF (mct(j1,kcs-3).EQ.mct(j2,6-kcs)) imatch=1
  120               CONTINUE
                    IF (imatch.EQ.1) goto 140
C...Check whether this colour tag belongs to the present junction
C...by seeing whether any parton with this colour tag has the same
C...mother as the junction.
                    jct=mct(j1,kcs-3)
                    imatch=0
                    DO 130 j2=mint(84)+1,n
                      imo2=k(j2,3)
C...First scattering partons have IMO1 = 3 and 4.
                      IF (imo2.EQ.mint(83)+3.OR.imo2.EQ.mint(83)+4)
     &                     imo2=imo2-2
                      IF (mct(j2,kcs-3).EQ.jct.AND.imo2.EQ.ijumo)
     &                     imatch=1
  130               CONTINUE
                    IF (imatch.EQ.0) goto 140
                    ia=j1
  140             CONTINUE
C...Check for junction-junction strings without intermediate final state
C...glue (not detected above).
                  IF (ia.EQ.0) THEN
                    DO 160 mju=1,njunc
                      iju2=ijunc(mju,0)
                      IF (iju2.EQ.i) goto 160
                      itju2=mod(k(iju2,4)/mstu(5),mstu(5))
C...Only opposite types of junctions can connect to each other.
                      IF (mod(itju2,2).EQ.mod(itjunc,2)) goto 160
                      is=0
                      DO 150 j=1,npiece
                        IF (ipiece(j,4).EQ.iju2) is=is+1
  150                 CONTINUE
                      IF (is.EQ.3) goto 160
                      ib=i
                      ia=iju2
  160               CONTINUE
                  ENDIF
C...Switch to other side of adjacent parton and step from there.
                  kcs=9-kcs
                  i1beg = i1
                  nstp = 0
                  goto 170
                ENDIF
              ELSE IF (ilc.NE.3) THEN
              ENDIF
            ENDIF
          ENDIF
 
C...Look for coloured string endpoint, or (later) leftover gluon.
          IF(k(i,1).NE.3) goto 240
          kc=pycomp(k(i,2))
          IF(kc.EQ.0) goto 240
          kq=kchg(kc,2)
          IF(kq.EQ.0.OR.(mqgst.LE.2.AND.kq.EQ.2)) goto 240
 
C...Pick up loose string end.
          kcs=4
          IF(kq*isign(1,k(i,2)).LT.0) kcs=5
          ia=i
          ib=i
          i1beg=i1
          nstp=0
  170     nstp=nstp+1
          IF(nstp.GT.4*n) THEN
            CALL pyerrm(14,'(PYPREP:) caught in infinite loop')
            mint(51)=1
            RETURN
          ENDIF
 
C...Copy undecayed parton. Finished if reached string endpoint.
          IF(k(ia,1).EQ.3) THEN
            IF(i1.GE.mstu(4)-mstu32-5) THEN
              CALL pyerrm(11,'(PYPREP:) no more memory left in PYJETS')
              mint(51)=1
              mstu(24)=1
              RETURN
            ENDIF
            i1=i1+1
            k(i1,1)=2
            IF(nstp.GE.2.AND.kchg(pycomp(k(ia,2)),2).NE.2) k(i1,1)=1
            k(i1,2)=k(ia,2)
            k(i1,3)=ia
            k(i1,4)=0
            k(i1,5)=0
            DO 180 j=1,5
              p(i1,j)=p(ia,j)
              v(i1,j)=v(ia,j)
  180       CONTINUE
            k(ia,1)=k(ia,1)+10
            IF(k(i1,1).EQ.1) goto 240
          ENDIF
 
C...Also finished (for now) if reached junction; then copy to end.
          IF(k(ia,1).EQ.42) THEN
            ncopy=i1-i1beg
            IF(i1.GE.mstu(4)-mstu32-ncopy-5) THEN
              CALL pyerrm(11,'(PYPREP:) no more memory left in PYJETS')
              mint(51)=1
              mstu(24)=1
              RETURN
            ENDIF
            IF (mqgst.LE.2.AND.ncopy.NE.0) THEN
              DO 200 icopy=1,ncopy
                DO 190 j=1,5
                  k(mstu(4)-mstu32-icopy,j)=k(i1beg+icopy,j)
                  p(mstu(4)-mstu32-icopy,j)=p(i1beg+icopy,j)
                  v(mstu(4)-mstu32-icopy,j)=v(i1beg+icopy,j)
  190           CONTINUE
  200         CONTINUE
            ENDIF
C...For junction-junction strings, find end leg and reorder junction
C...daughter pointers so MOD(K(I,4),MSTU(5)) always points to the
C...junction-junction string piece.
            IF (k(i,1).EQ.42.AND.mint(33).EQ.0) THEN
              itmp=mod(k(ia,4),mstu(5))
              IF (itmp.NE.ib) THEN
                IF (mod(k(ia,5),mstu(5)).EQ.ib) THEN
                  k(ia,5)=k(ia,5)+(itmp-ib)
                ELSE
                  k(ia,5)=k(ia,5)+(itmp-ib)*mstu(5)
                ENDIF
                k(ia,4)=k(ia,4)+(ib-itmp)
              ENDIF
            ENDIF
            npiece=npiece+1
C...IPIECE:
C...0: endpoint in original ER
C...1:
C...2:
C...3: Parton immediately next to junction
C...4: Junction
            ipiece(npiece,0)=i
            ipiece(npiece,1)=mstu32+1
            ipiece(npiece,2)=mstu32+ncopy
            ipiece(npiece,3)=ib
            ipiece(npiece,4)=ia
            mstu32=mstu32+ncopy
            i1=i1beg
            goto 240
          ENDIF
 
C...GOTO next parton in colour space.
          ib=ia
          IF (mint(33).EQ.0) THEN
            IF(mod(k(ib,kcs)/mstu(5)**2,2).EQ.0.AND.mod(k(ib,kcs),mstu(5
     &           )).NE.0) THEN
              ia=mod(k(ib,kcs),mstu(5))
              k(ib,kcs)=k(ib,kcs)+mstu(5)**2
              mrev=0
            ELSE
              IF(k(ib,kcs).GE.2*mstu(5)**2.OR.mod(k(ib,kcs)/mstu(5),
     &             mstu(5)).EQ.0) kcs=9-kcs
              ia=mod(k(ib,kcs)/mstu(5),mstu(5))
              k(ib,kcs)=k(ib,kcs)+2*mstu(5)**2
              mrev=1
            ENDIF
            IF(ia.LE.0.OR.ia.GT.n) THEN
              CALL pyerrm(12,'(PYPREP:) colour rearrangement failed')
              IF(nerrpr.LT.5) THEN
                nerrpr=nerrpr+1
                WRITE(mstu(11),*) 'started at:', i
                WRITE(mstu(11),*) 'ended going from',ib,' to',ia
                WRITE(mstu(11),*) 'MQGST =',mqgst
                CALL pylist(4)
              ENDIF
              mint(51)=1
              RETURN
            ENDIF
            IF(mod(k(ia,4)/mstu(5),mstu(5)).EQ.ib.OR.mod(k(ia,5)/mstu(5)
     &           ,mstu(5)).EQ.ib) THEN
              IF(mrev.EQ.1) kcs=9-kcs
              IF(mod(k(ia,kcs)/mstu(5),mstu(5)).NE.ib) kcs=9-kcs
              k(ia,kcs)=k(ia,kcs)+2*mstu(5)**2
            ELSE
              IF(mrev.EQ.0) kcs=9-kcs
              IF(mod(k(ia,kcs),mstu(5)).NE.ib) kcs=9-kcs
              k(ia,kcs)=k(ia,kcs)+mstu(5)**2
            ENDIF
            IF(ia.NE.i) goto 170
C...Use colour tag information
          ELSE
C...First create colour tags starting on IB if none already present.
            IF (mct(ib,kcs-3).EQ.0) THEN
              CALL pycttr(ib,kcs,ib)
              IF(mint(51).NE.0) RETURN
            ENDIF
            jct=mct(ib,kcs-3)
            ifound=0
C...Find final state tag partner
            DO 210 it=max(1,ip),n
              IF (it.EQ.ib) goto 210
              IF (mct(it,6-kcs).EQ.jct.AND.k(it,1).LT.10.AND.k(it,1).gt
     &             .0) THEN
                ifound=ifound+1
                ia=it
              ENDIF
  210       CONTINUE
C...Just copy and goto next if exactly one partner found.
            IF (ifound.EQ.1) THEN
              goto 170
C...When no match found, match is presumably junction.
            ELSEIF (ifound.EQ.0.AND.mqgst.LE.2) THEN
C...Check whether this colour tag matches a junction
C...by seeing whether any parton with this colour tag has the same
C...mother as a junction.
C...NB: Only type 1 and 2 junctions handled presently.
              DO 230 iju=1,njunc
                ijumo=k(ijunc(iju,0),3)
                itjunc=mod(k(ijunc(iju,0),4)/mstu(5),mstu(5))
C...Colours only connect to junctions, anti-colours to antijunctions:
                IF (mod(itjunc+1,2)+1.NE.kcs-3) goto 230
                imatch=0
                DO 220 j1=max(1,ip),n
                  IF (k(j1,1).LE.0) goto 220
C...First scattering partons have IMO1 = 3 and 4.
                  imo=k(j1,3)
                  IF (imo.EQ.mint(83)+3.OR.imo.EQ.mint(83)+4)
     &                 imo=imo-2
                  IF (mct(j1,kcs-3).EQ.jct.AND.imo.EQ.ijumo.AND.mod(k(j1
     &                 ,3+itjunc)/mstu(5),mstu(5)).EQ.ijunc(iju,0))
     &                 imatch=1
C...Attempt at handling type > 3 junctions also. Not tested.
                  IF (itjunc.GE.3.AND.mct(j1,6-kcs).EQ.jct.AND.imo.eq
     &                 .ijumo) imatch=1
  220           CONTINUE
                IF (imatch.EQ.0) goto 230
                ia=ijunc(iju,0)
                ifound=ifound+1
  230         CONTINUE
 
              IF (ifound.EQ.1) THEN
                goto 170
              ELSEIF (ifound.EQ.0) THEN
                WRITE(chtmp,*) jct
                CALL pyerrm(12,'(PYPREP:) no matching colour tag: '
     &               //chtmp)
                IF(nerrpr.LT.5) THEN
                  nerrpr=nerrpr+1
                  CALL pylist(4)
                ENDIF
                mint(51)=1
                RETURN
              ENDIF
            ELSEIF (ifound.GE.2) THEN
              WRITE(chtmp,*) jct
              CALL pyerrm(12
     &             ,'(PYPREP:) too many occurences of colour line: '//
     &             chtmp)
              IF(nerrpr.LT.5) THEN
                nerrpr=nerrpr+1
                CALL pylist(4)
              ENDIF
              mint(51)=1
              RETURN
            ENDIF
          ENDIF
          k(i1,1)=1
  240   CONTINUE
  250 CONTINUE
 
C...Junction systems remain.
      iju=0
      ijus=0
      ijucnt=0
      mrev=0
      ijjstr=0
  260 ijucnt=ijucnt+1
      IF (ijucnt.LE.njunc) THEN
C...If we are not processing a j-j string, treat this junction as new.
        IF (ijjstr.EQ.0) THEN
          iju=ijunc(ijucnt,0)
          mrev=0
C...If junction has already been read, ignore it.
          IF (ijunc(ijucnt,4).EQ.1) goto 260
C...If we are on a j-j string, goto second j-j junction.
        ELSE
          ijucnt=ijucnt-1
          iju=ijus
        ENDIF
C...Mark selected junction read.
        DO 270 j=1,njunc
          IF (ijunc(j,0).EQ.iju) ijunc(j,4)=1
  270   CONTINUE
C...Determine junction type
        itjunc = mod(k(iju,4)/mstu(5),mstu(5))
C...Type 1 and 2 junctions: ~chi -> q q q, ~chi -> qbar,qbar,qbar
C...Type 3 and 4 junctions: ~qbar -> q q , ~q -> qbar qbar
C...Type 5 and 6 junctions: ~g -> q q q, ~g -> qbar qbar qbar
        IF (itjunc.GE.1.AND.itjunc.LE.6) THEN
          ihk=0
  280     ihk=ihk+1
C...Find which quarks belong to given junction.
          ihf=0
          DO 290 ipc=1,npiece
            IF (ipiece(ipc,4).EQ.iju) THEN
              ihf=ihf+1
              IF (ihf.EQ.ihk) iend=ipiece(ipc,3)
            ENDIF
            IF (ihk.EQ.3.AND.ipiece(ipc,0).EQ.iju) iend=ipiece(ipc,3)
  290     CONTINUE
C...IHK = 3 is special. Either normal string piece, or j-j string.
          IF(ihk.EQ.3) THEN
            IF (mrev.NE.1) THEN
              DO 300 ipc=1,npiece
C...If there is a j-j string starting on the present junction which has
C...zero length, insert next junction immediately.
                IF (ipiece(ipc,0).EQ.iju.AND.k(ipiece(ipc,4),1)
     &          .EQ.42.AND.ipiece(ipc,1)-1-ipiece(ipc,2).EQ.0) THEN
                  ijjstr = 1
                  goto 340
                ENDIF
  300         CONTINUE
              mrev = 1
C...If MREV is 1 and IHK is 3 we are finished with this system.
            ELSE
              mrev=0
              goto 260
            ENDIF
          ENDIF
 
C...If we've gotten this far, then either IHK < 3, or
C...an interjunction string exists, or just a third normal string.
          ijunc(ijucnt,ihk)=0
          ijjstr = 0
C..Order pieces belonging to this junction. Also look for j-j.
          DO 310 ipc=1,npiece
            IF (ipiece(ipc,3).EQ.iend) ijunc(ijucnt,ihk)=ipc
            IF (ihk.EQ.3.AND.ipiece(ipc,0).EQ.ijunc(ijucnt,0)
     &      .AND.k(ipiece(ipc,4),1).EQ.42) THEN
              ijunc(ijucnt,ihk)=ipc
              ijjstr = 1
              mrev = 0
            ENDIF
  310     CONTINUE
C...Copy back chains in proper order. MREV=0/1 : descending/ascending
          ipc=ijunc(ijucnt,ihk)
C...Temporary solution to cover for bug.
          IF(ipc.LE.0) THEN
            CALL pyerrm(12,'(PYPREP:) fails to hook up junctions')
            mint(51)=1
            RETURN
          ENDIF
          DO 330 icp=ipiece(ipc,1+mrev),ipiece(ipc,2-mrev),1-2*mrev
            i1=i1+1
            DO 320 j=1,5
              k(i1,j)=k(mstu(4)-icp,j)
              p(i1,j)=p(mstu(4)-icp,j)
              v(i1,j)=v(mstu(4)-icp,j)
  320       CONTINUE
  330     CONTINUE
          k(i1,1)=2
C...Mark last quark.
          IF (mrev.EQ.1.AND.ihk.GE.2) k(i1,1)=1
C...Do not insert junctions at wrong places.
          IF(ihk.LT.2.OR.mrev.NE.0) goto 360
C...Insert junction.
  340     ijus = iju
          IF (ihk.EQ.3) THEN
C...Shift to end junction if a j-j string has been processed.
            IF (ijjstr.NE.0) ijus = ipiece(ipc,4)
            mrev= 1
          ENDIF
          i1=i1+1
          DO 350 j=1,5
            k(i1,j)=0
            p(i1,j)=0.
            v(i1,j)=0.
  350     CONTINUE
          k(i1,1)=41
          k(ijus,1)=k(ijus,1)+10
          k(i1,2)=k(ijus,2)
          k(i1,3)=ijus
  360     IF (ihk.LT.3) goto 280
        ELSE
          CALL pyerrm(12,'(PYPREP:) Unknown junction type')
          mint(51)=1
          RETURN
        ENDIF
        IF (ijucnt.NE.njunc) goto 260
      ENDIF
      n=i1
 
C...Rearrange three strings from junction, e.g. in case one has been
C...shortened by shower, so the last is the largest-energy one.
      IF(njunc.GE.1) THEN
C...Find systems with exactly one junction.
        mjun1=0
        nbeg=nold+1
        DO 470 i=nold+1,n
          IF(k(i,1).NE.1.AND.k(i,1).NE.41) THEN
          ELSEIF(k(i,1).EQ.41) THEN
            mjun1=mjun1+1
          ELSEIF(k(i,1).EQ.1.AND.mjun1.NE.1) THEN
            mjun1=0
            nbeg=i+1
          ELSE
            nend=i
C...Sum up energy-momentum in each junction string.
            DO 370 j=1,5
              pju(1,j)=0d0
              pju(2,j)=0d0
              pju(3,j)=0d0
  370       CONTINUE
            nju=0
            DO 390 i1=nbeg,nend
              IF(k(i1,2).NE.21) THEN
                nju=nju+1
                ijur(nju)=i1
              ENDIF
              DO 380 j=1,5
                pju(min(nju,3),j)=pju(min(nju,3),j)+p(i1,j)
  380         CONTINUE
  390       CONTINUE
C...Find which of them has highest energy (minus mass) in rest frame.
            DO 400 j=1,5
              pju(4,j)=pju(1,j)+pju(2,j)+pju(3,j)
  400       CONTINUE
            pmju=sqrt(max(0d0,pju(4,4)**2-pju(4,1)**2-pju(4,2)**2-
     &      pju(4,3)**2))
            DO 410 i2=1,3
              pju(i2,6)=(pju(4,4)*pju(i2,4)-pju(4,1)*pju(i2,1)-
     &        pju(4,2)*pju(i2,2)-pju(4,3)*pju(i2,3))/pmju-pju(i2,5)
  410       CONTINUE
            IF(pju(3,6).LT.min(pju(1,6),pju(2,6))) THEN
C...Decide how to rearrange so that new last has highest energy.
              IF(pju(1,6).LT.pju(2,6)) THEN
                irng(1,1)=ijur(1)
                irng(1,2)=ijur(2)-1
                irng(2,1)=ijur(4)
                irng(2,2)=ijur(3)+1
                irng(4,1)=ijur(3)-1
                irng(4,2)=ijur(2)
              ELSE
                irng(1,1)=ijur(4)
                irng(1,2)=ijur(3)+1
                irng(2,1)=ijur(2)
                irng(2,2)=ijur(3)-1
                irng(4,1)=ijur(2)-1
                irng(4,2)=ijur(1)
              ENDIF
              irng(3,1)=ijur(3)
              irng(3,2)=ijur(3)
C...Copy in correct order below bottom of current event record.
              i2=n
              DO 440 ii=1,4
                DO 430 i1=irng(ii,1),irng(ii,2),
     &          isign(1,irng(ii,2)-irng(ii,1))
                  i2=i2+1
                  IF(i2.GE.mstu(4)-mstu32-5) THEN
                    CALL pyerrm(11,
     &              '(PYPREP:) no more memory left in PYJETS')
                    mint(51)=1
                    mstu(24)=1
                    RETURN
                  ENDIF
                  DO 420 j=1,5
                    k(i2,j)=k(i1,j)
                    p(i2,j)=p(i1,j)
                    v(i2,j)=v(i1,j)
  420             CONTINUE
                  IF(k(i2,1).EQ.1) k(i2,1)=2
  430           CONTINUE
  440         CONTINUE
              k(i2,1)=1
C...Copy back up, overwriting but now in correct order.
              DO 460 i1=nbeg,nend
                i2=i1-nbeg+n+1
                DO 450 j=1,5
                  k(i1,j)=k(i2,j)
                  p(i1,j)=p(i2,j)
                  v(i1,j)=v(i2,j)
  450           CONTINUE
  460         CONTINUE
            ENDIF
            mjun1=0
            nbeg=i+1
          ENDIF
  470   CONTINUE
 
C...Check whether q-q-j-j-qbar-qbar systems should be collapsed
C...to two q-qbar systems.
C...(MSTJ(19)=1 forces q-q-j-j-qbar-qbar.)
        IF (mstj(19).NE.1) THEN
          mjun1  = 0
          jjglue = 0
          nbeg   = nold+1
C...Force collapse when MSTJ(19)=2.
          IF (mstj(19).EQ.2) THEN
            delmjj = 1d9
            delmqq = 0d0
          ENDIF
C...Find systems with exactly two junctions.
          DO 700 i=nold+1,n
C...Count junctions
            IF (k(i,1).EQ.41) THEN
              mjun1 = mjun1+1
C...Check for interjunction gluons
              IF (mjun1.EQ.2.AND.k(i-1,1).NE.41) THEN
                jjglue = 1
              ENDIF
            ELSEIF(k(i,1).EQ.1.AND.(mjun1.NE.2)) THEN
C...If end of system reached with either zero or one junction, restart
C...with next system.
              mjun1  = 0
              jjglue = 0
              nbeg   = i+1
            ELSEIF(k(i,1).EQ.1) THEN
C...If end of system reached with exactly two junctions, compute string
C...length measure for the (q-q-j-j-qbar-qbar) topology and compare with
C...length measure for the (q-qbar)(q-qbar) topology.
              nend=i
C...Loop down through chain.
              isid=0
              DO 480 i1=nbeg,nend
C...Store string piece division locations in event record
                IF (k(i1,2).NE.21) THEN
                  isid       = isid+1
                  ijcp(isid) = i1
                ENDIF
  480         CONTINUE
C...Randomly choose between (1,3)(2,4) and (1,4)(2,3) topologies.
              isw=0
              IF (pyr(0).LT.0.5d0) isw=1
C...Randomly choose which qqbar string gets the jj gluons.
              igs=1
              IF (pyr(0).GT.0.5d0) igs=2
C...Only compute string lengths when no topology forced.
              IF (mstj(19).EQ.0) THEN
C...Repeat following for each junction
                DO 570 iju=1,2
C...Initialize iterative procedure for finding JRF
                  ijrfit=0
                  DO 490 ix=1,3
                    tjuold(ix)=0d0
  490             CONTINUE
                  tjuold(4)=1d0
C...Start iteration. Sum up momenta in string pieces
  500             DO 540 ijs=1,3
C...JD=-1 for first junction, +1 for second junction.
C...Find out where piece starts and ends and which direction to go.
                    jd=2*iju-3
                    IF (ijs.LE.2) THEN
                      ia = ijcp((iju-1)*7 - jd*(ijs+1)) + jd
                      ib = ijcp((iju-1)*7 - jd*ijs)
                    ELSEIF (ijs.EQ.3) THEN
                      jd =-jd
                      ia = ijcp((iju-1)*7 + jd*(ijs)) + jd
                      ib = ijcp((iju-1)*7 + jd*(ijs+3))
                    ENDIF
C...Initialize junction pull 4-vector.
                    DO 510 j=1,5
                      pul(ijs,j)=0d0
  510               CONTINUE
C...Initialize weight
                    pwt = 0d0
                    pwtold = 0d0
C...Sum up (weighted) momenta along each string piece
                    DO 530 isp=ia,ib,jd
C...If present parton not last in chain
                      IF (isp.NE.ia.AND.isp.NE.ib) THEN
C...If last parton was a junction, store present weight
                        IF (k(isp-jd,2).EQ.88) THEN
                          pwtold = pwt
C...If last parton was a quark, reset to stored weight.
                        ELSEIF (k(isp-jd,2).NE.21) THEN
                          pwt = pwtold
                        ENDIF
                      ENDIF
C...Skip next parton if weight already large
                      IF (pwt.GT.10d0) goto 530
C...Compute momentum in TJUOLD frame:
                      tdp=tjuold(1)*p(isp,1)+tjuold(2)*p(isp,2)+tjuold(3
     &                     )*p(isp,3)
                      bfc=tdp/(1d0+tjuold(4))+p(isp,4)
                      DO 520 j=1,3
                        tmp=p(isp,j)+tjuold(j)*bfc
                        pul(ijs,j)=pul(ijs,j)+tmp*exp(-pwt)
  520                 CONTINUE
C...Boosted energy
                      tmp=tjuold(4)*p(isp,4)+tdp
                      pul(ijs,4)=pul(ijs,j)+tmp*exp(-pwt)
C...Update weight
                      pwt=pwt+tmp/parj(48)
C...Put |p| rather than m in 5th slot
                      pul(ijs,5)=sqrt(pul(ijs,1)**2+pul(ijs,2)**2
     &                     +pul(ijs,3)**2)
  530               CONTINUE
  540             CONTINUE
C...Compute boost
                  ijrfit=ijrfit+1
                  CALL pyjurf(pul,t)
C...Combine new boost (T) with old boost (TJUOLD)
                  tmp=t(1)*tjuold(1)+t(2)*tjuold(2)+t(3)*tjuold(3)
                  DO 550 ix=1,3
                    tjuold(ix)=t(ix)+tjuold(ix)*(tmp/(1d0+tjuold(4))+t(4
     &                   ))
  550             CONTINUE
                  tjuold(4)=sqrt(1d0+tjuold(1)**2+tjuold(2)**2+tjuold(3)
     &                 **2)
C...If last boost small, accept JRF, else iterate.
C...Also prevent possibility of infinite loop.
                  IF (abs((t(4)-1d0)/tjuold(4)).GT.0.01d0.AND.
     &                 ijrfit.LT.mstj(18))THEN
                    goto 500
                  ELSEIF (ijrfit.GE.mstj(18)) THEN
                    CALL pyerrm(1,'(PYPREP:) failed to converge on JRF')
                  ENDIF
C...Store final boost, with change of sign since TJJ motion vector.
                  DO 560 ix=1,3
                    tjj(iju,ix)=-tjuold(ix)
  560             CONTINUE
                  tjj(iju,4)=sqrt(1d0+tjj(iju,1)**2+tjj(iju,2)**2
     &                 +tjj(iju,3)**2)
  570           CONTINUE
C...String length measure for (q-qbar)(q-qbar) topology.
C...Note only momenta of nearest partons used (since rest of system
C...identical).
                IF (jjglue.EQ.0) THEN
                  delmqq=4d0*four(ijcp(2)-1,ijcp(4+isw)+1)*four(ijcp(3)
     &                 -1,ijcp(5-isw)+1)
                ELSE
C...Put jj gluons on selected string (IGS selected randomly above).
                  IF (igs.EQ.1) THEN
                    delmqq=8d0*four(ijcp(2)-1,ijcp(4)-1)*four(ijcp(3)+1
     &                   ,ijcp(4+isw)+1)*four(ijcp(3)-1,ijcp(5-isw)+1)
                  ELSE
                    delmqq=8d0*four(ijcp(2)-1,ijcp(4+isw)+1)
     &                   *four(ijcp(3)-1,ijcp(4)-1)*four(ijcp(3)+1
     &                   ,ijcp(5-isw)+1)
                  ENDIF
                ENDIF
C...String length measure for q-q-j-j-q-q topology.
                t1g1=0d0
                t2g2=0d0
                t1t2=0d0
                t1p1=0d0
                t1p2=0d0
                t2p3=0d0
                t2p4=0d0
                isgn=-1
C...Note only momenta of nearest partons used (since rest of system
C...identical).
                DO 580 ix=1,4
                  IF (ix.EQ.4) isgn=1
                  t1p1=t1p1+isgn*tjj(1,ix)*p(ijcp(2)-1,ix)
                  t1p2=t1p2+isgn*tjj(1,ix)*p(ijcp(3)-1,ix)
                  t2p3=t2p3+isgn*tjj(2,ix)*p(ijcp(4)+1,ix)
                  t2p4=t2p4+isgn*tjj(2,ix)*p(ijcp(5)+1,ix)
                  IF (jjglue.EQ.0) THEN
C...Junction motion vector dot product gives length when inter-junction
C...gluons absent.
                    t1t2=t1t2+isgn*tjj(1,ix)*tjj(2,ix)
                  ELSE
C...Junction motion vector dot products with gluon momenta give length
C...when inter-junction gluons present.
                    t1g1=t1g1+isgn*tjj(1,ix)*p(ijcp(3)+1,ix)
                    t2g2=t2g2+isgn*tjj(2,ix)*p(ijcp(4)-1,ix)
                  ENDIF
  580           CONTINUE
                delmjj=16d0*t1p1*t1p2*t2p3*t2p4
                IF (jjglue.EQ.0) THEN
                  delmjj=delmjj*(t1t2+sqrt(t1t2**2-1))
                ELSE
                  delmjj=delmjj*4d0*t1g1*t2g2
                ENDIF
              ENDIF
C...If delmjj > delmqq collapse string system to q-qbar q-qbar
C...(Always the case for MSTJ(19)=2 due to initialization above)
              IF (delmjj.GT.delmqq) THEN
C...Put new system at end of event record
                ncop=n
                DO 650 ist=1,2
                  DO 600 icop=ijcp(ist),ijcp(ist+1)-1
                    ncop=ncop+1
                    DO 590 ix=1,5
                      p(ncop,ix)=p(icop,ix)
                      k(ncop,ix)=k(icop,ix)
  590               CONTINUE
  600             CONTINUE
                  IF (jjglue.NE.0.AND.ist.EQ.igs) THEN
C...Insert inter-junction gluon string piece (reversed)
                    njjgl=0
                    DO 620 icop=ijcp(4)-1,ijcp(3)+1,-1
                      njjgl=njjgl+1
                      ncop=ncop+1
                      DO 610 ix=1,5
                        p(ncop,ix)=p(icop,ix)
                        k(ncop,ix)=k(icop,ix)
  610                 CONTINUE
  620               CONTINUE
                    ENDIF
                  ifc=-2*ist+3
                  DO 640 icop=ijcp(ist+ifc*isw+3)+1,ijcp(ist+ifc*isw+4)
                    ncop=ncop+1
                    DO 630 ix=1,5
                      p(ncop,ix)=p(icop,ix)
                      k(ncop,ix)=k(icop,ix)
  630               CONTINUE
  640             CONTINUE
                  k(ncop,1)=1
  650           CONTINUE
C...Copy system back in right order
                DO 670 icop=nbeg,nend-2
                  DO 660 ix=1,5
                    p(icop,ix)=p(n+icop-nbeg+1,ix)
                    k(icop,ix)=k(n+icop-nbeg+1,ix)
  660             CONTINUE
  670           CONTINUE
C...Shift down rest of event record
                DO 690 icop=nend+1,n
                  DO 680 ix=1,5
                    p(icop-2,ix)=p(icop,ix)
                    k(icop-2,ix)=k(icop,ix)
  680             CONTINUE
  690             CONTINUE
C...Update length of event record.
                n=n-2
              ENDIF
              mjun1=0
              nbeg=i+1
            ENDIF
  700     CONTINUE
        ENDIF
      ENDIF
 
C...Done if no checks on small-mass systems.
      IF(mstj(14).LT.0) RETURN
      IF(mstj(14).EQ.0) goto 1140
 
C...Find lowest-mass colour singlet jet system.
      ns=n
  710 nsin=n-ns
      pdmin=1d0+parj(32)
      ic=0
      DO 770 i=max(1,ip),n
        IF(k(i,1).NE.1.AND.k(i,1).NE.2) THEN
        ELSEIF(k(i,1).EQ.2.AND.ic.EQ.0) THEN
          nsin=nsin+1
          ic=i
          DO 720 j=1,4
            dps(j)=p(i,j)
  720     CONTINUE
          mstj(93)=1
          dps(5)=pymass(k(i,2))
        ELSEIF(k(i,1).EQ.2.AND.k(i,2).NE.21) THEN
          DO 730 j=1,4
            dps(j)=dps(j)+p(i,j)
  730     CONTINUE
          mstj(93)=1
          dps(5)=dps(5)+pymass(k(i,2))
        ELSEIF(k(i,1).EQ.2) THEN
          DO 740 j=1,4
            dps(j)=dps(j)+p(i,j)
  740     CONTINUE
        ELSEIF(ic.NE.0.AND.kchg(pycomp(k(i,2)),2).NE.0) THEN
          DO 750 j=1,4
            dps(j)=dps(j)+p(i,j)
  750     CONTINUE
          mstj(93)=1
          dps(5)=dps(5)+pymass(k(i,2))
          pd=sqrt(max(0d0,dps(4)**2-dps(1)**2-dps(2)**2-dps(3)**2))-
     &    dps(5)
          IF(pd.LT.pdmin) THEN
            pdmin=pd
            DO 760 j=1,5
              dpc(j)=dps(j)
  760       CONTINUE
            ic1=ic
            ic2=i
          ENDIF
          ic=0
        ELSE
          nsin=nsin+1
        ENDIF
  770 CONTINUE
 
C...Done if lowest-mass system above threshold for string frag.
      IF(pdmin.GE.parj(32)) goto 1140
 
C...Fill small-mass system as cluster.
      nsav=n
      pecm=sqrt(max(0d0,dpc(4)**2-dpc(1)**2-dpc(2)**2-dpc(3)**2))
      k(n+1,1)=11
      k(n+1,2)=91
      k(n+1,3)=ic1
      p(n+1,1)=dpc(1)
      p(n+1,2)=dpc(2)
      p(n+1,3)=dpc(3)
      p(n+1,4)=dpc(4)
      p(n+1,5)=pecm
 
C...Set up history, assuming cluster -> 2 hadrons.
      nbody=2
      k(n+1,4)=n+2
      k(n+1,5)=n+3
      k(n+2,1)=1
      k(n+3,1)=1
      IF(mstu(16).NE.2) THEN
        k(n+2,3)=n+1
        k(n+3,3)=n+1
      ELSE
        k(n+2,3)=ic1
        k(n+3,3)=ic2
      ENDIF
      k(n+2,4)=0
      k(n+3,4)=0
      k(n+2,5)=0
      k(n+3,5)=0
      v(n+1,5)=0d0
      v(n+2,5)=0d0
      v(n+3,5)=0d0
 
C...Find total flavour content - complicated by presence of junctions.
      nq=0
      ndiq=0
      DO 780 i=ic1,ic2
        IF((k(i,1).EQ.1.OR.k(i,1).EQ.2).AND.k(i,2).NE.21) THEN
          nq=nq+1
          kfq(nq)=k(i,2)
          IF(iabs(k(i,2)).GT.1000) ndiq=ndiq+1
        ENDIF
  780 CONTINUE
 
C...If several diquarks, split up one to give even number of flavours.
      IF(nq.EQ.3.AND.ndiq.GE.2) THEN
        i1=3
        IF(iabs(kfq(3)).LT.1000) i1=1
        kfq(4)=isign(mod(iabs(kfq(i1))/100,10),kfq(i1))
        kfq(i1)=kfq(i1)/1000
        nq=4
        ndiq=ndiq-1
      ENDIF
 
C...If four quark ends, join two to diquark.
      IF(nq.EQ.4.AND.ndiq.EQ.0) THEN
        i1=1
        i2=2
        IF(kfq(i1)*kfq(i2).LT.0) i2=3
        IF(i2.EQ.3.AND.kfq(i1)*kfq(i2).LT.0) i2=4
        kfls=2*int(pyr(0)+3d0*parj(4)/(1d0+3d0*parj(4)))+1
        IF(kfq(i1).EQ.kfq(i2)) kfls=3
        kfq(i1)=isign(1000*max(iabs(kfq(i1)),iabs(kfq(i2)))+
     &  100*min(iabs(kfq(i1)),iabs(kfq(i2)))+kfls,kfq(i1))
        kfq(i2)=kfq(4)
        nq=3
        ndiq=1
      ENDIF
 
C...If two quark ends, plus quark or diquark, join quarks to diquark.
      IF(nq.EQ.3) THEN
        i1=1
        i2=2
        IF(iabs(kfq(i1)).GT.1000) i1=3
        IF(iabs(kfq(i2)).GT.1000) i2=3
        kfls=2*int(pyr(0)+3d0*parj(4)/(1d0+3d0*parj(4)))+1
        IF(kfq(i1).EQ.kfq(i2)) kfls=3
        kfq(i1)=isign(1000*max(iabs(kfq(i1)),iabs(kfq(i2)))+
     &  100*min(iabs(kfq(i1)),iabs(kfq(i2)))+kfls,kfq(i1))
        kfq(i2)=kfq(3)
        nq=2
        ndiq=ndiq+1
      ENDIF
 
C...Form two particles from flavours of lowest-mass system, if feasible.
      ntry = 0
  790 ntry = ntry + 1
 
C...Open string with two specified endpoint flavours.
      IF(nq.EQ.2) THEN
        kc1=pycomp(kfq(1))
        kc2=pycomp(kfq(2))
        IF(kc1.EQ.0.OR.kc2.EQ.0) goto 1140
        kq1=kchg(kc1,2)*isign(1,kfq(1))
        kq2=kchg(kc2,2)*isign(1,kfq(2))
        IF(kq1+kq2.NE.0) goto 1140
C...Start with qq, if there is one. Only allow for rank 1 popcorn meson
  800   k1=kfq(1)
        IF(iabs(kfq(2)).GT.1000) k1=kfq(2)
        mstu(125)=0
        CALL pydcyk(k1,0,kfln,k(n+2,2))
        CALL pydcyk(kfq(1)+kfq(2)-k1,-kfln,kfldmp,k(n+3,2))
        IF(k(n+2,2).EQ.0.OR.k(n+3,2).EQ.0) goto 800
 
C...Open string with four specified flavours.
      ELSEIF(nq.EQ.4) THEN
        kc1=pycomp(kfq(1))
        kc2=pycomp(kfq(2))
        kc3=pycomp(kfq(3))
        kc4=pycomp(kfq(4))
        IF(kc1.EQ.0.OR.kc2.EQ.0.OR.kc3.EQ.0.OR.kc4.EQ.0) goto 1140
        kq1=kchg(kc1,2)*isign(1,kfq(1))
        kq2=kchg(kc2,2)*isign(1,kfq(2))
        kq3=kchg(kc3,2)*isign(1,kfq(3))
        kq4=kchg(kc4,2)*isign(1,kfq(4))
        IF(kq1+kq2+kq3+kq4.NE.0) goto 1140
C...Combine flavours pairwise to form two hadrons.
  810   i1=1
        i2=2
        IF(kq1*kq2.GT.0.OR.(iabs(kfq(1)).GT.1000.AND.
     &  iabs(kfq(2)).GT.1000)) i2=3
        IF(i2.EQ.3.AND.(kq1*kq3.GT.0.OR.(iabs(kfq(1)).GT.1000.AND.
     &  iabs(kfq(3)).GT.1000))) i2=4
        i3=3
        IF(i2.EQ.3) i3=2
        i4=10-i1-i2-i3
        CALL pydcyk(kfq(i1),kfq(i2),kfldmp,k(n+2,2))
        CALL pydcyk(kfq(i3),kfq(i4),kfldmp,k(n+3,2))
        IF(k(n+2,2).EQ.0.OR.k(n+3,2).EQ.0) goto 810
 
C...Closed string.
      ELSE
        IF(iabs(k(ic2,2)).NE.21) goto 1140
C...No room for popcorn mesons in closed string -> 2 hadrons.
        mstu(125)=0
  820   CALL pydcyk(1+int((2d0+parj(2))*pyr(0)),0,kfln,kfdmp)
        CALL pydcyk(kfln,0,kflm,k(n+2,2))
        CALL pydcyk(-kfln,-kflm,kfldmp,k(n+3,2))
        IF(k(n+2,2).EQ.0.OR.k(n+3,2).EQ.0) goto 820
      ENDIF
      p(n+2,5)=pymass(k(n+2,2))
      p(n+3,5)=pymass(k(n+3,2))
 
C...If it does not work: try again (a number of times), give up (if no
C...place to shuffle momentum or too many flavours), or form one hadron.
      IF(p(n+2,5)+p(n+3,5)+parj(64).GE.pecm) THEN
        IF(ntry.LT.mstj(17).OR.(nq.EQ.4.AND.ntry.LT.5*mstj(17))) THEN
          goto 790
        ELSEIF(nsin.EQ.1.OR.nq.EQ.4) THEN
          goto 1140
        ELSE
          goto 890
        END IF
      END IF
 
C...Perform two-particle decay of jet system.
C...First step: find reference axis in decaying system rest frame.
C...(Borrow slot N+2 for temporary direction.)
      DO 830 j=1,4
        p(n+2,j)=p(ic1,j)
  830 CONTINUE
      DO 850 i=ic1+1,ic2-1
        IF((k(i,1).EQ.1.OR.k(i,1).EQ.2).AND.
     &  kchg(pycomp(k(i,2)),2).NE.0) THEN
          frac1=four(ic2,i)/(four(ic1,i)+four(ic2,i))
          DO 840 j=1,4
            p(n+2,j)=p(n+2,j)+frac1*p(i,j)
  840     CONTINUE
        ENDIF
  850 CONTINUE
      CALL pyrobo(n+2,n+2,0d0,0d0,-dpc(1)/dpc(4),-dpc(2)/dpc(4),
     &-dpc(3)/dpc(4))
      the1=pyangl(p(n+2,3),sqrt(p(n+2,1)**2+p(n+2,2)**2))
      phi1=pyangl(p(n+2,1),p(n+2,2))
 
C...Second step: generate isotropic/anisotropic decay.
      pa=sqrt((pecm**2-(p(n+2,5)+p(n+3,5))**2)*(pecm**2-
     &(p(n+2,5)-p(n+3,5))**2))/(2d0*pecm)
  860 ue(3)=pyr(0)
      IF(parj(21).LE.0.01d0) ue(3)=1d0
      pt2=(1d0-ue(3)**2)*pa**2
      IF(mstj(16).LE.0) THEN
        prev=0.5d0
      ELSE
        IF(exp(-pt2/(2d0*max(0.01d0,parj(21))**2)).LT.pyr(0)) goto 860
        pr1=p(n+2,5)**2+pt2
        pr2=p(n+3,5)**2+pt2
        alambd=sqrt(max(0d0,(pecm**2-pr1-pr2)**2-4d0*pr1*pr2))
        prevcf=parj(42)
        IF(mstj(11).EQ.2) prevcf=parj(39)
        prev=1d0/(1d0+exp(min(50d0,prevcf*alambd*parj(40))))
      ENDIF
      IF(pyr(0).LT.prev) ue(3)=-ue(3)
      phi=paru(2)*pyr(0)
      ue(1)=sqrt(1d0-ue(3)**2)*cos(phi)
      ue(2)=sqrt(1d0-ue(3)**2)*sin(phi)
      DO 870 j=1,3
        p(n+2,j)=pa*ue(j)
        p(n+3,j)=-pa*ue(j)
  870 CONTINUE
      p(n+2,4)=sqrt(pa**2+p(n+2,5)**2)
      p(n+3,4)=sqrt(pa**2+p(n+3,5)**2)
 
C...Third step: move back to event frame and set production vertex.
      CALL pyrobo(n+2,n+3,the1,phi1,dpc(1)/dpc(4),dpc(2)/dpc(4),
     &dpc(3)/dpc(4))
      DO 880 j=1,4
        v(n+1,j)=v(ic1,j)
        v(n+2,j)=v(ic1,j)
        v(n+3,j)=v(ic2,j)
  880 CONTINUE
      n=n+3
      goto 1120
 
C...Else form one particle, if possible.
  890 nbody=1
      k(n+1,5)=n+2
      DO 900 j=1,4
        v(n+1,j)=v(ic1,j)
        v(n+2,j)=v(ic1,j)
  900 CONTINUE
 
C...Select hadron flavour from available quark flavours.
  910 IF(nq.EQ.2.AND.iabs(kfq(1)).GT.100.AND.iabs(kfq(2)).GT.100) THEN
        goto 1140
      ELSEIF(nq.EQ.2) THEN
        CALL pykfdi(kfq(1),kfq(2),kfldmp,k(n+2,2))
      ELSE
        kfln=1+int((2d0+parj(2))*pyr(0))
        CALL pykfdi(kfln,-kfln,kfldmp,k(n+2,2))
      ENDIF
      IF(k(n+2,2).EQ.0) goto 910
      p(n+2,5)=pymass(k(n+2,2))
 
C...Use old algorithm for E/p conservation? (EN)
      IF (mstj(16).LE.0) goto 1080
 
C...Find the string piece closest to the cluster by a loop
C...over the undecayed partons not in present cluster. (EN)
      dglomi=1d30
      ibeg=0
      i0=0
      njunc=0
      DO 940 i1=max(1,ip),n-1
        IF(k(i1,1).EQ.1) njunc=0
        IF(k(i1,1).EQ.41) njunc=njunc+1
        IF(k(i1,1).EQ.41) goto 940
        IF(i1.GE.ic1-1.AND.i1.LE.ic2) THEN
          i0=0
        ELSEIF(k(i1,1).EQ.2) THEN
          IF(i0.EQ.0) i0=i1
          i2=i1
  920     i2=i2+1
          IF(k(i2,1).EQ.41) goto 940
          IF(k(i2,1).GT.10) goto 920
          IF(kchg(pycomp(k(i2,2)),2).EQ.0) goto 920
          IF(k(i1,2).EQ.21.AND.k(i2,2).NE.21.AND.k(i2,1).NE.1.AND.
     &    njunc.EQ.0) goto 940
          IF(k(i1,2).NE.21.AND.k(i2,2).EQ.21.AND.njunc.NE.0) goto 940
          IF(k(i1,2).NE.21.AND.k(i2,2).NE.21.AND.(i1.GT.i0.OR.
     &    k(i2,1).NE.1)) goto 940
 
C...Define velocity vectors e1, e2, ecl and differences e3, e4.
          DO 930 j=1,3
            e1(j)=p(i1,j)/p(i1,4)
            e2(j)=p(i2,j)/p(i2,4)
            ecl(j)=p(n+1,j)/p(n+1,4)
            e3(j)=e2(j)-e1(j)
            e4(j)=ecl(j)-e1(j)
  930     CONTINUE
 
C...Calculate minimal D=(e4-alpha*e3)**2 for 0<alpha<1.
          e3s=e3(1)**2+e3(2)**2+e3(3)**2
          e4s=e4(1)**2+e4(2)**2+e4(3)**2
          e34=e3(1)*e4(1)+e3(2)*e4(2)+e3(3)*e4(3)
          IF(e34.LE.0d0) THEN
            ddmin=e4s
          ELSEIF(e34.LT.e3s) THEN
            ddmin=e4s-e34**2/e3s
          ELSE
            ddmin=e4s-2d0*e34+e3s
          ENDIF
 
C...Is this the smallest so far?
          IF(ddmin.LT.dglomi) THEN
            dglomi=ddmin
            ibeg=i0
            ipcs=i1
          ENDIF
        ELSEIF(k(i1,1).EQ.1.AND.kchg(pycomp(k(i1,2)),2).NE.0) THEN
          i0=0
        ENDIF
  940 CONTINUE
 
C... Check if there are any strings to connect to the new gluon. (EN)
      IF (ibeg.EQ.0) goto 1080
 
C...Delta_m = m_clus - m_had > 0: emit a 'gluon' (EN)
      IF (p(n+1,5).GE.p(n+2,5)) THEN
 
C...Construct 'gluon' that is needed to put hadron on the mass shell.
        frac=p(n+2,5)/p(n+1,5)
        DO 950 j=1,5
          p(n+2,j)=frac*p(n+1,j)
          pg(j)=(1d0-frac)*p(n+1,j)
  950   CONTINUE
 
C... Copy string with new gluon put in.
        n=n+2
        i=ibeg-1
  960   i=i+1
        IF(k(i,1).NE.1.AND.k(i,1).NE.2.AND.k(i,1).NE.41) goto 960
        IF(kchg(pycomp(k(i,2)),2).EQ.0.AND.k(i,1).NE.41) goto 960
        n=n+1
        DO 970 j=1,5
          k(n,j)=k(i,j)
          p(n,j)=p(i,j)
          v(n,j)=v(i,j)
  970   CONTINUE
        k(i,1)=k(i,1)+10
        k(i,4)=n
        k(i,5)=n
        k(n,3)=i
        IF(i.EQ.ipcs) THEN
          n=n+1
          DO 980 j=1,5
            k(n,j)=k(n-1,j)
            p(n,j)=pg(j)
            v(n,j)=v(n-1,j)
  980     CONTINUE
          k(n,2)=21
          k(n,3)=nsav+1
        ENDIF
        IF(k(i,1).EQ.12.OR.k(i,1).EQ.51) goto 960
        goto 1120
 
C...Delta_m = m_clus - m_had < 0: have to absorb a 'gluon' instead,
C...from string piece endpoints.
      ELSE
 
C...Begin by copying string that should give energy to cluster.
        n=n+2
        i=ibeg-1
  990   i=i+1
        IF(k(i,1).NE.1.AND.k(i,1).NE.2.AND.k(i,1).NE.41) goto 990
        IF(kchg(pycomp(k(i,2)),2).EQ.0.AND.k(i,1).NE.41) goto 990
        n=n+1
        DO 1000 j=1,5
          k(n,j)=k(i,j)
          p(n,j)=p(i,j)
          v(n,j)=v(i,j)
 1000   CONTINUE
        k(i,1)=k(i,1)+10
        k(i,4)=n
        k(i,5)=n
        k(n,3)=i
        IF(i.EQ.ipcs) i1=n
        IF(k(i,1).EQ.12.OR.k(i,1).EQ.51) goto 990
        i2=i1+1
 
C...Set initial Phad.
        DO 1010 j=1,4
          p(nsav+2,j)=p(nsav+1,j)
 1010   CONTINUE
 
C...Calculate Pg, a part of which will be added to Phad later. (EN)
 1020   IF(mstj(16).EQ.1) THEN
          alpha=1d0
          beta=1d0
        ELSE
          alpha=four(nsav+1,i2)/four(i1,i2)
          beta=four(nsav+1,i1)/four(i1,i2)
        ENDIF
        DO 1030 j=1,4
          pg(j)=alpha*p(i1,j)+beta*p(i2,j)
 1030   CONTINUE
        pg(5)=sqrt(max(1d-20,pg(4)**2-pg(1)**2-pg(2)**2-pg(3)**2))
 
C..Solve 2nd order equation, use the best (smallest) solution. (EN)
        pmscol=p(nsav+2,4)**2-p(nsav+2,1)**2-p(nsav+2,2)**2-
     &  p(nsav+2,3)**2
        pclpg=(p(nsav+2,4)*pg(4)-p(nsav+2,1)*pg(1)-
     &  p(nsav+2,2)*pg(2)-p(nsav+2,3)*pg(3))/pg(5)**2
        delta=sqrt(pclpg**2+(p(nsav+2,5)**2-pmscol)/pg(5)**2)-pclpg
 
C...If all gluon energy eaten, zero it and take a step back.
        iter=0
        IF(delta*alpha.GT.1d0.AND.i1.GT.nsav+3.AND.k(i1,2).EQ.21) THEN
          iter=1
          DO 1040 j=1,4
            p(nsav+2,j)=p(nsav+2,j)+p(i1,j)
            p(i1,j)=0d0
 1040     CONTINUE
          p(i1,5)=0d0
          k(i1,1)=k(i1,1)+10
          i1=i1-1
          IF(k(i1,1).EQ.41) iter=-1
        ENDIF
        IF(delta*beta.GT.1d0.AND.i2.LT.n.AND.k(i2,2).EQ.21) THEN
          iter=1
          DO 1050 j=1,4
            p(nsav+2,j)=p(nsav+2,j)+p(i2,j)
            p(i2,j)=0d0
 1050     CONTINUE
          p(i2,5)=0d0
          k(i2,1)=k(i2,1)+10
          i2=i2+1
          IF(k(i2,1).EQ.41) iter=-1
        ENDIF
        IF(iter.EQ.1) goto 1020
 
C...If also all endpoint energy eaten, revert to old procedure.
        IF((1d0-delta*alpha)*p(i1,4).LT.p(i1,5).OR.
     &  (1d0-delta*beta)*p(i2,4).LT.p(i2,5).OR.iter.EQ.-1) THEN
          DO 1060 i=nsav+3,n
            im=k(i,3)
            k(im,1)=k(im,1)-10
            k(im,4)=0
            k(im,5)=0
 1060     CONTINUE
          n=nsav
          goto 1080
        ENDIF
 
C... Construct the collapsed hadron and modified string partons.
        DO 1070 j=1,4
          p(nsav+2,j)=p(nsav+2,j)+delta*pg(j)
          p(i1,j)=(1d0-delta*alpha)*p(i1,j)
          p(i2,j)=(1d0-delta*beta)*p(i2,j)
 1070   CONTINUE
          p(i1,5)=(1d0-delta*alpha)*p(i1,5)
          p(i2,5)=(1d0-delta*beta)*p(i2,5)
 
C...Finished with string collapse in new scheme.
        goto 1120
      ENDIF
 
C... Use old algorithm; by choice or when in trouble.
 1080 CONTINUE
C...Find parton/particle which combines to largest extra mass.
      ir=0
      ha=0d0
      hsm=0d0
      DO 1100 mcomb=1,3
        IF(ir.NE.0) goto 1100
        DO 1090 i=max(1,ip),n
          IF(k(i,1).LE.0.OR.k(i,1).GT.10.OR.(i.GE.ic1.AND.i.LE.ic2
     &    .AND.k(i,1).GE.1.AND.k(i,1).LE.2)) goto 1090
          IF(mcomb.EQ.1) kci=pycomp(k(i,2))
          IF(mcomb.EQ.1.AND.kci.EQ.0) goto 1090
          IF(mcomb.EQ.1.AND.kchg(kci,2).EQ.0.AND.i.LE.ns) goto 1090
          IF(mcomb.EQ.2.AND.iabs(k(i,2)).GT.10.AND.iabs(k(i,2)).LE.100)
     &    goto 1090
          hcr=dpc(4)*p(i,4)-dpc(1)*p(i,1)-dpc(2)*p(i,2)-dpc(3)*p(i,3)
          hsr=2d0*hcr+pecm**2-p(n+2,5)**2-2d0*p(n+2,5)*p(i,5)
          IF(hsr.GT.hsm) THEN
            ir=i
            ha=hcr
            hsm=hsr
          ENDIF
 1090   CONTINUE
 1100 CONTINUE
 
C...Shuffle energy and momentum to put new particle on mass shell.
      IF(ir.NE.0) THEN
        hb=pecm**2+ha
        hc=p(n+2,5)**2+ha
        hd=p(ir,5)**2+ha
        hk2=0.5d0*(hb*sqrt(max(0d0,((hb+hc)**2-4d0*(hb+hd)*p(n+2,5)**2)/
     &  (ha**2-(pecm*p(ir,5))**2)))-(hb+hc))/(hb+hd)
        hk1=(0.5d0*(p(n+2,5)**2-pecm**2)+hd*hk2)/hb
        DO 1110 j=1,4
          p(n+2,j)=(1d0+hk1)*dpc(j)-hk2*p(ir,j)
          p(ir,j)=(1d0+hk2)*p(ir,j)-hk1*dpc(j)
 1110   CONTINUE
        n=n+2
      ELSE
        CALL pyerrm(3,'(PYPREP:) no match for collapsing cluster')
        RETURN
      ENDIF
 
C...Mark collapsed system and store daughter pointers. Iterate.
 1120 DO 1130 i=ic1,ic2
        IF((k(i,1).EQ.1.OR.k(i,1).EQ.2).AND.
     &  kchg(pycomp(k(i,2)),2).NE.0) THEN
          k(i,1)=k(i,1)+10
          IF(mstu(16).NE.2) THEN
            k(i,4)=nsav+1
            k(i,5)=nsav+1
          ELSE
            k(i,4)=nsav+2
            k(i,5)=nsav+1+nbody
          ENDIF
        ENDIF
        IF(k(i,1).EQ.41) k(i,1)=k(i,1)+10
 1130 CONTINUE
      IF(n.LT.mstu(4)-mstu(32)-5) goto 710
 
C...Check flavours and invariant masses in parton systems.
 1140 np=0
      kfn=0
      kqs=0
      nju=0
      DO 1150 j=1,5
        dps(j)=0d0
 1150 CONTINUE
      DO 1180 i=max(1,ip),n
        IF(k(i,1).EQ.41) nju=nju+1
        IF(k(i,1).LE.0.OR.k(i,1).GT.10) goto 1180
        kc=pycomp(k(i,2))
        IF(kc.EQ.0) goto 1180
        kq=kchg(kc,2)*isign(1,k(i,2))
        IF(kq.EQ.0) goto 1180
        np=np+1
        IF(kq.NE.2) THEN
          kfn=kfn+1
          kqs=kqs+kq
          mstj(93)=1
          dps(5)=dps(5)+pymass(k(i,2))
        ENDIF
        DO 1160 j=1,4
          dps(j)=dps(j)+p(i,j)
 1160   CONTINUE
        IF(k(i,1).EQ.1) THEN
          nferr=0
          IF(nju.EQ.0.AND.np.NE.1) THEN
            IF(kfn.EQ.1.OR.kfn.GE.3.OR.kqs.NE.0) nferr=1
          ELSEIF(nju.EQ.1) THEN
            IF(kfn.NE.3.OR.iabs(kqs).NE.3) nferr=1
          ELSEIF(nju.EQ.2) THEN
            IF(kfn.NE.4.OR.kqs.NE.0) nferr=1
          ELSEIF(nju.GE.3) THEN
            nferr=1
          ENDIF
          IF(nferr.EQ.1) THEN
            CALL pyerrm(2,'(PYPREP:) unphysical flavour combination')
            mint(51)=1
            RETURN
          ENDIF
          IF(np.NE.1.AND.dps(4)**2-dps(1)**2-dps(2)**2-dps(3)**2.LT.
     &    (0.9d0*parj(32)+dps(5))**2) CALL pyerrm(3,
     &    '(PYPREP:) too small mass in jet system')
          np=0
          kfn=0
          kqs=0
          nju=0
          DO 1170 j=1,5
            dps(j)=0d0
 1170     CONTINUE
        ENDIF
 1180 CONTINUE
 
      RETURN
      END