pymirm.f

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

 
C*********************************************************************
 
C...PYMIRM
C...Picks primordial kT and shares longitudinal momentum among
C...beam remnants.
 
      SUBROUTINE pymirm
 
C...Double precision and integer declarations.
      IMPLICIT DOUBLE PRECISION(a-h, o-z)
      IMPLICIT INTEGER(i-n)
      INTEGER pyk,pychge,pycomp
C...The event record
      common/pyjets/n,npad,k(4000,5),p(4000,5),v(4000,5)
C...Parameters
      common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
      common/pypars/mstp(200),parp(200),msti(200),pari(200)
      common/pyint1/mint(400),vint(400)
C...The common block of colour tags.
      common/pyctag/nct,mct(4000,2)
C...The common block of dangling ends
      common/pyintm/kfival(2,3),nmi(2),imi(2,800,2),nvc(2,-6:6),
     &     xassoc(2,-6:6,240),xpsvc(-6:6,-1:240),pvctot(2,-1:1),
     &     xmi(2,240),pt2mi(240),imisep(0:240)
      SAVE /pyjets/,/pydat1/,/pypars/,/pyint1/,/pyintm/,/pyctag/
C...Local variables
      dimension w(0:2,0:2),vb(3),nnxt(2),ivalq(2),icomq(2)
C...W(I,J)|  J=0    |   1   |   2   |
C...  I=0 | Wrem**2 |  W+   |  W-   |
C...    1 | W1**2   |  W1+  |  W1-  |
C...    2 | W2**2   |  W2+  |  W2-  |
C...4-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...Tentative parametrization of <kT> as a function of Q.
      sigpt(q)=max(parj(21),2.1d0*q/(7d0+q))
C      SIGPT(Q)=MAX(0.36D0,4D0*SQRT(Q)/(10D0+SQRT(Q))
C      SIGPT(Q)=MAX(PARJ(21),3D0*SQRT(Q)/(5D0+SQRT(Q))
      getpt(q,sigma)=min(sigma*sqrt(-log(pyr(0))),parp(93))
C...Lambda kinematic function.
      flam(a,b,c)=a**2+b**2+c**2-2d0*(a*b+b*c+c*a)
 
C...Beginning and end of beam remnant partons
      nout=mint(53)
      isub=mint(1)
 
C...Loopback point if kinematic choices gives impossible configuration.
      ntry=0
  100 ntry=ntry+1
 
C...Assign kT values on each side separately.
      DO 180 js=1,2
 
C...First zero all kT on this side. Skip if no kT to generate.
        DO 110 im=1,nmi(js)
          p(imi(js,im,1),1)=0d0
          p(imi(js,im,1),2)=0d0
  110   CONTINUE
        IF(mstp(91).LE.0) goto 180
 
C...Now assign kT to each (non-collapsed) parton in IMI.
        DO 170 im=1,nmi(js)
          i=imi(js,im,1)
C...Select kT according to truncated gaussian or 1/kt6 tails.
C...For first interaction, either use rms width = PARP(91) or fitted.
          IF (im.EQ.1) THEN
            sigma=parp(91)
            IF (mstp(91).GE.11.AND.mstp(91).LE.20) THEN
              q=sqrt(pt2mi(im))
              sigma=sigpt(q)
            ENDIF
          ELSE
C...For subsequent interactions and BR partons use fragmentation width.
            sigma=parj(21)
          ENDIF
          phi=paru(2)*pyr(0)
          pt=0d0
          IF(ntry.LE.100) THEN
 111        IF (mstp(91).EQ.1.OR.mstp(91).EQ.11) THEN
              pt=getpt(q,sigma)
              ptx=pt*cos(phi)
              pty=pt*sin(phi)
            ELSEIF (mstp(91).EQ.2) THEN
              CALL pyerrm(11,'(PYMIRM:) Sorry, MSTP(91)=2 not '//
     &          'available, using MSTP(91)=1.')
              CALL pygive('MSTP(91)=1')
              goto 111
            ELSEIF(mstp(91).EQ.3.OR.mstp(91).EQ.13) THEN
C...Use distribution with kt**6 tails, rms width = PARP(91).
              eps=sqrt(3d0/2d0)*sigma
C...Generate PTX and PTY separately, each propto 1/KT**6
              DO 119 ixy=1,2
C...Decide which interval to try
 112            p12=1d0/(1d0+27d0/40d0*sigma**6/eps**6)
                IF (pyr(0).LT.p12) THEN
C...Use flat approx with accept/reject up to EPS.
                  pt=pyr(0)*eps
                  wt=(3d0/2d0*sigma**2/(pt**2+3d0/2d0*sigma**2))**3
                  IF (pyr(0).GT.wt) goto 112
                ELSE
C...Above EPS, use 1/kt**6 approx with accept/reject.
                  pt=eps/(pyr(0)**(1d0/5d0))
                  wt=pt**6/(pt**2+3d0/2d0*sigma**2)**3
                  IF (pyr(0).GT.wt) goto 112
                ENDIF
                msign=1
                IF (pyr(0).GT.0.5d0) msign=-1
                IF (ixy.EQ.1) ptx=msign*pt
                IF (ixy.EQ.2) pty=msign*pt
 119          CONTINUE
            ELSEIF (mstp(91).EQ.4.OR.mstp(91).EQ.14) THEN
              ptx=sigma*(sqrt(6d0)*pyr(0)-sqrt(3d0/2d0))
              pty=sigma*(sqrt(6d0)*pyr(0)-sqrt(3d0/2d0))
            ENDIF
C...Adjust final PT. Impose upper cutoff, or zero for soft evts.
            pt=sqrt(ptx**2+pty**2)
            wt=1d0
            IF (pt.GT.parp(93)) wt=sqrt(parp(93)/pt)
            IF(isub.EQ.95.AND.im.EQ.1) wt=0d0
            ptx=ptx*wt
            pty=pty*wt
            pt=sqrt(ptx**2+pty**2)
          ENDIF
 
          p(i,1)=p(i,1)+ptx
          p(i,2)=p(i,2)+pty
 
C...Compensation kicks, with varying degree of local anticorrelations.
          mcorr=mstp(90)
          IF (mcorr.EQ.0.OR.isub.EQ.95) THEN
            ptcx=-ptx/(nmi(js)-1)
            ptcy=-pty/(nmi(js)-1)
            IF(isub.EQ.95) THEN
              ptcx=-ptx/(nmi(js)-2)
              ptcy=-pty/(nmi(js)-2)
            ENDIF
            DO 120 imc=1,nmi(js)
              IF (imc.EQ.im) goto 120
              IF(isub.EQ.95.AND.imc.EQ.1) goto 120
              p(imi(js,imc,1),1)=p(imi(js,imc,1),1)+ptcx
              p(imi(js,imc,1),2)=p(imi(js,imc,1),2)+ptcy
  120       CONTINUE
          ELSEIF (mcorr.GE.1) THEN
            DO 140 msid=4,5
              nnxt(msid-3)=0
C...Count up # of neighbours on either side
              imo=i
  130         imo=k(imo,msid)/mstu(5)
              IF (imo.EQ.0) goto 140
              nnxt(msid-3)=nnxt(msid-3)+1
C...Stop at quarks and junctions
              IF (mcorr.EQ.1.AND.k(imo,2).EQ.21) goto 130
  140       CONTINUE
C...How should compensation be shared when unequal numbers on the
C...two sides? 50/50 regardless? N1:N2? Assume latter for now.
            nsum=nnxt(1)+nnxt(2)
            t1=0
            DO 160 msid=4,5
C...Total momentum to be compensated on this side
              IF (nnxt(msid-3).EQ.0) goto 160
              ptcx=-(nnxt(msid-3)*ptx)/nsum
              ptcy=-(nnxt(msid-3)*pty)/nsum
C...RS: compensation supression factor as we go out from parton I.
C...Hardcoded behaviour RS=0.5, i.e. 1/2**n falloff,
C...since (for now) MSTP(90) provides enough variability.
              rs=0.5d0
              fac=(1d0-rs)/(rs*(1-rs**nnxt(msid-3)))
              imo=i
  150         ida=imo
              imo=k(imo,msid)/mstu(5)
              IF (imo.EQ.0) goto 160
              fac=fac*rs
              IF (k(imo,2).NE.88) THEN
                p(imo,1)=p(imo,1)+fac*ptcx
                p(imo,2)=p(imo,2)+fac*ptcy
                IF (mcorr.EQ.1.AND.k(imo,2).EQ.21) goto 150
C...If we reach junction, divide out the kT that would have been
C...assigned to the junction on each of its other legs.
              ELSE
                l1=mod(k(imo,4),mstu(5))
                l2=k(imo,5)/mstu(5)
                l3=mod(k(imo,5),mstu(5))
                p(l1,1)=p(l1,1)+0.5d0*fac*ptcx
                p(l1,2)=p(l1,2)+0.5d0*fac*ptcy
                p(l2,1)=p(l2,1)+0.5d0*fac*ptcx
                p(l2,2)=p(l2,2)+0.5d0*fac*ptcy
                p(l3,1)=p(l3,1)+0.5d0*fac*ptcx
                p(l3,2)=p(l3,2)+0.5d0*fac*ptcy
                p(ida,1)=p(ida,1)-0.5d0*fac*ptcx
                p(ida,2)=p(ida,2)-0.5d0*fac*ptcy
              ENDIF
 
  160       CONTINUE
          ENDIF
  170   CONTINUE
C...End assignment of kT values to initiators and remnants.
  180 CONTINUE
 
C...Check kinematics constraints for non-BR partons.
      DO 190 im=1,mint(31)
        shat=xmi(1,im)*xmi(2,im)*vint(2)
        pt1=sqrt(p(imi(1,im,1),1)**2+p(imi(1,im,1),2)**2)
        pt2=sqrt(p(imi(2,im,1),1)**2+p(imi(2,im,1),2)**2)
        pt1pt2=p(imi(1,im,1),1)*p(imi(2,im,1),1)
     &        +p(imi(1,im,1),2)*p(imi(2,im,1),2)
        IF (shat.LT.2d0*(pt1*pt2-pt1pt2).AND.ntry.LE.100) THEN
          IF(ntry.GE.100) THEN
C...Kill this event and start another.
            CALL pyerrm(11,
     &           '(PYMIRM:) No consistent (x,kT) sets found')
            mint(51)=1
            RETURN
          ENDIF
          goto 100
        ENDIF
  190 CONTINUE
 
C...Calculate W+ and W- available for combined remnant system.
      w(0,1)=vint(1)
      w(0,2)=vint(1)
      DO 200 im=1,mint(31)
        pt2 = (p(imi(1,im,1),1)+p(imi(2,im,1),1))**2
     &       +(p(imi(1,im,1),2)+p(imi(2,im,1),2))**2
        st=xmi(1,im)*xmi(2,im)*vint(2)+pt2
        w(0,1)=w(0,1)-sqrt(xmi(1,im)/xmi(2,im)*st)
        w(0,2)=w(0,2)-sqrt(xmi(2,im)/xmi(1,im)*st)
  200 CONTINUE
C...Also store Wrem**2 = W+ * W-
      w(0,0)=w(0,1)*w(0,2)
 
      IF (w(0,0).LT.0d0.AND.ntry.LE.100) THEN
          IF(ntry.GE.100) THEN
C...Kill this event and start another.
            CALL pyerrm(11,
     &    '(PYMIRM:) Negative beam remnant mass squared unavoidable')
            mint(51)=1
            RETURN
          ENDIF
          goto 100
      ENDIF
 
C...Assign unscaled x values to partons/hadrons in each of the
C...beam remnants and calculate unscaled W+ and W- from them.
      ntryx=0
  210 ntryx=ntryx+1
      DO 280 js=1,2
        w(js,1)=0d0
        w(js,2)=0d0
        DO 270 im=mint(31)+1,nmi(js)
          i=imi(js,im,1)
          kf=k(i,2)
          kfa=iabs(kf)
          icomp=imi(js,im,2)
 
C...Skip collapsed gluons and junctions. Reset.
          IF (kfa.EQ.21.AND.k(i,1).EQ.14) goto 270
          IF (kfa.EQ.88) goto 270
          x=0d0
          ivalq(1)=0
          ivalq(2)=0
          icomq(1)=0
          icomq(2)=0
 
C...If gluon then only beam remnant, so takes all.
          IF(kfa.EQ.21) THEN
            x=1d0
C...If valence quark then use parametrized valence distribution.
          ELSEIF(kfa.LE.6.AND.icomp.EQ.0) THEN
            ivalq(1)=kf
C...If companion quark then derive from companion x.
          ELSEIF(kfa.LE.6) THEN
            icomq(1)=icomp
C...If valence diquark then use two parametrized valence distributions.
          ELSEIF(kfa.GT.1000.AND.mod(kfa/10,10).EQ.0.AND.
     &    icomp.EQ.0) THEN
            ivalq(1)=isign(kfa/1000,kf)
            ivalq(2)=isign(mod(kfa/100,10),kf)
C...If valence+sea diquark then combine valence + companion choices.
          ELSEIF(kfa.GT.1000.AND.mod(kfa/10,10).EQ.0.AND.
     &    icomp.LT.mstu(5)) THEN
            IF(kfa/1000.EQ.iabs(k(icomp,2))) THEN
              ivalq(1)=isign(mod(kfa/100,10),kf)
            ELSE
              ivalq(1)=isign(kfa/1000,kf)
            ENDIF
            icomq(1)=icomp
C...Extra code: workaround for diquark made out of two sea
C...quarks, but where not (yet) ICOMP > MSTU(5).
            DO 220 im1=1,mint(31)
              IF(imi(js,im1,2).EQ.i.AND.imi(js,im1,1).NE.icomp) THEN
                icomq(2)=imi(js,im1,1)
                ivalq(1)=0
              ENDIF
  220       CONTINUE
C...If sea diquark then sum of two derived from companion x.
          ELSEIF(kfa.GT.1000.AND.mod(kfa/10,10).EQ.0) THEN
             icomq(1)=mod(icomp,mstu(5))
             icomq(2)=icomp/mstu(5)
C...If meson or baryon then use fragmentation function.
C...Somewhat arbitrary split into old and new flavour, but OK normally.
          ELSE
            kfl3=mod(kfa/10,10)
            IF(mod(kfa/1000,10).EQ.0) THEN
              kfl1=mod(kfa/100,10)
            ELSE
              kfl1=mod(kfa,10000)-10*kfl3-1
              IF(mod(kfa/1000,10).EQ.mod(kfa/100,10).AND.
     &        mod(kfa,10).EQ.2) kfl1=kfl1+2
            ENDIF
            pr=p(i,5)**2+p(i,1)**2+p(i,2)**2
            CALL pyzdis(kfl1,kfl3,pr,x)
          ENDIF
 
          DO 260 iq=1,2
C...Calculation of x of valence quark: assume form (1-x)^a/sqrt(x),
C...where a=3.5 for u in proton, =2 for d in proton and =0.8 for meson.
C...In other baryons combine u and d from proton appropriately.
            IF(ivalq(iq).NE.0) THEN
              nval=0
              IF(kfival(js,1).EQ.ivalq(iq)) nval=nval+1
              IF(kfival(js,2).EQ.ivalq(iq)) nval=nval+1
              IF(kfival(js,3).EQ.ivalq(iq)) nval=nval+1
C...Meson.
              IF(kfival(js,3).EQ.0) THEN
                mdu=0
C...Baryon with three identical quarks: mix u and d forms.
              ELSEIF(nval.EQ.3) THEN
                mdu=int(pyr(0)+5d0/3d0)
C...Baryon, one of two identical quarks: u form.
              ELSEIF(nval.EQ.2) THEN
                mdu=2
C...Baryon with two identical quarks, but not the one picked: d form.
              ELSEIF(kfival(js,1).EQ.kfival(js,2).OR.kfival(js,2).EQ.
     &        kfival(js,3).OR.kfival(js,1).EQ.kfival(js,3)) THEN
                mdu=1
C...Baryon with three nonidentical quarks: mix u and d forms.
              ELSE
                mdu=int(pyr(0)+5d0/3d0)
              ENDIF
              xpow=0.8d0
              IF(mdu.EQ.1) xpow=3.5d0
              IF(mdu.EQ.2) xpow=2d0
  230         xx=pyr(0)**2
              IF((1d0-xx)**xpow.LT.pyr(0)) goto 230
              x=x+xx
            ENDIF
 
C...Calculation of x of companion quark.
            IF(icomq(iq).NE.0) THEN
              xcomp=1d-4
              DO 240 im1=1,mint(31)
                IF(imi(js,im1,1).EQ.icomq(iq)) xcomp=xmi(js,im1)
  240         CONTINUE
              npow=max(0,min(4,mstp(87)))
  250         xx=xcomp*(1d0/(1d0-pyr(0)*(1d0-xcomp))-1d0)
              corr=((1d0-xcomp-xx)/(1d0-xcomp))**npow*
     &        (xcomp**2+xx**2)/(xcomp+xx)**2
              IF(corr.LT.pyr(0)) goto 250
              x=x+xx
            ENDIF
  260     CONTINUE
 
C...Optionally enchance x of composite systems (e.g. diquarks)
          IF (kfa.GT.100) x=parp(79)*x
 
C...Store x. Also calculate light cone energies of each system.
          xmi(js,im)=x
          w(js,js)=w(js,js)+x
          w(js,3-js)=w(js,3-js)+(p(i,5)**2+p(i,1)**2+p(i,2)**2)/x
  270   CONTINUE
        w(js,js)=w(js,js)*w(0,js)
        w(js,3-js)=w(js,3-js)/w(0,js)
        w(js,0)=w(js,1)*w(js,2)
  280 CONTINUE
 
C...Check W1 W2 < Wrem (can be done before rescaling, since W
C...insensitive to global rescalings of the BR x values).
      IF (sqrt(w(1,0))+sqrt(w(2,0)).GT.sqrt(w(0,0)).AND.ntryx.LE.100)
     &     THEN
        goto 210
      ELSEIF (ntryx.GT.100.AND.ntry.LE.100) THEN
        goto 100
      ELSEIF (ntryx.GT.100) THEN
        CALL pyerrm(11,'(PYMIRM:) No consistent (x,kT) sets found')
        mint(57)=mint(57)+1
        mint(51)=1
        RETURN
      ENDIF
 
C...Compute x rescaling factors
      comtrm=w(0,0)+sqrt(flam(w(0,0),w(1,0),w(2,0)))
      r1=(comtrm+w(1,0)-w(2,0))/(2d0*w(1,1)*w(0,2))
      r2=(comtrm+w(2,0)-w(1,0))/(2d0*w(2,2)*w(0,1))
 
      IF (r1.LT.0.OR.r2.LT.0) THEN
        CALL pyerrm(19,'(PYMIRM:) negative rescaling factors !')
        mint(57)=mint(57)+1
        mint(51)=1
      ENDIF
 
C...Rescale W(1,*) and W(2,*) (not really necessary, but consistent).
      w(1,1)=w(1,1)*r1
      w(1,2)=w(1,2)/r1
      w(2,1)=w(2,1)/r2
      w(2,2)=w(2,2)*r2
 
C...Rescale BR x values.
      DO 290 im=mint(31)+1,max(nmi(1),nmi(2))
        xmi(1,im)=xmi(1,im)*r1
        xmi(2,im)=xmi(2,im)*r2
  290 CONTINUE
 
C...Now we have a consistent set of x and kT values.
C...First set up the initiators and their daughters correctly.
      DO 300 im=1,mint(31)
        i1=imi(1,im,1)
        i2=imi(2,im,1)
        st=xmi(1,im)*xmi(2,im)*vint(2)+(p(i1,1)+p(i2,1))**2+
     &       (p(i1,2)+p(i2,2))**2
        pt12=p(i1,1)**2+p(i1,2)**2
        pt22=p(i2,1)**2+p(i2,2)**2
C...p_z
        p(i1,3)=sqrt(flam(st,pt12,pt22)/(4d0*st))
        p(i2,3)=-p(i1,3)
C...Energies (masses should be zero at this stage)
        p(i1,4)=sqrt(pt12+p(i1,3)**2)
        p(i2,4)=sqrt(pt22+p(i2,3)**2)
 
C...Transverse 12 system initiator velocity:
        vb(1)=(p(i1,1)+p(i2,1))/sqrt(st)
        vb(2)=(p(i1,2)+p(i2,2))/sqrt(st)
C...Boost to overall initiator system rest frame
        CALL pyrobo(i1,i1,0d0,0d0,-vb(1),-vb(2),0d0)
        CALL pyrobo(i2,i2,0d0,0d0,-vb(1),-vb(2),0d0)
C...Compute phi,theta coordinates of I1 and rotate z axis.
        phi=pyangl(p(i1,1),p(i1,2))
        the=pyangl(p(i1,3),sqrt(p(i1,1)**2+p(i1,2)**2))
        CALL pyrobo(i1,i1,0d0,-phi,0d0,0d0,0d0)
        CALL pyrobo(i2,i2,0d0,-phi,0d0,0d0,0d0)
        CALL pyrobo(i1,i1,-the,0d0,0d0,0d0,0d0)
        CALL pyrobo(i2,i2,-the,0d0,0d0,0d0,0d0)
 
C...Now boost initiators + daughters back to LAB system
C...(also update documentation lines for MI = 1.)
        vb(3)=(xmi(1,im)-xmi(2,im))/(xmi(1,im)+xmi(2,im))
        imin=imisep(im-1)+1
        IF (im.EQ.1) imin=mint(83)+5
        imax=imisep(im)
        CALL pyrobo(imin,imax,the,phi,vb(1),vb(2),0d0)
        CALL pyrobo(imin,imax,0d0,0d0,0d0,0d0,vb(3))
 
  300 CONTINUE
 
 
C...For the beam remnant partons/hadrons, we only need to set pz and E.
      DO 320 js=1,2
        DO 310 im=mint(31)+1,nmi(js)
          i=imi(js,im,1)
C...Skip collapsed gluons and junctions.
          IF (k(i,2).EQ.21.AND.k(i,1).EQ.14) goto 310
          IF (kfa.EQ.88) goto 310
          rmt2=p(i,5)**2+p(i,1)**2+p(i,2)**2
          p(i,4)=0.5d0*(xmi(js,im)*w(0,js)+rmt2/(xmi(js,im)*w(0,js)))
          p(i,3)=0.5d0*(xmi(js,im)*w(0,js)-rmt2/(xmi(js,im)*w(0,js)))
          IF (js.EQ.2) p(i,3)=-p(i,3)
  310   CONTINUE
  320 CONTINUE
 
 
C...Documentation lines
      DO 340 js=1,2
        in=mint(83)+js+2
        io=imi(js,1,1)
        k(in,1)=21
        k(in,2)=k(io,2)
        k(in,3)=mint(83)+js
        k(in,4)=0
        k(in,5)=0
        DO 330 j=1,5
          p(in,j)=p(io,j)
          v(in,j)=v(io,j)
  330   CONTINUE
        mct(in,1)=mct(io,1)
        mct(in,2)=mct(io,2)
  340 CONTINUE
 
C...Final state colour reconnections.
      IF (mstp(95).NE.1.OR.mint(31).LE.1) goto 380
 
C...Number of colour tags for which a recoupling will be tried.
      ntot=nct
C...Number of recouplings to try
      mint(34)=0
      nrecp=0
      niter=0
  350 nrecp=mint(34)
      niter=niter+1
      iiter=0
  360 iiter=iiter+1
      IF (iiter.LE.parp(78)*ntot) THEN
C...Select two colour tags at random
C...NB: jj strings do not have colour tags assigned to them,
C...thus they are as yet not affected by anything done here.
        jct=pyr(0)*nct+1
        kct=mod(int(jct+pyr(0)*nct),nct)+1
        ij1=0
        ij2=0
        ik1=0
        ik2=0
C...Find final state partons with this (anti)colour
        DO 370 i=mint(84)+1,n
          IF (k(i,1).EQ.3) THEN
            IF (mct(i,1).EQ.jct) ij1=i
            IF (mct(i,2).EQ.jct) ij2=i
            IF (mct(i,1).EQ.kct) ik1=i
            IF (mct(i,2).EQ.kct) ik2=i
          ENDIF
  370   CONTINUE
C...Only consider recouplings not involving junctions for now.
        IF (ij1.EQ.0.OR.ij2.EQ.0.OR.ik1.EQ.0.OR.ik2.EQ.0) goto 360
 
        rlo=2d0*four(ij1,ij2)*2d0*four(ik1,ik2)
        rln=2d0*four(ij1,ik2)*2d0*four(ik1,ij2)
        IF (rln.LT.rlo.AND.mct(ij2,1).NE.kct.AND.mct(ik2,1).NE.jct) THEN
          mct(ij2,2)=kct
          mct(ik2,2)=jct
C...Count up number of reconnections
          mint(34)=mint(34)+1
        ENDIF
        IF (mint(34).LE.1000) THEN
          goto 360
        ELSE
          CALL pyerrm(4,'(PYMIRM:) caught in infinite loop')
          goto 380
        ENDIF
      ENDIF
      IF (nrecp.LT.mint(34)) goto 350
 
C...Signal PYPREP to use /PYCTAG/ information rather than K(I,KCS).
  380 mint(33)=1
 
      RETURN
      END