pythiaeRHIC/pythia/pydecy.f

0001
0002 C*********************************************************************
0003
0004 C...PYDECY
0005 C...Handles the decay of unstable particles.
0006
0007       SUBROUTINE PYDECY(IP)
0008
0009 C...Double precision and integer declarations.
0010       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0011       IMPLICIT INTEGER(I-N)
0012       INTEGER PYK,PYCHGE,PYCOMP
0013 C...Commonblocks.
0014       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
0015       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0016       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0017       COMMON/PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
0018       SAVE /PYJETS/,/PYDAT1/,/PYDAT2/,/PYDAT3/
0019 C...Local arrays.
0020       DIMENSION VDCY(4),KFLO(4),KFL1(4),PV(10,5),RORD(10),UE(3),BE(3),
0021      &WTCOR(10),PTAU(4),PCMTAU(4),DBETAU(3)
0022       CHARACTER CIDC*4
0023       DATA WTCOR/2D0,5D0,15D0,60D0,250D0,1500D0,1.2D4,1.2D5,150D0,16D0/
0024
0025 C...Functions: momentum in two-particle decays and four-product.
0026       PAWT(A,B,C)=SQRT((A**2-(B+C)**2)*(A**2-(B-C)**2))/(2D0*A)
0027       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)
0028
0029 C...Initial values.
0030       NTRY=0
0031       NSAV=N
0032       KFA=IABS(K(IP,2))
0033       KFS=ISIGN(1,K(IP,2))
0034       KC=PYCOMP(KFA)
0035       MSTJ(92)=0
0036
0037 C...Choose lifetime and determine decay vertex.
0038       IF(K(IP,1).EQ.5) THEN
0039         V(IP,5)=0D0
0040       ELSEIF(K(IP,1).NE.4) THEN
0041         V(IP,5)=-PMAS(KC,4)*LOG(PYR(0))
0042       ENDIF
0043       DO 100 J=1,4
0044         VDCY(J)=V(IP,J)+V(IP,5)*P(IP,J)/P(IP,5)
0045   100 CONTINUE
0046
0047 C...Determine whether decay allowed or not.
0048       MOUT=0
0049       IF(MSTJ(22).EQ.2) THEN
0050         IF(PMAS(KC,4).GT.PARJ(71)) MOUT=1
0051       ELSEIF(MSTJ(22).EQ.3) THEN
0052         IF(VDCY(1)**2+VDCY(2)**2+VDCY(3)**2.GT.PARJ(72)**2) MOUT=1
0053       ELSEIF(MSTJ(22).EQ.4) THEN
0054         IF(VDCY(1)**2+VDCY(2)**2.GT.PARJ(73)**2) MOUT=1
0055         IF(ABS(VDCY(3)).GT.PARJ(74)) MOUT=1
0056       ENDIF
0057       IF(MOUT.EQ.1.AND.K(IP,1).NE.5) THEN
0058         K(IP,1)=4
0059         RETURN
0060       ENDIF
0061
0062 C...Interface to external tau decay library (for tau polarization).
0063       IF(KFA.EQ.15.AND.MSTJ(28).GE.1) THEN
0064
0065 C...Starting values for pointers and momenta.
0066         ITAU=IP
0067         DO 110 J=1,4
0068           PTAU(J)=P(ITAU,J)
0069           PCMTAU(J)=P(ITAU,J)
0070   110   CONTINUE
0071
0072 C...Iterate to find position and code of mother of tau.
0073         IMTAU=ITAU
0074   120   IMTAU=K(IMTAU,3)
0075
0076         IF(IMTAU.EQ.0) THEN
0077 C...If no known origin then impossible to do anything further.
0078           KFORIG=0
0079           IORIG=0
0080
0081         ELSEIF(K(IMTAU,2).EQ.K(ITAU,2)) THEN
0082 C...If tau -> tau + gamma then add gamma energy and loop.
0083           IF(K(K(IMTAU,4),2).EQ.22) THEN
0084             DO 130 J=1,4
0085               PCMTAU(J)=PCMTAU(J)+P(K(IMTAU,4),J)
0086   130       CONTINUE
0087           ELSEIF(K(K(IMTAU,5),2).EQ.22) THEN
0088             DO 140 J=1,4
0089               PCMTAU(J)=PCMTAU(J)+P(K(IMTAU,5),J)
0090   140       CONTINUE
0091           ENDIF
0092           GOTO 120
0093
0094         ELSEIF(IABS(K(IMTAU,2)).GT.100) THEN
0095 C...If coming from weak decay of hadron then W is not stored in record,
0096 C...but can be reconstructed by adding neutrino momentum.
0097           KFORIG=-ISIGN(24,K(ITAU,2))
0098           IORIG=0
0099           DO 160 II=K(IMTAU,4),K(IMTAU,5)
0100             IF(K(II,2)*ISIGN(1,K(ITAU,2)).EQ.-16) THEN
0101               DO 150 J=1,4
0102                 PCMTAU(J)=PCMTAU(J)+P(II,J)
0103   150         CONTINUE
0104             ENDIF
0105   160     CONTINUE
0106
0107         ELSE
0108 C...If coming from resonance decay then find latest copy of this
0109 C...resonance (may not completely agree).
0110           KFORIG=K(IMTAU,2)
0111           IORIG=IMTAU
0112           DO 170 II=IMTAU+1,IP-1
0113             IF(K(II,2).EQ.KFORIG.AND.K(II,3).EQ.IORIG.AND.
0114      &      ABS(P(II,5)-P(IORIG,5)).LT.1D-5*P(IORIG,5)) IORIG=II
0115   170     CONTINUE
0116           DO 180 J=1,4
0117             PCMTAU(J)=P(IORIG,J)
0118   180     CONTINUE
0119         ENDIF
0120
0121 C...Boost tau to rest frame of production process (where known)
0122 C...and rotate it to sit along +z axis.
0123         DO 190 J=1,3
0124           DBETAU(J)=PCMTAU(J)/PCMTAU(4)
0125   190   CONTINUE
0126         IF(KFORIG.NE.0) CALL PYROBO(ITAU,ITAU,0D0,0D0,-DBETAU(1),
0127      &  -DBETAU(2),-DBETAU(3))
0128         PHITAU=PYANGL(P(ITAU,1),P(ITAU,2))
0129         CALL PYROBO(ITAU,ITAU,0D0,-PHITAU,0D0,0D0,0D0)
0130         THETAU=PYANGL(P(ITAU,3),P(ITAU,1))
0131         CALL PYROBO(ITAU,ITAU,-THETAU,0D0,0D0,0D0,0D0)
0132
0133 C...Call tau decay routine (if meaningful) and fill extra info.
0134         IF(KFORIG.NE.0.OR.MSTJ(28).EQ.2) THEN
0135           CALL PYTAUD(ITAU,IORIG,KFORIG,NDECAY)
0136           DO 200 II=NSAV+1,NSAV+NDECAY
0137             K(II,1)=1
0138             K(II,3)=IP
0139             K(II,4)=0
0140             K(II,5)=0
0141   200     CONTINUE
0142           N=NSAV+NDECAY
0143         ENDIF
0144
0145 C...Boost back decay tau and decay products.
0146         DO 210 J=1,4
0147           P(ITAU,J)=PTAU(J)
0148   210   CONTINUE
0149         IF(KFORIG.NE.0.OR.MSTJ(28).EQ.2) THEN
0150           CALL PYROBO(NSAV+1,N,THETAU,PHITAU,0D0,0D0,0D0)
0151           IF(KFORIG.NE.0) CALL PYROBO(NSAV+1,N,0D0,0D0,DBETAU(1),
0152      &    DBETAU(2),DBETAU(3))
0153
0154 C...Skip past ordinary tau decay treatment.
0155           MMAT=0
0156           MBST=0
0157           ND=0
0158           GOTO 630
0159         ENDIF
0160       ENDIF
0161
0162 C...B-Bbar mixing: flip sign of meson appropriately.
0163       MMIX=0
0164       IF((KFA.EQ.511.OR.KFA.EQ.531).AND.MSTJ(26).GE.1) THEN
0165         XBBMIX=PARJ(76)
0166         IF(KFA.EQ.531) XBBMIX=PARJ(77)
0167         IF(SIN(0.5D0*XBBMIX*V(IP,5)/PMAS(KC,4))**2.GT.PYR(0)) MMIX=1
0168         IF(MMIX.EQ.1) KFS=-KFS
0169       ENDIF
0170
0171 C...Check existence of decay channels. Particle/antiparticle rules.
0172       KCA=KC
0173       IF(MDCY(KC,2).GT.0) THEN
0174         MDMDCY=MDME(MDCY(KC,2),2)
0175         IF(MDMDCY.GT.80.AND.MDMDCY.LE.90) KCA=MDMDCY
0176       ENDIF
0177       IF(MDCY(KCA,2).LE.0.OR.MDCY(KCA,3).LE.0) THEN
0178         CALL PYERRM(9,'(PYDECY:) no decay channel defined')
0179         RETURN
0180       ENDIF
0181       IF(MOD(KFA/1000,10).EQ.0.AND.KCA.EQ.85) KFS=-KFS
0182       IF(KCHG(KC,3).EQ.0) THEN
0183         KFSP=1
0184         KFSN=0
0185         IF(PYR(0).GT.0.5D0) KFS=-KFS
0186       ELSEIF(KFS.GT.0) THEN
0187         KFSP=1
0188         KFSN=0
0189       ELSE
0190         KFSP=0
0191         KFSN=1
0192       ENDIF
0193
0194 C...Sum branching ratios of allowed decay channels.
0195   220 NOPE=0
0196       BRSU=0D0
0197       DO 230 IDL=MDCY(KCA,2),MDCY(KCA,2)+MDCY(KCA,3)-1
0198         IF(MDME(IDL,1).NE.1.AND.KFSP*MDME(IDL,1).NE.2.AND.
0199      &  KFSN*MDME(IDL,1).NE.3) GOTO 230
0200         IF(MDME(IDL,2).GT.100) GOTO 230
0201         NOPE=NOPE+1
0202         BRSU=BRSU+BRAT(IDL)
0203   230 CONTINUE
0204       IF(NOPE.EQ.0) THEN
0205         CALL PYERRM(2,'(PYDECY:) all decay channels closed by user')
0206         RETURN
0207       ENDIF
0208
0209 C...Select decay channel among allowed ones.
0210   240 RBR=BRSU*PYR(0)
0211       IDL=MDCY(KCA,2)-1
0212   250 IDL=IDL+1
0213       IF(MDME(IDL,1).NE.1.AND.KFSP*MDME(IDL,1).NE.2.AND.
0214      &KFSN*MDME(IDL,1).NE.3) THEN
0215         IF(IDL.LT.MDCY(KCA,2)+MDCY(KCA,3)-1) GOTO 250
0216       ELSEIF(MDME(IDL,2).GT.100) THEN
0217         IF(IDL.LT.MDCY(KCA,2)+MDCY(KCA,3)-1) GOTO 250
0218       ELSE
0219         IDC=IDL
0220         RBR=RBR-BRAT(IDL)
0221         IF(IDL.LT.MDCY(KCA,2)+MDCY(KCA,3)-1.AND.RBR.GT.0D0) GOTO 250
0222       ENDIF
0223
0224 C...Start readout of decay channel: matrix element, reset counters.
0225       MMAT=MDME(IDC,2)
0226   260 NTRY=NTRY+1
0227       IF(MOD(NTRY,200).EQ.0) THEN
0228         WRITE(CIDC,'(I4)') IDC
0229 C...Do not print warning for some well-known special cases.
0230         IF(KFA.NE.113.AND.KFA.NE.115.AND.KFA.NE.215)
0231      &  CALL PYERRM(4,'(PYDECY:) caught in loop for decay channel'//
0232      &  CIDC)
0233         GOTO 240
0234       ENDIF
0235       IF(NTRY.GT.1000) THEN
0236         CALL PYERRM(14,'(PYDECY:) caught in infinite loop')
0237         IF(MSTU(21).GE.1) RETURN
0238       ENDIF
0239       I=N
0240       NP=0
0241       NQ=0
0242       MBST=0
0243       IF(MMAT.GE.11.AND.P(IP,4).GT.20D0*P(IP,5)) MBST=1
0244       DO 270 J=1,4
0245         PV(1,J)=0D0
0246         IF(MBST.EQ.0) PV(1,J)=P(IP,J)
0247   270 CONTINUE
0248       IF(MBST.EQ.1) PV(1,4)=P(IP,5)
0249       PV(1,5)=P(IP,5)
0250       PS=0D0
0251       PSQ=0D0
0252       MREM=0
0253       MHADDY=0
0254       IF(KFA.GT.80) MHADDY=1
0255 C.. Random flavour and popcorn system memory.
0256       IRNDMO=0
0257       JTMO=0
0258       MSTU(121)=0
0259       MSTU(125)=10
0260
0261 C...Read out decay products. Convert to standard flavour code.
0262       JTMAX=5
0263       IF(MDME(IDC+1,2).EQ.101) JTMAX=10
0264       DO 280 JT=1,JTMAX
0265         IF(JT.LE.5) KP=KFDP(IDC,JT)
0266         IF(JT.GE.6) KP=KFDP(IDC+1,JT-5)
0267         IF(KP.EQ.0) GOTO 280
0268         KPA=IABS(KP)
0269         KCP=PYCOMP(KPA)
0270         IF(KPA.GT.80) MHADDY=1
0271         IF(KCHG(KCP,3).EQ.0.AND.KPA.NE.81.AND.KPA.NE.82) THEN
0272           KFP=KP
0273         ELSEIF(KPA.NE.81.AND.KPA.NE.82) THEN
0274           KFP=KFS*KP
0275         ELSEIF(KPA.EQ.81.AND.MOD(KFA/1000,10).EQ.0) THEN
0276           KFP=-KFS*MOD(KFA/10,10)
0277         ELSEIF(KPA.EQ.81.AND.MOD(KFA/100,10).GE.MOD(KFA/10,10)) THEN
0278           KFP=KFS*(100*MOD(KFA/10,100)+3)
0279         ELSEIF(KPA.EQ.81) THEN
0280           KFP=KFS*(1000*MOD(KFA/10,10)+100*MOD(KFA/100,10)+1)
0281         ELSEIF(KP.EQ.82) THEN
0282           CALL PYDCYK(-KFS*INT(1D0+(2D0+PARJ(2))*PYR(0)),0,KFP,KDUMP)
0283           IF(KFP.EQ.0) GOTO 260
0284           KFP=-KFP
0285           IRNDMO=1
0286           MSTJ(93)=1
0287           IF(PV(1,5).LT.PARJ(32)+2D0*PYMASS(KFP)) GOTO 260
0288         ELSEIF(KP.EQ.-82) THEN
0289           KFP=MSTU(124)
0290         ENDIF
0291         IF(KPA.EQ.81.OR.KPA.EQ.82) KCP=PYCOMP(KFP)
0292
0293 C...Add decay product to event record or to quark flavour list.
0294         KFPA=IABS(KFP)
0295         KQP=KCHG(KCP,2)
0296         IF(MMAT.GE.11.AND.MMAT.LE.30.AND.KQP.NE.0) THEN
0297           NQ=NQ+1
0298           KFLO(NQ)=KFP
0299 C...set rndmflav popcorn system pointer
0300           IF(KP.EQ.82.AND.MSTU(121).GT.0) JTMO=NQ
0301           MSTJ(93)=2
0302           PSQ=PSQ+PYMASS(KFLO(NQ))
0303         ELSEIF((MMAT.EQ.42.OR.MMAT.EQ.43.OR.MMAT.EQ.48).AND.NP.EQ.3.AND.
0304      &    MOD(NQ,2).EQ.1) THEN
0305           NQ=NQ-1
0306           PS=PS-P(I,5)
0307           K(I,1)=1
0308           KFI=K(I,2)
0309           CALL PYKFDI(KFP,KFI,KFLDMP,K(I,2))
0310           IF(K(I,2).EQ.0) GOTO 260
0311           MSTJ(93)=1
0312           P(I,5)=PYMASS(K(I,2))
0313           PS=PS+P(I,5)
0314         ELSE
0315           I=I+1
0316           NP=NP+1
0317           IF(MMAT.NE.33.AND.KQP.NE.0) NQ=NQ+1
0318           IF(MMAT.EQ.33.AND.KQP.NE.0.AND.KQP.NE.2) NQ=NQ+1
0319           K(I,1)=1+MOD(NQ,2)
0320           IF(MMAT.EQ.4.AND.JT.LE.2.AND.KFP.EQ.21) K(I,1)=2
0321           IF(MMAT.EQ.4.AND.JT.EQ.3) K(I,1)=1
0322           K(I,2)=KFP
0323           K(I,3)=IP
0324           K(I,4)=0
0325           K(I,5)=0
0326           P(I,5)=PYMASS(KFP)
0327           PS=PS+P(I,5)
0328         ENDIF
0329   280 CONTINUE
0330
0331 C...Check masses for resonance decays.
0332       IF(MHADDY.EQ.0) THEN
0333         IF(PS+PARJ(64).GT.PV(1,5)) GOTO 240
0334       ENDIF
0335
0336 C...Choose decay multiplicity in phase space model.
0337   290 IF(MMAT.GE.11.AND.MMAT.LE.30) THEN
0338         PSP=PS
0339         CNDE=PARJ(61)*LOG(MAX((PV(1,5)-PS-PSQ)/PARJ(62),1.1D0))
0340         IF(MMAT.EQ.12) CNDE=CNDE+PARJ(63)
0341   300   NTRY=NTRY+1
0342 C...Reset popcorn flags if new attempt. Re-select rndmflav if failed.
0343         IF(IRNDMO.EQ.0) THEN
0344            MSTU(121)=0
0345            JTMO=0
0346         ELSEIF(IRNDMO.EQ.1) THEN
0347            IRNDMO=2
0348         ELSE
0349            GOTO 260
0350         ENDIF
0351         IF(NTRY.GT.1000) THEN
0352           CALL PYERRM(14,'(PYDECY:) caught in infinite loop')
0353           IF(MSTU(21).GE.1) RETURN
0354         ENDIF
0355         IF(MMAT.LE.20) THEN
0356           GAUSS=SQRT(-2D0*CNDE*LOG(MAX(1D-10,PYR(0))))*
0357      &    SIN(PARU(2)*PYR(0))
0358           ND=0.5D0+0.5D0*NP+0.25D0*NQ+CNDE+GAUSS
0359           IF(ND.LT.NP+NQ/2.OR.ND.LT.2.OR.ND.GT.10) GOTO 300
0360           IF(MMAT.EQ.13.AND.ND.EQ.2) GOTO 300
0361           IF(MMAT.EQ.14.AND.ND.LE.3) GOTO 300
0362           IF(MMAT.EQ.15.AND.ND.LE.4) GOTO 300
0363         ELSE
0364           ND=MMAT-20
0365         ENDIF
0366 C.. Set maximum popcorn meson number. Test rndmflav popcorn size.
0367         MSTU(125)=ND-NQ/2
0368         IF(MSTU(121).GT.MSTU(125)) GOTO 300
0369
0370 C...Form hadrons from flavour content.
0371         DO 310 JT=1,NQ
0372           KFL1(JT)=KFLO(JT)
0373   310   CONTINUE
0374         IF(ND.EQ.NP+NQ/2) GOTO 330
0375         DO 320 I=N+NP+1,N+ND-NQ/2
0376 C.. Stick to started popcorn system, else pick side at random
0377           JT=JTMO
0378           IF(JT.EQ.0) JT=1+INT((NQ-1)*PYR(0))
0379           CALL PYDCYK(KFL1(JT),0,KFL2,K(I,2))
0380           IF(K(I,2).EQ.0) GOTO 300
0381           MSTU(125)=MSTU(125)-1
0382           JTMO=0
0383           IF(MSTU(121).GT.0) JTMO=JT
0384           KFL1(JT)=-KFL2
0385   320   CONTINUE
0386   330   JT=2
0387         JT2=3
0388         JT3=4
0389         IF(NQ.EQ.4.AND.PYR(0).LT.PARJ(66)) JT=4
0390         IF(JT.EQ.4.AND.ISIGN(1,KFL1(1)*(10-IABS(KFL1(1))))*
0391      &  ISIGN(1,KFL1(JT)*(10-IABS(KFL1(JT)))).GT.0) JT=3
0392         IF(JT.EQ.3) JT2=2
0393         IF(JT.EQ.4) JT3=2
0394         CALL PYDCYK(KFL1(1),KFL1(JT),KFLDMP,K(N+ND-NQ/2+1,2))
0395         IF(K(N+ND-NQ/2+1,2).EQ.0) GOTO 300
0396         IF(NQ.EQ.4) CALL PYDCYK(KFL1(JT2),KFL1(JT3),KFLDMP,K(N+ND,2))
0397         IF(NQ.EQ.4.AND.K(N+ND,2).EQ.0) GOTO 300
0398
0399 C...Check that sum of decay product masses not too large.
0400         PS=PSP
0401         DO 340 I=N+NP+1,N+ND
0402           K(I,1)=1
0403           K(I,3)=IP
0404           K(I,4)=0
0405           K(I,5)=0
0406           P(I,5)=PYMASS(K(I,2))
0407           PS=PS+P(I,5)
0408   340   CONTINUE
0409         IF(PS+PARJ(64).GT.PV(1,5)) GOTO 300
0410
0411 C...Rescale energy to subtract off spectator quark mass.
0412       ELSEIF((MMAT.EQ.31.OR.MMAT.EQ.33.OR.MMAT.EQ.44)
0413      &  .AND.NP.GE.3) THEN
0414         PS=PS-P(N+NP,5)
0415         PQT=(P(N+NP,5)+PARJ(65))/PV(1,5)
0416         DO 350 J=1,5
0417           P(N+NP,J)=PQT*PV(1,J)
0418           PV(1,J)=(1D0-PQT)*PV(1,J)
0419   350   CONTINUE
0420         IF(PS+PARJ(64).GT.PV(1,5)) GOTO 260
0421         ND=NP-1
0422         MREM=1
0423
0424 C...Fully specified final state: check mass broadening effects.
0425       ELSE
0426         IF(NP.GE.2.AND.PS+PARJ(64).GT.PV(1,5)) GOTO 260
0427         ND=NP
0428       ENDIF
0429
0430 C...Determine position of grandmother, number of sisters.
0431       NM=0
0432       KFAS=0
0433       MSGN=0
0434       IF(MMAT.EQ.3) THEN
0435         IM=K(IP,3)
0436         IF(IM.LT.0.OR.IM.GE.IP) IM=0
0437         IF(IM.NE.0) KFAM=IABS(K(IM,2))
0438         IF(IM.NE.0) THEN
0439           DO 360 IL=MAX(IP-2,IM+1),MIN(IP+2,N)
0440             IF(K(IL,3).EQ.IM) NM=NM+1
0441             IF(K(IL,3).EQ.IM.AND.IL.NE.IP) ISIS=IL
0442   360     CONTINUE
0443           IF(NM.NE.2.OR.KFAM.LE.100.OR.MOD(KFAM,10).NE.1.OR.
0444      &    MOD(KFAM/1000,10).NE.0) NM=0
0445           IF(NM.EQ.2) THEN
0446             KFAS=IABS(K(ISIS,2))
0447             IF((KFAS.LE.100.OR.MOD(KFAS,10).NE.1.OR.
0448      &      MOD(KFAS/1000,10).NE.0).AND.KFAS.NE.22) NM=0
0449           ENDIF
0450         ENDIF
0451       ENDIF
0452
0453 C...Kinematics of one-particle decays.
0454       IF(ND.EQ.1) THEN
0455         DO 370 J=1,4
0456           P(N+1,J)=P(IP,J)
0457   370   CONTINUE
0458         GOTO 630
0459       ENDIF
0460
0461 C...Calculate maximum weight ND-particle decay.
0462       PV(ND,5)=P(N+ND,5)
0463       IF(ND.GE.3) THEN
0464         WTMAX=1D0/WTCOR(ND-2)
0465         PMAX=PV(1,5)-PS+P(N+ND,5)
0466         PMIN=0D0
0467         DO 380 IL=ND-1,1,-1
0468           PMAX=PMAX+P(N+IL,5)
0469           PMIN=PMIN+P(N+IL+1,5)
0470           WTMAX=WTMAX*PAWT(PMAX,PMIN,P(N+IL,5))
0471   380   CONTINUE
0472       ENDIF
0473
0474 C...Find virtual gamma mass in Dalitz decay.
0475   390 IF(ND.EQ.2) THEN
0476       ELSEIF(MMAT.EQ.2) THEN
0477         PMES=4D0*PMAS(11,1)**2
0478         PMRHO2=PMAS(131,1)**2
0479         PGRHO2=PMAS(131,2)**2
0480   400   PMST=PMES*(P(IP,5)**2/PMES)**PYR(0)
0481         WT=(1+0.5D0*PMES/PMST)*SQRT(MAX(0D0,1D0-PMES/PMST))*
0482      &  (1D0-PMST/P(IP,5)**2)**3*(1D0+PGRHO2/PMRHO2)/
0483      &  ((1D0-PMST/PMRHO2)**2+PGRHO2/PMRHO2)
0484         IF(WT.LT.PYR(0)) GOTO 400
0485         PV(2,5)=MAX(2.00001D0*PMAS(11,1),SQRT(PMST))
0486
0487 C...M-generator gives weight. If rejected, try again.
0488       ELSE
0489   410   RORD(1)=1D0
0490         DO 440 IL1=2,ND-1
0491           RSAV=PYR(0)
0492           DO 420 IL2=IL1-1,1,-1
0493             IF(RSAV.LE.RORD(IL2)) GOTO 430
0494             RORD(IL2+1)=RORD(IL2)
0495   420     CONTINUE
0496   430     RORD(IL2+1)=RSAV
0497   440   CONTINUE
0498         RORD(ND)=0D0
0499         WT=1D0
0500         DO 450 IL=ND-1,1,-1
0501           PV(IL,5)=PV(IL+1,5)+P(N+IL,5)+(RORD(IL)-RORD(IL+1))*
0502      &    (PV(1,5)-PS)
0503           WT=WT*PAWT(PV(IL,5),PV(IL+1,5),P(N+IL,5))
0504   450   CONTINUE
0505         IF(WT.LT.PYR(0)*WTMAX) GOTO 410
0506       ENDIF
0507
0508 C...Perform two-particle decays in respective CM frame.
0509   460 DO 480 IL=1,ND-1
0510         PA=PAWT(PV(IL,5),PV(IL+1,5),P(N+IL,5))
0511         UE(3)=2D0*PYR(0)-1D0
0512         PHI=PARU(2)*PYR(0)
0513         UE(1)=SQRT(1D0-UE(3)**2)*COS(PHI)
0514         UE(2)=SQRT(1D0-UE(3)**2)*SIN(PHI)
0515         DO 470 J=1,3
0516           P(N+IL,J)=PA*UE(J)
0517           PV(IL+1,J)=-PA*UE(J)
0518   470   CONTINUE
0519         P(N+IL,4)=SQRT(PA**2+P(N+IL,5)**2)
0520         PV(IL+1,4)=SQRT(PA**2+PV(IL+1,5)**2)
0521   480 CONTINUE
0522
0523 C...Lorentz transform decay products to lab frame.
0524       DO 490 J=1,4
0525         P(N+ND,J)=PV(ND,J)
0526   490 CONTINUE
0527       DO 530 IL=ND-1,1,-1
0528         DO 500 J=1,3
0529           BE(J)=PV(IL,J)/PV(IL,4)
0530   500   CONTINUE
0531         GA=PV(IL,4)/PV(IL,5)
0532         DO 520 I=N+IL,N+ND
0533           BEP=BE(1)*P(I,1)+BE(2)*P(I,2)+BE(3)*P(I,3)
0534           DO 510 J=1,3
0535             P(I,J)=P(I,J)+GA*(GA*BEP/(1D0+GA)+P(I,4))*BE(J)
0536   510     CONTINUE
0537           P(I,4)=GA*(P(I,4)+BEP)
0538   520   CONTINUE
0539   530 CONTINUE
0540
0541 C...Check that no infinite loop in matrix element weight.
0542       NTRY=NTRY+1
0543       IF(NTRY.GT.800) GOTO 560
0544
0545 C...Matrix elements for omega and phi decays.
0546       IF(MMAT.EQ.1) THEN
0547         WT=(P(N+1,5)*P(N+2,5)*P(N+3,5))**2-(P(N+1,5)*FOUR(N+2,N+3))**2
0548      &  -(P(N+2,5)*FOUR(N+1,N+3))**2-(P(N+3,5)*FOUR(N+1,N+2))**2
0549      &  +2D0*FOUR(N+1,N+2)*FOUR(N+1,N+3)*FOUR(N+2,N+3)
0550         IF(MAX(WT*WTCOR(9)/P(IP,5)**6,0.001D0).LT.PYR(0)) GOTO 390
0551
0552 C...Matrix elements for pi0 or eta Dalitz decay to gamma e+ e-.
0553       ELSEIF(MMAT.EQ.2) THEN
0554         FOUR12=FOUR(N+1,N+2)
0555         FOUR13=FOUR(N+1,N+3)
0556         WT=(PMST-0.5D0*PMES)*(FOUR12**2+FOUR13**2)+
0557      &  PMES*(FOUR12*FOUR13+FOUR12**2+FOUR13**2)
0558         IF(WT.LT.PYR(0)*0.25D0*PMST*(P(IP,5)**2-PMST)**2) GOTO 460
0559
0560 C...Matrix element for S0 -> S1 + V1 -> S1 + S2 + S3 (S scalar,
0561 C...V vector), of form cos**2(theta02) in V1 rest frame, and for
0562 C...S0 -> gamma + V1 -> gamma + S2 + S3, of form sin**2(theta02).
0563       ELSEIF(MMAT.EQ.3.AND.NM.EQ.2) THEN
0564         FOUR10=FOUR(IP,IM)
0565         FOUR12=FOUR(IP,N+1)
0566         FOUR02=FOUR(IM,N+1)
0567         PMS1=P(IP,5)**2
0568         PMS0=P(IM,5)**2
0569         PMS2=P(N+1,5)**2
0570         IF(KFAS.NE.22) HNUM=(FOUR10*FOUR12-PMS1*FOUR02)**2
0571         IF(KFAS.EQ.22) HNUM=PMS1*(2D0*FOUR10*FOUR12*FOUR02-
0572      &  PMS1*FOUR02**2-PMS0*FOUR12**2-PMS2*FOUR10**2+PMS1*PMS0*PMS2)
0573         HNUM=MAX(1D-6*PMS1**2*PMS0*PMS2,HNUM)
0574         HDEN=(FOUR10**2-PMS1*PMS0)*(FOUR12**2-PMS1*PMS2)
0575         IF(HNUM.LT.PYR(0)*HDEN) GOTO 460
0576
0577 C...Matrix element for "onium" -> g + g + g or gamma + g + g.
0578       ELSEIF(MMAT.EQ.4) THEN
0579         HX1=2D0*FOUR(IP,N+1)/P(IP,5)**2
0580         HX2=2D0*FOUR(IP,N+2)/P(IP,5)**2
0581         HX3=2D0*FOUR(IP,N+3)/P(IP,5)**2
0582         WT=((1D0-HX1)/(HX2*HX3))**2+((1D0-HX2)/(HX1*HX3))**2+
0583      &  ((1D0-HX3)/(HX1*HX2))**2
0584         IF(WT.LT.2D0*PYR(0)) GOTO 390
0585         IF(K(IP+1,2).EQ.22.AND.(1D0-HX1)*P(IP,5)**2.LT.4D0*PARJ(32)**2)
0586      &  GOTO 390
0587
0588 C...Effective matrix element for nu spectrum in tau -> nu + hadrons.
0589       ELSEIF(MMAT.EQ.41) THEN
0590         IF(MBST.EQ.0) HX1=2D0*FOUR(IP,N+1)/P(IP,5)**2
0591         IF(MBST.EQ.1) HX1=2D0*P(N+1,4)/P(IP,5)
0592         HXM=MIN(0.75D0,2D0*(1D0-PS/P(IP,5)))
0593         IF(HX1*(3D0-2D0*HX1).LT.PYR(0)*HXM*(3D0-2D0*HXM)) GOTO 390
0594
0595 C...Matrix elements for weak decays (only semileptonic for c and b)
0596       ELSEIF((MMAT.EQ.42.OR.MMAT.EQ.43.OR.MMAT.EQ.44.OR.MMAT.EQ.48)
0597      &  .AND.ND.EQ.3) THEN
0598         IF(MBST.EQ.0) WT=FOUR(IP,N+1)*FOUR(N+2,N+3)
0599         IF(MBST.EQ.1) WT=P(IP,5)*P(N+1,4)*FOUR(N+2,N+3)
0600         IF(WT.LT.PYR(0)*P(IP,5)*PV(1,5)**3/WTCOR(10)) GOTO 390
0601       ELSEIF(MMAT.EQ.42.OR.MMAT.EQ.43.OR.MMAT.EQ.44.OR.MMAT.EQ.48) THEN
0602         DO 550 J=1,4
0603           P(N+NP+1,J)=0D0
0604           DO 540 IS=N+3,N+NP
0605             P(N+NP+1,J)=P(N+NP+1,J)+P(IS,J)
0606   540     CONTINUE
0607   550   CONTINUE
0608         IF(MBST.EQ.0) WT=FOUR(IP,N+1)*FOUR(N+2,N+NP+1)
0609         IF(MBST.EQ.1) WT=P(IP,5)*P(N+1,4)*FOUR(N+2,N+NP+1)
0610         IF(WT.LT.PYR(0)*P(IP,5)*PV(1,5)**3/WTCOR(10)) GOTO 390
0611       ENDIF
0612
0613 C...Scale back energy and reattach spectator.
0614   560 IF(MREM.EQ.1) THEN
0615         DO 570 J=1,5
0616           PV(1,J)=PV(1,J)/(1D0-PQT)
0617   570   CONTINUE
0618         ND=ND+1
0619         MREM=0
0620       ENDIF
0621
0622 C...Low invariant mass for system with spectator quark gives particle,
0623 C...not two jets. Readjust momenta accordingly.
0624       IF(MMAT.EQ.31.AND.ND.EQ.3) THEN
0625         MSTJ(93)=1
0626         PM2=PYMASS(K(N+2,2))
0627         MSTJ(93)=1
0628         PM3=PYMASS(K(N+3,2))
0629         IF(P(N+2,5)**2+P(N+3,5)**2+2D0*FOUR(N+2,N+3).GE.
0630      &  (PARJ(32)+PM2+PM3)**2) GOTO 630
0631         K(N+2,1)=1
0632         KFTEMP=K(N+2,2)
0633         CALL PYKFDI(KFTEMP,K(N+3,2),KFLDMP,K(N+2,2))
0634         IF(K(N+2,2).EQ.0) GOTO 260
0635         P(N+2,5)=PYMASS(K(N+2,2))
0636         PS=P(N+1,5)+P(N+2,5)
0637         PV(2,5)=P(N+2,5)
0638         MMAT=0
0639         ND=2
0640         GOTO 460
0641       ELSEIF(MMAT.EQ.44) THEN
0642         MSTJ(93)=1
0643         PM3=PYMASS(K(N+3,2))
0644         MSTJ(93)=1
0645         PM4=PYMASS(K(N+4,2))
0646         IF(P(N+3,5)**2+P(N+4,5)**2+2D0*FOUR(N+3,N+4).GE.
0647      &  (PARJ(32)+PM3+PM4)**2) GOTO 600
0648         K(N+3,1)=1
0649         KFTEMP=K(N+3,2)
0650         CALL PYKFDI(KFTEMP,K(N+4,2),KFLDMP,K(N+3,2))
0651         IF(K(N+3,2).EQ.0) GOTO 260
0652         P(N+3,5)=PYMASS(K(N+3,2))
0653         DO 580 J=1,3
0654           P(N+3,J)=P(N+3,J)+P(N+4,J)
0655   580   CONTINUE
0656         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)
0657         HA=P(N+1,4)**2-P(N+2,4)**2
0658         HB=HA-(P(N+1,5)**2-P(N+2,5)**2)
0659         HC=(P(N+1,1)-P(N+2,1))**2+(P(N+1,2)-P(N+2,2))**2+
0660      &  (P(N+1,3)-P(N+2,3))**2
0661         HD=(PV(1,4)-P(N+3,4))**2
0662         HE=HA**2-2D0*HD*(P(N+1,4)**2+P(N+2,4)**2)+HD**2
0663         HF=HD*HC-HB**2
0664         HG=HD*HC-HA*HB
0665         HH=(SQRT(HG**2+HE*HF)-HG)/(2D0*HF)
0666         DO 590 J=1,3
0667           PCOR=HH*(P(N+1,J)-P(N+2,J))
0668           P(N+1,J)=P(N+1,J)+PCOR
0669           P(N+2,J)=P(N+2,J)-PCOR
0670   590   CONTINUE
0671         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)
0672         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)
0673         ND=ND-1
0674       ENDIF
0675
0676 C...Check invariant mass of W jets. May give one particle or start over.
0677   600 IF((MMAT.EQ.42.OR.MMAT.EQ.43.OR.MMAT.EQ.44.OR.MMAT.EQ.48)
0678      &.AND.IABS(K(N+1,2)).LT.10) THEN
0679         PMR=SQRT(MAX(0D0,P(N+1,5)**2+P(N+2,5)**2+2D0*FOUR(N+1,N+2)))
0680         MSTJ(93)=1
0681         PM1=PYMASS(K(N+1,2))
0682         MSTJ(93)=1
0683         PM2=PYMASS(K(N+2,2))
0684         IF(PMR.GT.PARJ(32)+PM1+PM2) GOTO 610
0685         KFLDUM=INT(1.5D0+PYR(0))
0686         CALL PYKFDI(K(N+1,2),-ISIGN(KFLDUM,K(N+1,2)),KFLDMP,KF1)
0687         CALL PYKFDI(K(N+2,2),-ISIGN(KFLDUM,K(N+2,2)),KFLDMP,KF2)
0688         IF(KF1.EQ.0.OR.KF2.EQ.0) GOTO 260
0689         PSM=PYMASS(KF1)+PYMASS(KF2)
0690         IF((MMAT.EQ.42.OR.MMAT.EQ.48).AND.PMR.GT.PARJ(64)+PSM) GOTO 610
0691         IF(MMAT.GE.43.AND.PMR.GT.0.2D0*PARJ(32)+PSM) GOTO 610
0692         IF(MMAT.EQ.48) GOTO 390
0693         IF(ND.EQ.4.OR.KFA.EQ.15) GOTO 260
0694         K(N+1,1)=1
0695         KFTEMP=K(N+1,2)
0696         CALL PYKFDI(KFTEMP,K(N+2,2),KFLDMP,K(N+1,2))
0697         IF(K(N+1,2).EQ.0) GOTO 260
0698         P(N+1,5)=PYMASS(K(N+1,2))
0699         K(N+2,2)=K(N+3,2)
0700         P(N+2,5)=P(N+3,5)
0701         PS=P(N+1,5)+P(N+2,5)
0702         IF(PS+PARJ(64).GT.PV(1,5)) GOTO 260
0703         PV(2,5)=P(N+3,5)
0704         MMAT=0
0705         ND=2
0706         GOTO 460
0707       ENDIF
0708
0709 C...Phase space decay of partons from W decay.
0710   610 IF((MMAT.EQ.42.OR.MMAT.EQ.48).AND.IABS(K(N+1,2)).LT.10) THEN
0711         KFLO(1)=K(N+1,2)
0712         KFLO(2)=K(N+2,2)
0713         K(N+1,1)=K(N+3,1)
0714         K(N+1,2)=K(N+3,2)
0715         DO 620 J=1,5
0716           PV(1,J)=P(N+1,J)+P(N+2,J)
0717           P(N+1,J)=P(N+3,J)
0718   620   CONTINUE
0719         PV(1,5)=PMR
0720         N=N+1
0721         NP=0
0722         NQ=2
0723         PS=0D0
0724         MSTJ(93)=2
0725         PSQ=PYMASS(KFLO(1))
0726         MSTJ(93)=2
0727         PSQ=PSQ+PYMASS(KFLO(2))
0728         MMAT=11
0729         GOTO 290
0730       ENDIF
0731
0732 C...Boost back for rapidly moving particle.
0733   630 N=N+ND
0734       IF(MBST.EQ.1) THEN
0735         DO 640 J=1,3
0736           BE(J)=P(IP,J)/P(IP,4)
0737   640   CONTINUE
0738         GA=P(IP,4)/P(IP,5)
0739         DO 660 I=NSAV+1,N
0740           BEP=BE(1)*P(I,1)+BE(2)*P(I,2)+BE(3)*P(I,3)
0741           DO 650 J=1,3
0742             P(I,J)=P(I,J)+GA*(GA*BEP/(1D0+GA)+P(I,4))*BE(J)
0743   650     CONTINUE
0744           P(I,4)=GA*(P(I,4)+BEP)
0745   660   CONTINUE
0746       ENDIF
0747
0748 C...Fill in position of decay vertex.
0749       DO 680 I=NSAV+1,N
0750         DO 670 J=1,4
0751           V(I,J)=VDCY(J)
0752   670   CONTINUE
0753         V(I,5)=0D0
0754   680 CONTINUE
0755
0756 C...Set up for parton shower evolution from jets.
0757       IF(MSTJ(23).GE.1.AND.MMAT.EQ.4.AND.K(NSAV+1,2).EQ.21) THEN
0758         K(NSAV+1,1)=3
0759         K(NSAV+2,1)=3
0760         K(NSAV+3,1)=3
0761         K(NSAV+1,4)=MSTU(5)*(NSAV+2)
0762         K(NSAV+1,5)=MSTU(5)*(NSAV+3)
0763         K(NSAV+2,4)=MSTU(5)*(NSAV+3)
0764         K(NSAV+2,5)=MSTU(5)*(NSAV+1)
0765         K(NSAV+3,4)=MSTU(5)*(NSAV+1)
0766         K(NSAV+3,5)=MSTU(5)*(NSAV+2)
0767         MSTJ(92)=-(NSAV+1)
0768       ELSEIF(MSTJ(23).GE.1.AND.MMAT.EQ.4) THEN
0769         K(NSAV+2,1)=3
0770         K(NSAV+3,1)=3
0771         K(NSAV+2,4)=MSTU(5)*(NSAV+3)
0772         K(NSAV+2,5)=MSTU(5)*(NSAV+3)
0773         K(NSAV+3,4)=MSTU(5)*(NSAV+2)
0774         K(NSAV+3,5)=MSTU(5)*(NSAV+2)
0775         MSTJ(92)=NSAV+2
0776       ELSEIF(MSTJ(23).GE.1.AND.(MMAT.EQ.32.OR.MMAT.EQ.44).AND.
0777      &  IABS(K(NSAV+1,2)).LE.10.AND.IABS(K(NSAV+2,2)).LE.10) THEN
0778         K(NSAV+1,1)=3
0779         K(NSAV+2,1)=3
0780         K(NSAV+1,4)=MSTU(5)*(NSAV+2)
0781         K(NSAV+1,5)=MSTU(5)*(NSAV+2)
0782         K(NSAV+2,4)=MSTU(5)*(NSAV+1)
0783         K(NSAV+2,5)=MSTU(5)*(NSAV+1)
0784         MSTJ(92)=NSAV+1
0785       ELSEIF(MSTJ(23).GE.1.AND.(MMAT.EQ.32.OR.MMAT.EQ.44).AND.
0786      &  IABS(K(NSAV+1,2)).LE.20.AND.IABS(K(NSAV+2,2)).LE.20) THEN
0787         MSTJ(92)=NSAV+1
0788       ELSEIF(MSTJ(23).GE.1.AND.MMAT.EQ.33.AND.IABS(K(NSAV+2,2)).EQ.21)
0789      &  THEN
0790         K(NSAV+1,1)=3
0791         K(NSAV+2,1)=3
0792         K(NSAV+3,1)=3
0793         KCP=PYCOMP(K(NSAV+1,2))
0794         KQP=KCHG(KCP,2)*ISIGN(1,K(NSAV+1,2))
0795         JCON=4
0796         IF(KQP.LT.0) JCON=5
0797         K(NSAV+1,JCON)=MSTU(5)*(NSAV+2)
0798         K(NSAV+2,9-JCON)=MSTU(5)*(NSAV+1)
0799         K(NSAV+2,JCON)=MSTU(5)*(NSAV+3)
0800         K(NSAV+3,9-JCON)=MSTU(5)*(NSAV+2)
0801         MSTJ(92)=NSAV+1
0802       ELSEIF(MSTJ(23).GE.1.AND.MMAT.EQ.33) THEN
0803         K(NSAV+1,1)=3
0804         K(NSAV+3,1)=3
0805         K(NSAV+1,4)=MSTU(5)*(NSAV+3)
0806         K(NSAV+1,5)=MSTU(5)*(NSAV+3)
0807         K(NSAV+3,4)=MSTU(5)*(NSAV+1)
0808         K(NSAV+3,5)=MSTU(5)*(NSAV+1)
0809         MSTJ(92)=NSAV+1
0810       ENDIF
0811
0812 C...Mark decayed particle; special option for B-Bbar mixing.
0813       IF(K(IP,1).EQ.5) K(IP,1)=15
0814       IF(K(IP,1).LE.10) K(IP,1)=11
0815       IF(MMIX.EQ.1.AND.MSTJ(26).EQ.2.AND.K(IP,1).EQ.11) K(IP,1)=12
0816       K(IP,4)=NSAV+1
0817       K(IP,5)=N
0818
0819       RETURN
0820       END