PYTHIA/jetset/ludecy.F

   1
   2 C*********************************************************************
   3
   4       SUBROUTINE LUDECY(IP)
   5
   6 C...Purpose: to handle the decay of unstable particles.
   7       COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5)
   8       COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
   9       COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
  10       COMMON/LUDAT3/MDCY(500,3),MDME(2000,2),BRAT(2000),KFDP(2000,5)
  11       SAVE /LUJETS/,/LUDAT1/,/LUDAT2/,/LUDAT3/
  12       DIMENSION VDCY(4),KFLO(4),KFL1(4),PV(10,5),RORD(10),UE(3),BE(3),
  13      &WTCOR(10),PTAU(4),PCMTAU(4)
  14       DOUBLE PRECISION DBETAU(3)
  15       DATA WTCOR/2.,5.,15.,60.,250.,1500.,1.2E4,1.2E5,150.,16./
  16
  17 C...Functions: momentum in two-particle decays, four-product and
  18 C...matrix element times phase space in weak decays.
  19       PAWT(A,B,C)=SQRT((A**2-(B+C)**2)*(A**2-(B-C)**2))/(2.*A)
  20       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)
  21       HMEPS(HA)=((1.-HRQ-HA)**2+3.*HA*(1.+HRQ-HA))*
  22      &SQRT((1.-HRQ-HA)**2-4.*HRQ*HA)
  23
  24 C...Initial values.
  25       NTRY=0
  26       NSAV=N
  27       KFA=IABS(K(IP,2))
  28       KFS=ISIGN(1,K(IP,2))
  29       KC=LUCOMP(KFA)
  30       MSTJ(92)=0
  31
  32 C...Choose lifetime and determine decay vertex.
  33       IF(K(IP,1).EQ.5) THEN
  34         V(IP,5)=0.
  35       ELSEIF(K(IP,1).NE.4) THEN
  36         V(IP,5)=-PMAS(KC,4)*LOG(RLU(0))
  37       ENDIF
  38       DO 100 J=1,4
  39       VDCY(J)=V(IP,J)+V(IP,5)*P(IP,J)/P(IP,5)
  40   100 CONTINUE
  41
  42 C...Determine whether decay allowed or not.
  43       MOUT=0
  44       IF(MSTJ(22).EQ.2) THEN
  45         IF(PMAS(KC,4).GT.PARJ(71)) MOUT=1
  46       ELSEIF(MSTJ(22).EQ.3) THEN
  47         IF(VDCY(1)**2+VDCY(2)**2+VDCY(3)**2.GT.PARJ(72)**2) MOUT=1
  48       ELSEIF(MSTJ(22).EQ.4) THEN
  49         IF(VDCY(1)**2+VDCY(2)**2.GT.PARJ(73)**2) MOUT=1
  50         IF(ABS(VDCY(3)).GT.PARJ(74)) MOUT=1
  51       ENDIF
  52       IF(MOUT.EQ.1.AND.K(IP,1).NE.5) THEN
  53         K(IP,1)=4
  54         RETURN
  55       ENDIF
  56
  57 C...Interface to external tau decay library (for tau polarization).
  58       IF(KFA.EQ.15.AND.MSTJ(28).GE.1) THEN
  59
  60 C...Starting values for pointers and momenta.
  61         ITAU=IP
  62         DO 110 J=1,4
  63         PTAU(J)=P(ITAU,J)
  64         PCMTAU(J)=P(ITAU,J)
  65   110   CONTINUE
  66
  67 C...Iterate to find position and code of mother of tau.
  68         IMTAU=ITAU
  69   120   IMTAU=K(IMTAU,3)
  70
  71         IF(IMTAU.EQ.0) THEN
  72 C...If no known origin then impossible to do anything further.
  73           KFORIG=0
  74           IORIG=0
  75
  76         ELSEIF(K(IMTAU,2).EQ.K(ITAU,2)) THEN
  77 C...If tau -> tau + gamma then add gamma energy and loop.
  78           IF(K(K(IMTAU,4),2).EQ.22) THEN
  79             DO 130 J=1,4
  80             PCMTAU(J)=PCMTAU(J)+P(K(IMTAU,4),J)
  81   130       CONTINUE
  82           ELSEIF(K(K(IMTAU,5),2).EQ.22) THEN
  83             DO 140 J=1,4
  84             PCMTAU(J)=PCMTAU(J)+P(K(IMTAU,5),J)
  85   140       CONTINUE
  86           ENDIF
  87           GOTO 120
  88
  89         ELSEIF(IABS(K(IMTAU,2)).GT.100) THEN
  90 C...If coming from weak decay of hadron then W is not stored in record,
  91 C...but can be reconstructed by adding neutrino momentum.
  92           KFORIG=-ISIGN(24,K(ITAU,2))
  93           IORIG=0
  94           DO 160 II=K(IMTAU,4),K(IMTAU,5)
  95           IF(K(II,2)*ISIGN(1,K(ITAU,2)).EQ.-16) THEN
  96             DO 150 J=1,4
  97             PCMTAU(J)=PCMTAU(J)+P(II,J)
  98   150       CONTINUE
  99           ENDIF
 100   160     CONTINUE
 101
 102         ELSE
 103 C...If coming from resonance decay then find latest copy of this
 104 C...resonance (may not completely agree).
 105           KFORIG=K(IMTAU,2)
 106           IORIG=IMTAU
 107           DO 170 II=IMTAU+1,IP-1
 108           IF(K(II,2).EQ.KFORIG.AND.K(II,3).EQ.IORIG.AND.
 109      &    ABS(P(II,5)-P(IORIG,5)).LT.1E-5*P(IORIG,5)) IORIG=II
 110   170     CONTINUE
 111           DO 180 J=1,4
 112           PCMTAU(J)=P(IORIG,J)
 113   180     CONTINUE
 114         ENDIF
 115
 116 C...Boost tau to rest frame of production process (where known)
 117 C...and rotate it to sit along +z axis.
 118         DO 190 J=1,3
 119         DBETAU(J)=PCMTAU(J)/PCMTAU(4)
 120   190   CONTINUE
 121         IF(KFORIG.NE.0) CALL LUDBRB(ITAU,ITAU,0.,0.,-DBETAU(1),
 122      &  -DBETAU(2),-DBETAU(3))
 123         PHITAU=ULANGL(P(ITAU,1),P(ITAU,2))
 124         CALL LUDBRB(ITAU,ITAU,0.,-PHITAU,0D0,0D0,0D0)
 125         THETAU=ULANGL(P(ITAU,3),P(ITAU,1))
 126         CALL LUDBRB(ITAU,ITAU,-THETAU,0.,0D0,0D0,0D0)
 127
 128 C...Call tau decay routine (if meaningful) and fill extra info.
 129         IF(KFORIG.NE.0.OR.MSTJ(28).EQ.2) THEN
 130           CALL LUTAUD(ITAU,IORIG,KFORIG,NDECAY)
 131           DO 200 II=NSAV+1,NSAV+NDECAY
 132           K(II,1)=1
 133           K(II,3)=IP
 134           K(II,4)=0
 135           K(II,5)=0
 136   200     CONTINUE
 137           N=NSAV+NDECAY
 138         ENDIF
 139
 140 C...Boost back decay tau and decay products.
 141         DO 210 J=1,4
 142         P(ITAU,J)=PTAU(J)
 143   210   CONTINUE
 144         IF(KFORIG.NE.0.OR.MSTJ(28).EQ.2) THEN
 145           CALL LUDBRB(NSAV+1,N,THETAU,PHITAU,0D0,0D0,0D0)
 146           IF(KFORIG.NE.0) CALL LUDBRB(NSAV+1,N,0.,0.,DBETAU(1),
 147      &    DBETAU(2),DBETAU(3))
 148
 149 C...Skip past ordinary tau decay treatment.
 150           MMAT=0
 151           MBST=0
 152           ND=0
 153           GOTO 660
 154         ENDIF
 155       ENDIF
 156
 157 C...B-B~ mixing: flip sign of meson appropriately.
 158       MMIX=0
 159       IF((KFA.EQ.511.OR.KFA.EQ.531).AND.MSTJ(26).GE.1) THEN
 160         XBBMIX=PARJ(76)
 161         IF(KFA.EQ.531) XBBMIX=PARJ(77)
 162         IF(SIN(0.5*XBBMIX*V(IP,5)/PMAS(KC,4))**2.GT.RLU(0)) MMIX=1
 163         IF(MMIX.EQ.1) KFS=-KFS
 164       ENDIF
 165
 166 C...Check existence of decay channels. Particle/antiparticle rules.
 167       KCA=KC
 168       IF(MDCY(KC,2).GT.0) THEN
 169         MDMDCY=MDME(MDCY(KC,2),2)
 170         IF(MDMDCY.GT.80.AND.MDMDCY.LE.90) KCA=MDMDCY
 171       ENDIF
 172       IF(MDCY(KCA,2).LE.0.OR.MDCY(KCA,3).LE.0) THEN
 173         CALL LUERRM(9,'(LUDECY:) no decay channel defined')
 174         RETURN
 175       ENDIF
 176       IF(MOD(KFA/1000,10).EQ.0.AND.(KCA.EQ.85.OR.KCA.EQ.87)) KFS=-KFS
 177       IF(KCHG(KC,3).EQ.0) THEN
 178         KFSP=1
 179         KFSN=0
 180         IF(RLU(0).GT.0.5) KFS=-KFS
 181       ELSEIF(KFS.GT.0) THEN
 182         KFSP=1
 183         KFSN=0
 184       ELSE
 185         KFSP=0
 186         KFSN=1
 187       ENDIF
 188
 189 C...Sum branching ratios of allowed decay channels.
 190   220 NOPE=0
 191       BRSU=0.
 192       DO 230 IDL=MDCY(KCA,2),MDCY(KCA,2)+MDCY(KCA,3)-1
 193       IF(MDME(IDL,1).NE.1.AND.KFSP*MDME(IDL,1).NE.2.AND.
 194      &KFSN*MDME(IDL,1).NE.3) GOTO 230
 195       IF(MDME(IDL,2).GT.100) GOTO 230
 196       NOPE=NOPE+1
 197       BRSU=BRSU+BRAT(IDL)
 198   230 CONTINUE
 199       IF(NOPE.EQ.0) THEN
 200         CALL LUERRM(2,'(LUDECY:) all decay channels closed by user')
 201         RETURN
 202       ENDIF
 203
 204 C...Select decay channel among allowed ones.
 205   240 RBR=BRSU*RLU(0)
 206       IDL=MDCY(KCA,2)-1
 207   250 IDL=IDL+1
 208       IF(MDME(IDL,1).NE.1.AND.KFSP*MDME(IDL,1).NE.2.AND.
 209      &KFSN*MDME(IDL,1).NE.3) THEN
 210         IF(IDL.LT.MDCY(KCA,2)+MDCY(KCA,3)-1) GOTO 250
 211       ELSEIF(MDME(IDL,2).GT.100) THEN
 212         IF(IDL.LT.MDCY(KCA,2)+MDCY(KCA,3)-1) GOTO 250
 213       ELSE
 214         IDC=IDL
 215         RBR=RBR-BRAT(IDL)
 216         IF(IDL.LT.MDCY(KCA,2)+MDCY(KCA,3)-1.AND.RBR.GT.0.) GOTO 250
 217       ENDIF
 218
 219 C...Start readout of decay channel: matrix element, reset counters.
 220       MMAT=MDME(IDC,2)
 221   260 NTRY=NTRY+1
 222       IF(NTRY.GT.1000) THEN
 223         CALL LUERRM(14,'(LUDECY:) caught in infinite loop')
 224         IF(MSTU(21).GE.1) RETURN
 225       ENDIF
 226       I=N
 227       NP=0
 228       NQ=0
 229       MBST=0
 230       IF(MMAT.GE.11.AND.MMAT.NE.46.AND.P(IP,4).GT.20.*P(IP,5)) MBST=1
 231       DO 270 J=1,4
 232       PV(1,J)=0.
 233       IF(MBST.EQ.0) PV(1,J)=P(IP,J)
 234   270 CONTINUE
 235       IF(MBST.EQ.1) PV(1,4)=P(IP,5)
 236       PV(1,5)=P(IP,5)
 237       PS=0.
 238       PSQ=0.
 239       MREM=0
 240       MHADDY=0
 241       IF(KFA.GT.80) MHADDY=1
 242
 243 C...Read out decay products. Convert to standard flavour code.
 244       JTMAX=5
 245       IF(MDME(IDC+1,2).EQ.101) JTMAX=10
 246       DO 280 JT=1,JTMAX
 247       IF(JT.LE.5) KP=KFDP(IDC,JT)
 248       IF(JT.GE.6) KP=KFDP(IDC+1,JT-5)
 249       IF(KP.EQ.0) GOTO 280
 250       KPA=IABS(KP)
 251       KCP=LUCOMP(KPA)
 252       IF(KPA.GT.80) MHADDY=1
 253       IF(KCHG(KCP,3).EQ.0.AND.KPA.NE.81.AND.KPA.NE.82) THEN
 254         KFP=KP
 255       ELSEIF(KPA.NE.81.AND.KPA.NE.82) THEN
 256         KFP=KFS*KP
 257       ELSEIF(KPA.EQ.81.AND.MOD(KFA/1000,10).EQ.0) THEN
 258         KFP=-KFS*MOD(KFA/10,10)
 259       ELSEIF(KPA.EQ.81.AND.MOD(KFA/100,10).GE.MOD(KFA/10,10)) THEN
 260         KFP=KFS*(100*MOD(KFA/10,100)+3)
 261       ELSEIF(KPA.EQ.81) THEN
 262         KFP=KFS*(1000*MOD(KFA/10,10)+100*MOD(KFA/100,10)+1)
 263       ELSEIF(KP.EQ.82) THEN
 264         CALL LUKFDI(-KFS*INT(1.+(2.+PARJ(2))*RLU(0)),0,KFP,KDUMP)
 265         IF(KFP.EQ.0) GOTO 260
 266         MSTJ(93)=1
 267         IF(PV(1,5).LT.PARJ(32)+2.*ULMASS(KFP)) GOTO 260
 268       ELSEIF(KP.EQ.-82) THEN
 269         KFP=-KFP
 270         IF(IABS(KFP).GT.10) KFP=KFP+ISIGN(10000,KFP)
 271       ENDIF
 272       IF(KPA.EQ.81.OR.KPA.EQ.82) KCP=LUCOMP(KFP)
 273
 274 C...Add decay product to event record or to quark flavour list.
 275       KFPA=IABS(KFP)
 276       KQP=KCHG(KCP,2)
 277       IF(MMAT.GE.11.AND.MMAT.LE.30.AND.KQP.NE.0) THEN
 278         NQ=NQ+1
 279         KFLO(NQ)=KFP
 280         MSTJ(93)=2
 281         PSQ=PSQ+ULMASS(KFLO(NQ))
 282       ELSEIF((MMAT.EQ.42.OR.MMAT.EQ.43.OR.MMAT.EQ.48).AND.NP.EQ.3.AND.
 283      &MOD(NQ,2).EQ.1) THEN
 284         NQ=NQ-1
 285         PS=PS-P(I,5)
 286         K(I,1)=1
 287         KFI=K(I,2)
 288         CALL LUKFDI(KFP,KFI,KFLDMP,K(I,2))
 289         IF(K(I,2).EQ.0) GOTO 260
 290         MSTJ(93)=1
 291         P(I,5)=ULMASS(K(I,2))
 292         PS=PS+P(I,5)
 293       ELSE
 294         I=I+1
 295         NP=NP+1
 296         IF(MMAT.NE.33.AND.KQP.NE.0) NQ=NQ+1
 297         IF(MMAT.EQ.33.AND.KQP.NE.0.AND.KQP.NE.2) NQ=NQ+1
 298         K(I,1)=1+MOD(NQ,2)
 299         IF(MMAT.EQ.4.AND.JT.LE.2.AND.KFP.EQ.21) K(I,1)=2
 300         IF(MMAT.EQ.4.AND.JT.EQ.3) K(I,1)=1
 301         K(I,2)=KFP
 302         K(I,3)=IP
 303         K(I,4)=0
 304         K(I,5)=0
 305         P(I,5)=ULMASS(KFP)
 306         IF(MMAT.EQ.45.AND.KFPA.EQ.89) P(I,5)=PARJ(32)
 307         PS=PS+P(I,5)
 308       ENDIF
 309   280 CONTINUE
 310
 311 C...Check masses for resonance decays.
 312       IF(MHADDY.EQ.0) THEN
 313         IF(PS+PARJ(64).GT.PV(1,5)) GOTO 240
 314       ENDIF
 315
 316 C...Choose decay multiplicity in phase space model.
 317   290 IF(MMAT.GE.11.AND.MMAT.LE.30) THEN
 318         PSP=PS
 319         CNDE=PARJ(61)*LOG(MAX((PV(1,5)-PS-PSQ)/PARJ(62),1.1))
 320         IF(MMAT.EQ.12) CNDE=CNDE+PARJ(63)
 321   300   NTRY=NTRY+1
 322         IF(NTRY.GT.1000) THEN
 323           CALL LUERRM(14,'(LUDECY:) caught in infinite loop')
 324           IF(MSTU(21).GE.1) RETURN
 325         ENDIF
 326         IF(MMAT.LE.20) THEN
 327           GAUSS=SQRT(-2.*CNDE*LOG(MAX(1E-10,RLU(0))))*
 328      &    SIN(PARU(2)*RLU(0))
 329           ND=0.5+0.5*NP+0.25*NQ+CNDE+GAUSS
 330           IF(ND.LT.NP+NQ/2.OR.ND.LT.2.OR.ND.GT.10) GOTO 300
 331           IF(MMAT.EQ.13.AND.ND.EQ.2) GOTO 300
 332           IF(MMAT.EQ.14.AND.ND.LE.3) GOTO 300
 333           IF(MMAT.EQ.15.AND.ND.LE.4) GOTO 300
 334         ELSE
 335           ND=MMAT-20
 336         ENDIF
 337
 338 C...Form hadrons from flavour content.
 339         DO 310 JT=1,4
 340         KFL1(JT)=KFLO(JT)
 341   310   CONTINUE
 342         IF(ND.EQ.NP+NQ/2) GOTO 330
 343         DO 320 I=N+NP+1,N+ND-NQ/2
 344         JT=1+INT((NQ-1)*RLU(0))
 345         CALL LUKFDI(KFL1(JT),0,KFL2,K(I,2))
 346         IF(K(I,2).EQ.0) GOTO 300
 347         KFL1(JT)=-KFL2
 348   320   CONTINUE
 349   330   JT=2
 350         JT2=3
 351         JT3=4
 352         IF(NQ.EQ.4.AND.RLU(0).LT.PARJ(66)) JT=4
 353         IF(JT.EQ.4.AND.ISIGN(1,KFL1(1)*(10-IABS(KFL1(1))))*
 354      &  ISIGN(1,KFL1(JT)*(10-IABS(KFL1(JT)))).GT.0) JT=3
 355         IF(JT.EQ.3) JT2=2
 356         IF(JT.EQ.4) JT3=2
 357         CALL LUKFDI(KFL1(1),KFL1(JT),KFLDMP,K(N+ND-NQ/2+1,2))
 358         IF(K(N+ND-NQ/2+1,2).EQ.0) GOTO 300
 359         IF(NQ.EQ.4) CALL LUKFDI(KFL1(JT2),KFL1(JT3),KFLDMP,K(N+ND,2))
 360         IF(NQ.EQ.4.AND.K(N+ND,2).EQ.0) GOTO 300
 361
 362 C...Check that sum of decay product masses not too large.
 363         PS=PSP
 364         DO 340 I=N+NP+1,N+ND
 365         K(I,1)=1
 366         K(I,3)=IP
 367         K(I,4)=0
 368         K(I,5)=0
 369         P(I,5)=ULMASS(K(I,2))
 370         PS=PS+P(I,5)
 371   340   CONTINUE
 372         IF(PS+PARJ(64).GT.PV(1,5)) GOTO 300
 373
 374 C...Rescale energy to subtract off spectator quark mass.
 375       ELSEIF((MMAT.EQ.31.OR.MMAT.EQ.33.OR.MMAT.EQ.44.OR.MMAT.EQ.45)
 376      &.AND.NP.GE.3) THEN
 377         PS=PS-P(N+NP,5)
 378         PQT=(P(N+NP,5)+PARJ(65))/PV(1,5)
 379         DO 350 J=1,5
 380         P(N+NP,J)=PQT*PV(1,J)
 381         PV(1,J)=(1.-PQT)*PV(1,J)
 382   350   CONTINUE
 383         IF(PS+PARJ(64).GT.PV(1,5)) GOTO 260
 384         ND=NP-1
 385         MREM=1
 386
 387 C...Phase space factors imposed in W decay.
 388       ELSEIF(MMAT.EQ.46) THEN
 389         MSTJ(93)=1
 390         PSMC=ULMASS(K(N+1,2))
 391         MSTJ(93)=1
 392         PSMC=PSMC+ULMASS(K(N+2,2))
 393         IF(MAX(PS,PSMC)+PARJ(32).GT.PV(1,5)) GOTO 240
 394         HR1=(P(N+1,5)/PV(1,5))**2
 395         HR2=(P(N+2,5)/PV(1,5))**2
 396         IF((1.-HR1-HR2)*(2.+HR1+HR2)*SQRT((1.-HR1-HR2)**2-4.*HR1*HR2)
 397      &  .LT.2.*RLU(0)) GOTO 240
 398         ND=NP
 399
 400 C...Fully specified final state: check mass broadening effects.
 401       ELSE
 402         IF(NP.GE.2.AND.PS+PARJ(64).GT.PV(1,5)) GOTO 260
 403         ND=NP
 404       ENDIF
 405
 406 C...Select W mass in decay Q -> W + q, without W propagator.
 407       IF(MMAT.EQ.45.AND.MSTJ(25).LE.0) THEN
 408         HLQ=(PARJ(32)/PV(1,5))**2
 409         HUQ=(1.-(P(N+2,5)+PARJ(64))/PV(1,5))**2
 410         HRQ=(P(N+2,5)/PV(1,5))**2
 411   360   HW=HLQ+RLU(0)*(HUQ-HLQ)
 412         IF(HMEPS(HW).LT.RLU(0)) GOTO 360
 413         P(N+1,5)=PV(1,5)*SQRT(HW)
 414
 415 C...Ditto, including W propagator. Divide mass range into three regions.
 416       ELSEIF(MMAT.EQ.45) THEN
 417         HQW=(PV(1,5)/PMAS(24,1))**2
 418         HLW=(PARJ(32)/PMAS(24,1))**2
 419         HUW=((PV(1,5)-P(N+2,5)-PARJ(64))/PMAS(24,1))**2
 420         HRQ=(P(N+2,5)/PV(1,5))**2
 421         HG=PMAS(24,2)/PMAS(24,1)
 422         HATL=ATAN((HLW-1.)/HG)
 423         HM=MIN(1.,HUW-0.001)
 424         HMV1=HMEPS(HM/HQW)/((HM-1.)**2+HG**2)
 425   370   HM=HM-HG
 426         HMV2=HMEPS(HM/HQW)/((HM-1.)**2+HG**2)
 427         IF(HMV2.GT.HMV1.AND.HM-HG.GT.HLW) THEN
 428           HMV1=HMV2
 429           GOTO 370
 430         ENDIF
 431         HMV=MIN(2.*HMV1,HMEPS(HM/HQW)/HG**2)
 432         HM1=1.-SQRT(1./HMV-HG**2)
 433         IF(HM1.GT.HLW.AND.HM1.LT.HM) THEN
 434           HM=HM1
 435         ELSEIF(HMV2.LE.HMV1) THEN
 436           HM=MAX(HLW,HM-MIN(0.1,1.-HM))
 437         ENDIF
 438         HATM=ATAN((HM-1.)/HG)
 439         HWT1=(HATM-HATL)/HG
 440         HWT2=HMV*(MIN(1.,HUW)-HM)
 441         HWT3=0.
 442         IF(HUW.GT.1.) THEN
 443           HATU=ATAN((HUW-1.)/HG)
 444           HMP1=HMEPS(1./HQW)
 445           HWT3=HMP1*HATU/HG
 446         ENDIF
 447
 448 C...Select mass region and W mass there. Accept according to weight.
 449   380   HREG=RLU(0)*(HWT1+HWT2+HWT3)
 450         IF(HREG.LE.HWT1) THEN
 451           HW=1.+HG*TAN(HATL+RLU(0)*(HATM-HATL))
 452           HACC=HMEPS(HW/HQW)
 453         ELSEIF(HREG.LE.HWT1+HWT2) THEN
 454           HW=HM+RLU(0)*(MIN(1.,HUW)-HM)
 455           HACC=HMEPS(HW/HQW)/((HW-1.)**2+HG**2)/HMV
 456         ELSE
 457           HW=1.+HG*TAN(RLU(0)*HATU)
 458           HACC=HMEPS(HW/HQW)/HMP1
 459         ENDIF
 460         IF(HACC.LT.RLU(0)) GOTO 380
 461         P(N+1,5)=PMAS(24,1)*SQRT(HW)
 462       ENDIF
 463
 464 C...Determine position of grandmother, number of sisters, Q -> W sign.
 465       NM=0
 466       KFAS=0
 467       MSGN=0
 468       IF(MMAT.EQ.3.OR.MMAT.EQ.46) THEN
 469         IM=K(IP,3)
 470         IF(IM.LT.0.OR.IM.GE.IP) IM=0
 471         IF(MMAT.EQ.46.AND.MSTJ(27).EQ.1) THEN
 472           IM=0
 473         ELSEIF(MMAT.EQ.46.AND.MSTJ(27).GE.2.AND.IM.NE.0) THEN
 474           IF(K(IM,2).EQ.94) THEN
 475             IM=K(K(IM,3),3)
 476             IF(IM.LT.0.OR.IM.GE.IP) IM=0
 477           ENDIF
 478         ENDIF
 479         IF(IM.NE.0) KFAM=IABS(K(IM,2))
 480         IF(IM.NE.0.AND.MMAT.EQ.3) THEN
 481           DO 390 IL=MAX(IP-2,IM+1),MIN(IP+2,N)
 482           IF(K(IL,3).EQ.IM) NM=NM+1
 483           IF(K(IL,3).EQ.IM.AND.IL.NE.IP) ISIS=IL
 484   390     CONTINUE
 485           IF(NM.NE.2.OR.KFAM.LE.100.OR.MOD(KFAM,10).NE.1.OR.
 486      &    MOD(KFAM/1000,10).NE.0) NM=0
 487           IF(NM.EQ.2) THEN
 488             KFAS=IABS(K(ISIS,2))
 489             IF((KFAS.LE.100.OR.MOD(KFAS,10).NE.1.OR.
 490      &      MOD(KFAS/1000,10).NE.0).AND.KFAS.NE.22) NM=0
 491           ENDIF
 492         ELSEIF(IM.NE.0.AND.MMAT.EQ.46) THEN
 493           MSGN=ISIGN(1,K(IM,2)*K(IP,2))
 494           IF(KFAM.GT.100.AND.MOD(KFAM/1000,10).EQ.0) MSGN=
 495      &    MSGN*(-1)**MOD(KFAM/100,10)
 496         ENDIF
 497       ENDIF
 498
 499 C...Kinematics of one-particle decays.
 500       IF(ND.EQ.1) THEN
 501         DO 400 J=1,4
 502         P(N+1,J)=P(IP,J)
 503   400   CONTINUE
 504         GOTO 660
 505       ENDIF
 506
 507 C...Calculate maximum weight ND-particle decay.
 508       PV(ND,5)=P(N+ND,5)
 509       IF(ND.GE.3) THEN
 510         WTMAX=1./WTCOR(ND-2)
 511         PMAX=PV(1,5)-PS+P(N+ND,5)
 512         PMIN=0.
 513         DO 410 IL=ND-1,1,-1
 514         PMAX=PMAX+P(N+IL,5)
 515         PMIN=PMIN+P(N+IL+1,5)
 516         WTMAX=WTMAX*PAWT(PMAX,PMIN,P(N+IL,5))
 517   410   CONTINUE
 518       ENDIF
 519
 520 C...Find virtual gamma mass in Dalitz decay.
 521   420 IF(ND.EQ.2) THEN
 522       ELSEIF(MMAT.EQ.2) THEN
 523         PMES=4.*PMAS(11,1)**2
 524         PMRHO2=PMAS(131,1)**2
 525         PGRHO2=PMAS(131,2)**2
 526   430   PMST=PMES*(P(IP,5)**2/PMES)**RLU(0)
 527         WT=(1+0.5*PMES/PMST)*SQRT(MAX(0.,1.-PMES/PMST))*
 528      &  (1.-PMST/P(IP,5)**2)**3*(1.+PGRHO2/PMRHO2)/
 529      &  ((1.-PMST/PMRHO2)**2+PGRHO2/PMRHO2)
 530         IF(WT.LT.RLU(0)) GOTO 430
 531         PV(2,5)=MAX(2.00001*PMAS(11,1),SQRT(PMST))
 532
 533 C...M-generator gives weight. If rejected, try again.
 534       ELSE
 535   440   RORD(1)=1.
 536         DO 470 IL1=2,ND-1
 537         RSAV=RLU(0)
 538         DO 450 IL2=IL1-1,1,-1
 539         IF(RSAV.LE.RORD(IL2)) GOTO 460
 540         RORD(IL2+1)=RORD(IL2)
 541   450   CONTINUE
 542   460   RORD(IL2+1)=RSAV
 543   470   CONTINUE
 544         RORD(ND)=0.
 545         WT=1.
 546         DO 480 IL=ND-1,1,-1
 547         PV(IL,5)=PV(IL+1,5)+P(N+IL,5)+(RORD(IL)-RORD(IL+1))*(PV(1,5)-PS)
 548         WT=WT*PAWT(PV(IL,5),PV(IL+1,5),P(N+IL,5))
 549   480   CONTINUE
 550         IF(WT.LT.RLU(0)*WTMAX) GOTO 440
 551       ENDIF
 552
 553 C...Perform two-particle decays in respective CM frame.
 554   490 DO 510 IL=1,ND-1
 555       PA=PAWT(PV(IL,5),PV(IL+1,5),P(N+IL,5))
 556       UE(3)=2.*RLU(0)-1.
 557       PHI=PARU(2)*RLU(0)
 558       UE(1)=SQRT(1.-UE(3)**2)*COS(PHI)
 559       UE(2)=SQRT(1.-UE(3)**2)*SIN(PHI)
 560       DO 500 J=1,3
 561       P(N+IL,J)=PA*UE(J)
 562       PV(IL+1,J)=-PA*UE(J)
 563   500 CONTINUE
 564       P(N+IL,4)=SQRT(PA**2+P(N+IL,5)**2)
 565       PV(IL+1,4)=SQRT(PA**2+PV(IL+1,5)**2)
 566   510 CONTINUE
 567
 568 C...Lorentz transform decay products to lab frame.
 569       DO 520 J=1,4
 570       P(N+ND,J)=PV(ND,J)
 571   520 CONTINUE
 572       DO 560 IL=ND-1,1,-1
 573       DO 530 J=1,3
 574       BE(J)=PV(IL,J)/PV(IL,4)
 575   530 CONTINUE
 576       GA=PV(IL,4)/PV(IL,5)
 577       DO 550 I=N+IL,N+ND
 578       BEP=BE(1)*P(I,1)+BE(2)*P(I,2)+BE(3)*P(I,3)
 579       DO 540 J=1,3
 580       P(I,J)=P(I,J)+GA*(GA*BEP/(1.+GA)+P(I,4))*BE(J)
 581   540 CONTINUE
 582       P(I,4)=GA*(P(I,4)+BEP)
 583   550 CONTINUE
 584   560 CONTINUE
 585
 586 C...Check that no infinite loop in matrix element weight.
 587       NTRY=NTRY+1
 588       IF(NTRY.GT.800) GOTO 590
 589
 590 C...Matrix elements for omega and phi decays.
 591       IF(MMAT.EQ.1) THEN
 592         WT=(P(N+1,5)*P(N+2,5)*P(N+3,5))**2-(P(N+1,5)*FOUR(N+2,N+3))**2
 593      &  -(P(N+2,5)*FOUR(N+1,N+3))**2-(P(N+3,5)*FOUR(N+1,N+2))**2
 594      &  +2.*FOUR(N+1,N+2)*FOUR(N+1,N+3)*FOUR(N+2,N+3)
 595         IF(MAX(WT*WTCOR(9)/P(IP,5)**6,0.001).LT.RLU(0)) GOTO 420
 596
 597 C...Matrix elements for pi0 or eta Dalitz decay to gamma e+ e-.
 598       ELSEIF(MMAT.EQ.2) THEN
 599         FOUR12=FOUR(N+1,N+2)
 600         FOUR13=FOUR(N+1,N+3)
 601         WT=(PMST-0.5*PMES)*(FOUR12**2+FOUR13**2)+
 602      &  PMES*(FOUR12*FOUR13+FOUR12**2+FOUR13**2)
 603         IF(WT.LT.RLU(0)*0.25*PMST*(P(IP,5)**2-PMST)**2) GOTO 490
 604
 605 C...Matrix element for S0 -> S1 + V1 -> S1 + S2 + S3 (S scalar,
 606 C...V vector), of form cos**2(theta02) in V1 rest frame, and for
 607 C...S0 -> gamma + V1 -> gamma + S2 + S3, of form sin**2(theta02).
 608       ELSEIF(MMAT.EQ.3.AND.NM.EQ.2) THEN
 609         FOUR10=FOUR(IP,IM)
 610         FOUR12=FOUR(IP,N+1)
 611         FOUR02=FOUR(IM,N+1)
 612         PMS1=P(IP,5)**2
 613         PMS0=P(IM,5)**2
 614         PMS2=P(N+1,5)**2
 615         IF(KFAS.NE.22) HNUM=(FOUR10*FOUR12-PMS1*FOUR02)**2
 616         IF(KFAS.EQ.22) HNUM=PMS1*(2.*FOUR10*FOUR12*FOUR02-
 617      &  PMS1*FOUR02**2-PMS0*FOUR12**2-PMS2*FOUR10**2+PMS1*PMS0*PMS2)
 618         HNUM=MAX(1E-6*PMS1**2*PMS0*PMS2,HNUM)
 619         HDEN=(FOUR10**2-PMS1*PMS0)*(FOUR12**2-PMS1*PMS2)
 620         IF(HNUM.LT.RLU(0)*HDEN) GOTO 490
 621
 622 C...Matrix element for "onium" -> g + g + g or gamma + g + g.
 623       ELSEIF(MMAT.EQ.4) THEN
 624         HX1=2.*FOUR(IP,N+1)/P(IP,5)**2
 625         HX2=2.*FOUR(IP,N+2)/P(IP,5)**2
 626         HX3=2.*FOUR(IP,N+3)/P(IP,5)**2
 627         WT=((1.-HX1)/(HX2*HX3))**2+((1.-HX2)/(HX1*HX3))**2+
 628      &  ((1.-HX3)/(HX1*HX2))**2
 629         IF(WT.LT.2.*RLU(0)) GOTO 420
 630         IF(K(IP+1,2).EQ.22.AND.(1.-HX1)*P(IP,5)**2.LT.4.*PARJ(32)**2)
 631      &  GOTO 420
 632
 633 C...Effective matrix element for nu spectrum in tau -> nu + hadrons.
 634       ELSEIF(MMAT.EQ.41) THEN
 635         HX1=2.*FOUR(IP,N+1)/P(IP,5)**2
 636         HXM=MIN(0.75,2.*(1.-PS/P(IP,5)))
 637         IF(HX1*(3.-2.*HX1).LT.RLU(0)*HXM*(3.-2.*HXM)) GOTO 420
 638
 639 C...Matrix elements for weak decays (only semileptonic for c and b)
 640       ELSEIF((MMAT.EQ.42.OR.MMAT.EQ.43.OR.MMAT.EQ.44.OR.MMAT.EQ.48)
 641      &.AND.ND.EQ.3) THEN
 642         IF(MBST.EQ.0) WT=FOUR(IP,N+1)*FOUR(N+2,N+3)
 643         IF(MBST.EQ.1) WT=P(IP,5)*P(N+1,4)*FOUR(N+2,N+3)
 644         IF(WT.LT.RLU(0)*P(IP,5)*PV(1,5)**3/WTCOR(10)) GOTO 420
 645       ELSEIF(MMAT.EQ.42.OR.MMAT.EQ.43.OR.MMAT.EQ.44.OR.MMAT.EQ.48) THEN
 646         DO 580 J=1,4
 647         P(N+NP+1,J)=0.
 648         DO 570 IS=N+3,N+NP
 649         P(N+NP+1,J)=P(N+NP+1,J)+P(IS,J)
 650   570   CONTINUE
 651   580   CONTINUE
 652         IF(MBST.EQ.0) WT=FOUR(IP,N+1)*FOUR(N+2,N+NP+1)
 653         IF(MBST.EQ.1) WT=P(IP,5)*P(N+1,4)*FOUR(N+2,N+NP+1)
 654         IF(WT.LT.RLU(0)*P(IP,5)*PV(1,5)**3/WTCOR(10)) GOTO 420
 655
 656 C...Angular distribution in W decay.
 657       ELSEIF(MMAT.EQ.46.AND.MSGN.NE.0) THEN
 658         IF(MSGN.GT.0) WT=FOUR(IM,N+1)*FOUR(N+2,IP+1)
 659         IF(MSGN.LT.0) WT=FOUR(IM,N+2)*FOUR(N+1,IP+1)
 660         IF(WT.LT.RLU(0)*P(IM,5)**4/WTCOR(10)) GOTO 490
 661       ENDIF
 662
 663 C...Scale back energy and reattach spectator.
 664   590 IF(MREM.EQ.1) THEN
 665         DO 600 J=1,5
 666         PV(1,J)=PV(1,J)/(1.-PQT)
 667   600   CONTINUE
 668         ND=ND+1
 669         MREM=0
 670       ENDIF
 671
 672 C...Low invariant mass for system with spectator quark gives particle,
 673 C...not two jets. Readjust momenta accordingly.
 674       IF((MMAT.EQ.31.OR.MMAT.EQ.45).AND.ND.EQ.3) THEN
 675         MSTJ(93)=1
 676         PM2=ULMASS(K(N+2,2))
 677         MSTJ(93)=1
 678         PM3=ULMASS(K(N+3,2))
 679         IF(P(N+2,5)**2+P(N+3,5)**2+2.*FOUR(N+2,N+3).GE.
 680      &  (PARJ(32)+PM2+PM3)**2) GOTO 660
 681         K(N+2,1)=1
 682         KFTEMP=K(N+2,2)
 683         CALL LUKFDI(KFTEMP,K(N+3,2),KFLDMP,K(N+2,2))
 684         IF(K(N+2,2).EQ.0) GOTO 260
 685         P(N+2,5)=ULMASS(K(N+2,2))
 686         PS=P(N+1,5)+P(N+2,5)
 687         PV(2,5)=P(N+2,5)
 688         MMAT=0
 689         ND=2
 690         GOTO 490
 691       ELSEIF(MMAT.EQ.44) THEN
 692         MSTJ(93)=1
 693         PM3=ULMASS(K(N+3,2))
 694         MSTJ(93)=1
 695         PM4=ULMASS(K(N+4,2))
 696         IF(P(N+3,5)**2+P(N+4,5)**2+2.*FOUR(N+3,N+4).GE.
 697      &  (PARJ(32)+PM3+PM4)**2) GOTO 630
 698         K(N+3,1)=1
 699         KFTEMP=K(N+3,2)
 700         CALL LUKFDI(KFTEMP,K(N+4,2),KFLDMP,K(N+3,2))
 701         IF(K(N+3,2).EQ.0) GOTO 260
 702         P(N+3,5)=ULMASS(K(N+3,2))
 703         DO 610 J=1,3
 704         P(N+3,J)=P(N+3,J)+P(N+4,J)
 705   610   CONTINUE
 706         P(N+3,4)=SQRT(P(N+3,1)**2+P(N+3,2)**2+P(N+3,3)**2+P(N+3,5)**2)
 707         HA=P(N+1,4)**2-P(N+2,4)**2
 708         HB=HA-(P(N+1,5)**2-P(N+2,5)**2)
 709         HC=(P(N+1,1)-P(N+2,1))**2+(P(N+1,2)-P(N+2,2))**2+
 710      &  (P(N+1,3)-P(N+2,3))**2
 711         HD=(PV(1,4)-P(N+3,4))**2
 712         HE=HA**2-2.*HD*(P(N+1,4)**2+P(N+2,4)**2)+HD**2
 713         HF=HD*HC-HB**2
 714         HG=HD*HC-HA*HB
 715         HH=(SQRT(HG**2+HE*HF)-HG)/(2.*HF)
 716         DO 620 J=1,3
 717         PCOR=HH*(P(N+1,J)-P(N+2,J))
 718         P(N+1,J)=P(N+1,J)+PCOR
 719         P(N+2,J)=P(N+2,J)-PCOR
 720   620   CONTINUE
 721         P(N+1,4)=SQRT(P(N+1,1)**2+P(N+1,2)**2+P(N+1,3)**2+P(N+1,5)**2)
 722         P(N+2,4)=SQRT(P(N+2,1)**2+P(N+2,2)**2+P(N+2,3)**2+P(N+2,5)**2)
 723         ND=ND-1
 724       ENDIF
 725
 726 C...Check invariant mass of W jets. May give one particle or start over.
 727   630 IF((MMAT.EQ.42.OR.MMAT.EQ.43.OR.MMAT.EQ.44.OR.MMAT.EQ.48)
 728      &.AND.IABS(K(N+1,2)).LT.10) THEN
 729         PMR=SQRT(MAX(0.,P(N+1,5)**2+P(N+2,5)**2+2.*FOUR(N+1,N+2)))
 730         MSTJ(93)=1
 731         PM1=ULMASS(K(N+1,2))
 732         MSTJ(93)=1
 733         PM2=ULMASS(K(N+2,2))
 734         IF(PMR.GT.PARJ(32)+PM1+PM2) GOTO 640
 735         KFLDUM=INT(1.5+RLU(0))
 736         CALL LUKFDI(K(N+1,2),-ISIGN(KFLDUM,K(N+1,2)),KFLDMP,KF1)
 737         CALL LUKFDI(K(N+2,2),-ISIGN(KFLDUM,K(N+2,2)),KFLDMP,KF2)
 738         IF(KF1.EQ.0.OR.KF2.EQ.0) GOTO 260
 739         PSM=ULMASS(KF1)+ULMASS(KF2)
 740         IF((MMAT.EQ.42.OR.MMAT.EQ.48).AND.PMR.GT.PARJ(64)+PSM) GOTO 640
 741         IF(MMAT.GE.43.AND.PMR.GT.0.2*PARJ(32)+PSM) GOTO 640
 742         IF(MMAT.EQ.48) GOTO 420
 743         IF(ND.EQ.4.OR.KFA.EQ.15) GOTO 260
 744         K(N+1,1)=1
 745         KFTEMP=K(N+1,2)
 746         CALL LUKFDI(KFTEMP,K(N+2,2),KFLDMP,K(N+1,2))
 747         IF(K(N+1,2).EQ.0) GOTO 260
 748         P(N+1,5)=ULMASS(K(N+1,2))
 749         K(N+2,2)=K(N+3,2)
 750         P(N+2,5)=P(N+3,5)
 751         PS=P(N+1,5)+P(N+2,5)
 752         IF(PS+PARJ(64).GT.PV(1,5)) GOTO 260
 753         PV(2,5)=P(N+3,5)
 754         MMAT=0
 755         ND=2
 756         GOTO 490
 757       ENDIF
 758
 759 C...Phase space decay of partons from W decay.
 760   640 IF((MMAT.EQ.42.OR.MMAT.EQ.48).AND.IABS(K(N+1,2)).LT.10) THEN
 761         KFLO(1)=K(N+1,2)
 762         KFLO(2)=K(N+2,2)
 763         K(N+1,1)=K(N+3,1)
 764         K(N+1,2)=K(N+3,2)
 765         DO 650 J=1,5
 766         PV(1,J)=P(N+1,J)+P(N+2,J)
 767         P(N+1,J)=P(N+3,J)
 768   650   CONTINUE
 769         PV(1,5)=PMR
 770         N=N+1
 771         NP=0
 772         NQ=2
 773         PS=0.
 774         MSTJ(93)=2
 775         PSQ=ULMASS(KFLO(1))
 776         MSTJ(93)=2
 777         PSQ=PSQ+ULMASS(KFLO(2))
 778         MMAT=11
 779         GOTO 290
 780       ENDIF
 781
 782 C...Boost back for rapidly moving particle.
 783   660 N=N+ND
 784       IF(MBST.EQ.1) THEN
 785         DO 670 J=1,3
 786         BE(J)=P(IP,J)/P(IP,4)
 787   670   CONTINUE
 788         GA=P(IP,4)/P(IP,5)
 789         DO 690 I=NSAV+1,N
 790         BEP=BE(1)*P(I,1)+BE(2)*P(I,2)+BE(3)*P(I,3)
 791         DO 680 J=1,3
 792         P(I,J)=P(I,J)+GA*(GA*BEP/(1.+GA)+P(I,4))*BE(J)
 793   680   CONTINUE
 794         P(I,4)=GA*(P(I,4)+BEP)
 795   690   CONTINUE
 796       ENDIF
 797
 798 C...Fill in position of decay vertex.
 799       DO 710 I=NSAV+1,N
 800       DO 700 J=1,4
 801       V(I,J)=VDCY(J)
 802   700 CONTINUE
 803       V(I,5)=0.
 804   710 CONTINUE
 805
 806 C...Set up for parton shower evolution from jets.
 807       IF(MSTJ(23).GE.1.AND.MMAT.EQ.4.AND.K(NSAV+1,2).EQ.21) THEN
 808         K(NSAV+1,1)=3
 809         K(NSAV+2,1)=3
 810         K(NSAV+3,1)=3
 811         K(NSAV+1,4)=MSTU(5)*(NSAV+2)
 812         K(NSAV+1,5)=MSTU(5)*(NSAV+3)
 813         K(NSAV+2,4)=MSTU(5)*(NSAV+3)
 814         K(NSAV+2,5)=MSTU(5)*(NSAV+1)
 815         K(NSAV+3,4)=MSTU(5)*(NSAV+1)
 816         K(NSAV+3,5)=MSTU(5)*(NSAV+2)
 817         MSTJ(92)=-(NSAV+1)
 818       ELSEIF(MSTJ(23).GE.1.AND.MMAT.EQ.4) THEN
 819         K(NSAV+2,1)=3
 820         K(NSAV+3,1)=3
 821         K(NSAV+2,4)=MSTU(5)*(NSAV+3)
 822         K(NSAV+2,5)=MSTU(5)*(NSAV+3)
 823         K(NSAV+3,4)=MSTU(5)*(NSAV+2)
 824         K(NSAV+3,5)=MSTU(5)*(NSAV+2)
 825         MSTJ(92)=NSAV+2
 826       ELSEIF(MSTJ(23).GE.1.AND.(MMAT.EQ.32.OR.MMAT.EQ.44.OR.MMAT.EQ.46)
 827      &.AND.IABS(K(NSAV+1,2)).LE.10.AND.IABS(K(NSAV+2,2)).LE.10) THEN
 828         K(NSAV+1,1)=3
 829         K(NSAV+2,1)=3
 830         K(NSAV+1,4)=MSTU(5)*(NSAV+2)
 831         K(NSAV+1,5)=MSTU(5)*(NSAV+2)
 832         K(NSAV+2,4)=MSTU(5)*(NSAV+1)
 833         K(NSAV+2,5)=MSTU(5)*(NSAV+1)
 834         MSTJ(92)=NSAV+1
 835       ELSEIF(MSTJ(23).GE.1.AND.(MMAT.EQ.32.OR.MMAT.EQ.44.OR.MMAT.EQ.46)
 836      &.AND.IABS(K(NSAV+1,2)).LE.20.AND.IABS(K(NSAV+2,2)).LE.20) THEN
 837         MSTJ(92)=NSAV+1
 838       ELSEIF(MSTJ(23).GE.1.AND.MMAT.EQ.33.AND.IABS(K(NSAV+2,2)).EQ.21)
 839      &THEN
 840         K(NSAV+1,1)=3
 841         K(NSAV+2,1)=3
 842         K(NSAV+3,1)=3
 843         KCP=LUCOMP(K(NSAV+1,2))
 844         KQP=KCHG(KCP,2)*ISIGN(1,K(NSAV+1,2))
 845         JCON=4
 846         IF(KQP.LT.0) JCON=5
 847         K(NSAV+1,JCON)=MSTU(5)*(NSAV+2)
 848         K(NSAV+2,9-JCON)=MSTU(5)*(NSAV+1)
 849         K(NSAV+2,JCON)=MSTU(5)*(NSAV+3)
 850         K(NSAV+3,9-JCON)=MSTU(5)*(NSAV+2)
 851         MSTJ(92)=NSAV+1
 852       ELSEIF(MSTJ(23).GE.1.AND.MMAT.EQ.33) THEN
 853         K(NSAV+1,1)=3
 854         K(NSAV+3,1)=3
 855         K(NSAV+1,4)=MSTU(5)*(NSAV+3)
 856         K(NSAV+1,5)=MSTU(5)*(NSAV+3)
 857         K(NSAV+3,4)=MSTU(5)*(NSAV+1)
 858         K(NSAV+3,5)=MSTU(5)*(NSAV+1)
 859         MSTJ(92)=NSAV+1
 860
 861 C...Set up for parton shower evolution in t -> W + b.
 862       ELSEIF(MSTJ(27).GE.1.AND.MMAT.EQ.45.AND.ND.EQ.3) THEN
 863         K(NSAV+2,1)=3
 864         K(NSAV+3,1)=3
 865         K(NSAV+2,4)=MSTU(5)*(NSAV+3)
 866         K(NSAV+2,5)=MSTU(5)*(NSAV+3)
 867         K(NSAV+3,4)=MSTU(5)*(NSAV+2)
 868         K(NSAV+3,5)=MSTU(5)*(NSAV+2)
 869         MSTJ(92)=NSAV+1
 870       ENDIF
 871
 872 C...Mark decayed particle; special option for B-B~ mixing.
 873       IF(K(IP,1).EQ.5) K(IP,1)=15
 874       IF(K(IP,1).LE.10) K(IP,1)=11
 875       IF(MMIX.EQ.1.AND.MSTJ(26).EQ.2.AND.K(IP,1).EQ.11) K(IP,1)=12
 876       K(IP,4)=NSAV+1
 877       K(IP,5)=N
 878
 879       RETURN
 880       END