pythiaeRHIC/pythia/pyprep.f

0001
0002 C*********************************************************************
0003
0004 C...PYPREP
0005 C...Rearranges partons along strings.
0006 C...Special considerations for systems with junctions, with
0007 C...possibility of junction-antijunction annihilation.
0008 C...Allows small systems to collapse into one or two particles.
0009 C...Checks flavours and colour singlet invariant masses.
0010
0011       SUBROUTINE PYPREP(IP)
0012
0013 C...Double precision and integer declarations.
0014       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0015       INTEGER PYK,PYCHGE,PYCOMP
0016 C...Commonblocks.
0017       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
0018       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0019       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0020       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0021       COMMON/PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
0022       COMMON/PYINT1/MINT(400),VINT(400)
0023 C...The common block of colour tags.
0024       COMMON/PYCTAG/NCT,MCT(4000,2)
0025       SAVE /PYJETS/,/PYDAT1/,/PYDAT2/,/PYDAT3/,/PYINT1/,/PYCTAG/,
0026      &/PYPARS/
0027       DATA NERRPR/0/
0028       SAVE NERRPR
0029 C...Local arrays.
0030       DIMENSION DPS(5),DPC(5),UE(3),PG(5),E1(3),E2(3),E3(3),E4(3),
0031      &ECL(3),IJUNC(10,0:4),IPIECE(30,0:4),KFEND(4),KFQ(4),
0032      &IJUR(4),PJU(4,6),IRNG(4,2),TJJ(2,5),T(5),PUL(3,5),
0033      &IJCP(0:6),TJUOLD(5)
0034       CHARACTER CHTMP*6
0035
0036 C...Function to give four-product.
0037       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)
0038
0039 C...Rearrange parton shower product listing along strings: begin loop.
0040       MSTU(24)=0
0041       NOLD=N
0042       I1=N
0043       NJUNC=0
0044       NPIECE=0
0045       NJJSTR=0
0046       MSTU32=MSTU(32)+1
0047       DO 100 I=MAX(1,IP),N
0048 C...First store junction positions.
0049         IF(K(I,1).EQ.42) THEN
0050           NJUNC=NJUNC+1
0051           IJUNC(NJUNC,0)=I
0052           IJUNC(NJUNC,4)=0
0053         ENDIF
0054   100 CONTINUE
0055
0056       DO 250 MQGST=1,3
0057         DO 240 I=MAX(1,IP),N
0058 C...Special treatment for junctions
0059           IF (K(I,1).LE.0) GOTO 240
0060           IF(K(I,1).EQ.42) THEN
0061 C...MQGST=2: Look for junction-junction strings (not detected in the
0062 C...main search below).
0063             IF (MQGST.EQ.2.AND.NPIECE.NE.3*NJUNC) THEN
0064               IF (NJJSTR.EQ.0) THEN
0065                 NJJSTR = (3*NJUNC-NPIECE)/2
0066               ENDIF
0067 C...Check how many already identified strings end on this junction
0068               ILC=0
0069               DO 110 J=1,NPIECE
0070                 IF (IPIECE(J,4).EQ.I) ILC=ILC+1
0071   110         CONTINUE
0072 C...If less than 3, remaining must be to another junction
0073               IF (ILC.LT.3) THEN
0074                 IF (ILC.NE.2) THEN
0075 C...Multiple j-j connections not handled yet.
0076                   CALL PYERRM(2,
0077      &            '(PYPREP:) Too many junction-junction strings.')
0078                   MINT(51)=1
0079                   RETURN
0080                 ENDIF
0081 C...The colour information in the junction is unreadable for the
0082 C...colour space search further down in this routine, so we must
0083 C...start on the colour mother of this junction and then "artificially"
0084 C...prevent the colour mother from connecting here again.
0085                 ITJUNC=MOD(K(I,4)/MSTU(5),MSTU(5))
0086                 KCS=4
0087                 IF (MOD(ITJUNC,2).EQ.0) KCS=5
0088 C...Switch colour if the junction-junction leg is presumably a
0089 C...junction mother leg rather than a junction daughter leg.
0090                 IF (ITJUNC.GE.3) KCS=9-KCS
0091                 IF (MINT(33).EQ.0) THEN
0092 C...Find the unconnected leg and reorder junction daughter pointers so
0093 C...MOD(K(I,4),MSTU(5)) always points to the junction-junction string
0094 C...piece.
0095                   IA=MOD(K(I,4),MSTU(5))
0096                   IF (K(IA,KCS)/MSTU(5)**2.GE.2) THEN
0097                     ITMP=MOD(K(I,5),MSTU(5))
0098                     IF (K(ITMP,KCS)/MSTU(5)**2.GE.2) THEN
0099                       ITMP=MOD(K(I,5)/MSTU(5),MSTU(5))
0100                       K(I,5)=K(I,5)+(IA-ITMP)*MSTU(5)
0101                     ELSE
0102                       K(I,5)=K(I,5)+(IA-ITMP)
0103                     ENDIF
0104                     K(I,4)=K(I,4)+(ITMP-IA)
0105                     IA=ITMP
0106                   ENDIF
0107                   IF (ITJUNC.LE.2) THEN
0108 C...Beam baryon junction
0109                     K(IA,KCS)   = K(IA,KCS) + 2*MSTU(5)**2
0110                     K(I,KCS)    = K(I,KCS) + 1*MSTU(5)**2
0111 C...Else 1 -> 2 decay junction
0112                   ELSE
0113                     K(IA,KCS)   = K(IA,KCS) + MSTU(5)**2
0114                     K(I,KCS)    = K(I,KCS) + 2*MSTU(5)**2
0115                   ENDIF
0116                   I1BEG = I1
0117                   NSTP = 0
0118                   GOTO 170
0119 C...Alternatively use colour tag information.
0120                 ELSE
0121 C...Find a final state parton with appropriate dangling colour tag.
0122                   JCT=0
0123                   IA=0
0124                   IJUMO=K(I,3)
0125                   DO 140 J1=MAX(1,IP),N
0126                     IF (K(J1,1).NE.3) GOTO 140
0127 C...Check for matching final-state colour tag
0128                     IMATCH=0
0129                     DO 120 J2=MAX(1,IP),N
0130                       IF (K(J2,1).NE.3) GOTO 120
0131                       IF (MCT(J1,KCS-3).EQ.MCT(J2,6-KCS)) IMATCH=1
0132   120               CONTINUE
0133                     IF (IMATCH.EQ.1) GOTO 140
0134 C...Check whether this colour tag belongs to the present junction
0135 C...by seeing whether any parton with this colour tag has the same
0136 C...mother as the junction.
0137                     JCT=MCT(J1,KCS-3)
0138                     IMATCH=0
0139                     DO 130 J2=MINT(84)+1,N
0140                       IMO2=K(J2,3)
0141 C...First scattering partons have IMO1 = 3 and 4.
0142                       IF (IMO2.EQ.MINT(83)+3.OR.IMO2.EQ.MINT(83)+4)
0143      &                     IMO2=IMO2-2
0144                       IF (MCT(J2,KCS-3).EQ.JCT.AND.IMO2.EQ.IJUMO)
0145      &                     IMATCH=1
0146   130               CONTINUE
0147                     IF (IMATCH.EQ.0) GOTO 140
0148                     IA=J1
0149   140             CONTINUE
0150 C...Check for junction-junction strings without intermediate final state
0151 C...glue (not detected above).
0152                   IF (IA.EQ.0) THEN
0153                     DO 160 MJU=1,NJUNC
0154                       IJU2=IJUNC(MJU,0)
0155                       IF (IJU2.EQ.I) GOTO 160
0156                       ITJU2=MOD(K(IJU2,4)/MSTU(5),MSTU(5))
0157 C...Only opposite types of junctions can connect to each other.
0158                       IF (MOD(ITJU2,2).EQ.MOD(ITJUNC,2)) GOTO 160
0159                       IS=0
0160                       DO 150 J=1,NPIECE
0161                         IF (IPIECE(J,4).EQ.IJU2) IS=IS+1
0162   150                 CONTINUE
0163                       IF (IS.EQ.3) GOTO 160
0164                       IB=I
0165                       IA=IJU2
0166   160               CONTINUE
0167                   ENDIF
0168 C...Switch to other side of adjacent parton and step from there.
0169                   KCS=9-KCS
0170                   I1BEG = I1
0171                   NSTP = 0
0172                   GOTO 170
0173                 ENDIF
0174               ELSE IF (ILC.NE.3) THEN
0175               ENDIF
0176             ENDIF
0177           ENDIF
0178
0179 C...Look for coloured string endpoint, or (later) leftover gluon.
0180           IF(K(I,1).NE.3) GOTO 240
0181           KC=PYCOMP(K(I,2))
0182           IF(KC.EQ.0) GOTO 240
0183           KQ=KCHG(KC,2)
0184           IF(KQ.EQ.0.OR.(MQGST.LE.2.AND.KQ.EQ.2)) GOTO 240
0185
0186 C...Pick up loose string end.
0187           KCS=4
0188           IF(KQ*ISIGN(1,K(I,2)).LT.0) KCS=5
0189           IA=I
0190           IB=I
0191           I1BEG=I1
0192           NSTP=0
0193   170     NSTP=NSTP+1
0194           IF(NSTP.GT.4*N) THEN
0195             CALL PYERRM(14,'(PYPREP:) caught in infinite loop')
0196             MINT(51)=1
0197             RETURN
0198           ENDIF
0199
0200 C...Copy undecayed parton. Finished if reached string endpoint.
0201           IF(K(IA,1).EQ.3) THEN
0202             IF(I1.GE.MSTU(4)-MSTU32-5) THEN
0203               CALL PYERRM(11,'(PYPREP:) no more memory left in PYJETS')
0204               MINT(51)=1
0205               MSTU(24)=1
0206               RETURN
0207             ENDIF
0208             I1=I1+1
0209             K(I1,1)=2
0210             IF(NSTP.GE.2.AND.KCHG(PYCOMP(K(IA,2)),2).NE.2) K(I1,1)=1
0211             K(I1,2)=K(IA,2)
0212             K(I1,3)=IA
0213             K(I1,4)=0
0214             K(I1,5)=0
0215             DO 180 J=1,5
0216               P(I1,J)=P(IA,J)
0217               V(I1,J)=V(IA,J)
0218   180       CONTINUE
0219             K(IA,1)=K(IA,1)+10
0220             IF(K(I1,1).EQ.1) GOTO 240
0221           ENDIF
0222
0223 C...Also finished (for now) if reached junction; then copy to end.
0224           IF(K(IA,1).EQ.42) THEN
0225             NCOPY=I1-I1BEG
0226             IF(I1.GE.MSTU(4)-MSTU32-NCOPY-5) THEN
0227               CALL PYERRM(11,'(PYPREP:) no more memory left in PYJETS')
0228               MINT(51)=1
0229               MSTU(24)=1
0230               RETURN
0231             ENDIF
0232             IF (MQGST.LE.2.AND.NCOPY.NE.0) THEN
0233               DO 200 ICOPY=1,NCOPY
0234                 DO 190 J=1,5
0235                   K(MSTU(4)-MSTU32-ICOPY,J)=K(I1BEG+ICOPY,J)
0236                   P(MSTU(4)-MSTU32-ICOPY,J)=P(I1BEG+ICOPY,J)
0237                   V(MSTU(4)-MSTU32-ICOPY,J)=V(I1BEG+ICOPY,J)
0238   190           CONTINUE
0239   200         CONTINUE
0240             ENDIF
0241 C...For junction-junction strings, find end leg and reorder junction
0242 C...daughter pointers so MOD(K(I,4),MSTU(5)) always points to the
0243 C...junction-junction string piece.
0244             IF (K(I,1).EQ.42.AND.MINT(33).EQ.0) THEN
0245               ITMP=MOD(K(IA,4),MSTU(5))
0246               IF (ITMP.NE.IB) THEN
0247                 IF (MOD(K(IA,5),MSTU(5)).EQ.IB) THEN
0248                   K(IA,5)=K(IA,5)+(ITMP-IB)
0249                 ELSE
0250                   K(IA,5)=K(IA,5)+(ITMP-IB)*MSTU(5)
0251                 ENDIF
0252                 K(IA,4)=K(IA,4)+(IB-ITMP)
0253               ENDIF
0254             ENDIF
0255             NPIECE=NPIECE+1
0256 C...IPIECE:
0257 C...0: endpoint in original ER
0258 C...1:
0259 C...2:
0260 C...3: Parton immediately next to junction
0261 C...4: Junction
0262             IPIECE(NPIECE,0)=I
0263             IPIECE(NPIECE,1)=MSTU32+1
0264             IPIECE(NPIECE,2)=MSTU32+NCOPY
0265             IPIECE(NPIECE,3)=IB
0266             IPIECE(NPIECE,4)=IA
0267             MSTU32=MSTU32+NCOPY
0268             I1=I1BEG
0269             GOTO 240
0270           ENDIF
0271
0272 C...GOTO next parton in colour space.
0273           IB=IA
0274           IF (MINT(33).EQ.0) THEN
0275             IF(MOD(K(IB,KCS)/MSTU(5)**2,2).EQ.0.AND.MOD(K(IB,KCS),MSTU(5
0276      &           )).NE.0) THEN
0277               IA=MOD(K(IB,KCS),MSTU(5))
0278               K(IB,KCS)=K(IB,KCS)+MSTU(5)**2
0279               MREV=0
0280             ELSE
0281               IF(K(IB,KCS).GE.2*MSTU(5)**2.OR.MOD(K(IB,KCS)/MSTU(5),
0282      &             MSTU(5)).EQ.0) KCS=9-KCS
0283               IA=MOD(K(IB,KCS)/MSTU(5),MSTU(5))
0284               K(IB,KCS)=K(IB,KCS)+2*MSTU(5)**2
0285               MREV=1
0286             ENDIF
0287             IF(IA.LE.0.OR.IA.GT.N) THEN
0288               CALL PYERRM(12,'(PYPREP:) colour rearrangement failed')
0289               IF(NERRPR.LT.5) THEN
0290                 NERRPR=NERRPR+1
0291                 WRITE(MSTU(11),*) 'started at:', I
0292                 WRITE(MSTU(11),*) 'ended going from',IB,' to',IA
0293                 WRITE(MSTU(11),*) 'MQGST =',MQGST
0294                 CALL PYLIST(4)
0295               ENDIF
0296               MINT(51)=1
0297               RETURN
0298             ENDIF
0299             IF(MOD(K(IA,4)/MSTU(5),MSTU(5)).EQ.IB.OR.MOD(K(IA,5)/MSTU(5)
0300      &           ,MSTU(5)).EQ.IB) THEN
0301               IF(MREV.EQ.1) KCS=9-KCS
0302               IF(MOD(K(IA,KCS)/MSTU(5),MSTU(5)).NE.IB) KCS=9-KCS
0303               K(IA,KCS)=K(IA,KCS)+2*MSTU(5)**2
0304             ELSE
0305               IF(MREV.EQ.0) KCS=9-KCS
0306               IF(MOD(K(IA,KCS),MSTU(5)).NE.IB) KCS=9-KCS
0307               K(IA,KCS)=K(IA,KCS)+MSTU(5)**2
0308             ENDIF
0309             IF(IA.NE.I) GOTO 170
0310 C...Use colour tag information
0311           ELSE
0312 C...First create colour tags starting on IB if none already present.
0313             IF (MCT(IB,KCS-3).EQ.0) THEN
0314               CALL PYCTTR(IB,KCS,IB)
0315               IF(MINT(51).NE.0) RETURN
0316             ENDIF
0317             JCT=MCT(IB,KCS-3)
0318             IFOUND=0
0319 C...Find final state tag partner
0320             DO 210 IT=MAX(1,IP),N
0321               IF (IT.EQ.IB) GOTO 210
0322               IF (MCT(IT,6-KCS).EQ.JCT.AND.K(IT,1).LT.10.AND.K(IT,1).GT
0323      &             .0) THEN
0324                 IFOUND=IFOUND+1
0325                 IA=IT
0326               ENDIF
0327   210       CONTINUE
0328 C...Just copy and goto next if exactly one partner found.
0329             IF (IFOUND.EQ.1) THEN
0330               GOTO 170
0331 C...When no match found, match is presumably junction.
0332             ELSEIF (IFOUND.EQ.0.AND.MQGST.LE.2) THEN
0333 C...Check whether this colour tag matches a junction
0334 C...by seeing whether any parton with this colour tag has the same
0335 C...mother as a junction.
0336 C...NB: Only type 1 and 2 junctions handled presently.
0337               DO 230 IJU=1,NJUNC
0338                 IJUMO=K(IJUNC(IJU,0),3)
0339                 ITJUNC=MOD(K(IJUNC(IJU,0),4)/MSTU(5),MSTU(5))
0340 C...Colours only connect to junctions, anti-colours to antijunctions:
0341                 IF (MOD(ITJUNC+1,2)+1.NE.KCS-3) GOTO 230
0342                 IMATCH=0
0343                 DO 220 J1=MAX(1,IP),N
0344                   IF (K(J1,1).LE.0) GOTO 220
0345 C...First scattering partons have IMO1 = 3 and 4.
0346                   IMO=K(J1,3)
0347                   IF (IMO.EQ.MINT(83)+3.OR.IMO.EQ.MINT(83)+4)
0348      &                 IMO=IMO-2
0349                   IF (MCT(J1,KCS-3).EQ.JCT.AND.IMO.EQ.IJUMO.AND.MOD(K(J1
0350      &                 ,3+ITJUNC)/MSTU(5),MSTU(5)).EQ.IJUNC(IJU,0))
0351      &                 IMATCH=1
0352 C...Attempt at handling type > 3 junctions also. Not tested.
0353                   IF (ITJUNC.GE.3.AND.MCT(J1,6-KCS).EQ.JCT.AND.IMO.EQ
0354      &                 .IJUMO) IMATCH=1
0355   220           CONTINUE
0356                 IF (IMATCH.EQ.0) GOTO 230
0357                 IA=IJUNC(IJU,0)
0358                 IFOUND=IFOUND+1
0359   230         CONTINUE
0360
0361               IF (IFOUND.EQ.1) THEN
0362                 GOTO 170
0363               ELSEIF (IFOUND.EQ.0) THEN
0364                 WRITE(CHTMP,*) JCT
0365                 CALL PYERRM(12,'(PYPREP:) no matching colour tag: '
0366      &               //CHTMP)
0367                 IF(NERRPR.LT.5) THEN
0368                   NERRPR=NERRPR+1
0369                   CALL PYLIST(4)
0370                 ENDIF
0371                 MINT(51)=1
0372                 RETURN
0373               ENDIF
0374             ELSEIF (IFOUND.GE.2) THEN
0375               WRITE(CHTMP,*) JCT
0376               CALL PYERRM(12
0377      &             ,'(PYPREP:) too many occurences of colour line: '//
0378      &             CHTMP)
0379               IF(NERRPR.LT.5) THEN
0380                 NERRPR=NERRPR+1
0381                 CALL PYLIST(4)
0382               ENDIF
0383               MINT(51)=1
0384               RETURN
0385             ENDIF
0386           ENDIF
0387           K(I1,1)=1
0388   240   CONTINUE
0389   250 CONTINUE
0390
0391 C...Junction systems remain.
0392       IJU=0
0393       IJUS=0
0394       IJUCNT=0
0395       MREV=0
0396       IJJSTR=0
0397   260 IJUCNT=IJUCNT+1
0398       IF (IJUCNT.LE.NJUNC) THEN
0399 C...If we are not processing a j-j string, treat this junction as new.
0400         IF (IJJSTR.EQ.0) THEN
0401           IJU=IJUNC(IJUCNT,0)
0402           MREV=0
0403 C...If junction has already been read, ignore it.
0404           IF (IJUNC(IJUCNT,4).EQ.1) GOTO 260
0405 C...If we are on a j-j string, goto second j-j junction.
0406         ELSE
0407           IJUCNT=IJUCNT-1
0408           IJU=IJUS
0409         ENDIF
0410 C...Mark selected junction read.
0411         DO 270 J=1,NJUNC
0412           IF (IJUNC(J,0).EQ.IJU) IJUNC(J,4)=1
0413   270   CONTINUE
0414 C...Determine junction type
0415         ITJUNC = MOD(K(IJU,4)/MSTU(5),MSTU(5))
0416 C...Type 1 and 2 junctions: ~chi -> q q q, ~chi -> qbar,qbar,qbar
0417 C...Type 3 and 4 junctions: ~qbar -> q q , ~q -> qbar qbar
0418 C...Type 5 and 6 junctions: ~g -> q q q, ~g -> qbar qbar qbar
0419         IF (ITJUNC.GE.1.AND.ITJUNC.LE.6) THEN
0420           IHK=0
0421   280     IHK=IHK+1
0422 C...Find which quarks belong to given junction.
0423           IHF=0
0424           DO 290 IPC=1,NPIECE
0425             IF (IPIECE(IPC,4).EQ.IJU) THEN
0426               IHF=IHF+1
0427               IF (IHF.EQ.IHK) IEND=IPIECE(IPC,3)
0428             ENDIF
0429             IF (IHK.EQ.3.AND.IPIECE(IPC,0).EQ.IJU) IEND=IPIECE(IPC,3)
0430   290     CONTINUE
0431 C...IHK = 3 is special. Either normal string piece, or j-j string.
0432           IF(IHK.EQ.3) THEN
0433             IF (MREV.NE.1) THEN
0434               DO 300 IPC=1,NPIECE
0435 C...If there is a j-j string starting on the present junction which has
0436 C...zero length, insert next junction immediately.
0437                 IF (IPIECE(IPC,0).EQ.IJU.AND.K(IPIECE(IPC,4),1)
0438      &          .EQ.42.AND.IPIECE(IPC,1)-1-IPIECE(IPC,2).EQ.0) THEN
0439                   IJJSTR = 1
0440                   GOTO 340
0441                 ENDIF
0442   300         CONTINUE
0443               MREV = 1
0444 C...If MREV is 1 and IHK is 3 we are finished with this system.
0445             ELSE
0446               MREV=0
0447               GOTO 260
0448             ENDIF
0449           ENDIF
0450
0451 C...If we've gotten this far, then either IHK < 3, or
0452 C...an interjunction string exists, or just a third normal string.
0453           IJUNC(IJUCNT,IHK)=0
0454           IJJSTR = 0
0455 C..Order pieces belonging to this junction. Also look for j-j.
0456           DO 310 IPC=1,NPIECE
0457             IF (IPIECE(IPC,3).EQ.IEND) IJUNC(IJUCNT,IHK)=IPC
0458             IF (IHK.EQ.3.AND.IPIECE(IPC,0).EQ.IJUNC(IJUCNT,0)
0459      &      .AND.K(IPIECE(IPC,4),1).EQ.42) THEN
0460               IJUNC(IJUCNT,IHK)=IPC
0461               IJJSTR = 1
0462               MREV = 0
0463             ENDIF
0464   310     CONTINUE
0465 C...Copy back chains in proper order. MREV=0/1 : descending/ascending
0466           IPC=IJUNC(IJUCNT,IHK)
0467 C...Temporary solution to cover for bug.
0468           IF(IPC.LE.0) THEN
0469             CALL PYERRM(12,'(PYPREP:) fails to hook up junctions')
0470             MINT(51)=1
0471             RETURN
0472           ENDIF
0473           DO 330 ICP=IPIECE(IPC,1+MREV),IPIECE(IPC,2-MREV),1-2*MREV
0474             I1=I1+1
0475             DO 320 J=1,5
0476               K(I1,J)=K(MSTU(4)-ICP,J)
0477               P(I1,J)=P(MSTU(4)-ICP,J)
0478               V(I1,J)=V(MSTU(4)-ICP,J)
0479   320       CONTINUE
0480   330     CONTINUE
0481           K(I1,1)=2
0482 C...Mark last quark.
0483           IF (MREV.EQ.1.AND.IHK.GE.2) K(I1,1)=1
0484 C...Do not insert junctions at wrong places.
0485           IF(IHK.LT.2.OR.MREV.NE.0) GOTO 360
0486 C...Insert junction.
0487   340     IJUS = IJU
0488           IF (IHK.EQ.3) THEN
0489 C...Shift to end junction if a j-j string has been processed.
0490             IF (IJJSTR.NE.0) IJUS = IPIECE(IPC,4)
0491             MREV= 1
0492           ENDIF
0493           I1=I1+1
0494           DO 350 J=1,5
0495             K(I1,J)=0
0496             P(I1,J)=0.
0497             V(I1,J)=0.
0498   350     CONTINUE
0499           K(I1,1)=41
0500           K(IJUS,1)=K(IJUS,1)+10
0501           K(I1,2)=K(IJUS,2)
0502           K(I1,3)=IJUS
0503   360     IF (IHK.LT.3) GOTO 280
0504         ELSE
0505           CALL PYERRM(12,'(PYPREP:) Unknown junction type')
0506           MINT(51)=1
0507           RETURN
0508         ENDIF
0509         IF (IJUCNT.NE.NJUNC) GOTO 260
0510       ENDIF
0511       N=I1
0512
0513 C...Rearrange three strings from junction, e.g. in case one has been
0514 C...shortened by shower, so the last is the largest-energy one.
0515       IF(NJUNC.GE.1) THEN
0516 C...Find systems with exactly one junction.
0517         MJUN1=0
0518         NBEG=NOLD+1
0519         DO 470 I=NOLD+1,N
0520           IF(K(I,1).NE.1.AND.K(I,1).NE.41) THEN
0521           ELSEIF(K(I,1).EQ.41) THEN
0522             MJUN1=MJUN1+1
0523           ELSEIF(K(I,1).EQ.1.AND.MJUN1.NE.1) THEN
0524             MJUN1=0
0525             NBEG=I+1
0526           ELSE
0527             NEND=I
0528 C...Sum up energy-momentum in each junction string.
0529             DO 370 J=1,5
0530               PJU(1,J)=0D0
0531               PJU(2,J)=0D0
0532               PJU(3,J)=0D0
0533   370       CONTINUE
0534             NJU=0
0535             DO 390 I1=NBEG,NEND
0536               IF(K(I1,2).NE.21) THEN
0537                 NJU=NJU+1
0538                 IJUR(NJU)=I1
0539               ENDIF
0540               DO 380 J=1,5
0541                 PJU(MIN(NJU,3),J)=PJU(MIN(NJU,3),J)+P(I1,J)
0542   380         CONTINUE
0543   390       CONTINUE
0544 C...Find which of them has highest energy (minus mass) in rest frame.
0545             DO 400 J=1,5
0546               PJU(4,J)=PJU(1,J)+PJU(2,J)+PJU(3,J)
0547   400       CONTINUE
0548             PMJU=SQRT(MAX(0D0,PJU(4,4)**2-PJU(4,1)**2-PJU(4,2)**2-
0549      &      PJU(4,3)**2))
0550             DO 410 I2=1,3
0551               PJU(I2,6)=(PJU(4,4)*PJU(I2,4)-PJU(4,1)*PJU(I2,1)-
0552      &        PJU(4,2)*PJU(I2,2)-PJU(4,3)*PJU(I2,3))/PMJU-PJU(I2,5)
0553   410       CONTINUE
0554             IF(PJU(3,6).LT.MIN(PJU(1,6),PJU(2,6))) THEN
0555 C...Decide how to rearrange so that new last has highest energy.
0556               IF(PJU(1,6).LT.PJU(2,6)) THEN
0557                 IRNG(1,1)=IJUR(1)
0558                 IRNG(1,2)=IJUR(2)-1
0559                 IRNG(2,1)=IJUR(4)
0560                 IRNG(2,2)=IJUR(3)+1
0561                 IRNG(4,1)=IJUR(3)-1
0562                 IRNG(4,2)=IJUR(2)
0563               ELSE
0564                 IRNG(1,1)=IJUR(4)
0565                 IRNG(1,2)=IJUR(3)+1
0566                 IRNG(2,1)=IJUR(2)
0567                 IRNG(2,2)=IJUR(3)-1
0568                 IRNG(4,1)=IJUR(2)-1
0569                 IRNG(4,2)=IJUR(1)
0570               ENDIF
0571               IRNG(3,1)=IJUR(3)
0572               IRNG(3,2)=IJUR(3)
0573 C...Copy in correct order below bottom of current event record.
0574               I2=N
0575               DO 440 II=1,4
0576                 DO 430 I1=IRNG(II,1),IRNG(II,2),
0577      &          ISIGN(1,IRNG(II,2)-IRNG(II,1))
0578                   I2=I2+1
0579                   IF(I2.GE.MSTU(4)-MSTU32-5) THEN
0580                     CALL PYERRM(11,
0581      &              '(PYPREP:) no more memory left in PYJETS')
0582                     MINT(51)=1
0583                     MSTU(24)=1
0584                     RETURN
0585                   ENDIF
0586                   DO 420 J=1,5
0587                     K(I2,J)=K(I1,J)
0588                     P(I2,J)=P(I1,J)
0589                     V(I2,J)=V(I1,J)
0590   420             CONTINUE
0591                   IF(K(I2,1).EQ.1) K(I2,1)=2
0592   430           CONTINUE
0593   440         CONTINUE
0594               K(I2,1)=1
0595 C...Copy back up, overwriting but now in correct order.
0596               DO 460 I1=NBEG,NEND
0597                 I2=I1-NBEG+N+1
0598                 DO 450 J=1,5
0599                   K(I1,J)=K(I2,J)
0600                   P(I1,J)=P(I2,J)
0601                   V(I1,J)=V(I2,J)
0602   450           CONTINUE
0603   460         CONTINUE
0604             ENDIF
0605             MJUN1=0
0606             NBEG=I+1
0607           ENDIF
0608   470   CONTINUE
0609
0610 C...Check whether q-q-j-j-qbar-qbar systems should be collapsed
0611 C...to two q-qbar systems.
0612 C...(MSTJ(19)=1 forces q-q-j-j-qbar-qbar.)
0613         IF (MSTJ(19).NE.1) THEN
0614           MJUN1  = 0
0615           JJGLUE = 0
0616           NBEG   = NOLD+1
0617 C...Force collapse when MSTJ(19)=2.
0618           IF (MSTJ(19).EQ.2) THEN
0619             DELMJJ = 1D9
0620             DELMQQ = 0D0
0621           ENDIF
0622 C...Find systems with exactly two junctions.
0623           DO 700 I=NOLD+1,N
0624 C...Count junctions
0625             IF (K(I,1).EQ.41) THEN
0626               MJUN1 = MJUN1+1
0627 C...Check for interjunction gluons
0628               IF (MJUN1.EQ.2.AND.K(I-1,1).NE.41) THEN
0629                 JJGLUE = 1
0630               ENDIF
0631             ELSEIF(K(I,1).EQ.1.AND.(MJUN1.NE.2)) THEN
0632 C...If end of system reached with either zero or one junction, restart
0633 C...with next system.
0634               MJUN1  = 0
0635               JJGLUE = 0
0636               NBEG   = I+1
0637             ELSEIF(K(I,1).EQ.1) THEN
0638 C...If end of system reached with exactly two junctions, compute string
0639 C...length measure for the (q-q-j-j-qbar-qbar) topology and compare with
0640 C...length measure for the (q-qbar)(q-qbar) topology.
0641               NEND=I
0642 C...Loop down through chain.
0643               ISID=0
0644               DO 480 I1=NBEG,NEND
0645 C...Store string piece division locations in event record
0646                 IF (K(I1,2).NE.21) THEN
0647                   ISID       = ISID+1
0648                   IJCP(ISID) = I1
0649                 ENDIF
0650   480         CONTINUE
0651 C...Randomly choose between (1,3)(2,4) and (1,4)(2,3) topologies.
0652               ISW=0
0653               IF (PYR(0).LT.0.5D0) ISW=1
0654 C...Randomly choose which qqbar string gets the jj gluons.
0655               IGS=1
0656               IF (PYR(0).GT.0.5D0) IGS=2
0657 C...Only compute string lengths when no topology forced.
0658               IF (MSTJ(19).EQ.0) THEN
0659 C...Repeat following for each junction
0660                 DO 570 IJU=1,2
0661 C...Initialize iterative procedure for finding JRF
0662                   IJRFIT=0
0663                   DO 490 IX=1,3
0664                     TJUOLD(IX)=0D0
0665   490             CONTINUE
0666                   TJUOLD(4)=1D0
0667 C...Start iteration. Sum up momenta in string pieces
0668   500             DO 540 IJS=1,3
0669 C...JD=-1 for first junction, +1 for second junction.
0670 C...Find out where piece starts and ends and which direction to go.
0671                     JD=2*IJU-3
0672                     IF (IJS.LE.2) THEN
0673                       IA = IJCP((IJU-1)*7 - JD*(IJS+1)) + JD
0674                       IB = IJCP((IJU-1)*7 - JD*IJS)
0675                     ELSEIF (IJS.EQ.3) THEN
0676                       JD =-JD
0677                       IA = IJCP((IJU-1)*7 + JD*(IJS)) + JD
0678                       IB = IJCP((IJU-1)*7 + JD*(IJS+3))
0679                     ENDIF
0680 C...Initialize junction pull 4-vector.
0681                     DO 510 J=1,5
0682                       PUL(IJS,J)=0D0
0683   510               CONTINUE
0684 C...Initialize weight
0685                     PWT = 0D0
0686                     PWTOLD = 0D0
0687 C...Sum up (weighted) momenta along each string piece
0688                     DO 530 ISP=IA,IB,JD
0689 C...If present parton not last in chain
0690                       IF (ISP.NE.IA.AND.ISP.NE.IB) THEN
0691 C...If last parton was a junction, store present weight
0692                         IF (K(ISP-JD,2).EQ.88) THEN
0693                           PWTOLD = PWT
0694 C...If last parton was a quark, reset to stored weight.
0695                         ELSEIF (K(ISP-JD,2).NE.21) THEN
0696                           PWT = PWTOLD
0697                         ENDIF
0698                       ENDIF
0699 C...Skip next parton if weight already large
0700                       IF (PWT.GT.10D0) GOTO 530
0701 C...Compute momentum in TJUOLD frame:
0702                       TDP=TJUOLD(1)*P(ISP,1)+TJUOLD(2)*P(ISP,2)+TJUOLD(3
0703      &                     )*P(ISP,3)
0704                       BFC=TDP/(1D0+TJUOLD(4))+P(ISP,4)
0705                       DO 520 J=1,3
0706                         TMP=P(ISP,J)+TJUOLD(J)*BFC
0707                         PUL(IJS,J)=PUL(IJS,J)+TMP*EXP(-PWT)
0708   520                 CONTINUE
0709 C...Boosted energy
0710                       TMP=TJUOLD(4)*P(ISP,4)+TDP
0711                       PUL(IJS,4)=PUL(IJS,J)+TMP*EXP(-PWT)
0712 C...Update weight
0713                       PWT=PWT+TMP/PARJ(48)
0714 C...Put |p| rather than m in 5th slot
0715                       PUL(IJS,5)=SQRT(PUL(IJS,1)**2+PUL(IJS,2)**2
0716      &                     +PUL(IJS,3)**2)
0717   530               CONTINUE
0718   540             CONTINUE
0719 C...Compute boost
0720                   IJRFIT=IJRFIT+1
0721                   CALL PYJURF(PUL,T)
0722 C...Combine new boost (T) with old boost (TJUOLD)
0723                   TMP=T(1)*TJUOLD(1)+T(2)*TJUOLD(2)+T(3)*TJUOLD(3)
0724                   DO 550 IX=1,3
0725                     TJUOLD(IX)=T(IX)+TJUOLD(IX)*(TMP/(1D0+TJUOLD(4))+T(4
0726      &                   ))
0727   550             CONTINUE
0728                   TJUOLD(4)=SQRT(1D0+TJUOLD(1)**2+TJUOLD(2)**2+TJUOLD(3)
0729      &                 **2)
0730 C...If last boost small, accept JRF, else iterate.
0731 C...Also prevent possibility of infinite loop.
0732                   IF (ABS((T(4)-1D0)/TJUOLD(4)).GT.0.01D0.AND.
0733      &                 IJRFIT.LT.MSTJ(18))THEN
0734                     GOTO 500
0735                   ELSEIF (IJRFIT.GE.MSTJ(18)) THEN
0736                     CALL PYERRM(1,'(PYPREP:) failed to converge on JRF')
0737                   ENDIF
0738 C...Store final boost, with change of sign since TJJ motion vector.
0739                   DO 560 IX=1,3
0740                     TJJ(IJU,IX)=-TJUOLD(IX)
0741   560             CONTINUE
0742                   TJJ(IJU,4)=SQRT(1D0+TJJ(IJU,1)**2+TJJ(IJU,2)**2
0743      &                 +TJJ(IJU,3)**2)
0744   570           CONTINUE
0745 C...String length measure for (q-qbar)(q-qbar) topology.
0746 C...Note only momenta of nearest partons used (since rest of system
0747 C...identical).
0748                 IF (JJGLUE.EQ.0) THEN
0749                   DELMQQ=4D0*FOUR(IJCP(2)-1,IJCP(4+ISW)+1)*FOUR(IJCP(3)
0750      &                 -1,IJCP(5-ISW)+1)
0751                 ELSE
0752 C...Put jj gluons on selected string (IGS selected randomly above).
0753                   IF (IGS.EQ.1) THEN
0754                     DELMQQ=8D0*FOUR(IJCP(2)-1,IJCP(4)-1)*FOUR(IJCP(3)+1
0755      &                   ,IJCP(4+ISW)+1)*FOUR(IJCP(3)-1,IJCP(5-ISW)+1)
0756                   ELSE
0757                     DELMQQ=8D0*FOUR(IJCP(2)-1,IJCP(4+ISW)+1)
0758      &                   *FOUR(IJCP(3)-1,IJCP(4)-1)*FOUR(IJCP(3)+1
0759      &                   ,IJCP(5-ISW)+1)
0760                   ENDIF
0761                 ENDIF
0762 C...String length measure for q-q-j-j-q-q topology.
0763                 T1G1=0D0
0764                 T2G2=0D0
0765                 T1T2=0D0
0766                 T1P1=0D0
0767                 T1P2=0D0
0768                 T2P3=0D0
0769                 T2P4=0D0
0770                 ISGN=-1
0771 C...Note only momenta of nearest partons used (since rest of system
0772 C...identical).
0773                 DO 580 IX=1,4
0774                   IF (IX.EQ.4) ISGN=1
0775                   T1P1=T1P1+ISGN*TJJ(1,IX)*P(IJCP(2)-1,IX)
0776                   T1P2=T1P2+ISGN*TJJ(1,IX)*P(IJCP(3)-1,IX)
0777                   T2P3=T2P3+ISGN*TJJ(2,IX)*P(IJCP(4)+1,IX)
0778                   T2P4=T2P4+ISGN*TJJ(2,IX)*P(IJCP(5)+1,IX)
0779                   IF (JJGLUE.EQ.0) THEN
0780 C...Junction motion vector dot product gives length when inter-junction
0781 C...gluons absent.
0782                     T1T2=T1T2+ISGN*TJJ(1,IX)*TJJ(2,IX)
0783                   ELSE
0784 C...Junction motion vector dot products with gluon momenta give length
0785 C...when inter-junction gluons present.
0786                     T1G1=T1G1+ISGN*TJJ(1,IX)*P(IJCP(3)+1,IX)
0787                     T2G2=T2G2+ISGN*TJJ(2,IX)*P(IJCP(4)-1,IX)
0788                   ENDIF
0789   580           CONTINUE
0790                 DELMJJ=16D0*T1P1*T1P2*T2P3*T2P4
0791                 IF (JJGLUE.EQ.0) THEN
0792                   DELMJJ=DELMJJ*(T1T2+SQRT(T1T2**2-1))
0793                 ELSE
0794                   DELMJJ=DELMJJ*4D0*T1G1*T2G2
0795                 ENDIF
0796               ENDIF
0797 C...If delmjj > delmqq collapse string system to q-qbar q-qbar
0798 C...(Always the case for MSTJ(19)=2 due to initialization above)
0799               IF (DELMJJ.GT.DELMQQ) THEN
0800 C...Put new system at end of event record
0801                 NCOP=N
0802                 DO 650 IST=1,2
0803                   DO 600 ICOP=IJCP(IST),IJCP(IST+1)-1
0804                     NCOP=NCOP+1
0805                     DO 590 IX=1,5
0806                       P(NCOP,IX)=P(ICOP,IX)
0807                       K(NCOP,IX)=K(ICOP,IX)
0808   590               CONTINUE
0809   600             CONTINUE
0810                   IF (JJGLUE.NE.0.AND.IST.EQ.IGS) THEN
0811 C...Insert inter-junction gluon string piece (reversed)
0812                     NJJGL=0
0813                     DO 620 ICOP=IJCP(4)-1,IJCP(3)+1,-1
0814                       NJJGL=NJJGL+1
0815                       NCOP=NCOP+1
0816                       DO 610 IX=1,5
0817                         P(NCOP,IX)=P(ICOP,IX)
0818                         K(NCOP,IX)=K(ICOP,IX)
0819   610                 CONTINUE
0820   620               CONTINUE
0821                     ENDIF
0822                   IFC=-2*IST+3
0823                   DO 640 ICOP=IJCP(IST+IFC*ISW+3)+1,IJCP(IST+IFC*ISW+4)
0824                     NCOP=NCOP+1
0825                     DO 630 IX=1,5
0826                       P(NCOP,IX)=P(ICOP,IX)
0827                       K(NCOP,IX)=K(ICOP,IX)
0828   630               CONTINUE
0829   640             CONTINUE
0830                   K(NCOP,1)=1
0831   650           CONTINUE
0832 C...Copy system back in right order
0833                 DO 670 ICOP=NBEG,NEND-2
0834                   DO 660 IX=1,5
0835                     P(ICOP,IX)=P(N+ICOP-NBEG+1,IX)
0836                     K(ICOP,IX)=K(N+ICOP-NBEG+1,IX)
0837   660             CONTINUE
0838   670           CONTINUE
0839 C...Shift down rest of event record
0840                 DO 690 ICOP=NEND+1,N
0841                   DO 680 IX=1,5
0842                     P(ICOP-2,IX)=P(ICOP,IX)
0843                     K(ICOP-2,IX)=K(ICOP,IX)
0844   680             CONTINUE
0845   690             CONTINUE
0846 C...Update length of event record.
0847                 N=N-2
0848               ENDIF
0849               MJUN1=0
0850               NBEG=I+1
0851             ENDIF
0852   700     CONTINUE
0853         ENDIF
0854       ENDIF
0855
0856 C...Done if no checks on small-mass systems.
0857       IF(MSTJ(14).LT.0) RETURN
0858       IF(MSTJ(14).EQ.0) GOTO 1140
0859
0860 C...Find lowest-mass colour singlet jet system.
0861       NS=N
0862   710 NSIN=N-NS
0863       PDMIN=1D0+PARJ(32)
0864       IC=0
0865       DO 770 I=MAX(1,IP),N
0866         IF(K(I,1).NE.1.AND.K(I,1).NE.2) THEN
0867         ELSEIF(K(I,1).EQ.2.AND.IC.EQ.0) THEN
0868           NSIN=NSIN+1
0869           IC=I
0870           DO 720 J=1,4
0871             DPS(J)=P(I,J)
0872   720     CONTINUE
0873           MSTJ(93)=1
0874           DPS(5)=PYMASS(K(I,2))
0875         ELSEIF(K(I,1).EQ.2.AND.K(I,2).NE.21) THEN
0876           DO 730 J=1,4
0877             DPS(J)=DPS(J)+P(I,J)
0878   730     CONTINUE
0879           MSTJ(93)=1
0880           DPS(5)=DPS(5)+PYMASS(K(I,2))
0881         ELSEIF(K(I,1).EQ.2) THEN
0882           DO 740 J=1,4
0883             DPS(J)=DPS(J)+P(I,J)
0884   740     CONTINUE
0885         ELSEIF(IC.NE.0.AND.KCHG(PYCOMP(K(I,2)),2).NE.0) THEN
0886           DO 750 J=1,4
0887             DPS(J)=DPS(J)+P(I,J)
0888   750     CONTINUE
0889           MSTJ(93)=1
0890           DPS(5)=DPS(5)+PYMASS(K(I,2))
0891           PD=SQRT(MAX(0D0,DPS(4)**2-DPS(1)**2-DPS(2)**2-DPS(3)**2))-
0892      &    DPS(5)
0893           IF(PD.LT.PDMIN) THEN
0894             PDMIN=PD
0895             DO 760 J=1,5
0896               DPC(J)=DPS(J)
0897   760       CONTINUE
0898             IC1=IC
0899             IC2=I
0900           ENDIF
0901           IC=0
0902         ELSE
0903           NSIN=NSIN+1
0904         ENDIF
0905   770 CONTINUE
0906
0907 C...Done if lowest-mass system above threshold for string frag.
0908       IF(PDMIN.GE.PARJ(32)) GOTO 1140
0909
0910 C...Fill small-mass system as cluster.
0911       NSAV=N
0912       PECM=SQRT(MAX(0D0,DPC(4)**2-DPC(1)**2-DPC(2)**2-DPC(3)**2))
0913       K(N+1,1)=11
0914       K(N+1,2)=91
0915       K(N+1,3)=IC1
0916       P(N+1,1)=DPC(1)
0917       P(N+1,2)=DPC(2)
0918       P(N+1,3)=DPC(3)
0919       P(N+1,4)=DPC(4)
0920       P(N+1,5)=PECM
0921
0922 C...Set up history, assuming cluster -> 2 hadrons.
0923       NBODY=2
0924       K(N+1,4)=N+2
0925       K(N+1,5)=N+3
0926       K(N+2,1)=1
0927       K(N+3,1)=1
0928       IF(MSTU(16).NE.2) THEN
0929         K(N+2,3)=N+1
0930         K(N+3,3)=N+1
0931       ELSE
0932         K(N+2,3)=IC1
0933         K(N+3,3)=IC2
0934       ENDIF
0935       K(N+2,4)=0
0936       K(N+3,4)=0
0937       K(N+2,5)=0
0938       K(N+3,5)=0
0939       V(N+1,5)=0D0
0940       V(N+2,5)=0D0
0941       V(N+3,5)=0D0
0942
0943 C...Find total flavour content - complicated by presence of junctions.
0944       NQ=0
0945       NDIQ=0
0946       DO 780 I=IC1,IC2
0947         IF((K(I,1).EQ.1.OR.K(I,1).EQ.2).AND.K(I,2).NE.21) THEN
0948           NQ=NQ+1
0949           KFQ(NQ)=K(I,2)
0950           IF(IABS(K(I,2)).GT.1000) NDIQ=NDIQ+1
0951         ENDIF
0952   780 CONTINUE
0953
0954 C...If several diquarks, split up one to give even number of flavours.
0955       IF(NQ.EQ.3.AND.NDIQ.GE.2) THEN
0956         I1=3
0957         IF(IABS(KFQ(3)).LT.1000) I1=1
0958         KFQ(4)=ISIGN(MOD(IABS(KFQ(I1))/100,10),KFQ(I1))
0959         KFQ(I1)=KFQ(I1)/1000
0960         NQ=4
0961         NDIQ=NDIQ-1
0962       ENDIF
0963
0964 C...If four quark ends, join two to diquark.
0965       IF(NQ.EQ.4.AND.NDIQ.EQ.0) THEN
0966         I1=1
0967         I2=2
0968         IF(KFQ(I1)*KFQ(I2).LT.0) I2=3
0969         IF(I2.EQ.3.AND.KFQ(I1)*KFQ(I2).LT.0) I2=4
0970         KFLS=2*INT(PYR(0)+3D0*PARJ(4)/(1D0+3D0*PARJ(4)))+1
0971         IF(KFQ(I1).EQ.KFQ(I2)) KFLS=3
0972         KFQ(I1)=ISIGN(1000*MAX(IABS(KFQ(I1)),IABS(KFQ(I2)))+
0973      &  100*MIN(IABS(KFQ(I1)),IABS(KFQ(I2)))+KFLS,KFQ(I1))
0974         KFQ(I2)=KFQ(4)
0975         NQ=3
0976         NDIQ=1
0977       ENDIF
0978
0979 C...If two quark ends, plus quark or diquark, join quarks to diquark.
0980       IF(NQ.EQ.3) THEN
0981         I1=1
0982         I2=2
0983         IF(IABS(KFQ(I1)).GT.1000) I1=3
0984         IF(IABS(KFQ(I2)).GT.1000) I2=3
0985         KFLS=2*INT(PYR(0)+3D0*PARJ(4)/(1D0+3D0*PARJ(4)))+1
0986         IF(KFQ(I1).EQ.KFQ(I2)) KFLS=3
0987         KFQ(I1)=ISIGN(1000*MAX(IABS(KFQ(I1)),IABS(KFQ(I2)))+
0988      &  100*MIN(IABS(KFQ(I1)),IABS(KFQ(I2)))+KFLS,KFQ(I1))
0989         KFQ(I2)=KFQ(3)
0990         NQ=2
0991         NDIQ=NDIQ+1
0992       ENDIF
0993
0994 C...Form two particles from flavours of lowest-mass system, if feasible.
0995       NTRY = 0
0996   790 NTRY = NTRY + 1
0997
0998 C...Open string with two specified endpoint flavours.
0999       IF(NQ.EQ.2) THEN
1000         KC1=PYCOMP(KFQ(1))
1001         KC2=PYCOMP(KFQ(2))
1002         IF(KC1.EQ.0.OR.KC2.EQ.0) GOTO 1140
1003         KQ1=KCHG(KC1,2)*ISIGN(1,KFQ(1))
1004         KQ2=KCHG(KC2,2)*ISIGN(1,KFQ(2))
1005         IF(KQ1+KQ2.NE.0) GOTO 1140
1006 C...Start with qq, if there is one. Only allow for rank 1 popcorn meson
1007   800   K1=KFQ(1)
1008         IF(IABS(KFQ(2)).GT.1000) K1=KFQ(2)
1009         MSTU(125)=0
1010         CALL PYDCYK(K1,0,KFLN,K(N+2,2))
1011         CALL PYDCYK(KFQ(1)+KFQ(2)-K1,-KFLN,KFLDMP,K(N+3,2))
1012         IF(K(N+2,2).EQ.0.OR.K(N+3,2).EQ.0) GOTO 800
1013
1014 C...Open string with four specified flavours.
1015       ELSEIF(NQ.EQ.4) THEN
1016         KC1=PYCOMP(KFQ(1))
1017         KC2=PYCOMP(KFQ(2))
1018         KC3=PYCOMP(KFQ(3))
1019         KC4=PYCOMP(KFQ(4))
1020         IF(KC1.EQ.0.OR.KC2.EQ.0.OR.KC3.EQ.0.OR.KC4.EQ.0) GOTO 1140
1021         KQ1=KCHG(KC1,2)*ISIGN(1,KFQ(1))
1022         KQ2=KCHG(KC2,2)*ISIGN(1,KFQ(2))
1023         KQ3=KCHG(KC3,2)*ISIGN(1,KFQ(3))
1024         KQ4=KCHG(KC4,2)*ISIGN(1,KFQ(4))
1025         IF(KQ1+KQ2+KQ3+KQ4.NE.0) GOTO 1140
1026 C...Combine flavours pairwise to form two hadrons.
1027   810   I1=1
1028         I2=2
1029         IF(KQ1*KQ2.GT.0.OR.(IABS(KFQ(1)).GT.1000.AND.
1030      &  IABS(KFQ(2)).GT.1000)) I2=3
1031         IF(I2.EQ.3.AND.(KQ1*KQ3.GT.0.OR.(IABS(KFQ(1)).GT.1000.AND.
1032      &  IABS(KFQ(3)).GT.1000))) I2=4
1033         I3=3
1034         IF(I2.EQ.3) I3=2
1035         I4=10-I1-I2-I3
1036         CALL PYDCYK(KFQ(I1),KFQ(I2),KFLDMP,K(N+2,2))
1037         CALL PYDCYK(KFQ(I3),KFQ(I4),KFLDMP,K(N+3,2))
1038         IF(K(N+2,2).EQ.0.OR.K(N+3,2).EQ.0) GOTO 810
1039
1040 C...Closed string.
1041       ELSE
1042         IF(IABS(K(IC2,2)).NE.21) GOTO 1140
1043 C...No room for popcorn mesons in closed string -> 2 hadrons.
1044         MSTU(125)=0
1045   820   CALL PYDCYK(1+INT((2D0+PARJ(2))*PYR(0)),0,KFLN,KFDMP)
1046         CALL PYDCYK(KFLN,0,KFLM,K(N+2,2))
1047         CALL PYDCYK(-KFLN,-KFLM,KFLDMP,K(N+3,2))
1048         IF(K(N+2,2).EQ.0.OR.K(N+3,2).EQ.0) GOTO 820
1049       ENDIF
1050       P(N+2,5)=PYMASS(K(N+2,2))
1051       P(N+3,5)=PYMASS(K(N+3,2))
1052
1053 C...If it does not work: try again (a number of times), give up (if no
1054 C...place to shuffle momentum or too many flavours), or form one hadron.
1055       IF(P(N+2,5)+P(N+3,5)+PARJ(64).GE.PECM) THEN
1056         IF(NTRY.LT.MSTJ(17).OR.(NQ.EQ.4.AND.NTRY.LT.5*MSTJ(17))) THEN
1057           GOTO 790
1058         ELSEIF(NSIN.EQ.1.OR.NQ.EQ.4) THEN
1059           GOTO 1140
1060         ELSE
1061           GOTO 890
1062         END IF
1063       END IF
1064
1065 C...Perform two-particle decay of jet system.
1066 C...First step: find reference axis in decaying system rest frame.
1067 C...(Borrow slot N+2 for temporary direction.)
1068       DO 830 J=1,4
1069         P(N+2,J)=P(IC1,J)
1070   830 CONTINUE
1071       DO 850 I=IC1+1,IC2-1
1072         IF((K(I,1).EQ.1.OR.K(I,1).EQ.2).AND.
1073      &  KCHG(PYCOMP(K(I,2)),2).NE.0) THEN
1074           FRAC1=FOUR(IC2,I)/(FOUR(IC1,I)+FOUR(IC2,I))
1075           DO 840 J=1,4
1076             P(N+2,J)=P(N+2,J)+FRAC1*P(I,J)
1077   840     CONTINUE
1078         ENDIF
1079   850 CONTINUE
1080       CALL PYROBO(N+2,N+2,0D0,0D0,-DPC(1)/DPC(4),-DPC(2)/DPC(4),
1081      &-DPC(3)/DPC(4))
1082       THE1=PYANGL(P(N+2,3),SQRT(P(N+2,1)**2+P(N+2,2)**2))
1083       PHI1=PYANGL(P(N+2,1),P(N+2,2))
1084
1085 C...Second step: generate isotropic/anisotropic decay.
1086       PA=SQRT((PECM**2-(P(N+2,5)+P(N+3,5))**2)*(PECM**2-
1087      &(P(N+2,5)-P(N+3,5))**2))/(2D0*PECM)
1088   860 UE(3)=PYR(0)
1089       IF(PARJ(21).LE.0.01D0) UE(3)=1D0
1090       PT2=(1D0-UE(3)**2)*PA**2
1091       IF(MSTJ(16).LE.0) THEN
1092         PREV=0.5D0
1093       ELSE
1094         IF(EXP(-PT2/(2D0*MAX(0.01D0,PARJ(21))**2)).LT.PYR(0)) GOTO 860
1095         PR1=P(N+2,5)**2+PT2
1096         PR2=P(N+3,5)**2+PT2
1097         ALAMBD=SQRT(MAX(0D0,(PECM**2-PR1-PR2)**2-4D0*PR1*PR2))
1098         PREVCF=PARJ(42)
1099         IF(MSTJ(11).EQ.2) PREVCF=PARJ(39)
1100         PREV=1D0/(1D0+EXP(MIN(50D0,PREVCF*ALAMBD*PARJ(40))))
1101       ENDIF
1102       IF(PYR(0).LT.PREV) UE(3)=-UE(3)
1103       PHI=PARU(2)*PYR(0)
1104       UE(1)=SQRT(1D0-UE(3)**2)*COS(PHI)
1105       UE(2)=SQRT(1D0-UE(3)**2)*SIN(PHI)
1106       DO 870 J=1,3
1107         P(N+2,J)=PA*UE(J)
1108         P(N+3,J)=-PA*UE(J)
1109   870 CONTINUE
1110       P(N+2,4)=SQRT(PA**2+P(N+2,5)**2)
1111       P(N+3,4)=SQRT(PA**2+P(N+3,5)**2)
1112
1113 C...Third step: move back to event frame and set production vertex.
1114       CALL PYROBO(N+2,N+3,THE1,PHI1,DPC(1)/DPC(4),DPC(2)/DPC(4),
1115      &DPC(3)/DPC(4))
1116       DO 880 J=1,4
1117         V(N+1,J)=V(IC1,J)
1118         V(N+2,J)=V(IC1,J)
1119         V(N+3,J)=V(IC2,J)
1120   880 CONTINUE
1121       N=N+3
1122       GOTO 1120
1123
1124 C...Else form one particle, if possible.
1125   890 NBODY=1
1126       K(N+1,5)=N+2
1127       DO 900 J=1,4
1128         V(N+1,J)=V(IC1,J)
1129         V(N+2,J)=V(IC1,J)
1130   900 CONTINUE
1131
1132 C...Select hadron flavour from available quark flavours.
1133   910 IF(NQ.EQ.2.AND.IABS(KFQ(1)).GT.100.AND.IABS(KFQ(2)).GT.100) THEN
1134         GOTO 1140
1135       ELSEIF(NQ.EQ.2) THEN
1136         CALL PYKFDI(KFQ(1),KFQ(2),KFLDMP,K(N+2,2))
1137       ELSE
1138         KFLN=1+INT((2D0+PARJ(2))*PYR(0))
1139         CALL PYKFDI(KFLN,-KFLN,KFLDMP,K(N+2,2))
1140       ENDIF
1141       IF(K(N+2,2).EQ.0) GOTO 910
1142       P(N+2,5)=PYMASS(K(N+2,2))
1143
1144 C...Use old algorithm for E/p conservation? (EN)
1145       IF (MSTJ(16).LE.0) GOTO 1080
1146
1147 C...Find the string piece closest to the cluster by a loop
1148 C...over the undecayed partons not in present cluster. (EN)
1149       DGLOMI=1D30
1150       IBEG=0
1151       I0=0
1152       NJUNC=0
1153       DO 940 I1=MAX(1,IP),N-1
1154         IF(K(I1,1).EQ.1) NJUNC=0
1155         IF(K(I1,1).EQ.41) NJUNC=NJUNC+1
1156         IF(K(I1,1).EQ.41) GOTO 940
1157         IF(I1.GE.IC1-1.AND.I1.LE.IC2) THEN
1158           I0=0
1159         ELSEIF(K(I1,1).EQ.2) THEN
1160           IF(I0.EQ.0) I0=I1
1161           I2=I1
1162   920     I2=I2+1
1163           IF(K(I2,1).EQ.41) GOTO 940
1164           IF(K(I2,1).GT.10) GOTO 920
1165           IF(KCHG(PYCOMP(K(I2,2)),2).EQ.0) GOTO 920
1166           IF(K(I1,2).EQ.21.AND.K(I2,2).NE.21.AND.K(I2,1).NE.1.AND.
1167      &    NJUNC.EQ.0) GOTO 940
1168           IF(K(I1,2).NE.21.AND.K(I2,2).EQ.21.AND.NJUNC.NE.0) GOTO 940
1169           IF(K(I1,2).NE.21.AND.K(I2,2).NE.21.AND.(I1.GT.I0.OR.
1170      &    K(I2,1).NE.1)) GOTO 940
1171
1172 C...Define velocity vectors e1, e2, ecl and differences e3, e4.
1173           DO 930 J=1,3
1174             E1(J)=P(I1,J)/P(I1,4)
1175             E2(J)=P(I2,J)/P(I2,4)
1176             ECL(J)=P(N+1,J)/P(N+1,4)
1177             E3(J)=E2(J)-E1(J)
1178             E4(J)=ECL(J)-E1(J)
1179   930     CONTINUE
1180
1181 C...Calculate minimal D=(e4-alpha*e3)**2 for 0<alpha<1.
1182           E3S=E3(1)**2+E3(2)**2+E3(3)**2
1183           E4S=E4(1)**2+E4(2)**2+E4(3)**2
1184           E34=E3(1)*E4(1)+E3(2)*E4(2)+E3(3)*E4(3)
1185           IF(E34.LE.0D0) THEN
1186             DDMIN=E4S
1187           ELSEIF(E34.LT.E3S) THEN
1188             DDMIN=E4S-E34**2/E3S
1189           ELSE
1190             DDMIN=E4S-2D0*E34+E3S
1191           ENDIF
1192
1193 C...Is this the smallest so far?
1194           IF(DDMIN.LT.DGLOMI) THEN
1195             DGLOMI=DDMIN
1196             IBEG=I0
1197             IPCS=I1
1198           ENDIF
1199         ELSEIF(K(I1,1).EQ.1.AND.KCHG(PYCOMP(K(I1,2)),2).NE.0) THEN
1200           I0=0
1201         ENDIF
1202   940 CONTINUE
1203
1204 C... Check if there are any strings to connect to the new gluon. (EN)
1205       IF (IBEG.EQ.0) GOTO 1080
1206
1207 C...Delta_m = m_clus - m_had > 0: emit a 'gluon' (EN)
1208       IF (P(N+1,5).GE.P(N+2,5)) THEN
1209
1210 C...Construct 'gluon' that is needed to put hadron on the mass shell.
1211         FRAC=P(N+2,5)/P(N+1,5)
1212         DO 950 J=1,5
1213           P(N+2,J)=FRAC*P(N+1,J)
1214           PG(J)=(1D0-FRAC)*P(N+1,J)
1215   950   CONTINUE
1216
1217 C... Copy string with new gluon put in.
1218         N=N+2
1219         I=IBEG-1
1220   960   I=I+1
1221         IF(K(I,1).NE.1.AND.K(I,1).NE.2.AND.K(I,1).NE.41) GOTO 960
1222         IF(KCHG(PYCOMP(K(I,2)),2).EQ.0.AND.K(I,1).NE.41) GOTO 960
1223         N=N+1
1224         DO 970 J=1,5
1225           K(N,J)=K(I,J)
1226           P(N,J)=P(I,J)
1227           V(N,J)=V(I,J)
1228   970   CONTINUE
1229         K(I,1)=K(I,1)+10
1230         K(I,4)=N
1231         K(I,5)=N
1232         K(N,3)=I
1233         IF(I.EQ.IPCS) THEN
1234           N=N+1
1235           DO 980 J=1,5
1236             K(N,J)=K(N-1,J)
1237             P(N,J)=PG(J)
1238             V(N,J)=V(N-1,J)
1239   980     CONTINUE
1240           K(N,2)=21
1241           K(N,3)=NSAV+1
1242         ENDIF
1243         IF(K(I,1).EQ.12.OR.K(I,1).EQ.51) GOTO 960
1244         GOTO 1120
1245
1246 C...Delta_m = m_clus - m_had < 0: have to absorb a 'gluon' instead,
1247 C...from string piece endpoints.
1248       ELSE
1249
1250 C...Begin by copying string that should give energy to cluster.
1251         N=N+2
1252         I=IBEG-1
1253   990   I=I+1
1254         IF(K(I,1).NE.1.AND.K(I,1).NE.2.AND.K(I,1).NE.41) GOTO 990
1255         IF(KCHG(PYCOMP(K(I,2)),2).EQ.0.AND.K(I,1).NE.41) GOTO 990
1256         N=N+1
1257         DO 1000 J=1,5
1258           K(N,J)=K(I,J)
1259           P(N,J)=P(I,J)
1260           V(N,J)=V(I,J)
1261  1000   CONTINUE
1262         K(I,1)=K(I,1)+10
1263         K(I,4)=N
1264         K(I,5)=N
1265         K(N,3)=I
1266         IF(I.EQ.IPCS) I1=N
1267         IF(K(I,1).EQ.12.OR.K(I,1).EQ.51) GOTO 990
1268         I2=I1+1
1269
1270 C...Set initial Phad.
1271         DO 1010 J=1,4
1272           P(NSAV+2,J)=P(NSAV+1,J)
1273  1010   CONTINUE
1274
1275 C...Calculate Pg, a part of which will be added to Phad later. (EN)
1276  1020   IF(MSTJ(16).EQ.1) THEN
1277           ALPHA=1D0
1278           BETA=1D0
1279         ELSE
1280           ALPHA=FOUR(NSAV+1,I2)/FOUR(I1,I2)
1281           BETA=FOUR(NSAV+1,I1)/FOUR(I1,I2)
1282         ENDIF
1283         DO 1030 J=1,4
1284           PG(J)=ALPHA*P(I1,J)+BETA*P(I2,J)
1285  1030   CONTINUE
1286         PG(5)=SQRT(MAX(1D-20,PG(4)**2-PG(1)**2-PG(2)**2-PG(3)**2))
1287
1288 C..Solve 2nd order equation, use the best (smallest) solution. (EN)
1289         PMSCOL=P(NSAV+2,4)**2-P(NSAV+2,1)**2-P(NSAV+2,2)**2-
1290      &  P(NSAV+2,3)**2
1291         PCLPG=(P(NSAV+2,4)*PG(4)-P(NSAV+2,1)*PG(1)-
1292      &  P(NSAV+2,2)*PG(2)-P(NSAV+2,3)*PG(3))/PG(5)**2
1293         DELTA=SQRT(PCLPG**2+(P(NSAV+2,5)**2-PMSCOL)/PG(5)**2)-PCLPG
1294
1295 C...If all gluon energy eaten, zero it and take a step back.
1296         ITER=0
1297         IF(DELTA*ALPHA.GT.1D0.AND.I1.GT.NSAV+3.AND.K(I1,2).EQ.21) THEN
1298           ITER=1
1299           DO 1040 J=1,4
1300             P(NSAV+2,J)=P(NSAV+2,J)+P(I1,J)
1301             P(I1,J)=0D0
1302  1040     CONTINUE
1303           P(I1,5)=0D0
1304           K(I1,1)=K(I1,1)+10
1305           I1=I1-1
1306           IF(K(I1,1).EQ.41) ITER=-1
1307         ENDIF
1308         IF(DELTA*BETA.GT.1D0.AND.I2.LT.N.AND.K(I2,2).EQ.21) THEN
1309           ITER=1
1310           DO 1050 J=1,4
1311             P(NSAV+2,J)=P(NSAV+2,J)+P(I2,J)
1312             P(I2,J)=0D0
1313  1050     CONTINUE
1314           P(I2,5)=0D0
1315           K(I2,1)=K(I2,1)+10
1316           I2=I2+1
1317           IF(K(I2,1).EQ.41) ITER=-1
1318         ENDIF
1319         IF(ITER.EQ.1) GOTO 1020
1320
1321 C...If also all endpoint energy eaten, revert to old procedure.
1322         IF((1D0-DELTA*ALPHA)*P(I1,4).LT.P(I1,5).OR.
1323      &  (1D0-DELTA*BETA)*P(I2,4).LT.P(I2,5).OR.ITER.EQ.-1) THEN
1324           DO 1060 I=NSAV+3,N
1325             IM=K(I,3)
1326             K(IM,1)=K(IM,1)-10
1327             K(IM,4)=0
1328             K(IM,5)=0
1329  1060     CONTINUE
1330           N=NSAV
1331           GOTO 1080
1332         ENDIF
1333
1334 C... Construct the collapsed hadron and modified string partons.
1335         DO 1070 J=1,4
1336           P(NSAV+2,J)=P(NSAV+2,J)+DELTA*PG(J)
1337           P(I1,J)=(1D0-DELTA*ALPHA)*P(I1,J)
1338           P(I2,J)=(1D0-DELTA*BETA)*P(I2,J)
1339  1070   CONTINUE
1340           P(I1,5)=(1D0-DELTA*ALPHA)*P(I1,5)
1341           P(I2,5)=(1D0-DELTA*BETA)*P(I2,5)
1342
1343 C...Finished with string collapse in new scheme.
1344         GOTO 1120
1345       ENDIF
1346
1347 C... Use old algorithm; by choice or when in trouble.
1348  1080 CONTINUE
1349 C...Find parton/particle which combines to largest extra mass.
1350       IR=0
1351       HA=0D0
1352       HSM=0D0
1353       DO 1100 MCOMB=1,3
1354         IF(IR.NE.0) GOTO 1100
1355         DO 1090 I=MAX(1,IP),N
1356           IF(K(I,1).LE.0.OR.K(I,1).GT.10.OR.(I.GE.IC1.AND.I.LE.IC2
1357      &    .AND.K(I,1).GE.1.AND.K(I,1).LE.2)) GOTO 1090
1358           IF(MCOMB.EQ.1) KCI=PYCOMP(K(I,2))
1359           IF(MCOMB.EQ.1.AND.KCI.EQ.0) GOTO 1090
1360           IF(MCOMB.EQ.1.AND.KCHG(KCI,2).EQ.0.AND.I.LE.NS) GOTO 1090
1361           IF(MCOMB.EQ.2.AND.IABS(K(I,2)).GT.10.AND.IABS(K(I,2)).LE.100)
1362      &    GOTO 1090
1363           HCR=DPC(4)*P(I,4)-DPC(1)*P(I,1)-DPC(2)*P(I,2)-DPC(3)*P(I,3)
1364           HSR=2D0*HCR+PECM**2-P(N+2,5)**2-2D0*P(N+2,5)*P(I,5)
1365           IF(HSR.GT.HSM) THEN
1366             IR=I
1367             HA=HCR
1368             HSM=HSR
1369           ENDIF
1370  1090   CONTINUE
1371  1100 CONTINUE
1372
1373 C...Shuffle energy and momentum to put new particle on mass shell.
1374       IF(IR.NE.0) THEN
1375         HB=PECM**2+HA
1376         HC=P(N+2,5)**2+HA
1377         HD=P(IR,5)**2+HA
1378         HK2=0.5D0*(HB*SQRT(MAX(0D0,((HB+HC)**2-4D0*(HB+HD)*P(N+2,5)**2)/
1379      &  (HA**2-(PECM*P(IR,5))**2)))-(HB+HC))/(HB+HD)
1380         HK1=(0.5D0*(P(N+2,5)**2-PECM**2)+HD*HK2)/HB
1381         DO 1110 J=1,4
1382           P(N+2,J)=(1D0+HK1)*DPC(J)-HK2*P(IR,J)
1383           P(IR,J)=(1D0+HK2)*P(IR,J)-HK1*DPC(J)
1384  1110   CONTINUE
1385         N=N+2
1386       ELSE
1387         CALL PYERRM(3,'(PYPREP:) no match for collapsing cluster')
1388         RETURN
1389       ENDIF
1390
1391 C...Mark collapsed system and store daughter pointers. Iterate.
1392  1120 DO 1130 I=IC1,IC2
1393         IF((K(I,1).EQ.1.OR.K(I,1).EQ.2).AND.
1394      &  KCHG(PYCOMP(K(I,2)),2).NE.0) THEN
1395           K(I,1)=K(I,1)+10
1396           IF(MSTU(16).NE.2) THEN
1397             K(I,4)=NSAV+1
1398             K(I,5)=NSAV+1
1399           ELSE
1400             K(I,4)=NSAV+2
1401             K(I,5)=NSAV+1+NBODY
1402           ENDIF
1403         ENDIF
1404         IF(K(I,1).EQ.41) K(I,1)=K(I,1)+10
1405  1130 CONTINUE
1406       IF(N.LT.MSTU(4)-MSTU(32)-5) GOTO 710
1407
1408 C...Check flavours and invariant masses in parton systems.
1409  1140 NP=0
1410       KFN=0
1411       KQS=0
1412       NJU=0
1413       DO 1150 J=1,5
1414         DPS(J)=0D0
1415  1150 CONTINUE
1416       DO 1180 I=MAX(1,IP),N
1417         IF(K(I,1).EQ.41) NJU=NJU+1
1418         IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 1180
1419         KC=PYCOMP(K(I,2))
1420         IF(KC.EQ.0) GOTO 1180
1421         KQ=KCHG(KC,2)*ISIGN(1,K(I,2))
1422         IF(KQ.EQ.0) GOTO 1180
1423         NP=NP+1
1424         IF(KQ.NE.2) THEN
1425           KFN=KFN+1
1426           KQS=KQS+KQ
1427           MSTJ(93)=1
1428           DPS(5)=DPS(5)+PYMASS(K(I,2))
1429         ENDIF
1430         DO 1160 J=1,4
1431           DPS(J)=DPS(J)+P(I,J)
1432  1160   CONTINUE
1433         IF(K(I,1).EQ.1) THEN
1434           NFERR=0
1435           IF(NJU.EQ.0.AND.NP.NE.1) THEN
1436             IF(KFN.EQ.1.OR.KFN.GE.3.OR.KQS.NE.0) NFERR=1
1437           ELSEIF(NJU.EQ.1) THEN
1438             IF(KFN.NE.3.OR.IABS(KQS).NE.3) NFERR=1
1439           ELSEIF(NJU.EQ.2) THEN
1440             IF(KFN.NE.4.OR.KQS.NE.0) NFERR=1
1441           ELSEIF(NJU.GE.3) THEN
1442             NFERR=1
1443           ENDIF
1444           IF(NFERR.EQ.1) THEN
1445             CALL PYERRM(2,'(PYPREP:) unphysical flavour combination')
1446             MINT(51)=1
1447             RETURN
1448           ENDIF
1449           IF(NP.NE.1.AND.DPS(4)**2-DPS(1)**2-DPS(2)**2-DPS(3)**2.LT.
1450      &    (0.9D0*PARJ(32)+DPS(5))**2) CALL PYERRM(3,
1451      &    '(PYPREP:) too small mass in jet system')
1452           NP=0
1453           KFN=0
1454           KQS=0
1455           NJU=0
1456           DO 1170 J=1,5
1457             DPS(J)=0D0
1458  1170     CONTINUE
1459         ENDIF
1460  1180 CONTINUE
1461
1462       RETURN
1463       END