pythiaeRHIC/pythia/pymihk.f

0001
0002 C*********************************************************************
0003
0004 C...PYMIHK
0005 C...Finds left-behind remnant flavour content and hooks up
0006 C...the colour flow between the hard scattering and remnants
0007
0008       SUBROUTINE PYMIHK
0009
0010 C...Double precision and integer declarations.
0011       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0012       IMPLICIT INTEGER(I-N)
0013       INTEGER PYK,PYCHGE,PYCOMP
0014 C...The event record
0015       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
0016 C...Parameters
0017       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0018       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0019       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0020       COMMON/PYINT1/MINT(400),VINT(400)
0021 C...The common block of dangling ends
0022       COMMON/PYINTM/KFIVAL(2,3),NMI(2),IMI(2,800,2),NVC(2,-6:6),
0023      &     XASSOC(2,-6:6,240),XPSVC(-6:6,-1:240),PVCTOT(2,-1:1),
0024      &     XMI(2,240),PT2MI(240),IMISEP(0:240)
0025       SAVE /PYJETS/,/PYDAT1/,/PYDAT2/,/PYPARS/,/PYINT1/,/PYINTM/
0026 C...Local variables
0027       PARAMETER (NERSIZ=4000)
0028       COMMON /PYCBLS/MCO(NERSIZ,2),NCC,JCCO(NERSIZ,2),JCCN(NERSIZ,2)
0029      &     ,MACCPT
0030       COMMON /PYCTAG/NCT,MCT(NERSIZ,2)
0031       SAVE /PYCBLS/,/PYCTAG/
0032       DIMENSION JST(2,3),IV(2,3),IDQ(3),NVSUM(2),NBRTOT(2),NG(2)
0033      &     ,ITJUNC(2),MOUT(2),INSR(1000,3),ISTR(6),YMI(240)
0034       DATA NERRPR/0/
0035       SAVE NERRPR
0036       FOUR(I,J)=P(I,4)*P(J,4)-P(I,3)*P(J,3)-P(I,2)*P(J,2)-P(I,1)*P(J,1)
0037
0038 C...Set up error checkers
0039       IBOOST=0
0040
0041 C...Initialize colour arrays: MCO (Original) and MCT (New)
0042       DO 110 I=MINT(84)+1,NERSIZ
0043         DO 100 JC=1,2
0044           MCT(I,JC)=0
0045           MCO(I,JC)=0
0046   100   CONTINUE
0047 C...Also zero colour tracing information, if existed.
0048         IF (I.LE.N) THEN
0049           K(I,4)=MOD(K(I,4),MSTU(5)**2)
0050           K(I,5)=MOD(K(I,5),MSTU(5)**2)
0051         ENDIF
0052   110 CONTINUE
0053
0054 C...Initialize colour tag collapse arrays:
0055 C...JCCO (Original) and JCCN (New).
0056       DO 130 MG=MINT(84)+1,NERSIZ
0057         DO 120 JC=1,2
0058           JCCO(MG,JC)=0
0059           JCCN(MG,JC)=0
0060   120   CONTINUE
0061   130 CONTINUE
0062
0063 C...Zero gluon insertion array
0064       DO 150 IM=1,1000
0065         DO 140 J=1,3
0066           INSR(IM,J)=0
0067   140   CONTINUE
0068   150 CONTINUE
0069
0070 C...Compute hard scattering system rapidities
0071       IF (MSTP(89).EQ.1) THEN
0072         DO 160 IM=1,240
0073           IF (IM.LE.MINT(31)) THEN
0074             YMI(IM)=LOG(XMI(1,IM)/XMI(2,IM))
0075           ELSE
0076 C...Set (unsigned) rapidity = 100 for beam remnant systems.
0077             YMI(IM)=100D0
0078           ENDIF
0079   160   CONTINUE
0080       ENDIF
0081
0082 C...Treat each side separately
0083       DO 290 JS=1,2
0084
0085 C...Initialize side.
0086         NG(JS)=0
0087         JV=0
0088         KFS=ISIGN(1,MINT(10+JS))
0089
0090 C...Set valence content of pi0, gamma, K0S, K0L if not yet done.
0091         IF(KFIVAL(JS,1).EQ.0) THEN
0092           IF(MINT(10+JS).EQ.111) THEN
0093             KFIVAL(JS,1)=INT(1.5D0+PYR(0))
0094             KFIVAL(JS,2)=-KFIVAL(JS,1)
0095           ELSEIF(MINT(10+JS).EQ.22) THEN
0096             PYRKF=PYR(0)
0097             KFIVAL(JS,1)=1
0098             IF(PYRKF.GT.0.1D0) KFIVAL(JS,1)=2
0099             IF(PYRKF.GT.0.5D0) KFIVAL(JS,1)=3
0100             IF(PYRKF.GT.0.6D0) KFIVAL(JS,1)=4
0101             KFIVAL(JS,2)=-KFIVAL(JS,1)
0102           ELSEIF(MINT(10+JS).EQ.130.OR.MINT(10+JS).EQ.310) THEN
0103             IF(PYR(0).GT.0.5D0) THEN
0104               KFIVAL(JS,1)=1
0105               KFIVAL(JS,2)=-3
0106             ELSE
0107               KFIVAL(JS,1)=3
0108               KFIVAL(JS,2)=-1
0109             ENDIF
0110           ENDIF
0111         ENDIF
0112
0113 C...Initialize beam remnant sea and valence content flavour by flavour.
0114         NVSUM(JS)=0
0115         NBRTOT(JS)=0
0116         DO 210 JFA=1,6
0117 C...Count up original number of JFA valence quarks and antiquarks.
0118           NVALQ=0
0119           NVALQB=0
0120           NSEA=0
0121           DO 170 J=1,3
0122             IF(KFIVAL(JS,J).EQ.JFA) NVALQ=NVALQ+1
0123             IF(KFIVAL(JS,J).EQ.-JFA) NVALQB=NVALQB+1
0124   170     CONTINUE
0125           NVSUM(JS)=NVSUM(JS)+NVALQ+NVALQB
0126 C...Subtract kicked out valence and determine sea from flavour cons.
0127           DO 180 IM=1,NMI(JS)
0128             IFL = K(IMI(JS,IM,1),2)
0129             IFA = IABS(IFL)
0130             IFS = ISIGN(1,IFL)
0131             IF (IFL.EQ.JFA.AND.IMI(JS,IM,2).EQ.0) THEN
0132 C...Subtract K.O. valence quark from remainder.
0133               NVALQ=NVALQ-1
0134               JV=NVSUM(JS)-NVALQ-NVALQB
0135               IV(JS,JV)=IMI(JS,IM,1)
0136             ELSEIF (IFL.EQ.-JFA.AND.IMI(JS,IM,2).EQ.0) THEN
0137 C...Subtract K.O. valence antiquark from remainder.
0138               NVALQB=NVALQB-1
0139               JV=NVSUM(JS)-NVALQ-NVALQB
0140               IV(JS,JV)=IMI(JS,IM,1)
0141             ELSEIF (IFA.EQ.JFA) THEN
0142 C...Outside sea without companion: add opposite sea flavour inside.
0143               IF (IMI(JS,IM,2).LT.0) NSEA=NSEA-IFS
0144             ENDIF
0145   180     CONTINUE
0146 C...Check if space left in PYJETS for additional BR flavours
0147           NFLSUM=IABS(NSEA)+NVALQ+NVALQB
0148           NBRTOT(JS)=NBRTOT(JS)+NFLSUM
0149           IF (N+NFLSUM+1.GT.MSTU(4)) THEN
0150             CALL PYERRM(11,'(PYMIHK:) no more memory left in PYJETS')
0151             MINT(51)=1
0152             RETURN
0153           ENDIF
0154 C...Add required val+sea content to beam remnant.
0155           IF (NFLSUM.GT.0) THEN
0156             DO 200 IA=1,NFLSUM
0157 C...Insert beam remnant quark as p.t. symbolic parton in ER.
0158               N=N+1
0159               DO 190 IX=1,5
0160                 K(N,IX)=0
0161                 P(N,IX)=0D0
0162                 V(N,IX)=0D0
0163   190         CONTINUE
0164               K(N,1)=3
0165               K(N,2)=ISIGN(JFA,NSEA)
0166               IF (IA.LE.NVALQ) K(N,2)=JFA
0167               IF (IA.GT.NVALQ.AND.IA.LE.NVALQ+NVALQB) K(N,2)=-JFA
0168               K(N,3)=MINT(83)+JS
0169 C...Also update NMI, IMI, and IV arrays.
0170               NMI(JS)=NMI(JS)+1
0171               IMI(JS,NMI(JS),1)=N
0172               IMI(JS,NMI(JS),2)=-1
0173               IF (IA.LE.NVALQ+NVALQB) THEN
0174                 IMI(JS,NMI(JS),2)=0
0175                 JV=JV+1
0176                 IV(JS,JV)=IMI(JS,NMI(JS),1)
0177               ENDIF
0178   200       CONTINUE
0179           ENDIF
0180   210   CONTINUE
0181
0182         IM=0
0183   220   IM=IM+1
0184         IF (IM.LE.NMI(JS)) THEN
0185           IF (K(IMI(JS,IM,1),2).EQ.21) THEN
0186             NG(JS)=NG(JS)+1
0187 C...Add fictitious parent gluons for companion pairs.
0188           ELSEIF (IMI(JS,IM,2).NE.0.AND.K(IMI(JS,IM,1),2).GT.0) THEN
0189 C...Randomly assign companions to sea quarks which have none.
0190             IF (IMI(JS,IM,2).LT.0) THEN
0191               IMC=PYR(0)*NMI(JS)
0192   230         IMC=MOD(IMC,NMI(JS))+1
0193               IF (K(IMI(JS,IMC,1),2).NE.-K(IMI(JS,IM,1),2)) GOTO 230
0194               IF (IMI(JS,IMC,2).GE.0) GOTO 230
0195               IMI(JS, IM,2) = IMI(JS,IMC,1)
0196               IMI(JS,IMC,2) = IMI(JS, IM,1)
0197             ENDIF
0198 C...Add fictitious parent gluon
0199             N=N+1
0200             DO 240 IX=1,5
0201               K(N,IX)=0
0202               P(N,IX)=0D0
0203               V(N,IX)=0D0
0204   240       CONTINUE
0205             K(N,1)=14
0206             K(N,2)=21
0207             K(N,3)=MINT(83)+JS
0208 C...Set gluon (anti-)colour daughter pointers
0209             K(N,4)=IMI(JS, IM,1)
0210             K(N,5)=IMI(JS, IM,2)
0211 C...Set quark (anti-)colour parent pointers
0212             K(IMI(JS, IM,2),5)=K(IMI(JS, IM,2),5)+MSTU(5)*N
0213             K(IMI(JS, IM,1),4)=K(IMI(JS, IM,1),4)+MSTU(5)*N
0214 C...Add gluon to IMI
0215             NMI(JS)=NMI(JS)+1
0216             IMI(JS,NMI(JS),1)=N
0217             IMI(JS,NMI(JS),2)=0
0218           ENDIF
0219           GOTO 220
0220         ENDIF
0221
0222 C...If incoming (anti-)baryon, insert inside (anti-)junction.
0223 C...Set up initial v-v-j-v configuration. Otherwise set up
0224 C...mesonic v-vbar configuration
0225         IF (IABS(MINT(10+JS)).GT.1000) THEN
0226 C...Determine junction type (1: B=1 2: B=-1)
0227           ITJUNC(JS) = (3-KFS)/2
0228 C...Insert junction.
0229           N=N+1
0230           DO 250 IX=1,5
0231             K(N,IX)=0
0232             P(N,IX)=0D0
0233             V(N,IX)=0D0
0234   250     CONTINUE
0235 C...Set special junction codes:
0236           K(N,1)=42
0237           K(N,2)=88
0238 C...Set parent to side.
0239           K(N,3)=MINT(83)+JS
0240           K(N,4)=ITJUNC(JS)*MSTU(5)
0241           K(N,5)=0
0242 C...Connect valence quarks to junction.
0243           MOUT(JS)=0
0244           MANTI=ITJUNC(JS)-1
0245 C...Set (anti)colour mother = junction.
0246           DO 260 JV=1,3
0247             K(IV(JS,JV),4+MANTI)=MOD(K(IV(JS,JV),4+MANTI),MSTU(5))
0248      &           +MSTU(5)*N
0249 C...Keep track of partons adjacent to junction:
0250             JST(JS,JV)=IV(JS,JV)
0251   260     CONTINUE
0252         ELSE
0253 C...Mesons: set up initial q-qbar topology
0254           ITJUNC(JS)=0
0255           IF (K(IV(JS,1),2).GT.0) THEN
0256             IQ=IV(JS,1)
0257             IQBAR=IV(JS,2)
0258           ELSE
0259             IQ=IV(JS,2)
0260             IQBAR=IV(JS,1)
0261           ENDIF
0262           IV(JS,3)=0
0263           JST(JS,1)=IQ
0264           JST(JS,2)=IQBAR
0265           JST(JS,3)=0
0266           K(IQ,4)=MOD(K(IQ,4),MSTU(5))+MSTU(5)*IQBAR
0267           K(IQBAR,5)=MOD(K(IQBAR,5),MSTU(5))+MSTU(5)*IQ
0268 C...Special for mesons. Insert gluon if BR empty.
0269           IF (NBRTOT(JS).EQ.0) THEN
0270             N=N+1
0271             DO 270 IX=1,5
0272               K(N,IX)=0
0273               P(N,IX)=0D0
0274               V(N,IX)=0D0
0275   270       CONTINUE
0276             K(N,1)=3
0277             K(N,2)=21
0278             K(N,3)=MINT(83)+JS
0279             K(N,4)=0
0280             K(N,5)=0
0281             NBRTOT(JS)=1
0282             NG(JS)=NG(JS)+1
0283 C...Add gluon to IMI
0284             NMI(JS)=NMI(JS)+1
0285             IMI(JS,NMI(JS),1)=N
0286             IMI(JS,NMI(JS),2)=0
0287           ENDIF
0288           MOUT(JS)=0
0289         ENDIF
0290
0291 C...Count up number of valence quarks outside BR.
0292         DO 280 JV=1,3
0293           IF (JST(JS,JV).LE.MINT(53).AND.JST(JS,JV).GT.0)
0294      &         MOUT(JS)=MOUT(JS)+1
0295   280   CONTINUE
0296
0297   290 CONTINUE
0298
0299 C...Now both sides have been prepared in an initial vvjv (baryonic) or
0300 C...v(g)vbar (mesonic) configuration.
0301
0302 C...Create colour line tags starting from initiators.
0303       NCT=0
0304       DO 320 IM=1,MINT(31)
0305 C...Consider each side in turn.
0306         DO 310 JS=1,2
0307           I1=IMI(JS,IM,1)
0308           I2=IMI(3-JS,IM,1)
0309           DO 300 JCS=4,5
0310             IF (K(I1,2).NE.21.AND.(9-2*JCS).NE.ISIGN(1,K(I1,2)))
0311      &           GOTO 300
0312             IF (K(I1,JCS)/MSTU(5)**2.NE.0) GOTO 300
0313
0314             KCS=JCS
0315             CALL PYCTTR(I1,KCS,I2)
0316             IF(MINT(51).NE.0) RETURN
0317
0318   300     CONTINUE
0319   310   CONTINUE
0320   320 CONTINUE
0321
0322       DO 340 JS=1,2
0323 C...Create colour tags for beam remnant partons.
0324         DO 330 IM=MINT(31)+1,NMI(JS)
0325           IP=IMI(JS,IM,1)
0326           IF (K(IP,2).NE.21) THEN
0327             JC=(3-ISIGN(1,K(IP,2)))/2
0328             IF (MCT(IP,JC).EQ.0) THEN
0329               NCT=NCT+1
0330               MCT(IP,JC)=NCT
0331             ENDIF
0332           ELSE
0333 C...Gluons
0334             ICD=K(IP,4)
0335             IAD=K(IP,5)
0336             IF (ICD.NE.0) THEN
0337 C...Fictituous gluons just inherit from their quark daughters.
0338               ICC=MCT(ICD,1)
0339               IAC=MCT(IAD,2)
0340             ELSE
0341 C...Real beam remnant gluons get their own colours
0342               ICC=NCT+1
0343               IAC=NCT+2
0344               NCT=NCT+2
0345             ENDIF
0346             MCT(IP,1)=ICC
0347             MCT(IP,2)=IAC
0348           ENDIF
0349   330   CONTINUE
0350   340 CONTINUE
0351
0352 C...Create colour tags for colour lines which are detached from the
0353 C...initial state.
0354
0355       DO 360 MQGST=1,2
0356         DO 350 I=MINT(84)+1,N
0357
0358 C...Look for coloured string endpoint, or (later) leftover gluon.
0359           IF (K(I,1).NE.3) GOTO 350
0360           KC=PYCOMP(K(I,2))
0361           IF(KC.EQ.0) GOTO 350
0362           KQ=KCHG(KC,2)
0363           IF(KQ.EQ.0.OR.(MQGST.EQ.1.AND.KQ.EQ.2)) GOTO 350
0364
0365 C...Pick up loose string end with no previous tag.
0366           KCS=4
0367           IF(KQ*ISIGN(1,K(I,2)).LT.0) KCS=5
0368           IF(MCT(I,KCS-3).NE.0) GOTO 350
0369
0370           CALL PYCTTR(I,KCS,I)
0371           IF(MINT(51).NE.0) RETURN
0372
0373   350   CONTINUE
0374   360 CONTINUE
0375
0376 C...Store original colour tags
0377       DO 370 I=MINT(84)+1,N
0378         MCO(I,1)=MCT(I,1)
0379         MCO(I,2)=MCT(I,2)
0380   370 CONTINUE
0381
0382 C...Iteratively add gluons to already existing string pieces, enforcing
0383 C...various possible orderings, and rejecting insertions that would give
0384 C...rise to singlet gluons.
0385 C...<kappa tau> normalization.
0386       RM0=1.5D0
0387       MRETRY=0
0388       PARP80=PARP(80)
0389
0390 C...Set up simplified kinematics.
0391 C...Boost hard interaction systems.
0392       IBOOST=IBOOST+1
0393       DO 380 IM=1,MINT(31)
0394         BETA=(XMI(1,IM)-XMI(2,IM))/(XMI(1,IM)+XMI(2,IM))
0395         CALL PYROBO(IMISEP(IM-1)+1,IMISEP(IM),0D0,0D0,0D0,0D0,BETA)
0396   380 CONTINUE
0397 C...Assign preliminary beam remnant momenta.
0398       DO 390 I=MINT(53)+1,N
0399         JS=K(I,3)
0400         P(I,1)=0D0
0401         P(I,2)=0D0
0402         IF (K(I,2).NE.88) THEN
0403           P(I,4)=0.5D0*VINT(142+JS)*VINT(1)/MAX(1,NMI(JS)-MINT(31))
0404           P(I,3)=P(I,4)
0405           IF (JS.EQ.2) P(I,3)=-P(I,3)
0406         ELSE
0407 C...Junctions are wildcards for the present.
0408           P(I,4)=0D0
0409           P(I,3)=0D0
0410         ENDIF
0411   390 CONTINUE
0412
0413 C...Reset colour processing information.
0414   400 DO 410 I=MINT(84)+1,N
0415         K(I,4)=MOD(K(I,4),MSTU(5)**2)
0416         K(I,5)=MOD(K(I,5),MSTU(5)**2)
0417   410 CONTINUE
0418
0419       NCC=0
0420       DO 430 JS=1,2
0421 C...If meson,  without gluon in BR, collapse q-qbar colour tags:
0422         IF (ITJUNC(JS).EQ.0) THEN
0423           JC1=MCT(JST(JS,1),1)
0424           JC2=MCT(JST(JS,2),2)
0425           NCC=NCC+1
0426           JCCO(NCC,1)=MAX(JC1,JC2)
0427           JCCO(NCC,2)=MIN(JC1,JC2)
0428 C...Collapse colour tags in event record
0429           DO 420 I=MINT(84)+1,N
0430             IF (MCT(I,1).EQ.JCCO(NCC,1)) MCT(I,1)=JCCO(NCC,2)
0431             IF (MCT(I,2).EQ.JCCO(NCC,1)) MCT(I,2)=JCCO(NCC,2)
0432   420     CONTINUE
0433         ENDIF
0434   430 CONTINUE
0435
0436   440 JS=1
0437       IF (PYR(0).GT.0.5D0.OR.NG(1).EQ.0) JS=2
0438       IF (NG(JS).GT.0) THEN
0439         NOPT=0
0440         RLOPT=1D9
0441 C...Start at random gluon (optimizes speed for random attachments)
0442         NMGL=0
0443         IMGL=PYR(0)*NMI(JS)+1
0444   450   IMGL=MOD(IMGL,NMI(JS))+1
0445         NMGL=NMGL+1
0446 C...Only loop through NMI once (with upper limit to save time)
0447         IF (NMGL.LE.NMI(JS).AND.NOPT.LE.3) THEN
0448           IGL  = IMI(JS,IMGL,1)
0449 C...If not gluon or if already connected, try next.
0450           IF (K(IGL,2).NE.21.OR.K(IGL,4)/MSTU(5).NE.0
0451      &         .OR.K(IGL,5)/MSTU(5).NE.0) GOTO 450
0452 C...Now loop through all possible insertions of this gluon.
0453           NMP1=0
0454           IMP1=PYR(0)*NMI(JS)+1
0455   460     IMP1=MOD(IMP1,NMI(JS))+1
0456           NMP1=NMP1+1
0457           IF (IMP1.EQ.IMGL) GOTO 460
0458 C...Only loop through NMI once (with upper limit to save time).
0459           IF (NMP1.LE.NMI(JS).AND.NOPT.LE.3) THEN
0460             IP1  = IMI(JS,IMP1,1)
0461 C...Try both colour mother and colour anti-mother.
0462 C...Randomly select which one to try first.
0463             NANTI=0
0464             MANTI=PYR(0)*2
0465   470       MANTI=MOD(MANTI+1,2)
0466             NANTI=NANTI+1
0467             IF (NANTI.LE.2) THEN
0468               IP2 =MOD(K(IP1,4+MANTI)/MSTU(5),MSTU(5))
0469 C...Reject if no appropriate mother (or if mother is fictitious
0470 C...parent gluon.)
0471               IF (IP2.LE.0) GOTO 470
0472               IF (K(IP2,2).EQ.21.AND.IP2.GT.MINT(53)) GOTO 470
0473 C...Also reject if this link has already been tried.
0474               IF (K(IP1,4+MANTI)/MSTU(5)**2.EQ.2) GOTO 470
0475               IF (K(IP2,5-MANTI)/MSTU(5)**2.EQ.2) GOTO 470
0476 C...Set flag to indicate that this link has now been tried for this
0477 C...gluon. IP2 may be junction, which has several mothers.
0478               K(IP1,4+MANTI)=K(IP1,4+MANTI)+2*MSTU(5)**2
0479               IF (K(IP2,2).NE.88) THEN
0480                 K(IP2,5-MANTI)=K(IP2,5-MANTI)+2*MSTU(5)**2
0481               ENDIF
0482
0483 C...JCG1: Original colour tag of gluon on IP1 side
0484 C...JCG2: Original colour tag of gluon on IP2 side
0485 C...JCP1: Original colour tag of IP1 on gluon side
0486 C...JCP2: Original colour tag of IP2 on gluon side.
0487               JCG1=MCO(IGL,2-MANTI)
0488               JCG2=MCO(IGL,1+MANTI)
0489               JCP1=MCO(IP1,1+MANTI)
0490               JCP2=MCO(IP2,2-MANTI)
0491
0492               CALL PYMIHG(JCP1,JCG1,JCP2,JCG2)
0493 C...Reject gluon attachments that give rise to singlet gluons.
0494               IF (MACCPT.EQ.0) GOTO 470
0495
0496 C...Update colours
0497               JCG1=MCT(IGL,2-MANTI)
0498               JCG2=MCT(IGL,1+MANTI)
0499               JCP1=MCT(IP1,1+MANTI)
0500               JCP2=MCT(IP2,2-MANTI)
0501
0502 C...Select whether to accept this insertion
0503               IF (MSTP(89).EQ.0) THEN
0504 C...Random insertions: no measure.
0505                 RL=1D0
0506 C...For random ordering, we want to suppress beam remnant breakups
0507 C...already at this point.
0508                 IF (IP1.GT.MINT(53).AND.IP2.GT.MINT(53)
0509      &               .AND.MOUT(JS).NE.0.AND.PYR(0).GT.PARP80) THEN
0510                   NMP1=0
0511                   NMGL=0
0512                   GOTO 470
0513                 ENDIF
0514               ELSEIF (MSTP(89).EQ.1) THEN
0515 C...Rapidity ordering:
0516 C...YGL = Rapidity of gluon.
0517                 YGL=YMI(IMGL)
0518 C...If fictitious gluon
0519                 IF (YGL.EQ.100D0) THEN
0520                   YGL=(3-2*JS)*100D0
0521                   IDA1=MOD(K(IGL,4),MSTU(5))
0522                   IDA2=MOD(K(IGL,5),MSTU(5))
0523                   DO 480 IMT=1,NMI(JS)
0524 C...Select (arbitrarily) the most central daughter.
0525                     IF (IMI(JS,IMT,1).EQ.IDA1.OR.IMI(JS,IMT,1).EQ.IDA2)
0526      &                   THEN
0527                       IF (ABS(YGL).GT.ABS(YMI(IMT))) YGL=YMI(IMT)
0528                     ENDIF
0529   480             CONTINUE
0530                 ENDIF
0531 C...YP1 = Rapidity IP1
0532                 YP1=YMI(IMP1)
0533 C...If fictitious gluon
0534                 IF (YP1.EQ.100D0) THEN
0535                   YP1=(3-2*JS)*YP1
0536                   IDA1=MOD(K(IP1,4),MSTU(5))
0537                   IDA2=MOD(K(IP1,5),MSTU(5))
0538                   DO 490 IMT=1,NMI(JS)
0539 C...Select (arbitrarily) the most central daughter.
0540                     IF (IMI(JS,IMT,1).EQ.IDA1.OR.IMI(JS,IMT,1).EQ.IDA2)
0541      &                   THEN
0542                       IF (ABS(YP1).GT.ABS(YMI(IMT))) YP1=YMI(IMT)
0543                     ENDIF
0544   490             CONTINUE
0545                 ENDIF
0546 C...YP2 = Rapidity of mother system
0547                 IF (K(IP2,2).NE.88) THEN
0548                   DO 500 IMT=1,NMI(JS)
0549                     IF (IMI(JS,IMT,1).EQ.IP2) YP2=YMI(IMT)
0550   500             CONTINUE
0551 C...If fictitious gluon
0552                   IF (YP2.EQ.100D0) THEN
0553                     YP2=(3-2*JS)*YP2
0554                     IDA1=MOD(K(IP2,4),MSTU(5))
0555                     IDA2=MOD(K(IP2,5),MSTU(5))
0556                     DO 510 IMT=1,NMI(JS)
0557 C...Select (arbitrarily) the most central daughter.
0558                       IF (IMI(JS,IMT,1).EQ.IDA1.OR.IMI(JS,IMT,1).EQ.IDA2
0559      &                     ) THEN
0560                         IF (ABS(YP2).GT.ABS(YMI(IMT))) YP2=YMI(IMT)
0561                       ENDIF
0562   510               CONTINUE
0563                   ENDIF
0564 C...Assign (arbitrarily) 100D0 to junction also
0565                 ELSE
0566                   YP2=(3-2*JS)*100D0
0567                 ENDIF
0568                 RL=ABS(YGL-YP1)+ABS(YGL-YP2)
0569               ELSEIF (MSTP(89).EQ.2) THEN
0570 C...Lambda ordering:
0571 C...Compute lambda measure for this insertion.
0572                 RL=1D0
0573                 DO 520 IST=1,6
0574                   ISTR(IST)=0
0575   520           CONTINUE
0576 C...If IP2 is junction, not caught below.
0577                 IF (JCP2.EQ.0) THEN
0578                   ITJU=MOD(K(IP2,4)/MSTU(5),MSTU(5))
0579 C...Anti-junction is colour endpoint et vv., always on JCG2.
0580                   ISTR(5-ITJU)=IP2
0581                 ENDIF
0582                 DO 530 I=MINT(84)+1,N
0583                   IF (K(I,1).LT.10) THEN
0584 C...The new string pieces
0585                     IF (MCT(I,1).EQ.JCG1) ISTR(1)=I
0586                     IF (MCT(I,2).EQ.JCG1) ISTR(2)=I
0587                     IF (MCT(I,1).EQ.JCG2) ISTR(3)=I
0588                     IF (MCT(I,2).EQ.JCG2) ISTR(4)=I
0589                   ENDIF
0590   530           CONTINUE
0591 C...Also identify junctions as string endpoints.
0592                 DO 540 I=MINT(84)+1,N
0593                   ICMO=MOD(K(I,4)/MSTU(5),MSTU(5))
0594                   IAMO=MOD(K(I,5)/MSTU(5),MSTU(5))
0595 C...Find partons adjacent to junctions.
0596                   IF (K(ICMO,1).EQ.42.AND.MCT(I,1).EQ.JCG1.AND.ISTR(2)
0597      &                 .EQ.0) ISTR(2) = ICMO
0598                   IF (K(IAMO,1).EQ.42.AND.MCT(I,2).EQ.JCG1.AND.ISTR(1)
0599      &                 .EQ.0) ISTR(1) = IAMO
0600                   IF (K(ICMO,1).EQ.42.AND.MCT(I,1).EQ.JCG2.AND.ISTR(4)
0601      &                 .EQ.0) ISTR(4) = ICMO
0602                   IF (K(IAMO,1).EQ.42.AND.MCT(I,2).EQ.JCG2.AND.ISTR(3)
0603      &                 .EQ.0) ISTR(3) = IAMO
0604   540           CONTINUE
0605 C...The old string piece
0606                 ISTR(5)=ISTR(1+2*MANTI)
0607                 ISTR(6)=ISTR(4-2*MANTI)
0608                 RL=MAX(1D0,FOUR(ISTR(1),ISTR(2)))*MAX(1D0,FOUR(ISTR(3)
0609      &               ,ISTR(4)))/MAX(1D0,FOUR(ISTR(5),ISTR(6)))
0610                 RL=LOG(RL)
0611               ENDIF
0612 C...Allow some breadth to speed things up.
0613               IF (ABS(1D0-RL/RLOPT).LT.0.05D0) THEN
0614                 NOPT=NOPT+1
0615               ELSEIF (RL.GT.RLOPT) THEN
0616                 GOTO 470
0617               ELSE
0618                 NOPT=1
0619                 RLOPT=RL
0620               ENDIF
0621 C...INSR(NOPT,1)=Gluon colour mother
0622 C...INSR(NOPT,2)=Gluon
0623 C...INSR(NOPT,3)=Gluon anticolour mother
0624               IF (NOPT.GT.1000) GOTO 470
0625               INSR(NOPT,1+2*MANTI)=IP2
0626               INSR(NOPT,2)=IGL
0627               INSR(NOPT,3-2*MANTI)=IP1
0628               IF (MSTP(89).GT.0.OR.NOPT.EQ.0) GOTO 470
0629             ENDIF
0630             IF (MSTP(89).GT.0.OR.NOPT.EQ.0) GOTO 460
0631           ENDIF
0632 C...Reset link test information.
0633           DO 550 I=MINT(84)+1,N
0634             K(I,4)=MOD(K(I,4),MSTU(5)**2)
0635             K(I,5)=MOD(K(I,5),MSTU(5)**2)
0636   550     CONTINUE
0637           IF (MSTP(89).GT.0.OR.NOPT.EQ.0) GOTO 450
0638         ENDIF
0639 C...Now we have a list of best gluon insertions, none of which cause
0640 C...singlets to arise. If list is empty, try again a few times. Note:
0641 C...this should never happen if we have a meson with a gluon inserted
0642 C...in the beam remnant, since that breaks up the colour line.
0643         IF (NOPT.EQ.0) THEN
0644 C...Abandon BR-g-BR suppression for retries. This is not serious, it
0645 C...just means we happened to start with trying a bad sequence.
0646           PARP80=1D0
0647           IF (MRETRY.LE.10.AND.(ITJUNC(1).NE.0.OR.JST(1,3).EQ.0).AND
0648      &         .(ITJUNC(2).NE.0.OR.JST(2,3).EQ.0)) THEN
0649             MRETRY=MRETRY+1
0650             DO 590 JS=1,2
0651               IF (ITJUNC(JS).NE.0) THEN
0652                 JST(JS,1)=IV(JS,1)
0653                 JST(JS,2)=IV(JS,2)
0654                 JST(JS,3)=IV(JS,3)
0655 C...Reset valence quark parent pointers
0656                 DO 560 I=MINT(53)+1,N
0657                   IF (K(I,2).EQ.88.AND.K(I,3).EQ.JS) IJU=I
0658   560           CONTINUE
0659                 MANTI=ITJUNC(JS)-1
0660 C...Set (anti)colour mother = junction.
0661                 DO 570 JV=1,3
0662                   K(IV(JS,JV),4+MANTI)=MOD(K(IV(JS,JV),4+MANTI),MSTU(5))
0663      &                 +MSTU(5)*IJU
0664   570           CONTINUE
0665               ELSE
0666 C...Same for mesons. JST unchanged, so needn't be restored.
0667                 IQ=JST(JS,1)
0668                 IQBAR=JST(JS,2)
0669                 K(IQ,4)=MOD(K(IQ,4),MSTU(5))+MSTU(5)*IQBAR
0670                 K(IQBAR,5)=MOD(K(IQBAR,5),MSTU(5))+MSTU(5)*IQ
0671               ENDIF
0672 C...Also reset gluon parent pointers.
0673               NG(JS)=0
0674               DO 580 IM=1,NMI(JS)
0675                 I=IMI(JS,IM,1)
0676                 IF (K(I,2).EQ.21) THEN
0677                   K(I,4)=MOD(K(I,4),MSTU(5))
0678                   K(I,5)=MOD(K(I,5),MSTU(5))
0679                   NG(JS)=NG(JS)+1
0680                 ENDIF
0681   580         CONTINUE
0682   590       CONTINUE
0683 C...Reset colour tags
0684             DO 600 I=MINT(84)+1,N
0685               MCT(I,1)=MCO(I,1)
0686               MCT(I,2)=MCO(I,2)
0687   600       CONTINUE
0688             GOTO 400
0689           ELSE
0690             IF(NERRPR.LT.5) THEN
0691               NERRPR=NERRPR+1
0692               CALL PYLIST(4)
0693               CALL PYERRM(19,'(PYMIHK:) No physical colour flow found!')
0694               WRITE(MSTU(11),*) 'NG:', NG,'   MOUT:', MOUT(JS)
0695             ENDIF
0696 C...Kill event and start another.
0697             MINT(51)=1
0698             RETURN
0699           ENDIF
0700         ELSE
0701 C...Select between insertions, suppressing insertions wholly in the BR.
0702           IIN=PYR(0)*NOPT+1
0703   610     IIN=MOD(IIN,NOPT)+1
0704           IF (INSR(IIN,1).GT.MINT(53).AND.INSR(IIN,3).GT.MINT(53)
0705      &         .AND.MOUT(JS).NE.0.AND.PYR(0).GT.PARP80) GOTO 610
0706         ENDIF
0707
0708 C...Now we know which gluon to insert where. Colour tags in JCCO and
0709 C...colour connection information should be updated, NG(JS) should be
0710 C...counted down, and a new loop performed if there are still gluons
0711 C...left on any side.
0712         ICM=INSR(IIN,1)
0713         IACM=INSR(IIN,3)
0714         IGL=INSR(IIN,2)
0715 C...JCG : Original gluon colour tag
0716 C...JCAG: Original gluon anticolour tag.
0717 C...JCM : Original anticolour tag of gluon colour mother
0718 C...JACM: Original colour tag of gluon anticolour mother
0719         JCG=MCO(IGL,1)
0720         JCM=MCO(ICM,2)
0721         JACG=MCO(IGL,2)
0722         JACM=MCO(IACM,1)
0723
0724         CALL PYMIHG(JACM,JACG,JCM,JCG)
0725         IF (MACCPT.EQ.0) THEN
0726           IF(NERRPR.LT.5) THEN
0727             NERRPR=NERRPR+1
0728             CALL PYLIST(4)
0729             CALL PYERRM(11,'(PYMIHK:) Unphysical colour flow!')
0730             WRITE(MSTU(11),*) 'attaching', IGL,' between', ICM, IACM
0731           ENDIF
0732 C...Kill event and start another.
0733           MINT(51)=1
0734           RETURN
0735         ELSE
0736 C...If everything went fine, store new JCCN in JCCO.
0737           NCC=NCC+1
0738           DO 620 ICC=1,NCC
0739             JCCO(ICC,1)=JCCN(ICC,1)
0740             JCCO(ICC,2)=JCCN(ICC,2)
0741   620     CONTINUE
0742         ENDIF
0743
0744 C...One gluon attached is counted as equivalent to one end outside.
0745         MOUT(JS)=1
0746 C...Set IGL colour mother = ICM.
0747         K(IGL,4)=MOD(K(IGL,4),MSTU(5))+MSTU(5)*ICM
0748 C...Set ICM anticolour mother = IGL colour.
0749         IF (K(ICM,2).NE.88) THEN
0750           K(ICM,5)=MOD(K(ICM,5),MSTU(5))+MSTU(5)*IGL
0751         ELSE
0752 C...If ICM is junction, just update JST array for now.
0753           DO 630 MSJ=1,3
0754             IF (JST(JS,MSJ).EQ.IACM) JST(JS,MSJ)=IGL
0755   630     CONTINUE
0756         ENDIF
0757 C...Set IGL anticolour mother = IACM.
0758         K(IGL,5)=MOD(K(IGL,5),MSTU(5))+MSTU(5)*IACM
0759 C...Set IACM anticolour mother = IGL anticolour.
0760         IF (K(IACM,2).NE.88) THEN
0761           K(IACM,4)=MOD(K(IACM,4),MSTU(5))+MSTU(5)*IGL
0762         ELSE
0763 C...If IACM is junction, just update JST array for now.
0764           DO 640 MSJ=1,3
0765             IF (JST(JS,MSJ).EQ.ICM) JST(JS,MSJ)=IGL
0766   640     CONTINUE
0767         ENDIF
0768 C...Count down # unconnected gluons.
0769         NG(JS)=NG(JS)-1
0770       ENDIF
0771       IF (NG(1).GT.0.OR.NG(2).GT.0) GOTO 440
0772
0773       DO 840 JS=1,2
0774 C...Collapse fictitious gluons.
0775         DO 670 IGL=MINT(53)+1,N
0776           IF (K(IGL,2).EQ.21.AND.K(IGL,3).EQ.MINT(83)+JS.AND.
0777      &         K(IGL,1).EQ.14) THEN
0778             ICM=K(IGL,4)/MSTU(5)
0779             IAM=K(IGL,5)/MSTU(5)
0780             ICD=MOD(K(IGL,4),MSTU(5))
0781             IAD=MOD(K(IGL,5),MSTU(5))
0782 C...Set gluon daughters pointing to gluon mothers
0783             K(IAD,5)=MOD(K(IAD,5),MSTU(5))+MSTU(5)*IAM
0784             K(ICD,4)=MOD(K(ICD,4),MSTU(5))+MSTU(5)*ICM
0785 C...Set gluon mothers pointing to gluon daughters.
0786             IF (K(ICM,2).NE.88) THEN
0787               K(ICM,5)=MOD(K(ICM,5),MSTU(5))+MSTU(5)*ICD
0788             ELSE
0789 C...Special case: mother=junction. Just update JST array for now.
0790               DO 650 MSJ=1,3
0791                 IF (JST(JS,MSJ).EQ.IGL) JST(JS,MSJ)=ICD
0792   650         CONTINUE
0793             ENDIF
0794             IF (K(IAM,2).NE.88) THEN
0795               K(IAM,4)=MOD(K(IAM,4),MSTU(5))+MSTU(5)*IAD
0796             ELSE
0797               DO 660 MSJ=1,3
0798                 IF (JST(JS,MSJ).EQ.IGL) JST(JS,MSJ)=IAD
0799   660         CONTINUE
0800             ENDIF
0801           ENDIF
0802   670   CONTINUE
0803
0804 C...Erase collapsed gluons from NMI and IMI (but keep them in ER)
0805         IM=NMI(JS)+1
0806   680   IM=IM-1
0807         IF (IM.GT.MINT(31).AND.K(IMI(JS,IM,1),2).NE.21) GOTO 680
0808         IF (IM.GT.MINT(31)) THEN
0809           NMI(JS)=NMI(JS)-1
0810           DO 690 IMR=IM,NMI(JS)
0811             IMI(JS,IMR,1)=IMI(JS,IMR+1,1)
0812             IMI(JS,IMR,2)=IMI(JS,IMR+1,2)
0813   690     CONTINUE
0814           GOTO 680
0815         ENDIF
0816
0817 C...Finally, connect junction.
0818         IF (ITJUNC(JS).NE.0) THEN
0819           DO 700 I=MINT(53)+1,N
0820             IF (K(I,2).EQ.88.AND.K(I,3).EQ.MINT(83)+JS) IJU=I
0821   700     CONTINUE
0822 C...NBRJQ counts # of jq, NBRVQ # of jv, inside BR.
0823           NBRJQ =0
0824           NBRVQ =0
0825           DO 720 MSJ=1,3
0826             IDQ(MSJ)=0
0827 C...Find jq with no glue inbetween inside beam remnant.
0828             IF (JST(JS,MSJ).GT.MINT(53).AND.IABS(K(JST(JS,MSJ),2)).LE.5)
0829      &           THEN
0830               NBRJQ=NBRJQ+1
0831 C...Set IDQ = -I if q non-valence and = +I if q valence.
0832               IDQ(NBRJQ)=-JST(JS,MSJ)
0833               DO 710 JV=1,3
0834                 IF (IV(JS,JV).EQ.JST(JS,MSJ)) THEN
0835                   IDQ(NBRJQ)=JST(JS,MSJ)
0836                   NBRVQ=NBRVQ+1
0837                 ENDIF
0838   710         CONTINUE
0839             ENDIF
0840             I12=MOD(MSJ+1,2)
0841             I45=5
0842             IF (MSJ.EQ.3) I45=4
0843             K(IJU,I45)=K(IJU,I45)+(MSTU(5)**I12)*JST(JS,MSJ)
0844   720     CONTINUE
0845
0846 C...Check if diquark can be formed.
0847           IF ((MSTP(88).GE.0.AND.NBRVQ.GE.2).OR.(NBRJQ.GE.2.AND.MSTP(88)
0848      &         .GE.1)) THEN
0849 C...If there is less than 2 valence quarks connected to junction
0850 C...and MSTP(88)>1, use random non-valence quarks to fill up.
0851             IF (NBRVQ.LE.1) THEN
0852               NDIQ=NBRVQ
0853   730         JFLIP=NBRJQ*PYR(0)+1
0854               IF (IDQ(JFLIP).LT.0) THEN
0855                 IDQ(JFLIP)=-IDQ(JFLIP)
0856                 NDIQ=NDIQ+1
0857               ENDIF
0858               IF (NDIQ.LE.1) GOTO 730
0859             ENDIF
0860 C...Place selected quarks first in IDQ, ordered in flavour.
0861             DO 740 JDQ=1,3
0862               IF (IDQ(JDQ).LE.0) THEN
0863                 ITEMP1  = IDQ(JDQ)
0864                 IDQ(JDQ)= IDQ(3)
0865                 IDQ(3)  = -ITEMP1
0866                 IF (IABS(K(IDQ(1),2)).LT.IABS(K(IDQ(2),2))) THEN
0867                   ITEMP1  = IDQ(1)
0868                   IDQ(1)  = IDQ(2)
0869                   IDQ(2)  = ITEMP1
0870                 ENDIF
0871               ENDIF
0872   740       CONTINUE
0873 C...Choose diquark spin.
0874             IF (NBRVQ.EQ.2) THEN
0875 C...If the selected quarks are both valence, we may use SU(6) rules
0876 C...to figure out which spin the diquark has, by a subdivision of the
0877 C...original beam hadron into the selected diquark system plus a kicked
0878 C...out quark, IKO.
0879               JKO=6
0880               DO 760 JDQ=1,2
0881                 DO 750 JV=1,3
0882                   IF (IDQ(JDQ).EQ.IV(JS,JV)) JKO=JKO-JV
0883   750           CONTINUE
0884   760         CONTINUE
0885               IKO=IV(JS,JKO)
0886               CALL PYSPLI(MINT(10+JS),K(IKO,2),KFDUM,KFDQ)
0887             ELSE
0888 C...If one or more of the selected quarks are not valence, we cannot use
0889 C...SU(6) subdivisions of the original beam hadron. Instead, with the
0890 C...flavours of the diquark already selected, we assume for now
0891 C...50:50 spin-1:spin-0 (where spin-0 possible).
0892               KFDQ=1000*K(IDQ(1),2)+100*K(IDQ(2),2)
0893               IS=3
0894               IF (K(IDQ(1),2).NE.K(IDQ(2),2).AND.
0895      &           (1D0+3D0*PARJ(4))*PYR(0).LT.1D0) IS=1
0896               KFDQ=KFDQ+ISIGN(IS,KFDQ)
0897             ENDIF
0898
0899 C...Collapse diquark-j-quark system to baryon, if allowed and possible.
0900 C...Note: third quark can per definition not also be valence,
0901 C...therefore we can only do this if we are allowed to use sea quarks.
0902   770       IF (IDQ(3).NE.0.AND.MSTP(88).GE.2) THEN
0903               NTRY=0
0904   780         NTRY=NTRY+1
0905               CALL PYKFDI(KFDQ,K(IABS(IDQ(3)),2),KFDUM,KFBAR)
0906               IF (KFBAR.EQ.0.AND.NTRY.LE.100) THEN
0907                 GOTO 780
0908               ELSEIF(NTRY.GT.100) THEN
0909 C...If no baryon can be found, give up and form diquark.
0910                 IDQ(3)=0
0911                 GOTO 770
0912               ELSE
0913 C...Replace junction by baryon.
0914                 K(IJU,1)=1
0915                 K(IJU,2)=KFBAR
0916                 K(IJU,3)=MINT(83)+JS
0917                 K(IJU,4)=0
0918                 K(IJU,5)=0
0919                 P(IJU,5)=PYMASS(KFBAR)
0920                 DO 790 MSJ=1,3
0921 C...Prepare removal of participating quarks from ER.
0922                   K(JST(JS,MSJ),1)=-1
0923   790           CONTINUE
0924               ENDIF
0925             ELSE
0926 C...If collapse to baryon not possible or not allowed, replace junction
0927 C...by diquark. This way, collapsed gluons that were pointing at the
0928 C...junction will now point (correctly) at diquark.
0929               MANTI=ITJUNC(JS)-1
0930               K(IJU,1)=3
0931               K(IJU,2)=KFDQ
0932               K(IJU,3)=MINT(83)+JS
0933               K(IJU,4)=0
0934               K(IJU,5)=0
0935               DO 800 MSJ=1,3
0936                 IP=JST(JS,MSJ)
0937                 IF (IP.NE.IDQ(1).AND.IP.NE.IDQ(2)) THEN
0938                   K(IJU,4+MANTI)=0
0939                   K(IJU,5-MANTI)=IP*MSTU(5)
0940                   K(IP,4+MANTI)=MOD(K(IP,4+MANTI),MSTU(5))+
0941      &                 MSTU(5)*IJU
0942                   MCT(IJU,2-MANTI)=MCT(IP,1+MANTI)
0943                 ELSE
0944 C...Prepare removal of participating quarks from ER.
0945                   K(IP,1)=-1
0946                 ENDIF
0947   800         CONTINUE
0948             ENDIF
0949
0950 C...Update so ER pointers to collapsed quarks
0951 C...now go to collapsed object.
0952             DO 820 I=MINT(84)+1,N
0953               IF ((K(I,3).EQ.MINT(83)+JS.OR.K(I,3).EQ.MINT(83)+2+JS).AND
0954      &             .K(I,1).GT.0) THEN
0955                 DO 810 ISID=4,5
0956                   IMO=K(I,ISID)/MSTU(5)
0957                   IDA=MOD(K(I,ISID),MSTU(5))
0958                   IF (IMO.GT.0) THEN
0959                     IF (K(IMO,1).EQ.-1) IMO=IJU
0960                   ENDIF
0961                   IF (IDA.GT.0) THEN
0962                     IF (K(IDA,1).EQ.-1) IDA=IJU
0963                   ENDIF
0964                   K(I,ISID)=IDA+MSTU(5)*IMO
0965   810           CONTINUE
0966               ENDIF
0967   820       CONTINUE
0968           ENDIF
0969         ENDIF
0970
0971 C...Finally, if beam remnant is empty, insert a gluon in beam remnant.
0972 C...(this only happens for baryons, where we want to force the gluon
0973 C...to sit next to the junction. Mesons handled above.)
0974         IF (NBRTOT(JS).EQ.0) THEN
0975           N=N+1
0976           DO 830 IX=1,5
0977             K(N,IX)=0
0978             P(N,IX)=0D0
0979             V(N,IX)=0D0
0980   830     CONTINUE
0981           IGL=N
0982           K(IGL,1)=3
0983           K(IGL,2)=21
0984           K(IGL,3)=MINT(83)+JS
0985           IF (ITJUNC(JS).NE.0) THEN
0986 C...Incoming baryons. Pick random leg in JST (NVSUM = 3 for baryons)
0987             JLEG=PYR(0)*NVSUM(JS)+1
0988             I1=JST(JS,JLEG)
0989             JST(JS,JLEG)=IGL
0990             JCT=MCT(I1,ITJUNC(JS))
0991             MCT(IGL,3-ITJUNC(JS))=JCT
0992             NCT=NCT+1
0993             MCT(IGL,ITJUNC(JS))=NCT
0994             MANTI=ITJUNC(JS)-1
0995           ELSE
0996 C...Meson. Should not happen.
0997             CALL PYERRM(19,'(PYMIHK:) Empty meson beam remnant')
0998             IF(NERRPR.LT.5) THEN
0999               WRITE(MSTU(11),*) 'This should not have been possible!'
1000               CALL PYLIST(4)
1001               NERRPR=NERRPR+1
1002             ENDIF
1003             MINT(51)=1
1004             RETURN
1005           ENDIF
1006           I2=MOD(K(I1,4+MANTI)/MSTU(5),MSTU(5))
1007           K(I1,4+MANTI)=MOD(K(I1,4+MANTI),MSTU(5))+MSTU(5)*IGL
1008           K(IGL,5-MANTI)=MOD(K(IGL,5-MANTI),MSTU(5))+MSTU(5)*I1
1009           K(IGL,4+MANTI)=MOD(K(IGL,4+MANTI),MSTU(5))+MSTU(5)*I2
1010           IF (K(I2,2).NE.88) THEN
1011             K(I2,5-MANTI)=MOD(K(I2,5-MANTI),MSTU(5))+MSTU(5)*IGL
1012           ELSE
1013             IF (MOD(K(I2,4),MSTU(5)).EQ.I1) THEN
1014               K(I2,4)=(K(I2,4)/MSTU(5))*MSTU(5)+IGL
1015             ELSEIF(MOD(K(I2,5)/MSTU(5),MSTU(5)).EQ.I1) THEN
1016               K(I2,5)=MOD(K(I2,5),MSTU(5))+MSTU(5)*IGL
1017             ELSE
1018               K(I2,5)=(K(I2,5)/MSTU(5))*MSTU(5)+IGL
1019             ENDIF
1020           ENDIF
1021         ENDIF
1022   840 CONTINUE
1023
1024 C...Remove collapsed quarks and junctions from ER and update IMI.
1025       CALL PYEDIT(11)
1026
1027 C...Also update beam remnant part of IMI.
1028       NMI(1)=MINT(31)
1029       NMI(2)=MINT(31)
1030       DO 850 I=MINT(53)+1,N
1031         IF (K(I,1).LE.0) GOTO 850
1032 C...Restore BR quark/diquark/baryon pointers in IMI.
1033         IF ((K(I,2).NE.21.OR.K(I,1).NE.14).AND.K(I,2).NE.88) THEN
1034           JS=K(I,3)-MINT(83)
1035           NMI(JS)=NMI(JS)+1
1036           IMI(JS,NMI(JS),1)=I
1037           IMI(JS,NMI(JS),2)=0
1038         ENDIF
1039   850 CONTINUE
1040
1041 C...Restore companion information from collapsed gluons.
1042       DO 870 I=MINT(53)+1,N
1043         IF (K(I,2).EQ.21.AND.K(I,1).EQ.14) THEN
1044           JS=K(I,3)-MINT(83)
1045           JCD=MOD(K(I,4),MSTU(5))
1046           JAD=MOD(K(I,5),MSTU(5))
1047           DO 860 IM=1,NMI(JS)
1048             IF (IMI(JS,IM,1).EQ.JCD) IMC=IM
1049             IF (IMI(JS,IM,1).EQ.JAD) IMA=IM
1050   860     CONTINUE
1051           IMI(JS,IMC,2)=IMI(JS,IMA,1)
1052           IMI(JS,IMA,2)=IMI(JS,IMC,1)
1053         ENDIF
1054   870 CONTINUE
1055
1056 C...Renumber colour lines (since some have disappeared)
1057       JCT=0
1058       JCD=0
1059   880 JCT=JCT+1
1060       MFOUND=0
1061       I=MINT(84)
1062   890 I=I+1
1063       IF (I.EQ.N+1) THEN
1064         IF (MFOUND.EQ.0) JCD=JCD+1
1065       ELSEIF (MCT(I,1).EQ.JCT.AND.K(I,1).GE.1) THEN
1066         MCT(I,1)=JCT-JCD
1067         MFOUND=1
1068       ELSEIF (MCT(I,2).EQ.JCT.AND.K(I,1).GE.1) THEN
1069         MCT(I,2)=JCT-JCD
1070         MFOUND=1
1071       ENDIF
1072       IF (I.LE.N) GOTO 890
1073       IF (JCT.LT.NCT) GOTO 880
1074       NCT=JCT-JCD
1075
1076 C...Reset hard interaction subsystems to their CM frames.
1077       IF (IBOOST.EQ.1) THEN
1078         DO 900 IM=1,MINT(31)
1079           BETA=-(XMI(1,IM)-XMI(2,IM))/(XMI(1,IM)+XMI(2,IM))
1080           CALL PYROBO(IMISEP(IM-1)+1,IMISEP(IM),0D0,0D0,0D0,0D0,BETA)
1081   900   CONTINUE
1082 C...Zero beam remnant longitudinal momenta and energies
1083         DO 910 I=MINT(53)+1,N
1084           P(I,3)=0D0
1085           P(I,4)=0D0
1086   910   CONTINUE
1087       ELSE
1088         CALL PYERRM(9
1089      &       ,'(PYMIHK:) Inconsistent kinematics. Too many boosts.')
1090 C...Kill event and start another.
1091         MINT(51)=1
1092         RETURN
1093       ENDIF
1094
1095  9999 RETURN
1096       END