pythiaeRHIC/pythia/pyindf.f

0001
0002 C*********************************************************************
0003
0004 C...PYINDF
0005 C...Handles the fragmentation of a jet system (or a single
0006 C...jet) according to independent fragmentation models.
0007
0008       SUBROUTINE PYINDF(IP)
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...Commonblocks.
0015       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
0016       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0017       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0018       SAVE /PYJETS/,/PYDAT1/,/PYDAT2/
0019 C...Local arrays.
0020       DIMENSION DPS(5),PSI(4),NFI(3),NFL(3),IFET(3),KFLF(3),
0021      &KFLO(2),PXO(2),PYO(2),WO(2)
0022
0023 C.. MOPS error message
0024       IF(MSTJ(12).GT.3) CALL PYERRM(9,'(PYINDF:) MSTJ(12)>3 options'//
0025      &' are not treated as expected in independent fragmentation')
0026
0027 C...Reset counters. Identify parton system and take copy. Check flavour.
0028       NSAV=N
0029       MSTU90=MSTU(90)
0030       NJET=0
0031       KQSUM=0
0032       DO 100 J=1,5
0033         DPS(J)=0D0
0034   100 CONTINUE
0035       I=IP-1
0036   110 I=I+1
0037       IF(I.GT.MIN(N,MSTU(4)-MSTU(32))) THEN
0038         CALL PYERRM(12,'(PYINDF:) failed to reconstruct jet system')
0039         IF(MSTU(21).GE.1) RETURN
0040       ENDIF
0041       IF(K(I,1).NE.1.AND.K(I,1).NE.2) GOTO 110
0042       KC=PYCOMP(K(I,2))
0043       IF(KC.EQ.0) GOTO 110
0044       KQ=KCHG(KC,2)*ISIGN(1,K(I,2))
0045       IF(KQ.EQ.0) GOTO 110
0046       NJET=NJET+1
0047       IF(KQ.NE.2) KQSUM=KQSUM+KQ
0048       DO 120 J=1,5
0049         K(NSAV+NJET,J)=K(I,J)
0050         P(NSAV+NJET,J)=P(I,J)
0051         DPS(J)=DPS(J)+P(I,J)
0052   120 CONTINUE
0053       K(NSAV+NJET,3)=I
0054       IF(K(I,1).EQ.2.OR.(MSTJ(3).LE.5.AND.N.GT.I.AND.
0055      &K(I+1,1).EQ.2)) GOTO 110
0056       IF(NJET.NE.1.AND.KQSUM.NE.0) THEN
0057         CALL PYERRM(12,'(PYINDF:) unphysical flavour combination')
0058         IF(MSTU(21).GE.1) RETURN
0059       ENDIF
0060
0061 C...Boost copied system to CM frame. Find CM energy and sum flavours.
0062       IF(NJET.NE.1) THEN
0063         MSTU(33)=1
0064         CALL PYROBO(NSAV+1,NSAV+NJET,0D0,0D0,-DPS(1)/DPS(4),
0065      &  -DPS(2)/DPS(4),-DPS(3)/DPS(4))
0066       ENDIF
0067       PECM=0D0
0068       DO 130 J=1,3
0069         NFI(J)=0
0070   130 CONTINUE
0071       DO 140 I=NSAV+1,NSAV+NJET
0072         PECM=PECM+P(I,4)
0073         KFA=IABS(K(I,2))
0074         IF(KFA.LE.3) THEN
0075           NFI(KFA)=NFI(KFA)+ISIGN(1,K(I,2))
0076         ELSEIF(KFA.GT.1000) THEN
0077           KFLA=MOD(KFA/1000,10)
0078           KFLB=MOD(KFA/100,10)
0079           IF(KFLA.LE.3) NFI(KFLA)=NFI(KFLA)+ISIGN(1,K(I,2))
0080           IF(KFLB.LE.3) NFI(KFLB)=NFI(KFLB)+ISIGN(1,K(I,2))
0081         ENDIF
0082   140 CONTINUE
0083
0084 C...Loop over attempts made. Reset counters.
0085       NTRY=0
0086   150 NTRY=NTRY+1
0087       IF(NTRY.GT.200) THEN
0088         CALL PYERRM(14,'(PYINDF:) caught in infinite loop')
0089         IF(MSTU(21).GE.1) RETURN
0090       ENDIF
0091       N=NSAV+NJET
0092       MSTU(90)=MSTU90
0093       DO 160 J=1,3
0094         NFL(J)=NFI(J)
0095         IFET(J)=0
0096         KFLF(J)=0
0097   160 CONTINUE
0098
0099 C...Loop over jets to be fragmented.
0100       DO 230 IP1=NSAV+1,NSAV+NJET
0101         MSTJ(91)=0
0102         NSAV1=N
0103         MSTU91=MSTU(90)
0104
0105 C...Initial flavour and momentum values. Jet along +z axis.
0106         KFLH=IABS(K(IP1,2))
0107         IF(KFLH.GT.10) KFLH=MOD(KFLH/1000,10)
0108         KFLO(2)=0
0109         WF=P(IP1,4)+SQRT(P(IP1,1)**2+P(IP1,2)**2+P(IP1,3)**2)
0110
0111 C...Initial values for quark or diquark jet.
0112   170   IF(IABS(K(IP1,2)).NE.21) THEN
0113           NSTR=1
0114           KFLO(1)=K(IP1,2)
0115           CALL PYPTDI(0,PXO(1),PYO(1))
0116           WO(1)=WF
0117
0118 C...Initial values for gluon treated like random quark jet.
0119         ELSEIF(MSTJ(2).LE.2) THEN
0120           NSTR=1
0121           IF(MSTJ(2).EQ.2) MSTJ(91)=1
0122           KFLO(1)=INT(1D0+(2D0+PARJ(2))*PYR(0))*(-1)**INT(PYR(0)+0.5D0)
0123           CALL PYPTDI(0,PXO(1),PYO(1))
0124           WO(1)=WF
0125
0126 C...Initial values for gluon treated like quark-antiquark jet pair,
0127 C...sharing energy according to Altarelli-Parisi splitting function.
0128         ELSE
0129           NSTR=2
0130           IF(MSTJ(2).EQ.4) MSTJ(91)=1
0131           KFLO(1)=INT(1D0+(2D0+PARJ(2))*PYR(0))*(-1)**INT(PYR(0)+0.5D0)
0132           KFLO(2)=-KFLO(1)
0133           CALL PYPTDI(0,PXO(1),PYO(1))
0134           PXO(2)=-PXO(1)
0135           PYO(2)=-PYO(1)
0136           WO(1)=WF*PYR(0)**(1D0/3D0)
0137           WO(2)=WF-WO(1)
0138         ENDIF
0139
0140 C...Initial values for rank, flavour, pT and W+.
0141         DO 220 ISTR=1,NSTR
0142   180     I=N
0143           MSTU(90)=MSTU91
0144           IRANK=0
0145           KFL1=KFLO(ISTR)
0146           PX1=PXO(ISTR)
0147           PY1=PYO(ISTR)
0148           W=WO(ISTR)
0149
0150 C...New hadron. Generate flavour and hadron species.
0151   190     I=I+1
0152           IF(I.GE.MSTU(4)-MSTU(32)-NJET-5) THEN
0153             CALL PYERRM(11,'(PYINDF:) no more memory left in PYJETS')
0154             IF(MSTU(21).GE.1) RETURN
0155           ENDIF
0156           IRANK=IRANK+1
0157           K(I,1)=1
0158           K(I,3)=IP1
0159           K(I,4)=0
0160           K(I,5)=0
0161   200     CALL PYKFDI(KFL1,0,KFL2,K(I,2))
0162           IF(K(I,2).EQ.0) GOTO 180
0163           IF(IRANK.EQ.1.AND.IABS(KFL1).LE.10.AND.IABS(KFL2).GT.10) THEN
0164             IF(PYR(0).GT.PARJ(19)) GOTO 200
0165           ENDIF
0166
0167 C...Find hadron mass. Generate four-momentum.
0168           P(I,5)=PYMASS(K(I,2))
0169           CALL PYPTDI(KFL1,PX2,PY2)
0170           P(I,1)=PX1+PX2
0171           P(I,2)=PY1+PY2
0172           PR=P(I,5)**2+P(I,1)**2+P(I,2)**2
0173           CALL PYZDIS(KFL1,KFL2,PR,Z)
0174           MZSAV=0
0175           IF(IABS(KFL1).GE.4.AND.IABS(KFL1).LE.8.AND.MSTU(90).LT.8) THEN
0176             MZSAV=1
0177             MSTU(90)=MSTU(90)+1
0178             MSTU(90+MSTU(90))=I
0179             PARU(90+MSTU(90))=Z
0180           ENDIF
0181           P(I,3)=0.5D0*(Z*W-PR/MAX(1D-4,Z*W))
0182           P(I,4)=0.5D0*(Z*W+PR/MAX(1D-4,Z*W))
0183           IF(MSTJ(3).GE.1.AND.IRANK.EQ.1.AND.KFLH.GE.4.AND.
0184      &    P(I,3).LE.0.001D0) THEN
0185             IF(W.GE.P(I,5)+0.5D0*PARJ(32)) GOTO 180
0186             P(I,3)=0.0001D0
0187             P(I,4)=SQRT(PR)
0188             Z=P(I,4)/W
0189           ENDIF
0190
0191 C...Remaining flavour and momentum.
0192           KFL1=-KFL2
0193           PX1=-PX2
0194           PY1=-PY2
0195           W=(1D0-Z)*W
0196           DO 210 J=1,5
0197             V(I,J)=0D0
0198   210     CONTINUE
0199
0200 C...Check if pL acceptable. Go back for new hadron if enough energy.
0201           IF(MSTJ(3).GE.0.AND.P(I,3).LT.0D0) THEN
0202             I=I-1
0203             IF(MZSAV.EQ.1) MSTU(90)=MSTU(90)-1
0204           ENDIF
0205           IF(W.GT.PARJ(31)) GOTO 190
0206           N=I
0207   220   CONTINUE
0208         IF(MOD(MSTJ(3),5).EQ.4.AND.N.EQ.NSAV1) WF=WF+0.1D0*PARJ(32)
0209         IF(MOD(MSTJ(3),5).EQ.4.AND.N.EQ.NSAV1) GOTO 170
0210
0211 C...Rotate jet to new direction.
0212         THE=PYANGL(P(IP1,3),SQRT(P(IP1,1)**2+P(IP1,2)**2))
0213         PHI=PYANGL(P(IP1,1),P(IP1,2))
0214         MSTU(33)=1
0215         CALL PYROBO(NSAV1+1,N,THE,PHI,0D0,0D0,0D0)
0216         K(K(IP1,3),4)=NSAV1+1
0217         K(K(IP1,3),5)=N
0218
0219 C...End of jet generation loop. Skip conservation in some cases.
0220   230 CONTINUE
0221       IF(NJET.EQ.1.OR.MSTJ(3).LE.0) GOTO 490
0222       IF(MOD(MSTJ(3),5).NE.0.AND.N-NSAV-NJET.LT.2) GOTO 150
0223
0224 C...Subtract off produced hadron flavours, finished if zero.
0225       DO 240 I=NSAV+NJET+1,N
0226         KFA=IABS(K(I,2))
0227         KFLA=MOD(KFA/1000,10)
0228         KFLB=MOD(KFA/100,10)
0229         KFLC=MOD(KFA/10,10)
0230         IF(KFLA.EQ.0) THEN
0231           IF(KFLB.LE.3) NFL(KFLB)=NFL(KFLB)-ISIGN(1,K(I,2))*(-1)**KFLB
0232           IF(KFLC.LE.3) NFL(KFLC)=NFL(KFLC)+ISIGN(1,K(I,2))*(-1)**KFLB
0233         ELSE
0234           IF(KFLA.LE.3) NFL(KFLA)=NFL(KFLA)-ISIGN(1,K(I,2))
0235           IF(KFLB.LE.3) NFL(KFLB)=NFL(KFLB)-ISIGN(1,K(I,2))
0236           IF(KFLC.LE.3) NFL(KFLC)=NFL(KFLC)-ISIGN(1,K(I,2))
0237         ENDIF
0238   240 CONTINUE
0239       NREQ=(IABS(NFL(1))+IABS(NFL(2))+IABS(NFL(3))-IABS(NFL(1)+
0240      &NFL(2)+NFL(3)))/2+IABS(NFL(1)+NFL(2)+NFL(3))/3
0241       IF(NREQ.EQ.0) GOTO 320
0242
0243 C...Take away flavour of low-momentum particles until enough freedom.
0244       NREM=0
0245   250 IREM=0
0246       P2MIN=PECM**2
0247       DO 260 I=NSAV+NJET+1,N
0248         P2=P(I,1)**2+P(I,2)**2+P(I,3)**2
0249         IF(K(I,1).EQ.1.AND.P2.LT.P2MIN) IREM=I
0250         IF(K(I,1).EQ.1.AND.P2.LT.P2MIN) P2MIN=P2
0251   260 CONTINUE
0252       IF(IREM.EQ.0) GOTO 150
0253       K(IREM,1)=7
0254       KFA=IABS(K(IREM,2))
0255       KFLA=MOD(KFA/1000,10)
0256       KFLB=MOD(KFA/100,10)
0257       KFLC=MOD(KFA/10,10)
0258       IF(KFLA.GE.4.OR.KFLB.GE.4) K(IREM,1)=8
0259       IF(K(IREM,1).EQ.8) GOTO 250
0260       IF(KFLA.EQ.0) THEN
0261         ISGN=ISIGN(1,K(IREM,2))*(-1)**KFLB
0262         IF(KFLB.LE.3) NFL(KFLB)=NFL(KFLB)+ISGN
0263         IF(KFLC.LE.3) NFL(KFLC)=NFL(KFLC)-ISGN
0264       ELSE
0265         IF(KFLA.LE.3) NFL(KFLA)=NFL(KFLA)+ISIGN(1,K(IREM,2))
0266         IF(KFLB.LE.3) NFL(KFLB)=NFL(KFLB)+ISIGN(1,K(IREM,2))
0267         IF(KFLC.LE.3) NFL(KFLC)=NFL(KFLC)+ISIGN(1,K(IREM,2))
0268       ENDIF
0269       NREM=NREM+1
0270       NREQ=(IABS(NFL(1))+IABS(NFL(2))+IABS(NFL(3))-IABS(NFL(1)+
0271      &NFL(2)+NFL(3)))/2+IABS(NFL(1)+NFL(2)+NFL(3))/3
0272       IF(NREQ.GT.NREM) GOTO 250
0273       DO 270 I=NSAV+NJET+1,N
0274         IF(K(I,1).EQ.8) K(I,1)=1
0275   270 CONTINUE
0276
0277 C...Find combination of existing and new flavours for hadron.
0278   280 NFET=2
0279       IF(NFL(1)+NFL(2)+NFL(3).NE.0) NFET=3
0280       IF(NREQ.LT.NREM) NFET=1
0281       IF(IABS(NFL(1))+IABS(NFL(2))+IABS(NFL(3)).EQ.0) NFET=0
0282       DO 290 J=1,NFET
0283         IFET(J)=1+(IABS(NFL(1))+IABS(NFL(2))+IABS(NFL(3)))*PYR(0)
0284         KFLF(J)=ISIGN(1,NFL(1))
0285         IF(IFET(J).GT.IABS(NFL(1))) KFLF(J)=ISIGN(2,NFL(2))
0286         IF(IFET(J).GT.IABS(NFL(1))+IABS(NFL(2))) KFLF(J)=ISIGN(3,NFL(3))
0287   290 CONTINUE
0288       IF(NFET.EQ.2.AND.(IFET(1).EQ.IFET(2).OR.KFLF(1)*KFLF(2).GT.0))
0289      &GOTO 280
0290       IF(NFET.EQ.3.AND.(IFET(1).EQ.IFET(2).OR.IFET(1).EQ.IFET(3).OR.
0291      &IFET(2).EQ.IFET(3).OR.KFLF(1)*KFLF(2).LT.0.OR.KFLF(1)*KFLF(3)
0292      &.LT.0.OR.KFLF(1)*(NFL(1)+NFL(2)+NFL(3)).LT.0)) GOTO 280
0293       IF(NFET.EQ.0) KFLF(1)=1+INT((2D0+PARJ(2))*PYR(0))
0294       IF(NFET.EQ.0) KFLF(2)=-KFLF(1)
0295       IF(NFET.EQ.1) KFLF(2)=ISIGN(1+INT((2D0+PARJ(2))*PYR(0)),-KFLF(1))
0296       IF(NFET.LE.2) KFLF(3)=0
0297       IF(KFLF(3).NE.0) THEN
0298         KFLFC=ISIGN(1000*MAX(IABS(KFLF(1)),IABS(KFLF(3)))+
0299      &  100*MIN(IABS(KFLF(1)),IABS(KFLF(3)))+1,KFLF(1))
0300         IF(KFLF(1).EQ.KFLF(3).OR.(1D0+3D0*PARJ(4))*PYR(0).GT.1D0)
0301      &  KFLFC=KFLFC+ISIGN(2,KFLFC)
0302       ELSE
0303         KFLFC=KFLF(1)
0304       ENDIF
0305       CALL PYKFDI(KFLFC,KFLF(2),KFLDMP,KF)
0306       IF(KF.EQ.0) GOTO 280
0307       DO 300 J=1,MAX(2,NFET)
0308         NFL(IABS(KFLF(J)))=NFL(IABS(KFLF(J)))-ISIGN(1,KFLF(J))
0309   300 CONTINUE
0310
0311 C...Store hadron at random among free positions.
0312       NPOS=MIN(1+INT(PYR(0)*NREM),NREM)
0313       DO 310 I=NSAV+NJET+1,N
0314         IF(K(I,1).EQ.7) NPOS=NPOS-1
0315         IF(K(I,1).EQ.1.OR.NPOS.NE.0) GOTO 310
0316         K(I,1)=1
0317         K(I,2)=KF
0318         P(I,5)=PYMASS(K(I,2))
0319         P(I,4)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2+P(I,5)**2)
0320   310 CONTINUE
0321       NREM=NREM-1
0322       NREQ=(IABS(NFL(1))+IABS(NFL(2))+IABS(NFL(3))-IABS(NFL(1)+
0323      &NFL(2)+NFL(3)))/2+IABS(NFL(1)+NFL(2)+NFL(3))/3
0324       IF(NREM.GT.0) GOTO 280
0325
0326 C...Compensate for missing momentum in global scheme (3 options).
0327   320 IF(MOD(MSTJ(3),5).NE.0.AND.MOD(MSTJ(3),5).NE.4) THEN
0328         DO 340 J=1,3
0329           PSI(J)=0D0
0330           DO 330 I=NSAV+NJET+1,N
0331             PSI(J)=PSI(J)+P(I,J)
0332   330     CONTINUE
0333   340   CONTINUE
0334         PSI(4)=PSI(1)**2+PSI(2)**2+PSI(3)**2
0335         PWS=0D0
0336         DO 350 I=NSAV+NJET+1,N
0337           IF(MOD(MSTJ(3),5).EQ.1) PWS=PWS+P(I,4)
0338           IF(MOD(MSTJ(3),5).EQ.2) PWS=PWS+SQRT(P(I,5)**2+(PSI(1)*P(I,1)+
0339      &    PSI(2)*P(I,2)+PSI(3)*P(I,3))**2/PSI(4))
0340           IF(MOD(MSTJ(3),5).EQ.3) PWS=PWS+1D0
0341   350   CONTINUE
0342         DO 370 I=NSAV+NJET+1,N
0343           IF(MOD(MSTJ(3),5).EQ.1) PW=P(I,4)
0344           IF(MOD(MSTJ(3),5).EQ.2) PW=SQRT(P(I,5)**2+(PSI(1)*P(I,1)+
0345      &    PSI(2)*P(I,2)+PSI(3)*P(I,3))**2/PSI(4))
0346           IF(MOD(MSTJ(3),5).EQ.3) PW=1D0
0347           DO 360 J=1,3
0348             P(I,J)=P(I,J)-PSI(J)*PW/PWS
0349   360     CONTINUE
0350           P(I,4)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2+P(I,5)**2)
0351   370   CONTINUE
0352
0353 C...Compensate for missing momentum withing each jet separately.
0354       ELSEIF(MOD(MSTJ(3),5).EQ.4) THEN
0355         DO 390 I=N+1,N+NJET
0356           K(I,1)=0
0357           DO 380 J=1,5
0358             P(I,J)=0D0
0359   380     CONTINUE
0360   390   CONTINUE
0361         DO 410 I=NSAV+NJET+1,N
0362           IR1=K(I,3)
0363           IR2=N+IR1-NSAV
0364           K(IR2,1)=K(IR2,1)+1
0365           PLS=(P(I,1)*P(IR1,1)+P(I,2)*P(IR1,2)+P(I,3)*P(IR1,3))/
0366      &    (P(IR1,1)**2+P(IR1,2)**2+P(IR1,3)**2)
0367           DO 400 J=1,3
0368             P(IR2,J)=P(IR2,J)+P(I,J)-PLS*P(IR1,J)
0369   400     CONTINUE
0370           P(IR2,4)=P(IR2,4)+P(I,4)
0371           P(IR2,5)=P(IR2,5)+PLS
0372   410   CONTINUE
0373         PSS=0D0
0374         DO 420 I=N+1,N+NJET
0375           IF(K(I,1).NE.0) PSS=PSS+P(I,4)/(PECM*(0.8D0*P(I,5)+0.2D0))
0376   420   CONTINUE
0377         DO 440 I=NSAV+NJET+1,N
0378           IR1=K(I,3)
0379           IR2=N+IR1-NSAV
0380           PLS=(P(I,1)*P(IR1,1)+P(I,2)*P(IR1,2)+P(I,3)*P(IR1,3))/
0381      &    (P(IR1,1)**2+P(IR1,2)**2+P(IR1,3)**2)
0382           DO 430 J=1,3
0383             P(I,J)=P(I,J)-P(IR2,J)/K(IR2,1)+(1D0/(P(IR2,5)*PSS)-1D0)*
0384      &      PLS*P(IR1,J)
0385   430     CONTINUE
0386           P(I,4)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2+P(I,5)**2)
0387   440   CONTINUE
0388       ENDIF
0389
0390 C...Scale momenta for energy conservation.
0391       IF(MOD(MSTJ(3),5).NE.0) THEN
0392         PMS=0D0
0393         PES=0D0
0394         PQS=0D0
0395         DO 450 I=NSAV+NJET+1,N
0396           PMS=PMS+P(I,5)
0397           PES=PES+P(I,4)
0398           PQS=PQS+P(I,5)**2/P(I,4)
0399   450   CONTINUE
0400         IF(PMS.GE.PECM) GOTO 150
0401         NECO=0
0402   460   NECO=NECO+1
0403         PFAC=(PECM-PQS)/(PES-PQS)
0404         PES=0D0
0405         PQS=0D0
0406         DO 480 I=NSAV+NJET+1,N
0407           DO 470 J=1,3
0408             P(I,J)=PFAC*P(I,J)
0409   470     CONTINUE
0410           P(I,4)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2+P(I,5)**2)
0411           PES=PES+P(I,4)
0412           PQS=PQS+P(I,5)**2/P(I,4)
0413   480   CONTINUE
0414         IF(NECO.LT.10.AND.ABS(PECM-PES).GT.2D-6*PECM) GOTO 460
0415       ENDIF
0416
0417 C...Origin of produced particles and parton daughter pointers.
0418   490 DO 500 I=NSAV+NJET+1,N
0419         IF(MSTU(16).NE.2) K(I,3)=NSAV+1
0420         IF(MSTU(16).EQ.2) K(I,3)=K(K(I,3),3)
0421   500 CONTINUE
0422       DO 510 I=NSAV+1,NSAV+NJET
0423         I1=K(I,3)
0424         K(I1,1)=K(I1,1)+10
0425         IF(MSTU(16).NE.2) THEN
0426           K(I1,4)=NSAV+1
0427           K(I1,5)=NSAV+1
0428         ELSE
0429           K(I1,4)=K(I1,4)-NJET+1
0430           K(I1,5)=K(I1,5)-NJET+1
0431           IF(K(I1,5).LT.K(I1,4)) THEN
0432             K(I1,4)=0
0433             K(I1,5)=0
0434           ENDIF
0435         ENDIF
0436   510 CONTINUE
0437
0438 C...Document independent fragmentation system. Remove copy of jets.
0439       NSAV=NSAV+1
0440       K(NSAV,1)=11
0441       K(NSAV,2)=93
0442       K(NSAV,3)=IP
0443       K(NSAV,4)=NSAV+1
0444       K(NSAV,5)=N-NJET+1
0445       DO 520 J=1,4
0446         P(NSAV,J)=DPS(J)
0447         V(NSAV,J)=V(IP,J)
0448   520 CONTINUE
0449       P(NSAV,5)=SQRT(MAX(0D0,DPS(4)**2-DPS(1)**2-DPS(2)**2-DPS(3)**2))
0450       V(NSAV,5)=0D0
0451       DO 540 I=NSAV+NJET,N
0452         DO 530 J=1,5
0453           K(I-NJET+1,J)=K(I,J)
0454           P(I-NJET+1,J)=P(I,J)
0455           V(I-NJET+1,J)=V(I,J)
0456   530   CONTINUE
0457   540 CONTINUE
0458       N=N-NJET+1
0459       DO 550 IZ=MSTU90+1,MSTU(90)
0460         MSTU(90+IZ)=MSTU(90+IZ)-NJET+1
0461   550 CONTINUE
0462
0463 C...Boost back particle system. Set production vertices.
0464       IF(NJET.NE.1) CALL PYROBO(NSAV+1,N,0D0,0D0,DPS(1)/DPS(4),
0465      &DPS(2)/DPS(4),DPS(3)/DPS(4))
0466       DO 570 I=NSAV+1,N
0467         DO 560 J=1,4
0468           V(I,J)=V(IP,J)
0469   560   CONTINUE
0470   570 CONTINUE
0471
0472       RETURN
0473       END