pythiaeRHIC/pythia/pymirm.f

0001
0002 C*********************************************************************
0003
0004 C...PYMIRM
0005 C...Picks primordial kT and shares longitudinal momentum among
0006 C...beam remnants.
0007
0008       SUBROUTINE PYMIRM
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/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0019       COMMON/PYINT1/MINT(400),VINT(400)
0020 C...The common block of colour tags.
0021       COMMON/PYCTAG/NCT,MCT(4000,2)
0022 C...The common block of dangling ends
0023       COMMON/PYINTM/KFIVAL(2,3),NMI(2),IMI(2,800,2),NVC(2,-6:6),
0024      &     XASSOC(2,-6:6,240),XPSVC(-6:6,-1:240),PVCTOT(2,-1:1),
0025      &     XMI(2,240),PT2MI(240),IMISEP(0:240)
0026       SAVE /PYJETS/,/PYDAT1/,/PYPARS/,/PYINT1/,/PYINTM/,/PYCTAG/
0027 C...Local variables
0028       DIMENSION W(0:2,0:2),VB(3),NNXT(2),IVALQ(2),ICOMQ(2)
0029 C...W(I,J)|  J=0    |   1   |   2   |
0030 C...  I=0 | Wrem**2 |  W+   |  W-   |
0031 C...    1 | W1**2   |  W1+  |  W1-  |
0032 C...    2 | W2**2   |  W2+  |  W2-  |
0033 C...4-product
0034       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)
0035 C...Tentative parametrization of <kT> as a function of Q.
0036       SIGPT(Q)=MAX(PARJ(21),2.1D0*Q/(7D0+Q))
0037 C      SIGPT(Q)=MAX(0.36D0,4D0*SQRT(Q)/(10D0+SQRT(Q))
0038 C      SIGPT(Q)=MAX(PARJ(21),3D0*SQRT(Q)/(5D0+SQRT(Q))
0039       GETPT(Q,SIGMA)=MIN(SIGMA*SQRT(-LOG(PYR(0))),PARP(93))
0040 C...Lambda kinematic function.
0041       FLAM(A,B,C)=A**2+B**2+C**2-2D0*(A*B+B*C+C*A)
0042
0043 C...Beginning and end of beam remnant partons
0044       NOUT=MINT(53)
0045       ISUB=MINT(1)
0046
0047 C...Loopback point if kinematic choices gives impossible configuration.
0048       NTRY=0
0049   100 NTRY=NTRY+1
0050
0051 C...Assign kT values on each side separately.
0052       DO 180 JS=1,2
0053
0054 C...First zero all kT on this side. Skip if no kT to generate.
0055         DO 110 IM=1,NMI(JS)
0056           P(IMI(JS,IM,1),1)=0D0
0057           P(IMI(JS,IM,1),2)=0D0
0058   110   CONTINUE
0059         IF(MSTP(91).LE.0) GOTO 180
0060
0061 C...Now assign kT to each (non-collapsed) parton in IMI.
0062         DO 170 IM=1,NMI(JS)
0063           I=IMI(JS,IM,1)
0064 C...Select kT according to truncated gaussian or 1/kt6 tails.
0065 C...For first interaction, either use rms width = PARP(91) or fitted.
0066           IF (IM.EQ.1) THEN
0067             SIGMA=PARP(91)
0068             IF (MSTP(91).GE.11.AND.MSTP(91).LE.20) THEN
0069               Q=SQRT(PT2MI(IM))
0070               SIGMA=SIGPT(Q)
0071             ENDIF
0072           ELSE
0073 C...For subsequent interactions and BR partons use fragmentation width.
0074             SIGMA=PARJ(21)
0075           ENDIF
0076           PHI=PARU(2)*PYR(0)
0077           PT=0D0
0078           IF(NTRY.LE.100) THEN
0079  111        IF (MSTP(91).EQ.1.OR.MSTP(91).EQ.11) THEN
0080               PT=GETPT(Q,SIGMA)
0081               PTX=PT*COS(PHI)
0082               PTY=PT*SIN(PHI)
0083             ELSEIF (MSTP(91).EQ.2) THEN
0084               CALL PYERRM(11,'(PYMIRM:) Sorry, MSTP(91)=2 not '//
0085      &          'available, using MSTP(91)=1.')
0086               CALL PYGIVE('MSTP(91)=1')
0087               GOTO 111
0088             ELSEIF(MSTP(91).EQ.3.OR.MSTP(91).EQ.13) THEN
0089 C...Use distribution with kt**6 tails, rms width = PARP(91).
0090               EPS=SQRT(3D0/2D0)*SIGMA
0091 C...Generate PTX and PTY separately, each propto 1/KT**6
0092               DO 119 IXY=1,2
0093 C...Decide which interval to try
0094  112            P12=1D0/(1D0+27D0/40D0*SIGMA**6/EPS**6)
0095                 IF (PYR(0).LT.P12) THEN
0096 C...Use flat approx with accept/reject up to EPS.
0097                   PT=PYR(0)*EPS
0098                   WT=(3D0/2D0*SIGMA**2/(PT**2+3D0/2D0*SIGMA**2))**3
0099                   IF (PYR(0).GT.WT) GOTO 112
0100                 ELSE
0101 C...Above EPS, use 1/kt**6 approx with accept/reject.
0102                   PT=EPS/(PYR(0)**(1D0/5D0))
0103                   WT=PT**6/(PT**2+3D0/2D0*SIGMA**2)**3
0104                   IF (PYR(0).GT.WT) GOTO 112
0105                 ENDIF
0106                 MSIGN=1
0107                 IF (PYR(0).GT.0.5D0) MSIGN=-1
0108                 IF (IXY.EQ.1) PTX=MSIGN*PT
0109                 IF (IXY.EQ.2) PTY=MSIGN*PT
0110  119          CONTINUE
0111             ELSEIF (MSTP(91).EQ.4.OR.MSTP(91).EQ.14) THEN
0112               PTX=SIGMA*(SQRT(6D0)*PYR(0)-SQRT(3D0/2D0))
0113               PTY=SIGMA*(SQRT(6D0)*PYR(0)-SQRT(3D0/2D0))
0114             ENDIF
0115 C...Adjust final PT. Impose upper cutoff, or zero for soft evts.
0116             PT=SQRT(PTX**2+PTY**2)
0117             WT=1D0
0118             IF (PT.GT.PARP(93)) WT=SQRT(PARP(93)/PT)
0119             IF(ISUB.EQ.95.AND.IM.EQ.1) WT=0D0
0120             PTX=PTX*WT
0121             PTY=PTY*WT
0122             PT=SQRT(PTX**2+PTY**2)
0123           ENDIF
0124
0125           P(I,1)=P(I,1)+PTX
0126           P(I,2)=P(I,2)+PTY
0127
0128 C...Compensation kicks, with varying degree of local anticorrelations.
0129           MCORR=MSTP(90)
0130           IF (MCORR.EQ.0.OR.ISUB.EQ.95) THEN
0131             PTCX=-PTX/(NMI(JS)-1)
0132             PTCY=-PTY/(NMI(JS)-1)
0133             IF(ISUB.EQ.95) THEN
0134               PTCX=-PTX/(NMI(JS)-2)
0135               PTCY=-PTY/(NMI(JS)-2)
0136             ENDIF
0137             DO 120 IMC=1,NMI(JS)
0138               IF (IMC.EQ.IM) GOTO 120
0139               IF(ISUB.EQ.95.AND.IMC.EQ.1) GOTO 120
0140               P(IMI(JS,IMC,1),1)=P(IMI(JS,IMC,1),1)+PTCX
0141               P(IMI(JS,IMC,1),2)=P(IMI(JS,IMC,1),2)+PTCY
0142   120       CONTINUE
0143           ELSEIF (MCORR.GE.1) THEN
0144             DO 140 MSID=4,5
0145               NNXT(MSID-3)=0
0146 C...Count up # of neighbours on either side
0147               IMO=I
0148   130         IMO=K(IMO,MSID)/MSTU(5)
0149               IF (IMO.EQ.0) GOTO 140
0150               NNXT(MSID-3)=NNXT(MSID-3)+1
0151 C...Stop at quarks and junctions
0152               IF (MCORR.EQ.1.AND.K(IMO,2).EQ.21) GOTO 130
0153   140       CONTINUE
0154 C...How should compensation be shared when unequal numbers on the
0155 C...two sides? 50/50 regardless? N1:N2? Assume latter for now.
0156             NSUM=NNXT(1)+NNXT(2)
0157             T1=0
0158             DO 160 MSID=4,5
0159 C...Total momentum to be compensated on this side
0160               IF (NNXT(MSID-3).EQ.0) GOTO 160
0161               PTCX=-(NNXT(MSID-3)*PTX)/NSUM
0162               PTCY=-(NNXT(MSID-3)*PTY)/NSUM
0163 C...RS: compensation supression factor as we go out from parton I.
0164 C...Hardcoded behaviour RS=0.5, i.e. 1/2**n falloff,
0165 C...since (for now) MSTP(90) provides enough variability.
0166               RS=0.5D0
0167               FAC=(1D0-RS)/(RS*(1-RS**NNXT(MSID-3)))
0168               IMO=I
0169   150         IDA=IMO
0170               IMO=K(IMO,MSID)/MSTU(5)
0171               IF (IMO.EQ.0) GOTO 160
0172               FAC=FAC*RS
0173               IF (K(IMO,2).NE.88) THEN
0174                 P(IMO,1)=P(IMO,1)+FAC*PTCX
0175                 P(IMO,2)=P(IMO,2)+FAC*PTCY
0176                 IF (MCORR.EQ.1.AND.K(IMO,2).EQ.21) GOTO 150
0177 C...If we reach junction, divide out the kT that would have been
0178 C...assigned to the junction on each of its other legs.
0179               ELSE
0180                 L1=MOD(K(IMO,4),MSTU(5))
0181                 L2=K(IMO,5)/MSTU(5)
0182                 L3=MOD(K(IMO,5),MSTU(5))
0183                 P(L1,1)=P(L1,1)+0.5D0*FAC*PTCX
0184                 P(L1,2)=P(L1,2)+0.5D0*FAC*PTCY
0185                 P(L2,1)=P(L2,1)+0.5D0*FAC*PTCX
0186                 P(L2,2)=P(L2,2)+0.5D0*FAC*PTCY
0187                 P(L3,1)=P(L3,1)+0.5D0*FAC*PTCX
0188                 P(L3,2)=P(L3,2)+0.5D0*FAC*PTCY
0189                 P(IDA,1)=P(IDA,1)-0.5D0*FAC*PTCX
0190                 P(IDA,2)=P(IDA,2)-0.5D0*FAC*PTCY
0191               ENDIF
0192
0193   160       CONTINUE
0194           ENDIF
0195   170   CONTINUE
0196 C...End assignment of kT values to initiators and remnants.
0197   180 CONTINUE
0198
0199 C...Check kinematics constraints for non-BR partons.
0200       DO 190 IM=1,MINT(31)
0201         SHAT=XMI(1,IM)*XMI(2,IM)*VINT(2)
0202         PT1=SQRT(P(IMI(1,IM,1),1)**2+P(IMI(1,IM,1),2)**2)
0203         PT2=SQRT(P(IMI(2,IM,1),1)**2+P(IMI(2,IM,1),2)**2)
0204         PT1PT2=P(IMI(1,IM,1),1)*P(IMI(2,IM,1),1)
0205      &        +P(IMI(1,IM,1),2)*P(IMI(2,IM,1),2)
0206         IF (SHAT.LT.2D0*(PT1*PT2-PT1PT2).AND.NTRY.LE.100) THEN
0207           IF(NTRY.GE.100) THEN
0208 C...Kill this event and start another.
0209             CALL PYERRM(11,
0210      &           '(PYMIRM:) No consistent (x,kT) sets found')
0211             MINT(51)=1
0212             RETURN
0213           ENDIF
0214           GOTO 100
0215         ENDIF
0216   190 CONTINUE
0217
0218 C...Calculate W+ and W- available for combined remnant system.
0219       W(0,1)=VINT(1)
0220       W(0,2)=VINT(1)
0221       DO 200 IM=1,MINT(31)
0222         PT2 = (P(IMI(1,IM,1),1)+P(IMI(2,IM,1),1))**2
0223      &       +(P(IMI(1,IM,1),2)+P(IMI(2,IM,1),2))**2
0224         ST=XMI(1,IM)*XMI(2,IM)*VINT(2)+PT2
0225         W(0,1)=W(0,1)-SQRT(XMI(1,IM)/XMI(2,IM)*ST)
0226         W(0,2)=W(0,2)-SQRT(XMI(2,IM)/XMI(1,IM)*ST)
0227   200 CONTINUE
0228 C...Also store Wrem**2 = W+ * W-
0229       W(0,0)=W(0,1)*W(0,2)
0230
0231       IF (W(0,0).LT.0D0.AND.NTRY.LE.100) THEN
0232           IF(NTRY.GE.100) THEN
0233 C...Kill this event and start another.
0234             CALL PYERRM(11,
0235      &    '(PYMIRM:) Negative beam remnant mass squared unavoidable')
0236             MINT(51)=1
0237             RETURN
0238           ENDIF
0239           GOTO 100
0240       ENDIF
0241
0242 C...Assign unscaled x values to partons/hadrons in each of the
0243 C...beam remnants and calculate unscaled W+ and W- from them.
0244       NTRYX=0
0245   210 NTRYX=NTRYX+1
0246       DO 280 JS=1,2
0247         W(JS,1)=0D0
0248         W(JS,2)=0D0
0249         DO 270 IM=MINT(31)+1,NMI(JS)
0250           I=IMI(JS,IM,1)
0251           KF=K(I,2)
0252           KFA=IABS(KF)
0253           ICOMP=IMI(JS,IM,2)
0254
0255 C...Skip collapsed gluons and junctions. Reset.
0256           IF (KFA.EQ.21.AND.K(I,1).EQ.14) GOTO 270
0257           IF (KFA.EQ.88) GOTO 270
0258           X=0D0
0259           IVALQ(1)=0
0260           IVALQ(2)=0
0261           ICOMQ(1)=0
0262           ICOMQ(2)=0
0263
0264 C...If gluon then only beam remnant, so takes all.
0265           IF(KFA.EQ.21) THEN
0266             X=1D0
0267 C...If valence quark then use parametrized valence distribution.
0268           ELSEIF(KFA.LE.6.AND.ICOMP.EQ.0) THEN
0269             IVALQ(1)=KF
0270 C...If companion quark then derive from companion x.
0271           ELSEIF(KFA.LE.6) THEN
0272             ICOMQ(1)=ICOMP
0273 C...If valence diquark then use two parametrized valence distributions.
0274           ELSEIF(KFA.GT.1000.AND.MOD(KFA/10,10).EQ.0.AND.
0275      &    ICOMP.EQ.0) THEN
0276             IVALQ(1)=ISIGN(KFA/1000,KF)
0277             IVALQ(2)=ISIGN(MOD(KFA/100,10),KF)
0278 C...If valence+sea diquark then combine valence + companion choices.
0279           ELSEIF(KFA.GT.1000.AND.MOD(KFA/10,10).EQ.0.AND.
0280      &    ICOMP.LT.MSTU(5)) THEN
0281             IF(KFA/1000.EQ.IABS(K(ICOMP,2))) THEN
0282               IVALQ(1)=ISIGN(MOD(KFA/100,10),KF)
0283             ELSE
0284               IVALQ(1)=ISIGN(KFA/1000,KF)
0285             ENDIF
0286             ICOMQ(1)=ICOMP
0287 C...Extra code: workaround for diquark made out of two sea
0288 C...quarks, but where not (yet) ICOMP > MSTU(5).
0289             DO 220 IM1=1,MINT(31)
0290               IF(IMI(JS,IM1,2).EQ.I.AND.IMI(JS,IM1,1).NE.ICOMP) THEN
0291                 ICOMQ(2)=IMI(JS,IM1,1)
0292                 IVALQ(1)=0
0293               ENDIF
0294   220       CONTINUE
0295 C...If sea diquark then sum of two derived from companion x.
0296           ELSEIF(KFA.GT.1000.AND.MOD(KFA/10,10).EQ.0) THEN
0297              ICOMQ(1)=MOD(ICOMP,MSTU(5))
0298              ICOMQ(2)=ICOMP/MSTU(5)
0299 C...If meson or baryon then use fragmentation function.
0300 C...Somewhat arbitrary split into old and new flavour, but OK normally.
0301           ELSE
0302             KFL3=MOD(KFA/10,10)
0303             IF(MOD(KFA/1000,10).EQ.0) THEN
0304               KFL1=MOD(KFA/100,10)
0305             ELSE
0306               KFL1=MOD(KFA,10000)-10*KFL3-1
0307               IF(MOD(KFA/1000,10).EQ.MOD(KFA/100,10).AND.
0308      &        MOD(KFA,10).EQ.2) KFL1=KFL1+2
0309             ENDIF
0310             PR=P(I,5)**2+P(I,1)**2+P(I,2)**2
0311             CALL PYZDIS(KFL1,KFL3,PR,X)
0312           ENDIF
0313
0314           DO 260 IQ=1,2
0315 C...Calculation of x of valence quark: assume form (1-x)^a/sqrt(x),
0316 C...where a=3.5 for u in proton, =2 for d in proton and =0.8 for meson.
0317 C...In other baryons combine u and d from proton appropriately.
0318             IF(IVALQ(IQ).NE.0) THEN
0319               NVAL=0
0320               IF(KFIVAL(JS,1).EQ.IVALQ(IQ)) NVAL=NVAL+1
0321               IF(KFIVAL(JS,2).EQ.IVALQ(IQ)) NVAL=NVAL+1
0322               IF(KFIVAL(JS,3).EQ.IVALQ(IQ)) NVAL=NVAL+1
0323 C...Meson.
0324               IF(KFIVAL(JS,3).EQ.0) THEN
0325                 MDU=0
0326 C...Baryon with three identical quarks: mix u and d forms.
0327               ELSEIF(NVAL.EQ.3) THEN
0328                 MDU=INT(PYR(0)+5D0/3D0)
0329 C...Baryon, one of two identical quarks: u form.
0330               ELSEIF(NVAL.EQ.2) THEN
0331                 MDU=2
0332 C...Baryon with two identical quarks, but not the one picked: d form.
0333               ELSEIF(KFIVAL(JS,1).EQ.KFIVAL(JS,2).OR.KFIVAL(JS,2).EQ.
0334      &        KFIVAL(JS,3).OR.KFIVAL(JS,1).EQ.KFIVAL(JS,3)) THEN
0335                 MDU=1
0336 C...Baryon with three nonidentical quarks: mix u and d forms.
0337               ELSE
0338                 MDU=INT(PYR(0)+5D0/3D0)
0339               ENDIF
0340               XPOW=0.8D0
0341               IF(MDU.EQ.1) XPOW=3.5D0
0342               IF(MDU.EQ.2) XPOW=2D0
0343   230         XX=PYR(0)**2
0344               IF((1D0-XX)**XPOW.LT.PYR(0)) GOTO 230
0345               X=X+XX
0346             ENDIF
0347
0348 C...Calculation of x of companion quark.
0349             IF(ICOMQ(IQ).NE.0) THEN
0350               XCOMP=1D-4
0351               DO 240 IM1=1,MINT(31)
0352                 IF(IMI(JS,IM1,1).EQ.ICOMQ(IQ)) XCOMP=XMI(JS,IM1)
0353   240         CONTINUE
0354               NPOW=MAX(0,MIN(4,MSTP(87)))
0355   250         XX=XCOMP*(1D0/(1D0-PYR(0)*(1D0-XCOMP))-1D0)
0356               CORR=((1D0-XCOMP-XX)/(1D0-XCOMP))**NPOW*
0357      &        (XCOMP**2+XX**2)/(XCOMP+XX)**2
0358               IF(CORR.LT.PYR(0)) GOTO 250
0359               X=X+XX
0360             ENDIF
0361   260     CONTINUE
0362
0363 C...Optionally enchance x of composite systems (e.g. diquarks)
0364           IF (KFA.GT.100) X=PARP(79)*X
0365
0366 C...Store x. Also calculate light cone energies of each system.
0367           XMI(JS,IM)=X
0368           W(JS,JS)=W(JS,JS)+X
0369           W(JS,3-JS)=W(JS,3-JS)+(P(I,5)**2+P(I,1)**2+P(I,2)**2)/X
0370   270   CONTINUE
0371         W(JS,JS)=W(JS,JS)*W(0,JS)
0372         W(JS,3-JS)=W(JS,3-JS)/W(0,JS)
0373         W(JS,0)=W(JS,1)*W(JS,2)
0374   280 CONTINUE
0375
0376 C...Check W1 W2 < Wrem (can be done before rescaling, since W
0377 C...insensitive to global rescalings of the BR x values).
0378       IF (SQRT(W(1,0))+SQRT(W(2,0)).GT.SQRT(W(0,0)).AND.NTRYX.LE.100)
0379      &     THEN
0380         GOTO 210
0381       ELSEIF (NTRYX.GT.100.AND.NTRY.LE.100) THEN
0382         GOTO 100
0383       ELSEIF (NTRYX.GT.100) THEN
0384         CALL PYERRM(11,'(PYMIRM:) No consistent (x,kT) sets found')
0385         MINT(57)=MINT(57)+1
0386         MINT(51)=1
0387         RETURN
0388       ENDIF
0389
0390 C...Compute x rescaling factors
0391       COMTRM=W(0,0)+SQRT(FLAM(W(0,0),W(1,0),W(2,0)))
0392       R1=(COMTRM+W(1,0)-W(2,0))/(2D0*W(1,1)*W(0,2))
0393       R2=(COMTRM+W(2,0)-W(1,0))/(2D0*W(2,2)*W(0,1))
0394
0395       IF (R1.LT.0.OR.R2.LT.0) THEN
0396         CALL PYERRM(19,'(PYMIRM:) negative rescaling factors !')
0397         MINT(57)=MINT(57)+1
0398         MINT(51)=1
0399       ENDIF
0400
0401 C...Rescale W(1,*) and W(2,*) (not really necessary, but consistent).
0402       W(1,1)=W(1,1)*R1
0403       W(1,2)=W(1,2)/R1
0404       W(2,1)=W(2,1)/R2
0405       W(2,2)=W(2,2)*R2
0406
0407 C...Rescale BR x values.
0408       DO 290 IM=MINT(31)+1,MAX(NMI(1),NMI(2))
0409         XMI(1,IM)=XMI(1,IM)*R1
0410         XMI(2,IM)=XMI(2,IM)*R2
0411   290 CONTINUE
0412
0413 C...Now we have a consistent set of x and kT values.
0414 C...First set up the initiators and their daughters correctly.
0415       DO 300 IM=1,MINT(31)
0416         I1=IMI(1,IM,1)
0417         I2=IMI(2,IM,1)
0418         ST=XMI(1,IM)*XMI(2,IM)*VINT(2)+(P(I1,1)+P(I2,1))**2+
0419      &       (P(I1,2)+P(I2,2))**2
0420         PT12=P(I1,1)**2+P(I1,2)**2
0421         PT22=P(I2,1)**2+P(I2,2)**2
0422 C...p_z
0423         P(I1,3)=SQRT(FLAM(ST,PT12,PT22)/(4D0*ST))
0424         P(I2,3)=-P(I1,3)
0425 C...Energies (masses should be zero at this stage)
0426         P(I1,4)=SQRT(PT12+P(I1,3)**2)
0427         P(I2,4)=SQRT(PT22+P(I2,3)**2)
0428
0429 C...Transverse 12 system initiator velocity:
0430         VB(1)=(P(I1,1)+P(I2,1))/SQRT(ST)
0431         VB(2)=(P(I1,2)+P(I2,2))/SQRT(ST)
0432 C...Boost to overall initiator system rest frame
0433         CALL PYROBO(I1,I1,0D0,0D0,-VB(1),-VB(2),0D0)
0434         CALL PYROBO(I2,I2,0D0,0D0,-VB(1),-VB(2),0D0)
0435 C...Compute phi,theta coordinates of I1 and rotate z axis.
0436         PHI=PYANGL(P(I1,1),P(I1,2))
0437         THE=PYANGL(P(I1,3),SQRT(P(I1,1)**2+P(I1,2)**2))
0438         CALL PYROBO(I1,I1,0D0,-PHI,0D0,0D0,0D0)
0439         CALL PYROBO(I2,I2,0D0,-PHI,0D0,0D0,0D0)
0440         CALL PYROBO(I1,I1,-THE,0D0,0D0,0D0,0D0)
0441         CALL PYROBO(I2,I2,-THE,0D0,0D0,0D0,0D0)
0442
0443 C...Now boost initiators + daughters back to LAB system
0444 C...(also update documentation lines for MI = 1.)
0445         VB(3)=(XMI(1,IM)-XMI(2,IM))/(XMI(1,IM)+XMI(2,IM))
0446         IMIN=IMISEP(IM-1)+1
0447         IF (IM.EQ.1) IMIN=MINT(83)+5
0448         IMAX=IMISEP(IM)
0449         CALL PYROBO(IMIN,IMAX,THE,PHI,VB(1),VB(2),0D0)
0450         CALL PYROBO(IMIN,IMAX,0D0,0D0,0D0,0D0,VB(3))
0451
0452   300 CONTINUE
0453
0454
0455 C...For the beam remnant partons/hadrons, we only need to set pz and E.
0456       DO 320 JS=1,2
0457         DO 310 IM=MINT(31)+1,NMI(JS)
0458           I=IMI(JS,IM,1)
0459 C...Skip collapsed gluons and junctions.
0460           IF (K(I,2).EQ.21.AND.K(I,1).EQ.14) GOTO 310
0461           IF (KFA.EQ.88) GOTO 310
0462           RMT2=P(I,5)**2+P(I,1)**2+P(I,2)**2
0463           P(I,4)=0.5D0*(XMI(JS,IM)*W(0,JS)+RMT2/(XMI(JS,IM)*W(0,JS)))
0464           P(I,3)=0.5D0*(XMI(JS,IM)*W(0,JS)-RMT2/(XMI(JS,IM)*W(0,JS)))
0465           IF (JS.EQ.2) P(I,3)=-P(I,3)
0466   310   CONTINUE
0467   320 CONTINUE
0468
0469
0470 C...Documentation lines
0471       DO 340 JS=1,2
0472         IN=MINT(83)+JS+2
0473         IO=IMI(JS,1,1)
0474         K(IN,1)=21
0475         K(IN,2)=K(IO,2)
0476         K(IN,3)=MINT(83)+JS
0477         K(IN,4)=0
0478         K(IN,5)=0
0479         DO 330 J=1,5
0480           P(IN,J)=P(IO,J)
0481           V(IN,J)=V(IO,J)
0482   330   CONTINUE
0483         MCT(IN,1)=MCT(IO,1)
0484         MCT(IN,2)=MCT(IO,2)
0485   340 CONTINUE
0486
0487 C...Final state colour reconnections.
0488       IF (MSTP(95).NE.1.OR.MINT(31).LE.1) GOTO 380
0489
0490 C...Number of colour tags for which a recoupling will be tried.
0491       NTOT=NCT
0492 C...Number of recouplings to try
0493       MINT(34)=0
0494       NRECP=0
0495       NITER=0
0496   350 NRECP=MINT(34)
0497       NITER=NITER+1
0498       IITER=0
0499   360 IITER=IITER+1
0500       IF (IITER.LE.PARP(78)*NTOT) THEN
0501 C...Select two colour tags at random
0502 C...NB: jj strings do not have colour tags assigned to them,
0503 C...thus they are as yet not affected by anything done here.
0504         JCT=PYR(0)*NCT+1
0505         KCT=MOD(INT(JCT+PYR(0)*NCT),NCT)+1
0506         IJ1=0
0507         IJ2=0
0508         IK1=0
0509         IK2=0
0510 C...Find final state partons with this (anti)colour
0511         DO 370 I=MINT(84)+1,N
0512           IF (K(I,1).EQ.3) THEN
0513             IF (MCT(I,1).EQ.JCT) IJ1=I
0514             IF (MCT(I,2).EQ.JCT) IJ2=I
0515             IF (MCT(I,1).EQ.KCT) IK1=I
0516             IF (MCT(I,2).EQ.KCT) IK2=I
0517           ENDIF
0518   370   CONTINUE
0519 C...Only consider recouplings not involving junctions for now.
0520         IF (IJ1.EQ.0.OR.IJ2.EQ.0.OR.IK1.EQ.0.OR.IK2.EQ.0) GOTO 360
0521
0522         RLO=2D0*FOUR(IJ1,IJ2)*2D0*FOUR(IK1,IK2)
0523         RLN=2D0*FOUR(IJ1,IK2)*2D0*FOUR(IK1,IJ2)
0524         IF (RLN.LT.RLO.AND.MCT(IJ2,1).NE.KCT.AND.MCT(IK2,1).NE.JCT) THEN
0525           MCT(IJ2,2)=KCT
0526           MCT(IK2,2)=JCT
0527 C...Count up number of reconnections
0528           MINT(34)=MINT(34)+1
0529         ENDIF
0530         IF (MINT(34).LE.1000) THEN
0531           GOTO 360
0532         ELSE
0533           CALL PYERRM(4,'(PYMIRM:) caught in infinite loop')
0534           GOTO 380
0535         ENDIF
0536       ENDIF
0537       IF (NRECP.LT.MINT(34)) GOTO 350
0538
0539 C...Signal PYPREP to use /PYCTAG/ information rather than K(I,KCS).
0540   380 MINT(33)=1
0541
0542       RETURN
0543       END