PHOS/shaker/ludecy.f

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