pythiaeRHIC/pythia/pyreco.f

0001
0002 C***********************************************************************
0003
0004 C...PYRECO
0005 C...Handles the possibility of colour reconnection in W+W- events,
0006 C...Based on the main scenarios of the Sjostrand and Khoze study:
0007 C...I, II, II', intermediate and instantaneous; plus one model
0008 C...along the lines of the Gustafson and Hakkinen: GH.
0009 C...Note: also handles Z0 Z0 and W-W+ events, but notation below
0010 C...is as if first resonance is W+ and second W-.
0011
0012       SUBROUTINE PYRECO(IW1,IW2,NSD1,NAFT1)
0013
0014 C...Double precision and integer declarations.
0015       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0016       IMPLICIT INTEGER(I-N)
0017       INTEGER PYK,PYCHGE,PYCOMP
0018 C...Parameter value; number of points in MC integration.
0019       PARAMETER (NPT=100)
0020 C...Commonblocks.
0021       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
0022       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0023       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0024       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0025       COMMON/PYINT1/MINT(400),VINT(400)
0026       SAVE /PYJETS/,/PYDAT1/,/PYDAT2/,/PYPARS/,/PYINT1/
0027 C...Local arrays.
0028       DIMENSION NBEG(2),NEND(2),INP(50),INM(50),BEWW(3),XP(3),XM(3),
0029      &V1(3),V2(3),BETP(50,4),DIRP(50,3),BETM(50,4),DIRM(50,3),
0030      &XD(4),XB(4),IAP(NPT),IAM(NPT),WTA(NPT),V1P(3),V2P(3),V1M(3),
0031      &V2M(3),Q(4,3),XPP(3),XMM(3),IPC(20),IMC(20),TC(0:20),TPC(20),
0032      &TMC(20),IJOIN(100)
0033
0034 C...Functions to give four-product and to do determinants.
0035       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)
0036       DETER(I,J,L)=Q(I,1)*Q(J,2)*Q(L,3)-Q(I,1)*Q(L,2)*Q(J,3)+
0037      &Q(J,1)*Q(L,2)*Q(I,3)-Q(J,1)*Q(I,2)*Q(L,3)+
0038      &Q(L,1)*Q(I,2)*Q(J,3)-Q(L,1)*Q(J,2)*Q(I,3)
0039
0040 C...Only allow fraction of recoupling for GH, intermediate and
0041 C...instantaneous.
0042       IF(MSTP(115).EQ.5.OR.MSTP(115).EQ.11.OR.MSTP(115).EQ.12) THEN
0043         IF(PYR(0).GT.PARP(120)) RETURN
0044       ENDIF
0045       ISUB=MINT(1)
0046
0047 C...Common part for scenarios I, II, II', and GH.
0048       IF(MSTP(115).EQ.1.OR.MSTP(115).EQ.2.OR.MSTP(115).EQ.3.OR.
0049      &MSTP(115).EQ.5) THEN
0050
0051 C...Read out frequently-used parameters.
0052         PI=PARU(1)
0053         HBAR=PARU(3)
0054         PMW=PMAS(24,1)
0055         IF(ISUB.EQ.22) PMW=PMAS(23,1)
0056         PGW=PMAS(24,2)
0057         IF(ISUB.EQ.22) PGW=PMAS(23,2)
0058         TFRAG=PARP(115)
0059         RHAD=PARP(116)
0060         FACT=PARP(117)
0061         BLOWR=PARP(118)
0062         BLOWT=PARP(119)
0063
0064 C...Find range of decay products of the W's.
0065 C...Background: the W's are stored in IW1 and IW2.
0066 C...Their direct decay products in NSD1+1 through NSD1+4.
0067 C...Products after shower (if any) in NSD1+5 through NAFT1
0068 C...for first W and in NAFT1+1 through N for the second.
0069         IF(NAFT1.GT.NSD1+4) THEN
0070           NBEG(1)=NSD1+5
0071           NEND(1)=NAFT1
0072         ELSE
0073           NBEG(1)=NSD1+1
0074           NEND(1)=NSD1+2
0075         ENDIF
0076         IF(N.GT.NAFT1) THEN
0077           NBEG(2)=NAFT1+1
0078           NEND(2)=N
0079         ELSE
0080           NBEG(2)=NSD1+3
0081           NEND(2)=NSD1+4
0082         ENDIF
0083
0084 C...Rearrange parton shower products along strings.
0085         NOLD=N
0086         CALL PYPREP(NSD1+1)
0087         IF(MINT(51).NE.0) RETURN
0088
0089 C...Find partons pointing back to W+ and W-; store them with quark
0090 C...end of string first.
0091         NNP=0
0092         NNM=0
0093         ISGP=0
0094         ISGM=0
0095         DO 120 I=NOLD+1,N
0096           IF(K(I,1).NE.1.AND.K(I,1).NE.2) GOTO 120
0097           IF(IABS(K(I,2)).GE.22) GOTO 120
0098           IF(K(I,3).GE.NBEG(1).AND.K(I,3).LE.NEND(1)) THEN
0099             IF(ISGP.EQ.0) ISGP=ISIGN(1,K(I,2))
0100             NNP=NNP+1
0101             IF(ISGP.EQ.1) THEN
0102               INP(NNP)=I
0103             ELSE
0104               DO 100 I1=NNP,2,-1
0105                 INP(I1)=INP(I1-1)
0106   100         CONTINUE
0107               INP(1)=I
0108             ENDIF
0109             IF(K(I,1).EQ.1) ISGP=0
0110           ELSEIF(K(I,3).GE.NBEG(2).AND.K(I,3).LE.NEND(2)) THEN
0111             IF(ISGM.EQ.0) ISGM=ISIGN(1,K(I,2))
0112             NNM=NNM+1
0113             IF(ISGM.EQ.1) THEN
0114               INM(NNM)=I
0115             ELSE
0116               DO 110 I1=NNM,2,-1
0117                 INM(I1)=INM(I1-1)
0118   110         CONTINUE
0119               INM(1)=I
0120             ENDIF
0121             IF(K(I,1).EQ.1) ISGM=0
0122           ENDIF
0123   120   CONTINUE
0124
0125 C...Boost to W+W- rest frame (not strictly needed).
0126         DO 130 J=1,3
0127           BEWW(J)=(P(IW1,J)+P(IW2,J))/(P(IW1,4)+P(IW2,4))
0128   130   CONTINUE
0129         CALL PYROBO(IW1,IW1,0D0,0D0,-BEWW(1),-BEWW(2),-BEWW(3))
0130         CALL PYROBO(IW2,IW2,0D0,0D0,-BEWW(1),-BEWW(2),-BEWW(3))
0131         CALL PYROBO(NOLD+1,N,0D0,0D0,-BEWW(1),-BEWW(2),-BEWW(3))
0132
0133 C...Select decay vertices of W+ and W-.
0134         TP=HBAR*(-LOG(PYR(0)))*P(IW1,4)/
0135      &  SQRT((P(IW1,5)**2-PMW**2)**2+(P(IW1,5)**2*PGW/PMW)**2)
0136         TM=HBAR*(-LOG(PYR(0)))*P(IW2,4)/
0137      &  SQRT((P(IW2,5)**2-PMW**2)**2+(P(IW2,5)**2*PGW/PMW)**2)
0138         GTMAX=MAX(TP,TM)
0139         DO 140 J=1,3
0140           XP(J)=TP*P(IW1,J)/P(IW1,4)
0141           XM(J)=TM*P(IW2,J)/P(IW2,4)
0142   140   CONTINUE
0143
0144 C...Begin scenario I specifics.
0145         IF(MSTP(115).EQ.1) THEN
0146
0147 C...Reconstruct velocity and direction of W+ string pieces.
0148           DO 170 IIP=1,NNP-1
0149             IF(K(INP(IIP),2).LT.0) GOTO 170
0150             I1=INP(IIP)
0151             I2=INP(IIP+1)
0152             P1A=SQRT(P(I1,1)**2+P(I1,2)**2+P(I1,3)**2)
0153             P2A=SQRT(P(I2,1)**2+P(I2,2)**2+P(I2,3)**2)
0154             DO 150 J=1,3
0155               V1(J)=P(I1,J)/P1A
0156               V2(J)=P(I2,J)/P2A
0157               BETP(IIP,J)=0.5D0*(V1(J)+V2(J))
0158               DIRP(IIP,J)=V1(J)-V2(J)
0159   150       CONTINUE
0160             BETP(IIP,4)=1D0/SQRT(1D0-BETP(IIP,1)**2-BETP(IIP,2)**2-
0161      &      BETP(IIP,3)**2)
0162             DIRL=SQRT(DIRP(IIP,1)**2+DIRP(IIP,2)**2+DIRP(IIP,3)**2)
0163             DO 160 J=1,3
0164               DIRP(IIP,J)=DIRP(IIP,J)/DIRL
0165   160       CONTINUE
0166   170     CONTINUE
0167
0168 C...Reconstruct velocity and direction of W- string pieces.
0169           DO 200 IIM=1,NNM-1
0170             IF(K(INM(IIM),2).LT.0) GOTO 200
0171             I1=INM(IIM)
0172             I2=INM(IIM+1)
0173             P1A=SQRT(P(I1,1)**2+P(I1,2)**2+P(I1,3)**2)
0174             P2A=SQRT(P(I2,1)**2+P(I2,2)**2+P(I2,3)**2)
0175             DO 180 J=1,3
0176               V1(J)=P(I1,J)/P1A
0177               V2(J)=P(I2,J)/P2A
0178               BETM(IIM,J)=0.5D0*(V1(J)+V2(J))
0179               DIRM(IIM,J)=V1(J)-V2(J)
0180   180       CONTINUE
0181             BETM(IIM,4)=1D0/SQRT(1D0-BETM(IIM,1)**2-BETM(IIM,2)**2-
0182      &      BETM(IIM,3)**2)
0183             DIRL=SQRT(DIRM(IIM,1)**2+DIRM(IIM,2)**2+DIRM(IIM,3)**2)
0184             DO 190 J=1,3
0185               DIRM(IIM,J)=DIRM(IIM,J)/DIRL
0186   190       CONTINUE
0187   200     CONTINUE
0188
0189 C...Loop over number of space-time points.
0190           NACC=0
0191           SUM=0D0
0192           DO 250 IPT=1,NPT
0193
0194 C...Pick x,y,z,t Gaussian (width RHAD and TFRAG, respectively).
0195             R=SQRT(-LOG(PYR(0)))
0196             PHI=2D0*PI*PYR(0)
0197             X=BLOWR*RHAD*R*COS(PHI)
0198             Y=BLOWR*RHAD*R*SIN(PHI)
0199             R=SQRT(-LOG(PYR(0)))
0200             PHI=2D0*PI*PYR(0)
0201             Z=BLOWR*RHAD*R*COS(PHI)
0202             T=GTMAX+BLOWT*SQRT(0.5D0)*TFRAG*R*ABS(SIN(PHI))
0203
0204 C...Reject impossible points. Weight for sample distribution.
0205             IF(T**2-X**2-Y**2-Z**2.LT.0D0) GOTO 250
0206             WTSMP=EXP(-(X**2+Y**2+Z**2)/(BLOWR*RHAD)**2)*
0207      &      EXP(-2D0*(T-GTMAX)**2/(BLOWT*TFRAG)**2)
0208
0209 C...Loop over W+ string pieces and find one with largest weight.
0210             IMAXP=0
0211             WTMAXP=1D-10
0212             XD(1)=X-XP(1)
0213             XD(2)=Y-XP(2)
0214             XD(3)=Z-XP(3)
0215             XD(4)=T-TP
0216             DO 220 IIP=1,NNP-1
0217               IF(K(INP(IIP),2).LT.0) GOTO 220
0218               BED=BETP(IIP,1)*XD(1)+BETP(IIP,2)*XD(2)+BETP(IIP,3)*XD(3)
0219               BEDG=BETP(IIP,4)*(BETP(IIP,4)*BED/(1D0+BETP(IIP,4))-XD(4))
0220               DO 210 J=1,3
0221                 XB(J)=XD(J)+BEDG*BETP(IIP,J)
0222   210         CONTINUE
0223               XB(4)=BETP(IIP,4)*(XD(4)-BED)
0224               SR2=XB(1)**2+XB(2)**2+XB(3)**2
0225               SZ2=(DIRP(IIP,1)*XB(1)+DIRP(IIP,2)*XB(2)+
0226      &        DIRP(IIP,3)*XB(3))**2
0227               WTP=EXP(-(SR2-SZ2)/(2D0*RHAD**2))*EXP(-(XB(4)**2-SZ2)/
0228      &        TFRAG**2)
0229               IF(XB(4)-SQRT(SR2).LT.0D0) WTP=0D0
0230               IF(WTP.GT.WTMAXP) THEN
0231                 IMAXP=IIP
0232                 WTMAXP=WTP
0233               ENDIF
0234   220       CONTINUE
0235
0236 C...Loop over W- string pieces and find one with largest weight.
0237             IMAXM=0
0238             WTMAXM=1D-10
0239             XD(1)=X-XM(1)
0240             XD(2)=Y-XM(2)
0241             XD(3)=Z-XM(3)
0242             XD(4)=T-TM
0243             DO 240 IIM=1,NNM-1
0244               IF(K(INM(IIM),2).LT.0) GOTO 240
0245               BED=BETM(IIM,1)*XD(1)+BETM(IIM,2)*XD(2)+BETM(IIM,3)*XD(3)
0246               BEDG=BETM(IIM,4)*(BETM(IIM,4)*BED/(1D0+BETM(IIM,4))-XD(4))
0247               DO 230 J=1,3
0248                 XB(J)=XD(J)+BEDG*BETM(IIM,J)
0249   230         CONTINUE
0250               XB(4)=BETM(IIM,4)*(XD(4)-BED)
0251               SR2=XB(1)**2+XB(2)**2+XB(3)**2
0252               SZ2=(DIRM(IIM,1)*XB(1)+DIRM(IIM,2)*XB(2)+
0253      &        DIRM(IIM,3)*XB(3))**2
0254               WTM=EXP(-(SR2-SZ2)/(2D0*RHAD**2))*EXP(-(XB(4)**2-SZ2)/
0255      &        TFRAG**2)
0256               IF(XB(4)-SQRT(SR2).LT.0D0) WTM=0D0
0257               IF(WTM.GT.WTMAXM) THEN
0258                 IMAXM=IIM
0259                 WTMAXM=WTM
0260               ENDIF
0261   240       CONTINUE
0262
0263 C...Result of integration.
0264             WT=0D0
0265             IF(IMAXP.NE.0.AND.IMAXM.NE.0) THEN
0266               WT=WTMAXP*WTMAXM/WTSMP
0267               SUM=SUM+WT
0268               NACC=NACC+1
0269               IAP(NACC)=IMAXP
0270               IAM(NACC)=IMAXM
0271               WTA(NACC)=WT
0272             ENDIF
0273   250     CONTINUE
0274           RES=BLOWR**3*BLOWT*SUM/NPT
0275
0276 C...Decide whether to reconnect and, if so, where.
0277           IACC=0
0278           PREC=1D0-EXP(-FACT*RES)
0279           IF(PREC.GT.PYR(0)) THEN
0280             RSUM=PYR(0)*SUM
0281             DO 260 IA=1,NACC
0282               IACC=IA
0283               RSUM=RSUM-WTA(IA)
0284               IF(RSUM.LE.0D0) GOTO 270
0285   260       CONTINUE
0286   270       IIP=IAP(IACC)
0287             IIM=IAM(IACC)
0288           ENDIF
0289
0290 C...Begin scenario II and II' specifics.
0291         ELSEIF(MSTP(115).EQ.2.OR.MSTP(115).EQ.3) THEN
0292
0293 C...Loop through all string pieces, one from W+ and one from W-.
0294           NCROSS=0
0295           TC(0)=0D0
0296           DO 340 IIP=1,NNP-1
0297             IF(K(INP(IIP),2).LT.0) GOTO 340
0298             I1P=INP(IIP)
0299             I2P=INP(IIP+1)
0300             DO 330 IIM=1,NNM-1
0301               IF(K(INM(IIM),2).LT.0) GOTO 330
0302               I1M=INM(IIM)
0303               I2M=INM(IIM+1)
0304
0305 C...Find endpoint velocity vectors.
0306               DO 280 J=1,3
0307                 V1P(J)=P(I1P,J)/P(I1P,4)
0308                 V2P(J)=P(I2P,J)/P(I2P,4)
0309                 V1M(J)=P(I1M,J)/P(I1M,4)
0310                 V2M(J)=P(I2M,J)/P(I2M,4)
0311   280         CONTINUE
0312
0313 C...Define q matrix and find t.
0314               DO 290 J=1,3
0315                 Q(1,J)=V2P(J)-V1P(J)
0316                 Q(2,J)=-(V2M(J)-V1M(J))
0317                 Q(3,J)=XP(J)-XM(J)-TP*V1P(J)+TM*V1M(J)
0318                 Q(4,J)=V1P(J)-V1M(J)
0319   290         CONTINUE
0320               T=-DETER(1,2,3)/DETER(1,2,4)
0321
0322 C...Find alpha and beta; i.e. coordinates of crossing point.
0323               S11=Q(1,1)*(T-TP)
0324               S12=Q(2,1)*(T-TM)
0325               S13=Q(3,1)+Q(4,1)*T
0326               S21=Q(1,2)*(T-TP)
0327               S22=Q(2,2)*(T-TM)
0328               S23=Q(3,2)+Q(4,2)*T
0329               DEN=S11*S22-S12*S21
0330               ALP=(S12*S23-S22*S13)/DEN
0331               BET=(S21*S13-S11*S23)/DEN
0332
0333 C...Check if solution acceptable.
0334               IANSW=1
0335               IF(T.LT.GTMAX) IANSW=0
0336               IF(ALP.LT.0D0.OR.ALP.GT.1D0) IANSW=0
0337               IF(BET.LT.0D0.OR.BET.GT.1D0) IANSW=0
0338
0339 C...Find point of crossing and check that not inconsistent.
0340               DO 300 J=1,3
0341                 XPP(J)=XP(J)+(V1P(J)+ALP*(V2P(J)-V1P(J)))*(T-TP)
0342                 XMM(J)=XM(J)+(V1M(J)+BET*(V2M(J)-V1M(J)))*(T-TM)
0343   300         CONTINUE
0344               D2PM=(XPP(1)-XMM(1))**2+(XPP(2)-XMM(2))**2+
0345      &        (XPP(3)-XMM(3))**2
0346               D2P=XPP(1)**2+XPP(2)**2+XPP(3)**2
0347               D2M=XMM(1)**2+XMM(2)**2+XMM(3)**2
0348               IF(D2PM.GT.1D-4*(D2P+D2M)) IANSW=-1
0349
0350 C...Find string eigentimes at crossing.
0351               IF(IANSW.EQ.1) THEN
0352                 TAUP=SQRT(MAX(0D0,(T-TP)**2-(XPP(1)-XP(1))**2-
0353      &          (XPP(2)-XP(2))**2-(XPP(3)-XP(3))**2))
0354                 TAUM=SQRT(MAX(0D0,(T-TM)**2-(XMM(1)-XM(1))**2-
0355      &          (XMM(2)-XM(2))**2-(XMM(3)-XM(3))**2))
0356               ELSE
0357                 TAUP=0D0
0358                 TAUM=0D0
0359               ENDIF
0360
0361 C...Order crossings by time. End loop over crossings.
0362               IF(IANSW.EQ.1.AND.NCROSS.LT.20) THEN
0363                 NCROSS=NCROSS+1
0364                 DO 310 I1=NCROSS,1,-1
0365                   IF(T.GT.TC(I1-1).OR.I1.EQ.1) THEN
0366                     IPC(I1)=IIP
0367                     IMC(I1)=IIM
0368                     TC(I1)=T
0369                     TPC(I1)=TAUP
0370                     TMC(I1)=TAUM
0371                     GOTO 320
0372                   ELSE
0373                     IPC(I1)=IPC(I1-1)
0374                     IMC(I1)=IMC(I1-1)
0375                     TC(I1)=TC(I1-1)
0376                     TPC(I1)=TPC(I1-1)
0377                     TMC(I1)=TMC(I1-1)
0378                   ENDIF
0379   310           CONTINUE
0380   320           CONTINUE
0381               ENDIF
0382   330       CONTINUE
0383   340     CONTINUE
0384
0385 C...Loop over crossings; find first (if any) acceptable one.
0386           IACC=0
0387           IF(NCROSS.GE.1) THEN
0388             DO 350 IC=1,NCROSS
0389               PNFRAG=EXP(-(TPC(IC)**2+TMC(IC)**2)/TFRAG**2)
0390               IF(PNFRAG.GT.PYR(0)) THEN
0391 C...Scenario II: only compare with fragmentation time.
0392                 IF(MSTP(115).EQ.2) THEN
0393                   IACC=IC
0394                   IIP=IPC(IACC)
0395                   IIM=IMC(IACC)
0396                   GOTO 360
0397 C...Scenario II': also require that string length decreases.
0398                 ELSE
0399                   IIP=IPC(IC)
0400                   IIM=IMC(IC)
0401                   I1P=INP(IIP)
0402                   I2P=INP(IIP+1)
0403                   I1M=INM(IIM)
0404                   I2M=INM(IIM+1)
0405                   ELOLD=FOUR(I1P,I2P)*FOUR(I1M,I2M)
0406                   ELNEW=FOUR(I1P,I2M)*FOUR(I1M,I2P)
0407                   IF(ELNEW.LT.ELOLD) THEN
0408                     IACC=IC
0409                     IIP=IPC(IACC)
0410                     IIM=IMC(IACC)
0411                     GOTO 360
0412                   ENDIF
0413                 ENDIF
0414               ENDIF
0415   350       CONTINUE
0416   360       CONTINUE
0417           ENDIF
0418
0419 C...Begin scenario GH specifics.
0420         ELSEIF(MSTP(115).EQ.5) THEN
0421
0422 C...Loop through all string pieces, one from W+ and one from W-.
0423           IACC=0
0424           ELMIN=1D0
0425           DO 380 IIP=1,NNP-1
0426             IF(K(INP(IIP),2).LT.0) GOTO 380
0427             I1P=INP(IIP)
0428             I2P=INP(IIP+1)
0429             DO 370 IIM=1,NNM-1
0430               IF(K(INM(IIM),2).LT.0) GOTO 370
0431               I1M=INM(IIM)
0432               I2M=INM(IIM+1)
0433
0434 C...Look for largest decrease of (exponent of) Lambda measure.
0435               ELOLD=FOUR(I1P,I2P)*FOUR(I1M,I2M)
0436               ELNEW=FOUR(I1P,I2M)*FOUR(I1M,I2P)
0437               ELDIF=ELNEW/MAX(1D-10,ELOLD)
0438               IF(ELDIF.LT.ELMIN) THEN
0439                 IACC=IIP+IIM
0440                 ELMIN=ELDIF
0441                 IPC(1)=IIP
0442                 IMC(1)=IIM
0443               ENDIF
0444   370       CONTINUE
0445   380     CONTINUE
0446           IIP=IPC(1)
0447           IIM=IMC(1)
0448         ENDIF
0449
0450 C...Common for scenarios I, II, II' and GH: reconnect strings.
0451         IF(IACC.NE.0) THEN
0452           MINT(32)=1
0453           NJOIN=0
0454           DO 390 IS=1,NNP+NNM
0455             NJOIN=NJOIN+1
0456             IF(IS.LE.IIP) THEN
0457               I=INP(IS)
0458             ELSEIF(IS.LE.IIP+NNM-IIM) THEN
0459               I=INM(IS-IIP+IIM)
0460             ELSEIF(IS.LE.IIP+NNM) THEN
0461               I=INM(IS-IIP-NNM+IIM)
0462             ELSE
0463               I=INP(IS-NNM)
0464             ENDIF
0465             IJOIN(NJOIN)=I
0466             IF(K(I,2).LT.0) THEN
0467               CALL PYJOIN(NJOIN,IJOIN)
0468               NJOIN=0
0469             ENDIF
0470   390     CONTINUE
0471
0472 C...Restore original event record if no reconnection.
0473         ELSE
0474           DO 400 I=NSD1+1,NOLD
0475             IF(K(I,1).EQ.13.OR.K(I,1).EQ.14) THEN
0476               K(I,4)=MOD(K(I,4),MSTU(5)**2)
0477               K(I,5)=MOD(K(I,5),MSTU(5)**2)
0478             ENDIF
0479   400     CONTINUE
0480           DO 410 I=NOLD+1,N
0481             K(K(I,3),1)=3
0482   410     CONTINUE
0483           N=NOLD
0484         ENDIF
0485
0486 C...Boost back system.
0487         CALL PYROBO(IW1,IW1,0D0,0D0,BEWW(1),BEWW(2),BEWW(3))
0488         CALL PYROBO(IW2,IW2,0D0,0D0,BEWW(1),BEWW(2),BEWW(3))
0489         IF(N.GT.NOLD) CALL PYROBO(NOLD+1,N,0D0,0D0,
0490      &  BEWW(1),BEWW(2),BEWW(3))
0491
0492 C...Common part for intermediate and instantaneous scenarios.
0493       ELSEIF(MSTP(115).EQ.11.OR.MSTP(115).EQ.12) THEN
0494         MINT(32)=1
0495
0496 C...Remove old shower products and reset showering ones.
0497         N=NSD1+4
0498         DO 420 I=NSD1+1,NSD1+4
0499           K(I,1)=3
0500           K(I,4)=MOD(K(I,4),MSTU(5)**2)
0501           K(I,5)=MOD(K(I,5),MSTU(5)**2)
0502   420   CONTINUE
0503
0504 C...Identify quark-antiquark pairs.
0505         IQ1=NSD1+1
0506         IQ2=NSD1+2
0507         IQ3=NSD1+3
0508         IF(K(IQ1,2)*K(IQ3,2).LT.0) IQ3=NSD1+4
0509         IQ4=2*NSD1+7-IQ3
0510
0511 C...Reconnect strings.
0512         IJOIN(1)=IQ1
0513         IJOIN(2)=IQ4
0514         CALL PYJOIN(2,IJOIN)
0515         IJOIN(1)=IQ3
0516         IJOIN(2)=IQ2
0517         CALL PYJOIN(2,IJOIN)
0518
0519 C...Do new parton showers in intermediate scenario.
0520         IF(MSTP(71).GE.1.AND.MSTP(115).EQ.11) THEN
0521           MSTJ50=MSTJ(50)
0522           MSTJ(50)=0
0523           CALL PYSHOW(IQ1,IQ2,P(IW1,5))
0524           CALL PYSHOW(IQ3,IQ4,P(IW2,5))
0525           MSTJ(50)=MSTJ50
0526
0527 C...Do new parton showers in instantaneous scenario.
0528         ELSEIF(MSTP(71).GE.1.AND.MSTP(115).EQ.12) THEN
0529           PPM2=(P(IQ1,4)+P(IQ4,4))**2-(P(IQ1,1)+P(IQ4,1))**2-
0530      &    (P(IQ1,2)+P(IQ4,2))**2-(P(IQ1,3)+P(IQ4,3))**2
0531           PPM=SQRT(MAX(0D0,PPM2))
0532           CALL PYSHOW(IQ1,IQ4,PPM)
0533           PPM2=(P(IQ3,4)+P(IQ2,4))**2-(P(IQ3,1)+P(IQ2,1))**2-
0534      &    (P(IQ3,2)+P(IQ2,2))**2-(P(IQ3,3)+P(IQ2,3))**2
0535           PPM=SQRT(MAX(0D0,PPM2))
0536           CALL PYSHOW(IQ3,IQ2,PPM)
0537         ENDIF
0538       ENDIF
0539
0540       RETURN
0541       END