PYTHIA/jetset/luindf.F

   1
   2 C*********************************************************************
   3
   4       SUBROUTINE LUINDF(IP)
   5
   6 C...Purpose: to handle the fragmentation of a jet system (or a single
   7 C...jet) according to independent fragmentation models.
   8       IMPLICIT DOUBLE PRECISION(D)
   9       COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5)
  10       COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
  11       COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
  12       SAVE /LUJETS/,/LUDAT1/,/LUDAT2/
  13       DIMENSION DPS(5),PSI(4),NFI(3),NFL(3),IFET(3),KFLF(3),
  14      &KFLO(2),PXO(2),PYO(2),WO(2)
  15
  16 C...Reset counters. Identify parton system and take copy. Check flavour.
  17       NSAV=N
  18       MSTU90=MSTU(90)
  19       NJET=0
  20       KQSUM=0
  21       DO 100 J=1,5
  22       DPS(J)=0.
  23   100 CONTINUE
  24       I=IP-1
  25   110 I=I+1
  26       IF(I.GT.MIN(N,MSTU(4)-MSTU(32))) THEN
  27         CALL LUERRM(12,'(LUINDF:) failed to reconstruct jet system')
  28         IF(MSTU(21).GE.1) RETURN
  29       ENDIF
  30       IF(K(I,1).NE.1.AND.K(I,1).NE.2) GOTO 110
  31       KC=LUCOMP(K(I,2))
  32       IF(KC.EQ.0) GOTO 110
  33       KQ=KCHG(KC,2)*ISIGN(1,K(I,2))
  34       IF(KQ.EQ.0) GOTO 110
  35       NJET=NJET+1
  36       IF(KQ.NE.2) KQSUM=KQSUM+KQ
  37       DO 120 J=1,5
  38       K(NSAV+NJET,J)=K(I,J)
  39       P(NSAV+NJET,J)=P(I,J)
  40       DPS(J)=DPS(J)+P(I,J)
  41   120 CONTINUE
  42       K(NSAV+NJET,3)=I
  43       IF(K(I,1).EQ.2.OR.(MSTJ(3).LE.5.AND.N.GT.I.AND.
  44      &K(I+1,1).EQ.2)) GOTO 110
  45       IF(NJET.NE.1.AND.KQSUM.NE.0) THEN
  46         CALL LUERRM(12,'(LUINDF:) unphysical flavour combination')
  47         IF(MSTU(21).GE.1) RETURN
  48       ENDIF
  49
  50 C...Boost copied system to CM frame. Find CM energy and sum flavours.
  51       IF(NJET.NE.1) THEN
  52         MSTU(33)=1
  53         CALL LUDBRB(NSAV+1,NSAV+NJET,0.,0.,-DPS(1)/DPS(4),
  54      &  -DPS(2)/DPS(4),-DPS(3)/DPS(4))
  55       ENDIF
  56       PECM=0.
  57       DO 130 J=1,3
  58       NFI(J)=0
  59   130 CONTINUE
  60       DO 140 I=NSAV+1,NSAV+NJET
  61       PECM=PECM+P(I,4)
  62       KFA=IABS(K(I,2))
  63       IF(KFA.LE.3) THEN
  64         NFI(KFA)=NFI(KFA)+ISIGN(1,K(I,2))
  65       ELSEIF(KFA.GT.1000) THEN
  66         KFLA=MOD(KFA/1000,10)
  67         KFLB=MOD(KFA/100,10)
  68         IF(KFLA.LE.3) NFI(KFLA)=NFI(KFLA)+ISIGN(1,K(I,2))
  69         IF(KFLB.LE.3) NFI(KFLB)=NFI(KFLB)+ISIGN(1,K(I,2))
  70       ENDIF
  71   140 CONTINUE
  72
  73 C...Loop over attempts made. Reset counters.
  74       NTRY=0
  75   150 NTRY=NTRY+1
  76       IF(NTRY.GT.200) THEN
  77         CALL LUERRM(14,'(LUINDF:) caught in infinite loop')
  78         IF(MSTU(21).GE.1) RETURN
  79       ENDIF
  80       N=NSAV+NJET
  81       MSTU(90)=MSTU90
  82       DO 160 J=1,3
  83       NFL(J)=NFI(J)
  84       IFET(J)=0
  85       KFLF(J)=0
  86   160 CONTINUE
  87
  88 C...Loop over jets to be fragmented.
  89       DO 230 IP1=NSAV+1,NSAV+NJET
  90       MSTJ(91)=0
  91       NSAV1=N
  92       MSTU91=MSTU(90)
  93
  94 C...Initial flavour and momentum values. Jet along +z axis.
  95       KFLH=IABS(K(IP1,2))
  96       IF(KFLH.GT.10) KFLH=MOD(KFLH/1000,10)
  97       KFLO(2)=0
  98       WF=P(IP1,4)+SQRT(P(IP1,1)**2+P(IP1,2)**2+P(IP1,3)**2)
  99
 100 C...Initial values for quark or diquark jet.
 101   170 IF(IABS(K(IP1,2)).NE.21) THEN
 102         NSTR=1
 103         KFLO(1)=K(IP1,2)
 104         CALL LUPTDI(0,PXO(1),PYO(1))
 105         WO(1)=WF
 106
 107 C...Initial values for gluon treated like random quark jet.
 108       ELSEIF(MSTJ(2).LE.2) THEN
 109         NSTR=1
 110         IF(MSTJ(2).EQ.2) MSTJ(91)=1
 111         KFLO(1)=INT(1.+(2.+PARJ(2))*RLU(0))*(-1)**INT(RLU(0)+0.5)
 112         CALL LUPTDI(0,PXO(1),PYO(1))
 113         WO(1)=WF
 114
 115 C...Initial values for gluon treated like quark-antiquark jet pair,
 116 C...sharing energy according to Altarelli-Parisi splitting function.
 117       ELSE
 118         NSTR=2
 119         IF(MSTJ(2).EQ.4) MSTJ(91)=1
 120         KFLO(1)=INT(1.+(2.+PARJ(2))*RLU(0))*(-1)**INT(RLU(0)+0.5)
 121         KFLO(2)=-KFLO(1)
 122         CALL LUPTDI(0,PXO(1),PYO(1))
 123         PXO(2)=-PXO(1)
 124         PYO(2)=-PYO(1)
 125         WO(1)=WF*RLU(0)**(1./3.)
 126         WO(2)=WF-WO(1)
 127       ENDIF
 128
 129 C...Initial values for rank, flavour, pT and W+.
 130       DO 220 ISTR=1,NSTR
 131   180 I=N
 132       MSTU(90)=MSTU91
 133       IRANK=0
 134       KFL1=KFLO(ISTR)
 135       PX1=PXO(ISTR)
 136       PY1=PYO(ISTR)
 137       W=WO(ISTR)
 138
 139 C...New hadron. Generate flavour and hadron species.
 140   190 I=I+1
 141       IF(I.GE.MSTU(4)-MSTU(32)-NJET-5) THEN
 142         CALL LUERRM(11,'(LUINDF:) no more memory left in LUJETS')
 143         IF(MSTU(21).GE.1) RETURN
 144       ENDIF
 145       IRANK=IRANK+1
 146       K(I,1)=1
 147       K(I,3)=IP1
 148       K(I,4)=0
 149       K(I,5)=0
 150   200 CALL LUKFDI(KFL1,0,KFL2,K(I,2))
 151       IF(K(I,2).EQ.0) GOTO 180
 152       IF(MSTJ(12).GE.3.AND.IRANK.EQ.1.AND.IABS(KFL1).LE.10.AND.
 153      &IABS(KFL2).GT.10) THEN
 154         IF(RLU(0).GT.PARJ(19)) GOTO 200
 155       ENDIF
 156
 157 C...Find hadron mass. Generate four-momentum.
 158       P(I,5)=ULMASS(K(I,2))
 159       CALL LUPTDI(KFL1,PX2,PY2)
 160       P(I,1)=PX1+PX2
 161       P(I,2)=PY1+PY2
 162       PR=P(I,5)**2+P(I,1)**2+P(I,2)**2
 163       CALL LUZDIS(KFL1,KFL2,PR,Z)
 164       MZSAV=0
 165       IF(IABS(KFL1).GE.4.AND.IABS(KFL1).LE.8.AND.MSTU(90).LT.8) THEN
 166         MZSAV=1
 167         MSTU(90)=MSTU(90)+1
 168         MSTU(90+MSTU(90))=I
 169         PARU(90+MSTU(90))=Z
 170       ENDIF
 171       P(I,3)=0.5*(Z*W-PR/MAX(1E-4,Z*W))
 172       P(I,4)=0.5*(Z*W+PR/MAX(1E-4,Z*W))
 173       IF(MSTJ(3).GE.1.AND.IRANK.EQ.1.AND.KFLH.GE.4.AND.
 174      &P(I,3).LE.0.001) THEN
 175         IF(W.GE.P(I,5)+0.5*PARJ(32)) GOTO 180
 176         P(I,3)=0.0001
 177         P(I,4)=SQRT(PR)
 178         Z=P(I,4)/W
 179       ENDIF
 180
 181 C...Remaining flavour and momentum.
 182       KFL1=-KFL2
 183       PX1=-PX2
 184       PY1=-PY2
 185       W=(1.-Z)*W
 186       DO 210 J=1,5
 187       V(I,J)=0.
 188   210 CONTINUE
 189
 190 C...Check if pL acceptable. Go back for new hadron if enough energy.
 191       IF(MSTJ(3).GE.0.AND.P(I,3).LT.0.) THEN
 192         I=I-1
 193         IF(MZSAV.EQ.1) MSTU(90)=MSTU(90)-1
 194       ENDIF
 195       IF(W.GT.PARJ(31)) GOTO 190
 196       N=I
 197   220 CONTINUE
 198       IF(MOD(MSTJ(3),5).EQ.4.AND.N.EQ.NSAV1) WF=WF+0.1*PARJ(32)
 199       IF(MOD(MSTJ(3),5).EQ.4.AND.N.EQ.NSAV1) GOTO 170
 200
 201 C...Rotate jet to new direction.
 202       THE=ULANGL(P(IP1,3),SQRT(P(IP1,1)**2+P(IP1,2)**2))
 203       PHI=ULANGL(P(IP1,1),P(IP1,2))
 204       MSTU(33)=1
 205       CALL LUDBRB(NSAV1+1,N,THE,PHI,0D0,0D0,0D0)
 206       K(K(IP1,3),4)=NSAV1+1
 207       K(K(IP1,3),5)=N
 208
 209 C...End of jet generation loop. Skip conservation in some cases.
 210   230 CONTINUE
 211       IF(NJET.EQ.1.OR.MSTJ(3).LE.0) GOTO 490
 212       IF(MOD(MSTJ(3),5).NE.0.AND.N-NSAV-NJET.LT.2) GOTO 150
 213
 214 C...Subtract off produced hadron flavours, finished if zero.
 215       DO 240 I=NSAV+NJET+1,N
 216       KFA=IABS(K(I,2))
 217       KFLA=MOD(KFA/1000,10)
 218       KFLB=MOD(KFA/100,10)
 219       KFLC=MOD(KFA/10,10)
 220       IF(KFLA.EQ.0) THEN
 221         IF(KFLB.LE.3) NFL(KFLB)=NFL(KFLB)-ISIGN(1,K(I,2))*(-1)**KFLB
 222         IF(KFLC.LE.3) NFL(KFLC)=NFL(KFLC)+ISIGN(1,K(I,2))*(-1)**KFLB
 223       ELSE
 224         IF(KFLA.LE.3) NFL(KFLA)=NFL(KFLA)-ISIGN(1,K(I,2))
 225         IF(KFLB.LE.3) NFL(KFLB)=NFL(KFLB)-ISIGN(1,K(I,2))
 226         IF(KFLC.LE.3) NFL(KFLC)=NFL(KFLC)-ISIGN(1,K(I,2))
 227       ENDIF
 228   240 CONTINUE
 229       NREQ=(IABS(NFL(1))+IABS(NFL(2))+IABS(NFL(3))-IABS(NFL(1)+
 230      &NFL(2)+NFL(3)))/2+IABS(NFL(1)+NFL(2)+NFL(3))/3
 231       IF(NREQ.EQ.0) GOTO 320
 232
 233 C...Take away flavour of low-momentum particles until enough freedom.
 234       NREM=0
 235   250 IREM=0
 236       P2MIN=PECM**2
 237       DO 260 I=NSAV+NJET+1,N
 238       P2=P(I,1)**2+P(I,2)**2+P(I,3)**2
 239       IF(K(I,1).EQ.1.AND.P2.LT.P2MIN) IREM=I
 240       IF(K(I,1).EQ.1.AND.P2.LT.P2MIN) P2MIN=P2
 241   260 CONTINUE
 242       IF(IREM.EQ.0) GOTO 150
 243       K(IREM,1)=7
 244       KFA=IABS(K(IREM,2))
 245       KFLA=MOD(KFA/1000,10)
 246       KFLB=MOD(KFA/100,10)
 247       KFLC=MOD(KFA/10,10)
 248       IF(KFLA.GE.4.OR.KFLB.GE.4) K(IREM,1)=8
 249       IF(K(IREM,1).EQ.8) GOTO 250
 250       IF(KFLA.EQ.0) THEN
 251         ISGN=ISIGN(1,K(IREM,2))*(-1)**KFLB
 252         IF(KFLB.LE.3) NFL(KFLB)=NFL(KFLB)+ISGN
 253         IF(KFLC.LE.3) NFL(KFLC)=NFL(KFLC)-ISGN
 254       ELSE
 255         IF(KFLA.LE.3) NFL(KFLA)=NFL(KFLA)+ISIGN(1,K(IREM,2))
 256         IF(KFLB.LE.3) NFL(KFLB)=NFL(KFLB)+ISIGN(1,K(IREM,2))
 257         IF(KFLC.LE.3) NFL(KFLC)=NFL(KFLC)+ISIGN(1,K(IREM,2))
 258       ENDIF
 259       NREM=NREM+1
 260       NREQ=(IABS(NFL(1))+IABS(NFL(2))+IABS(NFL(3))-IABS(NFL(1)+
 261      &NFL(2)+NFL(3)))/2+IABS(NFL(1)+NFL(2)+NFL(3))/3
 262       IF(NREQ.GT.NREM) GOTO 250
 263       DO 270 I=NSAV+NJET+1,N
 264       IF(K(I,1).EQ.8) K(I,1)=1
 265   270 CONTINUE
 266
 267 C...Find combination of existing and new flavours for hadron.
 268   280 NFET=2
 269       IF(NFL(1)+NFL(2)+NFL(3).NE.0) NFET=3
 270       IF(NREQ.LT.NREM) NFET=1
 271       IF(IABS(NFL(1))+IABS(NFL(2))+IABS(NFL(3)).EQ.0) NFET=0
 272       DO 290 J=1,NFET
 273       IFET(J)=1+(IABS(NFL(1))+IABS(NFL(2))+IABS(NFL(3)))*RLU(0)
 274       KFLF(J)=ISIGN(1,NFL(1))
 275       IF(IFET(J).GT.IABS(NFL(1))) KFLF(J)=ISIGN(2,NFL(2))
 276       IF(IFET(J).GT.IABS(NFL(1))+IABS(NFL(2))) KFLF(J)=ISIGN(3,NFL(3))
 277   290 CONTINUE
 278       IF(NFET.EQ.2.AND.(IFET(1).EQ.IFET(2).OR.KFLF(1)*KFLF(2).GT.0))
 279      &GOTO 280
 280       IF(NFET.EQ.3.AND.(IFET(1).EQ.IFET(2).OR.IFET(1).EQ.IFET(3).OR.
 281      &IFET(2).EQ.IFET(3).OR.KFLF(1)*KFLF(2).LT.0.OR.KFLF(1)*KFLF(3)
 282      &.LT.0.OR.KFLF(1)*(NFL(1)+NFL(2)+NFL(3)).LT.0)) GOTO 280
 283       IF(NFET.EQ.0) KFLF(1)=1+INT((2.+PARJ(2))*RLU(0))
 284       IF(NFET.EQ.0) KFLF(2)=-KFLF(1)
 285       IF(NFET.EQ.1) KFLF(2)=ISIGN(1+INT((2.+PARJ(2))*RLU(0)),-KFLF(1))
 286       IF(NFET.LE.2) KFLF(3)=0
 287       IF(KFLF(3).NE.0) THEN
 288         KFLFC=ISIGN(1000*MAX(IABS(KFLF(1)),IABS(KFLF(3)))+
 289      &  100*MIN(IABS(KFLF(1)),IABS(KFLF(3)))+1,KFLF(1))
 290         IF(KFLF(1).EQ.KFLF(3).OR.(1.+3.*PARJ(4))*RLU(0).GT.1.)
 291      &  KFLFC=KFLFC+ISIGN(2,KFLFC)
 292       ELSE
 293         KFLFC=KFLF(1)
 294       ENDIF
 295       CALL LUKFDI(KFLFC,KFLF(2),KFLDMP,KF)
 296       IF(KF.EQ.0) GOTO 280
 297       DO 300 J=1,MAX(2,NFET)
 298       NFL(IABS(KFLF(J)))=NFL(IABS(KFLF(J)))-ISIGN(1,KFLF(J))
 299   300 CONTINUE
 300
 301 C...Store hadron at random among free positions.
 302       NPOS=MIN(1+INT(RLU(0)*NREM),NREM)
 303       DO 310 I=NSAV+NJET+1,N
 304       IF(K(I,1).EQ.7) NPOS=NPOS-1
 305       IF(K(I,1).EQ.1.OR.NPOS.NE.0) GOTO 310
 306       K(I,1)=1
 307       K(I,2)=KF
 308       P(I,5)=ULMASS(K(I,2))
 309       P(I,4)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2+P(I,5)**2)
 310   310 CONTINUE
 311       NREM=NREM-1
 312       NREQ=(IABS(NFL(1))+IABS(NFL(2))+IABS(NFL(3))-IABS(NFL(1)+
 313      &NFL(2)+NFL(3)))/2+IABS(NFL(1)+NFL(2)+NFL(3))/3
 314       IF(NREM.GT.0) GOTO 280
 315
 316 C...Compensate for missing momentum in global scheme (3 options).
 317   320 IF(MOD(MSTJ(3),5).NE.0.AND.MOD(MSTJ(3),5).NE.4) THEN
 318         DO 340 J=1,3
 319         PSI(J)=0.
 320         DO 330 I=NSAV+NJET+1,N
 321         PSI(J)=PSI(J)+P(I,J)
 322   330   CONTINUE
 323   340   CONTINUE
 324         PSI(4)=PSI(1)**2+PSI(2)**2+PSI(3)**2
 325         PWS=0.
 326         DO 350 I=NSAV+NJET+1,N
 327         IF(MOD(MSTJ(3),5).EQ.1) PWS=PWS+P(I,4)
 328         IF(MOD(MSTJ(3),5).EQ.2) PWS=PWS+SQRT(P(I,5)**2+(PSI(1)*P(I,1)+
 329      &  PSI(2)*P(I,2)+PSI(3)*P(I,3))**2/PSI(4))
 330         IF(MOD(MSTJ(3),5).EQ.3) PWS=PWS+1.
 331   350   CONTINUE
 332         DO 370 I=NSAV+NJET+1,N
 333         IF(MOD(MSTJ(3),5).EQ.1) PW=P(I,4)
 334         IF(MOD(MSTJ(3),5).EQ.2) PW=SQRT(P(I,5)**2+(PSI(1)*P(I,1)+
 335      &  PSI(2)*P(I,2)+PSI(3)*P(I,3))**2/PSI(4))
 336         IF(MOD(MSTJ(3),5).EQ.3) PW=1.
 337         DO 360 J=1,3
 338         P(I,J)=P(I,J)-PSI(J)*PW/PWS
 339   360   CONTINUE
 340         P(I,4)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2+P(I,5)**2)
 341   370   CONTINUE
 342
 343 C...Compensate for missing momentum withing each jet separately.
 344       ELSEIF(MOD(MSTJ(3),5).EQ.4) THEN
 345         DO 390 I=N+1,N+NJET
 346         K(I,1)=0
 347         DO 380 J=1,5
 348         P(I,J)=0.
 349   380   CONTINUE
 350   390   CONTINUE
 351         DO 410 I=NSAV+NJET+1,N
 352         IR1=K(I,3)
 353         IR2=N+IR1-NSAV
 354         K(IR2,1)=K(IR2,1)+1
 355         PLS=(P(I,1)*P(IR1,1)+P(I,2)*P(IR1,2)+P(I,3)*P(IR1,3))/
 356      &  (P(IR1,1)**2+P(IR1,2)**2+P(IR1,3)**2)
 357         DO 400 J=1,3
 358         P(IR2,J)=P(IR2,J)+P(I,J)-PLS*P(IR1,J)
 359   400   CONTINUE
 360         P(IR2,4)=P(IR2,4)+P(I,4)
 361         P(IR2,5)=P(IR2,5)+PLS
 362   410   CONTINUE
 363         PSS=0.
 364         DO 420 I=N+1,N+NJET
 365         IF(K(I,1).NE.0) PSS=PSS+P(I,4)/(PECM*(0.8*P(I,5)+0.2))
 366   420   CONTINUE
 367         DO 440 I=NSAV+NJET+1,N
 368         IR1=K(I,3)
 369         IR2=N+IR1-NSAV
 370         PLS=(P(I,1)*P(IR1,1)+P(I,2)*P(IR1,2)+P(I,3)*P(IR1,3))/
 371      &  (P(IR1,1)**2+P(IR1,2)**2+P(IR1,3)**2)
 372         DO 430 J=1,3
 373         P(I,J)=P(I,J)-P(IR2,J)/K(IR2,1)+(1./(P(IR2,5)*PSS)-1.)*PLS*
 374      &  P(IR1,J)
 375   430   CONTINUE
 376         P(I,4)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2+P(I,5)**2)
 377   440   CONTINUE
 378       ENDIF
 379
 380 C...Scale momenta for energy conservation.
 381       IF(MOD(MSTJ(3),5).NE.0) THEN
 382         PMS=0.
 383         PES=0.
 384         PQS=0.
 385         DO 450 I=NSAV+NJET+1,N
 386         PMS=PMS+P(I,5)
 387         PES=PES+P(I,4)
 388         PQS=PQS+P(I,5)**2/P(I,4)
 389   450   CONTINUE
 390         IF(PMS.GE.PECM) GOTO 150
 391         NECO=0
 392   460   NECO=NECO+1
 393         PFAC=(PECM-PQS)/(PES-PQS)
 394         PES=0.
 395         PQS=0.
 396         DO 480 I=NSAV+NJET+1,N
 397         DO 470 J=1,3
 398         P(I,J)=PFAC*P(I,J)
 399   470   CONTINUE
 400         P(I,4)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2+P(I,5)**2)
 401         PES=PES+P(I,4)
 402         PQS=PQS+P(I,5)**2/P(I,4)
 403   480   CONTINUE
 404         IF(NECO.LT.10.AND.ABS(PECM-PES).GT.2E-6*PECM) GOTO 460
 405       ENDIF
 406
 407 C...Origin of produced particles and parton daughter pointers.
 408   490 DO 500 I=NSAV+NJET+1,N
 409       IF(MSTU(16).NE.2) K(I,3)=NSAV+1
 410       IF(MSTU(16).EQ.2) K(I,3)=K(K(I,3),3)
 411   500 CONTINUE
 412       DO 510 I=NSAV+1,NSAV+NJET
 413       I1=K(I,3)
 414       K(I1,1)=K(I1,1)+10
 415       IF(MSTU(16).NE.2) THEN
 416         K(I1,4)=NSAV+1
 417         K(I1,5)=NSAV+1
 418       ELSE
 419         K(I1,4)=K(I1,4)-NJET+1
 420         K(I1,5)=K(I1,5)-NJET+1
 421         IF(K(I1,5).LT.K(I1,4)) THEN
 422           K(I1,4)=0
 423           K(I1,5)=0
 424         ENDIF
 425       ENDIF
 426   510 CONTINUE
 427
 428 C...Document independent fragmentation system. Remove copy of jets.
 429       NSAV=NSAV+1
 430       K(NSAV,1)=11
 431       K(NSAV,2)=93
 432       K(NSAV,3)=IP
 433       K(NSAV,4)=NSAV+1
 434       K(NSAV,5)=N-NJET+1
 435       DO 520 J=1,4
 436       P(NSAV,J)=DPS(J)
 437       V(NSAV,J)=V(IP,J)
 438   520 CONTINUE
 439       P(NSAV,5)=SQRT(MAX(0D0,DPS(4)**2-DPS(1)**2-DPS(2)**2-DPS(3)**2))
 440       V(NSAV,5)=0.
 441       DO 540 I=NSAV+NJET,N
 442       DO 530 J=1,5
 443       K(I-NJET+1,J)=K(I,J)
 444       P(I-NJET+1,J)=P(I,J)
 445       V(I-NJET+1,J)=V(I,J)
 446   530 CONTINUE
 447   540 CONTINUE
 448       N=N-NJET+1
 449       DO 550 IZ=MSTU90+1,MSTU(90)
 450       MSTU(90+IZ)=MSTU(90+IZ)-NJET+1
 451   550 CONTINUE
 452
 453 C...Boost back particle system. Set production vertices.
 454       IF(NJET.NE.1) CALL LUDBRB(NSAV+1,N,0.,0.,DPS(1)/DPS(4),
 455      &DPS(2)/DPS(4),DPS(3)/DPS(4))
 456       DO 570 I=NSAV+1,N
 457       DO 560 J=1,4
 458       V(I,J)=V(IP,J)
 459   560 CONTINUE
 460   570 CONTINUE
 461
 462       RETURN
 463       END