hijing/src/ludecy.f

0001
0002 C*********************************************************************
0003
0004       SUBROUTINE LUDECY(IP)
0005
0006 C...Purpose: to handle the decay of unstable particles.
0007       COMMON/LUJETS/N,K(9000,5),P(9000,5),V(9000,5)
0008       SAVE /LUJETS/
0009       COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0010       SAVE /LUDAT1/
0011       COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
0012       SAVE /LUDAT2/
0013       COMMON/LUDAT3/MDCY(500,3),MDME(2000,2),BRAT(2000),KFDP(2000,5)
0014       SAVE /LUDAT3/
0015       DIMENSION VDCY(4),KFLO(4),KFL1(4),PV(10,5),RORD(10),UE(3),BE(3),
0016      &WTCOR(10)
0017       DATA WTCOR/2.,5.,15.,60.,250.,1500.,1.2E4,1.2E5,150.,16./
0018
0019 C...Functions: momentum in two-particle decays, four-product and
0020 C...matrix element times phase space in weak decays.
0021       PAWT(A,B,C)=SQRT((A**2-(B+C)**2)*(A**2-(B-C)**2))/(2.*A)
0022       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)
0023       HMEPS(HA)=((1.-HRQ-HA)**2+3.*HA*(1.+HRQ-HA))*
0024      &SQRT((1.-HRQ-HA)**2-4.*HRQ*HA)
0025
0026 C...Initial values.
0027       NTRY=0
0028       NSAV=N
0029       KFA=IABS(K(IP,2))
0030       KFS=ISIGN(1,K(IP,2))
0031       KC=LUCOMP(KFA)
0032       MSTJ(92)=0
0033
0034 C...Choose lifetime and determine decay vertex.
0035       IF(K(IP,1).EQ.5) THEN
0036         V(IP,5)=0.
0037       ELSEIF(K(IP,1).NE.4) THEN
0038         V(IP,5)=-PMAS(KC,4)*LOG(RLU(0))
0039       ENDIF
0040       DO 100 J=1,4
0041   100 VDCY(J)=V(IP,J)+V(IP,5)*P(IP,J)/P(IP,5)
0042
0043 C...Determine whether decay allowed or not.
0044       MOUT=0
0045       IF(MSTJ(22).EQ.2) THEN
0046         IF(PMAS(KC,4).GT.PARJ(71)) MOUT=1
0047       ELSEIF(MSTJ(22).EQ.3) THEN
0048         IF(VDCY(1)**2+VDCY(2)**2+VDCY(3)**2.GT.PARJ(72)**2) MOUT=1
0049       ELSEIF(MSTJ(22).EQ.4) THEN
0050         IF(VDCY(1)**2+VDCY(2)**2.GT.PARJ(73)**2) MOUT=1
0051         IF(ABS(VDCY(3)).GT.PARJ(74)) MOUT=1
0052       ENDIF
0053       IF(MOUT.EQ.1.AND.K(IP,1).NE.5) THEN
0054         K(IP,1)=4
0055         RETURN
0056       ENDIF
0057
0058 C...Check existence of decay channels. Particle/antiparticle rules.
0059       KCA=KC
0060       IF(MDCY(KC,2).GT.0) THEN
0061         MDMDCY=MDME(MDCY(KC,2),2)
0062         IF(MDMDCY.GT.80.AND.MDMDCY.LE.90) KCA=MDMDCY
0063       ENDIF
0064       IF(MDCY(KCA,2).LE.0.OR.MDCY(KCA,3).LE.0) THEN
0065         CALL LUERRM(9,'(LUDECY:) no decay channel defined')
0066         RETURN
0067       ENDIF
0068       IF(MOD(KFA/1000,10).EQ.0.AND.(KCA.EQ.85.OR.KCA.EQ.87)) KFS=-KFS
0069       IF(KCHG(KC,3).EQ.0) THEN
0070         KFSP=1
0071         KFSN=0
0072         IF(RLU(0).GT.0.5) KFS=-KFS
0073       ELSEIF(KFS.GT.0) THEN
0074         KFSP=1
0075         KFSN=0
0076       ELSE
0077         KFSP=0
0078         KFSN=1
0079       ENDIF
0080
0081 C...Sum branching ratios of allowed decay channels.
0082   110 NOPE=0
0083       BRSU=0.
0084       DO 120 IDL=MDCY(KCA,2),MDCY(KCA,2)+MDCY(KCA,3)-1
0085       IF(MDME(IDL,1).NE.1.AND.KFSP*MDME(IDL,1).NE.2.AND.
0086      &KFSN*MDME(IDL,1).NE.3) GOTO 120
0087       IF(MDME(IDL,2).GT.100) GOTO 120
0088       NOPE=NOPE+1
0089       BRSU=BRSU+BRAT(IDL)
0090   120 CONTINUE
0091       IF(NOPE.EQ.0) THEN
0092         CALL LUERRM(2,'(LUDECY:) all decay channels closed by user')
0093         RETURN
0094       ENDIF
0095
0096 C...Select decay channel among allowed ones.
0097   130 RBR=BRSU*RLU(0)
0098       IDL=MDCY(KCA,2)-1
0099   140 IDL=IDL+1
0100       IF(MDME(IDL,1).NE.1.AND.KFSP*MDME(IDL,1).NE.2.AND.
0101      &KFSN*MDME(IDL,1).NE.3) THEN
0102         IF(IDL.LT.MDCY(KCA,2)+MDCY(KCA,3)-1) GOTO 140
0103       ELSEIF(MDME(IDL,2).GT.100) THEN
0104         IF(IDL.LT.MDCY(KCA,2)+MDCY(KCA,3)-1) GOTO 140
0105       ELSE
0106         IDC=IDL
0107         RBR=RBR-BRAT(IDL)
0108         IF(IDL.LT.MDCY(KCA,2)+MDCY(KCA,3)-1.AND.RBR.GT.0.) GOTO 140
0109       ENDIF
0110
0111 C...Start readout of decay channel: matrix element, reset counters.
0112       MMAT=MDME(IDC,2)
0113   150 NTRY=NTRY+1
0114       IF(NTRY.GT.1000) THEN
0115         CALL LUERRM(14,'(LUDECY:) caught in infinite loop')
0116         IF(MSTU(21).GE.1) RETURN
0117       ENDIF
0118       I=N
0119       NP=0
0120       NQ=0
0121       MBST=0
0122       IF(MMAT.GE.11.AND.MMAT.NE.46.AND.P(IP,4).GT.20.*P(IP,5)) MBST=1
0123       DO 160 J=1,4
0124       PV(1,J)=0.
0125   160 IF(MBST.EQ.0) PV(1,J)=P(IP,J)
0126       IF(MBST.EQ.1) PV(1,4)=P(IP,5)
0127       PV(1,5)=P(IP,5)
0128       PS=0.
0129       PSQ=0.
0130       MREM=0
0131
0132 C...Read out decay products. Convert to standard flavour code.
0133       JTMAX=5
0134       IF(MDME(IDC+1,2).EQ.101) JTMAX=10
0135       DO 170 JT=1,JTMAX
0136       IF(JT.LE.5) KP=KFDP(IDC,JT)
0137       IF(JT.GE.6) KP=KFDP(IDC+1,JT-5)
0138       IF(KP.EQ.0) GOTO 170
0139       KPA=IABS(KP)
0140       KCP=LUCOMP(KPA)
0141       IF(KCHG(KCP,3).EQ.0.AND.KPA.NE.81.AND.KPA.NE.82) THEN
0142         KFP=KP
0143       ELSEIF(KPA.NE.81.AND.KPA.NE.82) THEN
0144         KFP=KFS*KP
0145       ELSEIF(KPA.EQ.81.AND.MOD(KFA/1000,10).EQ.0) THEN
0146         KFP=-KFS*MOD(KFA/10,10)
0147       ELSEIF(KPA.EQ.81.AND.MOD(KFA/100,10).GE.MOD(KFA/10,10)) THEN
0148         KFP=KFS*(100*MOD(KFA/10,100)+3)
0149       ELSEIF(KPA.EQ.81) THEN
0150         KFP=KFS*(1000*MOD(KFA/10,10)+100*MOD(KFA/100,10)+1)
0151       ELSEIF(KP.EQ.82) THEN
0152         CALL LUKFDI(-KFS*INT(1.+(2.+PARJ(2))*RLU(0)),0,KFP,KDUMP)
0153         IF(KFP.EQ.0) GOTO 150
0154         MSTJ(93)=1
0155         IF(PV(1,5).LT.PARJ(32)+2.*ULMASS(KFP)) GOTO 150
0156       ELSEIF(KP.EQ.-82) THEN
0157         KFP=-KFP
0158         IF(IABS(KFP).GT.10) KFP=KFP+ISIGN(10000,KFP)
0159       ENDIF
0160       IF(KPA.EQ.81.OR.KPA.EQ.82) KCP=LUCOMP(KFP)
0161
0162 C...Add decay product to event record or to quark flavour list.
0163       KFPA=IABS(KFP)
0164       KQP=KCHG(KCP,2)
0165       IF(MMAT.GE.11.AND.MMAT.LE.30.AND.KQP.NE.0) THEN
0166         NQ=NQ+1
0167         KFLO(NQ)=KFP
0168         MSTJ(93)=2
0169         PSQ=PSQ+ULMASS(KFLO(NQ))
0170       ELSEIF(MMAT.GE.42.AND.MMAT.LE.43.AND.NP.EQ.3.AND.MOD(NQ,2).EQ.1)
0171      &THEN
0172         NQ=NQ-1
0173         PS=PS-P(I,5)
0174         K(I,1)=1
0175         KFI=K(I,2)
0176         CALL LUKFDI(KFP,KFI,KFLDMP,K(I,2))
0177         IF(K(I,2).EQ.0) GOTO 150
0178         MSTJ(93)=1
0179         P(I,5)=ULMASS(K(I,2))
0180         PS=PS+P(I,5)
0181       ELSE
0182         I=I+1
0183         NP=NP+1
0184         IF(MMAT.NE.33.AND.KQP.NE.0) NQ=NQ+1
0185         IF(MMAT.EQ.33.AND.KQP.NE.0.AND.KQP.NE.2) NQ=NQ+1
0186         K(I,1)=1+MOD(NQ,2)
0187         IF(MMAT.EQ.4.AND.JT.LE.2.AND.KFP.EQ.21) K(I,1)=2
0188         IF(MMAT.EQ.4.AND.JT.EQ.3) K(I,1)=1
0189         K(I,2)=KFP
0190         K(I,3)=IP
0191         K(I,4)=0
0192         K(I,5)=0
0193         P(I,5)=ULMASS(KFP)
0194         IF(MMAT.EQ.45.AND.KFPA.EQ.89) P(I,5)=PARJ(32)
0195         PS=PS+P(I,5)
0196       ENDIF
0197   170 CONTINUE
0198
0199 C...Choose decay multiplicity in phase space model.
0200   180 IF(MMAT.GE.11.AND.MMAT.LE.30) THEN
0201         PSP=PS
0202         CNDE=PARJ(61)*LOG(MAX((PV(1,5)-PS-PSQ)/PARJ(62),1.1))
0203         IF(MMAT.EQ.12) CNDE=CNDE+PARJ(63)
0204   190   NTRY=NTRY+1
0205         IF(NTRY.GT.1000) THEN
0206           CALL LUERRM(14,'(LUDECY:) caught in infinite loop')
0207           IF(MSTU(21).GE.1) RETURN
0208         ENDIF
0209         IF(MMAT.LE.20) THEN
0210           GAUSS=SQRT(-2.*CNDE*LOG(MAX(1E-10,RLU(0))))*
0211      &    SIN(PARU(2)*RLU(0))
0212           ND=0.5+0.5*NP+0.25*NQ+CNDE+GAUSS
0213           IF(ND.LT.NP+NQ/2.OR.ND.LT.2.OR.ND.GT.10) GOTO 190
0214           IF(MMAT.EQ.13.AND.ND.EQ.2) GOTO 190
0215           IF(MMAT.EQ.14.AND.ND.LE.3) GOTO 190
0216           IF(MMAT.EQ.15.AND.ND.LE.4) GOTO 190
0217         ELSE
0218           ND=MMAT-20
0219         ENDIF
0220
0221 C...Form hadrons from flavour content.
0222         DO 200 JT=1,4
0223   200   KFL1(JT)=KFLO(JT)
0224         IF(ND.EQ.NP+NQ/2) GOTO 220
0225         DO 210 I=N+NP+1,N+ND-NQ/2
0226         JT=1+INT((NQ-1)*RLU(0))
0227         CALL LUKFDI(KFL1(JT),0,KFL2,K(I,2))
0228         IF(K(I,2).EQ.0) GOTO 190
0229   210   KFL1(JT)=-KFL2
0230   220   JT=2
0231         JT2=3
0232         JT3=4
0233         IF(NQ.EQ.4.AND.RLU(0).LT.PARJ(66)) JT=4
0234         IF(JT.EQ.4.AND.ISIGN(1,KFL1(1)*(10-IABS(KFL1(1))))*
0235      &  ISIGN(1,KFL1(JT)*(10-IABS(KFL1(JT)))).GT.0) JT=3
0236         IF(JT.EQ.3) JT2=2
0237         IF(JT.EQ.4) JT3=2
0238         CALL LUKFDI(KFL1(1),KFL1(JT),KFLDMP,K(N+ND-NQ/2+1,2))
0239         IF(K(N+ND-NQ/2+1,2).EQ.0) GOTO 190
0240         IF(NQ.EQ.4) CALL LUKFDI(KFL1(JT2),KFL1(JT3),KFLDMP,K(N+ND,2))
0241         IF(NQ.EQ.4.AND.K(N+ND,2).EQ.0) GOTO 190
0242
0243 C...Check that sum of decay product masses not too large.
0244         PS=PSP
0245         DO 230 I=N+NP+1,N+ND
0246         K(I,1)=1
0247         K(I,3)=IP
0248         K(I,4)=0
0249         K(I,5)=0
0250         P(I,5)=ULMASS(K(I,2))
0251   230   PS=PS+P(I,5)
0252         IF(PS+PARJ(64).GT.PV(1,5)) GOTO 190
0253
0254 C...Rescale energy to subtract off spectator quark mass.
0255       ELSEIF((MMAT.EQ.31.OR.MMAT.EQ.33.OR.MMAT.EQ.44.OR.MMAT.EQ.45).
0256      &AND.NP.GE.3) THEN
0257         PS=PS-P(N+NP,5)
0258         PQT=(P(N+NP,5)+PARJ(65))/PV(1,5)
0259         DO 240 J=1,5
0260         P(N+NP,J)=PQT*PV(1,J)
0261   240   PV(1,J)=(1.-PQT)*PV(1,J)
0262         IF(PS+PARJ(64).GT.PV(1,5)) GOTO 150
0263         ND=NP-1
0264         MREM=1
0265
0266 C...Phase space factors imposed in W decay.
0267       ELSEIF(MMAT.EQ.46) THEN
0268         MSTJ(93)=1
0269         PSMC=ULMASS(K(N+1,2))
0270         MSTJ(93)=1
0271         PSMC=PSMC+ULMASS(K(N+2,2))
0272         IF(MAX(PS,PSMC)+PARJ(32).GT.PV(1,5)) GOTO 130
0273         HR1=(P(N+1,5)/PV(1,5))**2
0274         HR2=(P(N+2,5)/PV(1,5))**2
0275         IF((1.-HR1-HR2)*(2.+HR1+HR2)*SQRT((1.-HR1-HR2)**2-4.*HR1*HR2).
0276      &  LT.2.*RLU(0)) GOTO 130
0277         ND=NP
0278
0279 C...Fully specified final state: check mass broadening effects.
0280       ELSE
0281         IF(NP.GE.2.AND.PS+PARJ(64).GT.PV(1,5)) GOTO 150
0282         ND=NP
0283       ENDIF
0284
0285 C...Select W mass in decay Q -> W + q, without W propagator.
0286       IF(MMAT.EQ.45.AND.MSTJ(25).LE.0) THEN
0287         HLQ=(PARJ(32)/PV(1,5))**2
0288         HUQ=(1.-(P(N+2,5)+PARJ(64))/PV(1,5))**2
0289         HRQ=(P(N+2,5)/PV(1,5))**2
0290   250   HW=HLQ+RLU(0)*(HUQ-HLQ)
0291         IF(HMEPS(HW).LT.RLU(0)) GOTO 250
0292         P(N+1,5)=PV(1,5)*SQRT(HW)
0293
0294 C...Ditto, including W propagator. Divide mass range into three regions.
0295       ELSEIF(MMAT.EQ.45) THEN
0296         HQW=(PV(1,5)/PMAS(24,1))**2
0297         HLW=(PARJ(32)/PMAS(24,1))**2
0298         HUW=((PV(1,5)-P(N+2,5)-PARJ(64))/PMAS(24,1))**2
0299         HRQ=(P(N+2,5)/PV(1,5))**2
0300         HG=PMAS(24,2)/PMAS(24,1)
0301         HATL=ATAN((HLW-1.)/HG)
0302         HM=MIN(1.,HUW-0.001)
0303         HMV1=HMEPS(HM/HQW)/((HM-1.)**2+HG**2)
0304   260   HM=HM-HG
0305         HMV2=HMEPS(HM/HQW)/((HM-1.)**2+HG**2)
0306         HSAV1=HMEPS(HM/HQW)
0307         HSAV2=1./((HM-1.)**2+HG**2)
0308         IF(HMV2.GT.HMV1.AND.HM-HG.GT.HLW) THEN
0309           HMV1=HMV2
0310           GOTO 260
0311         ENDIF
0312         HMV=MIN(2.*HMV1,HMEPS(HM/HQW)/HG**2)
0313         HM1=1.-SQRT(1./HMV-HG**2)
0314         IF(HM1.GT.HLW.AND.HM1.LT.HM) THEN
0315           HM=HM1
0316         ELSEIF(HMV2.LE.HMV1) THEN
0317           HM=MAX(HLW,HM-MIN(0.1,1.-HM))
0318         ENDIF
0319         HATM=ATAN((HM-1.)/HG)
0320         HWT1=(HATM-HATL)/HG
0321         HWT2=HMV*(MIN(1.,HUW)-HM)
0322         HWT3=0.
0323         IF(HUW.GT.1.) THEN
0324           HATU=ATAN((HUW-1.)/HG)
0325           HMP1=HMEPS(1./HQW)
0326           HWT3=HMP1*HATU/HG
0327         ENDIF
0328
0329 C...Select mass region and W mass there. Accept according to weight.
0330   270   HREG=RLU(0)*(HWT1+HWT2+HWT3)
0331         IF(HREG.LE.HWT1) THEN
0332           HW=1.+HG*TAN(HATL+RLU(0)*(HATM-HATL))
0333           HACC=HMEPS(HW/HQW)
0334         ELSEIF(HREG.LE.HWT1+HWT2) THEN
0335           HW=HM+RLU(0)*(MIN(1.,HUW)-HM)
0336           HACC=HMEPS(HW/HQW)/((HW-1.)**2+HG**2)/HMV
0337         ELSE
0338           HW=1.+HG*TAN(RLU(0)*HATU)
0339           HACC=HMEPS(HW/HQW)/HMP1
0340         ENDIF
0341         IF(HACC.LT.RLU(0)) GOTO 270
0342         P(N+1,5)=PMAS(24,1)*SQRT(HW)
0343       ENDIF
0344
0345 C...Determine position of grandmother, number of sisters, Q -> W sign.
0346       NM=0
0347       MSGN=0
0348       IF(MMAT.EQ.3.OR.MMAT.EQ.46) THEN
0349         IM=K(IP,3)
0350         IF(IM.LT.0.OR.IM.GE.IP) IM=0
0351         IF(IM.NE.0) KFAM=IABS(K(IM,2))
0352         IF(IM.NE.0.AND.MMAT.EQ.3) THEN
0353           DO 280 IL=MAX(IP-2,IM+1),MIN(IP+2,N)
0354   280     IF(K(IL,3).EQ.IM) NM=NM+1
0355           IF(NM.NE.2.OR.KFAM.LE.100.OR.MOD(KFAM,10).NE.1.OR.
0356      &    MOD(KFAM/1000,10).NE.0) NM=0
0357         ELSEIF(IM.NE.0.AND.MMAT.EQ.46) THEN
0358           MSGN=ISIGN(1,K(IM,2)*K(IP,2))
0359           IF(KFAM.GT.100.AND.MOD(KFAM/1000,10).EQ.0) MSGN=
0360      &    MSGN*(-1)**MOD(KFAM/100,10)
0361         ENDIF
0362       ENDIF
0363
0364 C...Kinematics of one-particle decays.
0365       IF(ND.EQ.1) THEN
0366         DO 290 J=1,4
0367   290   P(N+1,J)=P(IP,J)
0368         GOTO 510
0369       ENDIF
0370
0371 C...Calculate maximum weight ND-particle decay.
0372       PV(ND,5)=P(N+ND,5)
0373       IF(ND.GE.3) THEN
0374         WTMAX=1./WTCOR(ND-2)
0375         PMAX=PV(1,5)-PS+P(N+ND,5)
0376         PMIN=0.
0377         DO 300 IL=ND-1,1,-1
0378         PMAX=PMAX+P(N+IL,5)
0379         PMIN=PMIN+P(N+IL+1,5)
0380   300   WTMAX=WTMAX*PAWT(PMAX,PMIN,P(N+IL,5))
0381       ENDIF
0382
0383 C...Find virtual gamma mass in Dalitz decay.
0384   310 IF(ND.EQ.2) THEN
0385       ELSEIF(MMAT.EQ.2) THEN
0386         PMES=4.*PMAS(11,1)**2
0387         PMRHO2=PMAS(131,1)**2
0388         PGRHO2=PMAS(131,2)**2
0389   320   PMST=PMES*(P(IP,5)**2/PMES)**RLU(0)
0390         WT=(1+0.5*PMES/PMST)*SQRT(MAX(0.,1.-PMES/PMST))*
0391      &  (1.-PMST/P(IP,5)**2)**3*(1.+PGRHO2/PMRHO2)/
0392      &  ((1.-PMST/PMRHO2)**2+PGRHO2/PMRHO2)
0393         IF(WT.LT.RLU(0)) GOTO 320
0394         PV(2,5)=MAX(2.00001*PMAS(11,1),SQRT(PMST))
0395
0396 C...M-generator gives weight. If rejected, try again.
0397       ELSE
0398   330   RORD(1)=1.
0399         DO 350 IL1=2,ND-1
0400         RSAV=RLU(0)
0401         DO 340 IL2=IL1-1,1,-1
0402         IF(RSAV.LE.RORD(IL2)) GOTO 350
0403   340   RORD(IL2+1)=RORD(IL2)
0404   350   RORD(IL2+1)=RSAV
0405         RORD(ND)=0.
0406         WT=1.
0407         DO 360 IL=ND-1,1,-1
0408         PV(IL,5)=PV(IL+1,5)+P(N+IL,5)+(RORD(IL)-RORD(IL+1))*(PV(1,5)-PS)
0409   360   WT=WT*PAWT(PV(IL,5),PV(IL+1,5),P(N+IL,5))
0410         IF(WT.LT.RLU(0)*WTMAX) GOTO 330
0411       ENDIF
0412
0413 C...Perform two-particle decays in respective CM frame.
0414   370 DO 390 IL=1,ND-1
0415       PA=PAWT(PV(IL,5),PV(IL+1,5),P(N+IL,5))
0416       UE(3)=2.*RLU(0)-1.
0417       PHI=PARU(2)*RLU(0)
0418       UE(1)=SQRT(1.-UE(3)**2)*COS(PHI)
0419       UE(2)=SQRT(1.-UE(3)**2)*SIN(PHI)
0420       DO 380 J=1,3
0421       P(N+IL,J)=PA*UE(J)
0422   380 PV(IL+1,J)=-PA*UE(J)
0423       P(N+IL,4)=SQRT(PA**2+P(N+IL,5)**2)
0424   390 PV(IL+1,4)=SQRT(PA**2+PV(IL+1,5)**2)
0425
0426 C...Lorentz transform decay products to lab frame.
0427       DO 400 J=1,4
0428   400 P(N+ND,J)=PV(ND,J)
0429       DO 430 IL=ND-1,1,-1
0430       DO 410 J=1,3
0431   410 BE(J)=PV(IL,J)/PV(IL,4)
0432       GA=PV(IL,4)/PV(IL,5)
0433       DO 430 I=N+IL,N+ND
0434       BEP=BE(1)*P(I,1)+BE(2)*P(I,2)+BE(3)*P(I,3)
0435       DO 420 J=1,3
0436   420 P(I,J)=P(I,J)+GA*(GA*BEP/(1.+GA)+P(I,4))*BE(J)
0437   430 P(I,4)=GA*(P(I,4)+BEP)
0438
0439 C...Matrix elements for omega and phi decays.
0440       IF(MMAT.EQ.1) THEN
0441         WT=(P(N+1,5)*P(N+2,5)*P(N+3,5))**2-(P(N+1,5)*FOUR(N+2,N+3))**2
0442      &  -(P(N+2,5)*FOUR(N+1,N+3))**2-(P(N+3,5)*FOUR(N+1,N+2))**2
0443      &  +2.*FOUR(N+1,N+2)*FOUR(N+1,N+3)*FOUR(N+2,N+3)
0444         IF(MAX(WT*WTCOR(9)/P(IP,5)**6,0.001).LT.RLU(0)) GOTO 310
0445
0446 C...Matrix elements for pi0 or eta Dalitz decay to gamma e+ e-.
0447       ELSEIF(MMAT.EQ.2) THEN
0448         FOUR12=FOUR(N+1,N+2)
0449         FOUR13=FOUR(N+1,N+3)
0450         FOUR23=0.5*PMST-0.25*PMES
0451         WT=(PMST-0.5*PMES)*(FOUR12**2+FOUR13**2)+
0452      &  PMES*(FOUR12*FOUR13+FOUR12**2+FOUR13**2)
0453         IF(WT.LT.RLU(0)*0.25*PMST*(P(IP,5)**2-PMST)**2) GOTO 370
0454
0455 C...Matrix element for S0 -> S1 + V1 -> S1 + S2 + S3 (S scalar,
0456 C...V vector), of form cos**2(theta02) in V1 rest frame.
0457       ELSEIF(MMAT.EQ.3.AND.NM.EQ.2) THEN
0458         IF((P(IP,5)**2*FOUR(IM,N+1)-FOUR(IP,IM)*FOUR(IP,N+1))**2.LE.
0459      &  RLU(0)*(FOUR(IP,IM)**2-(P(IP,5)*P(IM,5))**2)*(FOUR(IP,N+1)**2-
0460      &  (P(IP,5)*P(N+1,5))**2)) GOTO 370
0461
0462 C...Matrix element for "onium" -> g + g + g or gamma + g + g.
0463       ELSEIF(MMAT.EQ.4) THEN
0464         HX1=2.*FOUR(IP,N+1)/P(IP,5)**2
0465         HX2=2.*FOUR(IP,N+2)/P(IP,5)**2
0466         HX3=2.*FOUR(IP,N+3)/P(IP,5)**2
0467         WT=((1.-HX1)/(HX2*HX3))**2+((1.-HX2)/(HX1*HX3))**2+
0468      &  ((1.-HX3)/(HX1*HX2))**2
0469         IF(WT.LT.2.*RLU(0)) GOTO 310
0470         IF(K(IP+1,2).EQ.22.AND.(1.-HX1)*P(IP,5)**2.LT.4.*PARJ(32)**2)
0471      &  GOTO 310
0472
0473 C...Effective matrix element for nu spectrum in tau -> nu + hadrons.
0474       ELSEIF(MMAT.EQ.41) THEN
0475         HX1=2.*FOUR(IP,N+1)/P(IP,5)**2
0476         IF(8.*HX1*(3.-2.*HX1)/9..LT.RLU(0)) GOTO 310
0477
0478 C...Matrix elements for weak decays (only semileptonic for c and b)
0479       ELSEIF(MMAT.GE.42.AND.MMAT.LE.44.AND.ND.EQ.3) THEN
0480         IF(MBST.EQ.0) WT=FOUR(IP,N+1)*FOUR(N+2,N+3)
0481         IF(MBST.EQ.1) WT=P(IP,5)*P(N+1,4)*FOUR(N+2,N+3)
0482         IF(WT.LT.RLU(0)*P(IP,5)*PV(1,5)**3/WTCOR(10)) GOTO 310
0483       ELSEIF(MMAT.GE.42.AND.MMAT.LE.44) THEN
0484         DO 440 J=1,4
0485         P(N+NP+1,J)=0.
0486         DO 440 IS=N+3,N+NP
0487   440   P(N+NP+1,J)=P(N+NP+1,J)+P(IS,J)
0488         IF(MBST.EQ.0) WT=FOUR(IP,N+1)*FOUR(N+2,N+NP+1)
0489         IF(MBST.EQ.1) WT=P(IP,5)*P(N+1,4)*FOUR(N+2,N+NP+1)
0490         IF(WT.LT.RLU(0)*P(IP,5)*PV(1,5)**3/WTCOR(10)) GOTO 310
0491
0492 C...Angular distribution in W decay.
0493       ELSEIF(MMAT.EQ.46.AND.MSGN.NE.0) THEN
0494         IF(MSGN.GT.0) WT=FOUR(IM,N+1)*FOUR(N+2,IP+1)
0495         IF(MSGN.LT.0) WT=FOUR(IM,N+2)*FOUR(N+1,IP+1)
0496         IF(WT.LT.RLU(0)*P(IM,5)**4/WTCOR(10)) GOTO 370
0497       ENDIF
0498
0499 C...Scale back energy and reattach spectator.
0500       IF(MREM.EQ.1) THEN
0501         DO 450 J=1,5
0502   450   PV(1,J)=PV(1,J)/(1.-PQT)
0503         ND=ND+1
0504         MREM=0
0505       ENDIF
0506
0507 C...Low invariant mass for system with spectator quark gives particle,
0508 C...not two jets. Readjust momenta accordingly.
0509       IF((MMAT.EQ.31.OR.MMAT.EQ.45).AND.ND.EQ.3) THEN
0510         MSTJ(93)=1
0511         PM2=ULMASS(K(N+2,2))
0512         MSTJ(93)=1
0513         PM3=ULMASS(K(N+3,2))
0514         IF(P(N+2,5)**2+P(N+3,5)**2+2.*FOUR(N+2,N+3).GE.
0515      &  (PARJ(32)+PM2+PM3)**2) GOTO 510
0516         K(N+2,1)=1
0517         KFTEMP=K(N+2,2)
0518         CALL LUKFDI(KFTEMP,K(N+3,2),KFLDMP,K(N+2,2))
0519         IF(K(N+2,2).EQ.0) GOTO 150
0520         P(N+2,5)=ULMASS(K(N+2,2))
0521         PS=P(N+1,5)+P(N+2,5)
0522         PV(2,5)=P(N+2,5)
0523         MMAT=0
0524         ND=2
0525         GOTO 370
0526       ELSEIF(MMAT.EQ.44) THEN
0527         MSTJ(93)=1
0528         PM3=ULMASS(K(N+3,2))
0529         MSTJ(93)=1
0530         PM4=ULMASS(K(N+4,2))
0531         IF(P(N+3,5)**2+P(N+4,5)**2+2.*FOUR(N+3,N+4).GE.
0532      &  (PARJ(32)+PM3+PM4)**2) GOTO 480
0533         K(N+3,1)=1
0534         KFTEMP=K(N+3,2)
0535         CALL LUKFDI(KFTEMP,K(N+4,2),KFLDMP,K(N+3,2))
0536         IF(K(N+3,2).EQ.0) GOTO 150
0537         P(N+3,5)=ULMASS(K(N+3,2))
0538         DO 460 J=1,3
0539   460   P(N+3,J)=P(N+3,J)+P(N+4,J)
0540         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)
0541         HA=P(N+1,4)**2-P(N+2,4)**2
0542         HB=HA-(P(N+1,5)**2-P(N+2,5)**2)
0543         HC=(P(N+1,1)-P(N+2,1))**2+(P(N+1,2)-P(N+2,2))**2+
0544      &  (P(N+1,3)-P(N+2,3))**2
0545         HD=(PV(1,4)-P(N+3,4))**2
0546         HE=HA**2-2.*HD*(P(N+1,4)**2+P(N+2,4)**2)+HD**2
0547         HF=HD*HC-HB**2
0548         HG=HD*HC-HA*HB
0549         HH=(SQRT(HG**2+HE*HF)-HG)/(2.*HF)
0550         DO 470 J=1,3
0551         PCOR=HH*(P(N+1,J)-P(N+2,J))
0552         P(N+1,J)=P(N+1,J)+PCOR
0553   470   P(N+2,J)=P(N+2,J)-PCOR
0554         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)
0555         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)
0556         ND=ND-1
0557       ENDIF
0558
0559 C...Check invariant mass of W jets. May give one particle or start over.
0560   480 IF(MMAT.GE.42.AND.MMAT.LE.44.AND.IABS(K(N+1,2)).LT.10) THEN
0561         PMR=SQRT(MAX(0.,P(N+1,5)**2+P(N+2,5)**2+2.*FOUR(N+1,N+2)))
0562         MSTJ(93)=1
0563         PM1=ULMASS(K(N+1,2))
0564         MSTJ(93)=1
0565         PM2=ULMASS(K(N+2,2))
0566         IF(PMR.GT.PARJ(32)+PM1+PM2) GOTO 490
0567         KFLDUM=INT(1.5+RLU(0))
0568         CALL LUKFDI(K(N+1,2),-ISIGN(KFLDUM,K(N+1,2)),KFLDMP,KF1)
0569         CALL LUKFDI(K(N+2,2),-ISIGN(KFLDUM,K(N+2,2)),KFLDMP,KF2)
0570         IF(KF1.EQ.0.OR.KF2.EQ.0) GOTO 150
0571         PSM=ULMASS(KF1)+ULMASS(KF2)
0572         IF(MMAT.EQ.42.AND.PMR.GT.PARJ(64)+PSM) GOTO 490
0573         IF(MMAT.GE.43.AND.PMR.GT.0.2*PARJ(32)+PSM) GOTO 490
0574         IF(ND.EQ.4.OR.KFA.EQ.15) GOTO 150
0575         K(N+1,1)=1
0576         KFTEMP=K(N+1,2)
0577         CALL LUKFDI(KFTEMP,K(N+2,2),KFLDMP,K(N+1,2))
0578         IF(K(N+1,2).EQ.0) GOTO 150
0579         P(N+1,5)=ULMASS(K(N+1,2))
0580         K(N+2,2)=K(N+3,2)
0581         P(N+2,5)=P(N+3,5)
0582         PS=P(N+1,5)+P(N+2,5)
0583         PV(2,5)=P(N+3,5)
0584         MMAT=0
0585         ND=2
0586         GOTO 370
0587       ENDIF
0588
0589 C...Phase space decay of partons from W decay.
0590   490 IF(MMAT.EQ.42.AND.IABS(K(N+1,2)).LT.10) THEN
0591         KFLO(1)=K(N+1,2)
0592         KFLO(2)=K(N+2,2)
0593         K(N+1,1)=K(N+3,1)
0594         K(N+1,2)=K(N+3,2)
0595         DO 500 J=1,5
0596         PV(1,J)=P(N+1,J)+P(N+2,J)
0597   500   P(N+1,J)=P(N+3,J)
0598         PV(1,5)=PMR
0599         N=N+1
0600         NP=0
0601         NQ=2
0602         PS=0.
0603         MSTJ(93)=2
0604         PSQ=ULMASS(KFLO(1))
0605         MSTJ(93)=2
0606         PSQ=PSQ+ULMASS(KFLO(2))
0607         MMAT=11
0608         GOTO 180
0609       ENDIF
0610
0611 C...Boost back for rapidly moving particle.
0612   510 N=N+ND
0613       IF(MBST.EQ.1) THEN
0614         DO 520 J=1,3
0615   520   BE(J)=P(IP,J)/P(IP,4)
0616         GA=P(IP,4)/P(IP,5)
0617         DO 540 I=NSAV+1,N
0618         BEP=BE(1)*P(I,1)+BE(2)*P(I,2)+BE(3)*P(I,3)
0619         DO 530 J=1,3
0620   530   P(I,J)=P(I,J)+GA*(GA*BEP/(1.+GA)+P(I,4))*BE(J)
0621   540   P(I,4)=GA*(P(I,4)+BEP)
0622       ENDIF
0623
0624 C...Fill in position of decay vertex.
0625       DO 560 I=NSAV+1,N
0626       DO 550 J=1,4
0627   550 V(I,J)=VDCY(J)
0628   560 V(I,5)=0.
0629
0630 C...Set up for parton shower evolution from jets.
0631       IF(MSTJ(23).GE.1.AND.MMAT.EQ.4.AND.K(NSAV+1,2).EQ.21) THEN
0632         K(NSAV+1,1)=3
0633         K(NSAV+2,1)=3
0634         K(NSAV+3,1)=3
0635         K(NSAV+1,4)=MSTU(5)*(NSAV+2)
0636         K(NSAV+1,5)=MSTU(5)*(NSAV+3)
0637         K(NSAV+2,4)=MSTU(5)*(NSAV+3)
0638         K(NSAV+2,5)=MSTU(5)*(NSAV+1)
0639         K(NSAV+3,4)=MSTU(5)*(NSAV+1)
0640         K(NSAV+3,5)=MSTU(5)*(NSAV+2)
0641         MSTJ(92)=-(NSAV+1)
0642       ELSEIF(MSTJ(23).GE.1.AND.MMAT.EQ.4) THEN
0643         K(NSAV+2,1)=3
0644         K(NSAV+3,1)=3
0645         K(NSAV+2,4)=MSTU(5)*(NSAV+3)
0646         K(NSAV+2,5)=MSTU(5)*(NSAV+3)
0647         K(NSAV+3,4)=MSTU(5)*(NSAV+2)
0648         K(NSAV+3,5)=MSTU(5)*(NSAV+2)
0649         MSTJ(92)=NSAV+2
0650       ELSEIF(MSTJ(23).GE.1.AND.(MMAT.EQ.32.OR.MMAT.EQ.44.OR.MMAT.EQ.46).
0651      &AND.IABS(K(NSAV+1,2)).LE.10.AND.IABS(K(NSAV+2,2)).LE.10) THEN
0652         K(NSAV+1,1)=3
0653         K(NSAV+2,1)=3
0654         K(NSAV+1,4)=MSTU(5)*(NSAV+2)
0655         K(NSAV+1,5)=MSTU(5)*(NSAV+2)
0656         K(NSAV+2,4)=MSTU(5)*(NSAV+1)
0657         K(NSAV+2,5)=MSTU(5)*(NSAV+1)
0658         MSTJ(92)=NSAV+1
0659       ELSEIF(MSTJ(23).GE.1.AND.MMAT.EQ.33.AND.IABS(K(NSAV+2,2)).EQ.21)
0660      &THEN
0661         K(NSAV+1,1)=3
0662         K(NSAV+2,1)=3
0663         K(NSAV+3,1)=3
0664         KCP=LUCOMP(K(NSAV+1,2))
0665         KQP=KCHG(KCP,2)*ISIGN(1,K(NSAV+1,2))
0666         JCON=4
0667         IF(KQP.LT.0) JCON=5
0668         K(NSAV+1,JCON)=MSTU(5)*(NSAV+2)
0669         K(NSAV+2,9-JCON)=MSTU(5)*(NSAV+1)
0670         K(NSAV+2,JCON)=MSTU(5)*(NSAV+3)
0671         K(NSAV+3,9-JCON)=MSTU(5)*(NSAV+2)
0672         MSTJ(92)=NSAV+1
0673       ELSEIF(MSTJ(23).GE.1.AND.MMAT.EQ.33) THEN
0674         K(NSAV+1,1)=3
0675         K(NSAV+3,1)=3
0676         K(NSAV+1,4)=MSTU(5)*(NSAV+3)
0677         K(NSAV+1,5)=MSTU(5)*(NSAV+3)
0678         K(NSAV+3,4)=MSTU(5)*(NSAV+1)
0679         K(NSAV+3,5)=MSTU(5)*(NSAV+1)
0680         MSTJ(92)=NSAV+1
0681       ENDIF
0682
0683 C...Mark decayed particle.
0684       IF(K(IP,1).EQ.5) K(IP,1)=15
0685       IF(K(IP,1).LE.10) K(IP,1)=11
0686       K(IP,4)=NSAV+1
0687       K(IP,5)=N
0688
0689       RETURN
0690       END