hijing/src/lushow.f

0001
0002 C*********************************************************************
0003
0004       SUBROUTINE LUSHOW(IP1,IP2,QMAX)
0005
0006 C...Purpose: to generate timelike parton showers from given partons.
0007       IMPLICIT DOUBLE PRECISION(D)
0008       COMMON/LUJETS/N,K(9000,5),P(9000,5),V(9000,5)
0009       SAVE /LUJETS/
0010       COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0011       SAVE /LUDAT1/
0012       COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
0013       SAVE /LUDAT2/
0014       DIMENSION PMTH(5,40),PS(5),PMA(4),PMSD(4),IEP(4),IPA(4),
0015      &KFLA(4),KFLD(4),KFL(4),ITRY(4),ISI(4),ISL(4),DP(4),DPT(5,4)
0016
0017 C...Initialization of cutoff masses etc.
0018       IF(MSTJ(41).LE.0.OR.(MSTJ(41).EQ.1.AND.QMAX.LE.PARJ(82)).OR.
0019      &QMAX.LE.MIN(PARJ(82),PARJ(83)).OR.MSTJ(41).GE.3) RETURN
0020       PMTH(1,21)=ULMASS(21)
0021       PMTH(2,21)=SQRT(PMTH(1,21)**2+0.25*PARJ(82)**2)
0022       PMTH(3,21)=2.*PMTH(2,21)
0023       PMTH(4,21)=PMTH(3,21)
0024       PMTH(5,21)=PMTH(3,21)
0025       PMTH(1,22)=ULMASS(22)
0026       PMTH(2,22)=SQRT(PMTH(1,22)**2+0.25*PARJ(83)**2)
0027       PMTH(3,22)=2.*PMTH(2,22)
0028       PMTH(4,22)=PMTH(3,22)
0029       PMTH(5,22)=PMTH(3,22)
0030       PMQTH1=PARJ(82)
0031       IF(MSTJ(41).EQ.2) PMQTH1=MIN(PARJ(82),PARJ(83))
0032       PMQTH2=PMTH(2,21)
0033       IF(MSTJ(41).EQ.2) PMQTH2=MIN(PMTH(2,21),PMTH(2,22))
0034       DO 100 IF=1,8
0035       PMTH(1,IF)=ULMASS(IF)
0036       PMTH(2,IF)=SQRT(PMTH(1,IF)**2+0.25*PMQTH1**2)
0037       PMTH(3,IF)=PMTH(2,IF)+PMQTH2
0038       PMTH(4,IF)=SQRT(PMTH(1,IF)**2+0.25*PARJ(82)**2)+PMTH(2,21)
0039   100 PMTH(5,IF)=SQRT(PMTH(1,IF)**2+0.25*PARJ(83)**2)+PMTH(2,22)
0040       PT2MIN=MAX(0.5*PARJ(82),1.1*PARJ(81))**2
0041       ALAMS=PARJ(81)**2
0042       ALFM=LOG(PT2MIN/ALAMS)
0043
0044 C...Store positions of shower initiating partons.
0045       M3JC=0
0046       IF(IP1.GT.0.AND.IP1.LE.MIN(N,MSTU(4)-MSTU(32)).AND.IP2.EQ.0) THEN
0047         NPA=1
0048         IPA(1)=IP1
0049       ELSEIF(MIN(IP1,IP2).GT.0.AND.MAX(IP1,IP2).LE.MIN(N,MSTU(4)-
0050      &MSTU(32))) THEN
0051         NPA=2
0052         IPA(1)=IP1
0053         IPA(2)=IP2
0054       ELSEIF(IP1.GT.0.AND.IP1.LE.MIN(N,MSTU(4)-MSTU(32)).AND.IP2.LT.0.
0055      &AND.IP2.GE.-3) THEN
0056         NPA=IABS(IP2)
0057         DO 110 I=1,NPA
0058   110   IPA(I)=IP1+I-1
0059       ELSE
0060         CALL LUERRM(12,
0061      &  '(LUSHOW:) failed to reconstruct showering system')
0062         IF(MSTU(21).GE.1) RETURN
0063       ENDIF
0064
0065 C...Check on phase space available for emission.
0066       IREJ=0
0067       DO 120 J=1,5
0068   120 PS(J)=0.
0069       PM=0.
0070       DO 130 I=1,NPA
0071       KFLA(I)=IABS(K(IPA(I),2))
0072       PMA(I)=P(IPA(I),5)
0073       IF(KFLA(I).NE.0.AND.(KFLA(I).LE.8.OR.KFLA(I).EQ.21))
0074      &PMA(I)=PMTH(3,KFLA(I))
0075       PM=PM+PMA(I)
0076       IF(KFLA(I).EQ.0.OR.(KFLA(I).GT.8.AND.KFLA(I).NE.21).OR.
0077      &PMA(I).GT.QMAX) IREJ=IREJ+1
0078       DO 130 J=1,4
0079   130 PS(J)=PS(J)+P(IPA(I),J)
0080       IF(IREJ.EQ.NPA) RETURN
0081       PS(5)=SQRT(MAX(0.,PS(4)**2-PS(1)**2-PS(2)**2-PS(3)**2))
0082       IF(NPA.EQ.1) PS(5)=PS(4)
0083       IF(PS(5).LE.PM+PMQTH1) RETURN
0084       IF(NPA.EQ.2.AND.MSTJ(47).GE.1) THEN
0085         IF(KFLA(1).GE.1.AND.KFLA(1).LE.8.AND.KFLA(2).GE.1.AND.
0086      &  KFLA(2).LE.8) M3JC=1
0087         IF(MSTJ(47).GE.2) M3JC=1
0088       ENDIF
0089
0090 C...Define imagined single initiator of shower for parton system.
0091       NS=N
0092       IF(N.GT.MSTU(4)-MSTU(32)-5) THEN
0093         CALL LUERRM(11,'(LUSHOW:) no more memory left in LUJETS')
0094         IF(MSTU(21).GE.1) RETURN
0095       ENDIF
0096       IF(NPA.GE.2) THEN
0097         K(N+1,1)=11
0098         K(N+1,2)=21
0099         K(N+1,3)=0
0100         K(N+1,4)=0
0101         K(N+1,5)=0
0102         P(N+1,1)=0.
0103         P(N+1,2)=0.
0104         P(N+1,3)=0.
0105         P(N+1,4)=PS(5)
0106         P(N+1,5)=PS(5)
0107         V(N+1,5)=PS(5)**2
0108         N=N+1
0109       ENDIF
0110
0111 C...Loop over partons that may branch.
0112       NEP=NPA
0113       IM=NS
0114       IF(NPA.EQ.1) IM=NS-1
0115   140 IM=IM+1
0116       IF(N.GT.NS) THEN
0117         IF(IM.GT.N) GOTO 380
0118         KFLM=IABS(K(IM,2))
0119         IF(KFLM.EQ.0.OR.(KFLM.GT.8.AND.KFLM.NE.21)) GOTO 140
0120         IF(P(IM,5).LT.PMTH(2,KFLM)) GOTO 140
0121         IGM=K(IM,3)
0122       ELSE
0123         IGM=-1
0124       ENDIF
0125       IF(N+NEP.GT.MSTU(4)-MSTU(32)-5) THEN
0126         CALL LUERRM(11,'(LUSHOW:) no more memory left in LUJETS')
0127         IF(MSTU(21).GE.1) RETURN
0128       ENDIF
0129
0130 C...Position of aunt (sister to branching parton).
0131 C...Origin and flavour of daughters.
0132       IAU=0
0133       IF(IGM.GT.0) THEN
0134         IF(K(IM-1,3).EQ.IGM) IAU=IM-1
0135         IF(N.GE.IM+1.AND.K(IM+1,3).EQ.IGM) IAU=IM+1
0136       ENDIF
0137       IF(IGM.GE.0) THEN
0138         K(IM,4)=N+1
0139         DO 150 I=1,NEP
0140   150   K(N+I,3)=IM
0141       ELSE
0142         K(N+1,3)=IPA(1)
0143       ENDIF
0144       IF(IGM.LE.0) THEN
0145         DO 160 I=1,NEP
0146   160   K(N+I,2)=K(IPA(I),2)
0147       ELSEIF(KFLM.NE.21) THEN
0148         K(N+1,2)=K(IM,2)
0149         K(N+2,2)=K(IM,5)
0150       ELSEIF(K(IM,5).EQ.21) THEN
0151         K(N+1,2)=21
0152         K(N+2,2)=21
0153       ELSE
0154         K(N+1,2)=K(IM,5)
0155         K(N+2,2)=-K(IM,5)
0156       ENDIF
0157
0158 C...Reset flags on daughers and tries made.
0159       DO 170 IP=1,NEP
0160       K(N+IP,1)=3
0161       K(N+IP,4)=0
0162       K(N+IP,5)=0
0163       KFLD(IP)=IABS(K(N+IP,2))
0164       ITRY(IP)=0
0165       ISL(IP)=0
0166       ISI(IP)=0
0167   170 IF(KFLD(IP).GT.0.AND.(KFLD(IP).LE.8.OR.KFLD(IP).EQ.21)) ISI(IP)=1
0168       ISLM=0
0169
0170 C...Maximum virtuality of daughters.
0171       IF(IGM.LE.0) THEN
0172         DO 180 I=1,NPA
0173         IF(NPA.GE.3) P(N+I,4)=(PS(4)*P(IPA(I),4)-PS(1)*P(IPA(I),1)-
0174      &  PS(2)*P(IPA(I),2)-PS(3)*P(IPA(I),3))/PS(5)
0175         P(N+I,5)=MIN(QMAX,PS(5))
0176         IF(NPA.GE.3) P(N+I,5)=MIN(P(N+I,5),P(N+I,4))
0177   180   IF(ISI(I).EQ.0) P(N+I,5)=P(IPA(I),5)
0178       ELSE
0179         IF(MSTJ(43).LE.2) PEM=V(IM,2)
0180         IF(MSTJ(43).GE.3) PEM=P(IM,4)
0181         P(N+1,5)=MIN(P(IM,5),V(IM,1)*PEM)
0182         P(N+2,5)=MIN(P(IM,5),(1.-V(IM,1))*PEM)
0183         IF(K(N+2,2).EQ.22) P(N+2,5)=PMTH(1,22)
0184       ENDIF
0185       DO 190 I=1,NEP
0186       PMSD(I)=P(N+I,5)
0187       IF(ISI(I).EQ.1) THEN
0188         IF(P(N+I,5).LE.PMTH(3,KFLD(I))) P(N+I,5)=PMTH(1,KFLD(I))
0189       ENDIF
0190   190 V(N+I,5)=P(N+I,5)**2
0191
0192 C...Choose one of the daughters for evolution.
0193   200 INUM=0
0194       IF(NEP.EQ.1) INUM=1
0195       DO 210 I=1,NEP
0196   210 IF(INUM.EQ.0.AND.ISL(I).EQ.1) INUM=I
0197       DO 220 I=1,NEP
0198       IF(INUM.EQ.0.AND.ITRY(I).EQ.0.AND.ISI(I).EQ.1) THEN
0199         IF(P(N+I,5).GE.PMTH(2,KFLD(I))) INUM=I
0200       ENDIF
0201   220 CONTINUE
0202       IF(INUM.EQ.0) THEN
0203         RMAX=0.
0204         DO 230 I=1,NEP
0205         IF(ISI(I).EQ.1.AND.PMSD(I).GE.PMQTH2) THEN
0206           RPM=P(N+I,5)/PMSD(I)
0207           IF(RPM.GT.RMAX.AND.P(N+I,5).GE.PMTH(2,KFLD(I))) THEN
0208             RMAX=RPM
0209             INUM=I
0210           ENDIF
0211         ENDIF
0212   230   CONTINUE
0213       ENDIF
0214
0215 C...Store information on choice of evolving daughter.
0216       INUM=MAX(1,INUM)
0217       IEP(1)=N+INUM
0218       DO 240 I=2,NEP
0219       IEP(I)=IEP(I-1)+1
0220   240 IF(IEP(I).GT.N+NEP) IEP(I)=N+1
0221       DO 250 I=1,NEP
0222   250 KFL(I)=IABS(K(IEP(I),2))
0223       ITRY(INUM)=ITRY(INUM)+1
0224       IF(ITRY(INUM).GT.200) THEN
0225         CALL LUERRM(14,'(LUSHOW:) caught in infinite loop')
0226         IF(MSTU(21).GE.1) RETURN
0227       ENDIF
0228       Z=0.5
0229       IF(KFL(1).EQ.0.OR.(KFL(1).GT.8.AND.KFL(1).NE.21)) GOTO 300
0230       IF(P(IEP(1),5).LT.PMTH(2,KFL(1))) GOTO 300
0231
0232 C...Calculate allowed z range.
0233       IF(NEP.EQ.1) THEN
0234         PMED=PS(4)
0235       ELSEIF(IGM.EQ.0.OR.MSTJ(43).LE.2) THEN
0236         PMED=P(IM,5)
0237       ELSE
0238         IF(INUM.EQ.1) PMED=V(IM,1)*PEM
0239         IF(INUM.EQ.2) PMED=(1.-V(IM,1))*PEM
0240       ENDIF
0241       IF(MOD(MSTJ(43),2).EQ.1) THEN
0242         ZC=PMTH(2,21)/PMED
0243         ZCE=PMTH(2,22)/PMED
0244       ELSE
0245         ZC=0.5*(1.-SQRT(MAX(0.,1.-(2.*PMTH(2,21)/PMED)**2)))
0246         IF(ZC.LT.1E-4) ZC=(PMTH(2,21)/PMED)**2
0247         ZCE=0.5*(1.-SQRT(MAX(0.,1.-(2.*PMTH(2,22)/PMED)**2)))
0248         IF(ZCE.LT.1E-4) ZCE=(PMTH(2,22)/PMED)**2
0249       ENDIF
0250       ZC=MIN(ZC,0.491)
0251       ZCE=MIN(ZCE,0.491)
0252       IF((MSTJ(41).EQ.1.AND.ZC.GT.0.49).OR.(MSTJ(41).EQ.2.AND.
0253      &MIN(ZC,ZCE).GT.0.49)) THEN
0254         P(IEP(1),5)=PMTH(1,KFL(1))
0255         V(IEP(1),5)=P(IEP(1),5)**2
0256         GOTO 300
0257       ENDIF
0258
0259 C...Integral of Altarelli-Parisi z kernel for QCD.
0260       IF(MSTJ(49).EQ.0.AND.KFL(1).EQ.21) THEN
0261         FBR=6.*LOG((1.-ZC)/ZC)+MSTJ(45)*(0.5-ZC)
0262       ELSEIF(MSTJ(49).EQ.0) THEN
0263         FBR=(8./3.)*LOG((1.-ZC)/ZC)
0264
0265 C...Integral of Altarelli-Parisi z kernel for scalar gluon.
0266       ELSEIF(MSTJ(49).EQ.1.AND.KFL(1).EQ.21) THEN
0267         FBR=(PARJ(87)+MSTJ(45)*PARJ(88))*(1.-2.*ZC)
0268       ELSEIF(MSTJ(49).EQ.1) THEN
0269         FBR=(1.-2.*ZC)/3.
0270         IF(IGM.EQ.0.AND.M3JC.EQ.1) FBR=4.*FBR
0271
0272 C...Integral of Altarelli-Parisi z kernel for Abelian vector gluon.
0273       ELSEIF(KFL(1).EQ.21) THEN
0274         FBR=6.*MSTJ(45)*(0.5-ZC)
0275       ELSE
0276         FBR=2.*LOG((1.-ZC)/ZC)
0277       ENDIF
0278
0279 C...Integral of Altarelli-Parisi kernel for photon emission.
0280       IF(MSTJ(41).EQ.2.AND.KFL(1).GE.1.AND.KFL(1).LE.8)
0281      &FBRE=(KCHG(KFL(1),1)/3.)**2*2.*LOG((1.-ZCE)/ZCE)
0282
0283 C...Inner veto algorithm starts. Find maximum mass for evolution.
0284   260 PMS=V(IEP(1),5)
0285       IF(IGM.GE.0) THEN
0286         PM2=0.
0287         DO 270 I=2,NEP
0288         PM=P(IEP(I),5)
0289         IF(KFL(I).GT.0.AND.(KFL(I).LE.8.OR.KFL(I).EQ.21)) PM=
0290      &  PMTH(2,KFL(I))
0291   270   PM2=PM2+PM
0292         PMS=MIN(PMS,(P(IM,5)-PM2)**2)
0293       ENDIF
0294
0295 C...Select mass for daughter in QCD evolution.
0296       B0=27./6.
0297       DO 280 IF=4,MSTJ(45)
0298   280 IF(PMS.GT.4.*PMTH(2,IF)**2) B0=(33.-2.*IF)/6.
0299       IF(MSTJ(44).LE.0) THEN
0300         PMSQCD=PMS*EXP(MAX(-100.,LOG(RLU(0))*PARU(2)/(PARU(111)*FBR)))
0301       ELSEIF(MSTJ(44).EQ.1) THEN
0302         PMSQCD=4.*ALAMS*(0.25*PMS/ALAMS)**(RLU(0)**(B0/FBR))
0303       ELSE
0304         PMSQCD=PMS*RLU(0)**(ALFM*B0/FBR)
0305       ENDIF
0306       IF(ZC.GT.0.49.OR.PMSQCD.LE.PMTH(4,KFL(1))**2) PMSQCD=
0307      &PMTH(2,KFL(1))**2
0308       V(IEP(1),5)=PMSQCD
0309       MCE=1
0310
0311 C...Select mass for daughter in QED evolution.
0312       IF(MSTJ(41).EQ.2.AND.KFL(1).GE.1.AND.KFL(1).LE.8) THEN
0313         PMSQED=PMS*EXP(MAX(-100.,LOG(RLU(0))*PARU(2)/(PARU(101)*FBRE)))
0314         IF(ZCE.GT.0.49.OR.PMSQED.LE.PMTH(5,KFL(1))**2) PMSQED=
0315      &  PMTH(2,KFL(1))**2
0316         IF(PMSQED.GT.PMSQCD) THEN
0317           V(IEP(1),5)=PMSQED
0318           MCE=2
0319         ENDIF
0320       ENDIF
0321
0322 C...Check whether daughter mass below cutoff.
0323       P(IEP(1),5)=SQRT(V(IEP(1),5))
0324       IF(P(IEP(1),5).LE.PMTH(3,KFL(1))) THEN
0325         P(IEP(1),5)=PMTH(1,KFL(1))
0326         V(IEP(1),5)=P(IEP(1),5)**2
0327         GOTO 300
0328       ENDIF
0329
0330 C...Select z value of branching: q -> qgamma.
0331       IF(MCE.EQ.2) THEN
0332         Z=1.-(1.-ZCE)*(ZCE/(1.-ZCE))**RLU(0)
0333         IF(1.+Z**2.LT.2.*RLU(0)) GOTO 260
0334         K(IEP(1),5)=22
0335
0336 C...Select z value of branching: q -> qg, g -> gg, g -> qqbar.
0337       ELSEIF(MSTJ(49).NE.1.AND.KFL(1).NE.21) THEN
0338         Z=1.-(1.-ZC)*(ZC/(1.-ZC))**RLU(0)
0339         IF(1.+Z**2.LT.2.*RLU(0)) GOTO 260
0340         K(IEP(1),5)=21
0341       ELSEIF(MSTJ(49).EQ.0.AND.MSTJ(45)*(0.5-ZC).LT.RLU(0)*FBR) THEN
0342         Z=(1.-ZC)*(ZC/(1.-ZC))**RLU(0)
0343         IF(RLU(0).GT.0.5) Z=1.-Z
0344         IF((1.-Z*(1.-Z))**2.LT.RLU(0)) GOTO 260
0345         K(IEP(1),5)=21
0346       ELSEIF(MSTJ(49).NE.1) THEN
0347         Z=ZC+(1.-2.*ZC)*RLU(0)
0348         IF(Z**2+(1.-Z)**2.LT.RLU(0)) GOTO 260
0349         KFLB=1+INT(MSTJ(45)*RLU(0))
0350         PMQ=4.*PMTH(2,KFLB)**2/V(IEP(1),5)
0351         IF(PMQ.GE.1.) GOTO 260
0352         PMQ0=4.*PMTH(2,21)**2/V(IEP(1),5)
0353         IF(MOD(MSTJ(43),2).EQ.0.AND.(1.+0.5*PMQ)*SQRT(1.-PMQ).LT.
0354      &  RLU(0)*(1.+0.5*PMQ0)*SQRT(1.-PMQ0)) GOTO 260
0355         K(IEP(1),5)=KFLB
0356
0357 C...Ditto for scalar gluon model.
0358       ELSEIF(KFL(1).NE.21) THEN
0359         Z=1.-SQRT(ZC**2+RLU(0)*(1.-2.*ZC))
0360         K(IEP(1),5)=21
0361       ELSEIF(RLU(0)*(PARJ(87)+MSTJ(45)*PARJ(88)).LE.PARJ(87)) THEN
0362         Z=ZC+(1.-2.*ZC)*RLU(0)
0363         K(IEP(1),5)=21
0364       ELSE
0365         Z=ZC+(1.-2.*ZC)*RLU(0)
0366         KFLB=1+INT(MSTJ(45)*RLU(0))
0367         PMQ=4.*PMTH(2,KFLB)**2/V(IEP(1),5)
0368         IF(PMQ.GE.1.) GOTO 260
0369         K(IEP(1),5)=KFLB
0370       ENDIF
0371       IF(MCE.EQ.1.AND.MSTJ(44).GE.2) THEN
0372         IF(Z*(1.-Z)*V(IEP(1),5).LT.PT2MIN) GOTO 260
0373         IF(ALFM/LOG(V(IEP(1),5)*Z*(1.-Z)/ALAMS).LT.RLU(0)) GOTO 260
0374       ENDIF
0375
0376 C...Check if z consistent with chosen m.
0377       IF(KFL(1).EQ.21) THEN
0378         KFLGD1=IABS(K(IEP(1),5))
0379         KFLGD2=KFLGD1
0380       ELSE
0381         KFLGD1=KFL(1)
0382         KFLGD2=IABS(K(IEP(1),5))
0383       ENDIF
0384       IF(NEP.EQ.1) THEN
0385         PED=PS(4)
0386       ELSEIF(NEP.GE.3) THEN
0387         PED=P(IEP(1),4)
0388       ELSEIF(IGM.EQ.0.OR.MSTJ(43).LE.2) THEN
0389         PED=0.5*(V(IM,5)+V(IEP(1),5)-PM2**2)/P(IM,5)
0390       ELSE
0391         IF(IEP(1).EQ.N+1) PED=V(IM,1)*PEM
0392         IF(IEP(1).EQ.N+2) PED=(1.-V(IM,1))*PEM
0393       ENDIF
0394       IF(MOD(MSTJ(43),2).EQ.1) THEN
0395         PMQTH3=0.5*PARJ(82)
0396         IF(KFLGD2.EQ.22) PMQTH3=0.5*PARJ(83)
0397         PMQ1=(PMTH(1,KFLGD1)**2+PMQTH3**2)/V(IEP(1),5)
0398         PMQ2=(PMTH(1,KFLGD2)**2+PMQTH3**2)/V(IEP(1),5)
0399         ZD=SQRT(MAX(0.,(1.-V(IEP(1),5)/PED**2)*((1.-PMQ1-PMQ2)**2-
0400      &  4.*PMQ1*PMQ2)))
0401         ZH=1.+PMQ1-PMQ2
0402       ELSE
0403         ZD=SQRT(MAX(0.,1.-V(IEP(1),5)/PED**2))
0404         ZH=1.
0405       ENDIF
0406       ZL=0.5*(ZH-ZD)
0407       ZU=0.5*(ZH+ZD)
0408       IF(Z.LT.ZL.OR.Z.GT.ZU) GOTO 260
0409       IF(KFL(1).EQ.21) V(IEP(1),3)=LOG(ZU*(1.-ZL)/MAX(1E-20,ZL*
0410      &(1.-ZU)))
0411       IF(KFL(1).NE.21) V(IEP(1),3)=LOG((1.-ZL)/MAX(1E-10,1.-ZU))
0412
0413 C...Three-jet matrix element correction.
0414       IF(IGM.EQ.0.AND.M3JC.EQ.1) THEN
0415         X1=Z*(1.+V(IEP(1),5)/V(NS+1,5))
0416         X2=1.-V(IEP(1),5)/V(NS+1,5)
0417         X3=(1.-X1)+(1.-X2)
0418         IF(MCE.EQ.2) THEN
0419           KI1=K(IPA(INUM),2)
0420           KI2=K(IPA(3-INUM),2)
0421           QF1=KCHG(IABS(KI1),1)*ISIGN(1,KI1)/3.
0422           QF2=KCHG(IABS(KI2),1)*ISIGN(1,KI2)/3.
0423           WSHOW=QF1**2*(1.-X1)/X3*(1.+(X1/(2.-X2))**2)+
0424      &    QF2**2*(1.-X2)/X3*(1.+(X2/(2.-X1))**2)
0425           WME=(QF1*(1.-X1)/X3-QF2*(1.-X2)/X3)**2*(X1**2+X2**2)
0426         ELSEIF(MSTJ(49).NE.1) THEN
0427           WSHOW=1.+(1.-X1)/X3*(X1/(2.-X2))**2+
0428      &    (1.-X2)/X3*(X2/(2.-X1))**2
0429           WME=X1**2+X2**2
0430         ELSE
0431           WSHOW=4.*X3*((1.-X1)/(2.-X2)**2+(1.-X2)/(2.-X1)**2)
0432           WME=X3**2
0433         ENDIF
0434         IF(WME.LT.RLU(0)*WSHOW) GOTO 260
0435
0436 C...Impose angular ordering by rejection of nonordered emission.
0437       ELSEIF(MCE.EQ.1.AND.IGM.GT.0.AND.MSTJ(42).GE.2) THEN
0438         MAOM=1
0439         ZM=V(IM,1)
0440         IF(IEP(1).EQ.N+2) ZM=1.-V(IM,1)
0441         THE2ID=Z*(1.-Z)*(ZM*P(IM,4))**2/V(IEP(1),5)
0442         IAOM=IM
0443   290   IF(K(IAOM,5).EQ.22) THEN
0444           IAOM=K(IAOM,3)
0445           IF(K(IAOM,3).LE.NS) MAOM=0
0446           IF(MAOM.EQ.1) GOTO 290
0447         ENDIF
0448         IF(MAOM.EQ.1) THEN
0449           THE2IM=V(IAOM,1)*(1.-V(IAOM,1))*P(IAOM,4)**2/V(IAOM,5)
0450           IF(THE2ID.LT.THE2IM) GOTO 260
0451         ENDIF
0452       ENDIF
0453
0454 C...Impose user-defined maximum angle at first branching.
0455       IF(MSTJ(48).EQ.1) THEN
0456         IF(NEP.EQ.1.AND.IM.EQ.NS) THEN
0457           THE2ID=Z*(1.-Z)*PS(4)**2/V(IEP(1),5)
0458           IF(THE2ID.LT.1./PARJ(85)**2) GOTO 260
0459         ELSEIF(NEP.EQ.2.AND.IEP(1).EQ.NS+2) THEN
0460           THE2ID=Z*(1.-Z)*(0.5*P(IM,4))**2/V(IEP(1),5)
0461           IF(THE2ID.LT.1./PARJ(85)**2) GOTO 260
0462         ELSEIF(NEP.EQ.2.AND.IEP(1).EQ.NS+3) THEN
0463           THE2ID=Z*(1.-Z)*(0.5*P(IM,4))**2/V(IEP(1),5)
0464           IF(THE2ID.LT.1./PARJ(86)**2) GOTO 260
0465         ENDIF
0466       ENDIF
0467
0468 C...End of inner veto algorithm. Check if only one leg evolved so far.
0469   300 V(IEP(1),1)=Z
0470       ISL(1)=0
0471       ISL(2)=0
0472       IF(NEP.EQ.1) GOTO 330
0473       IF(NEP.EQ.2.AND.P(IEP(1),5)+P(IEP(2),5).GE.P(IM,5)) GOTO 200
0474       DO 310 I=1,NEP
0475       IF(ITRY(I).EQ.0.AND.KFLD(I).GT.0.AND.(KFLD(I).LE.8.OR.KFLD(I).EQ.
0476      &21)) THEN
0477         IF(P(N+I,5).GE.PMTH(2,KFLD(I))) GOTO 200
0478       ENDIF
0479   310 CONTINUE
0480
0481 C...Check if chosen multiplet m1,m2,z1,z2 is physical.
0482       IF(NEP.EQ.3) THEN
0483         PA1S=(P(N+1,4)+P(N+1,5))*(P(N+1,4)-P(N+1,5))
0484         PA2S=(P(N+2,4)+P(N+2,5))*(P(N+2,4)-P(N+2,5))
0485         PA3S=(P(N+3,4)+P(N+3,5))*(P(N+3,4)-P(N+3,5))
0486         PTS=0.25*(2.*PA1S*PA2S+2.*PA1S*PA3S+2.*PA2S*PA3S-
0487      &  PA1S**2-PA2S**2-PA3S**2)/PA1S
0488         IF(PTS.LE.0.) GOTO 200
0489       ELSEIF(IGM.EQ.0.OR.MSTJ(43).LE.2.OR.MOD(MSTJ(43),2).EQ.0) THEN
0490         DO 320 I1=N+1,N+2
0491         KFLDA=IABS(K(I1,2))
0492         IF(KFLDA.EQ.0.OR.(KFLDA.GT.8.AND.KFLDA.NE.21)) GOTO 320
0493         IF(P(I1,5).LT.PMTH(2,KFLDA)) GOTO 320
0494         IF(KFLDA.EQ.21) THEN
0495           KFLGD1=IABS(K(I1,5))
0496           KFLGD2=KFLGD1
0497         ELSE
0498           KFLGD1=KFLDA
0499           KFLGD2=IABS(K(I1,5))
0500         ENDIF
0501         I2=2*N+3-I1
0502         IF(IGM.EQ.0.OR.MSTJ(43).LE.2) THEN
0503           PED=0.5*(V(IM,5)+V(I1,5)-V(I2,5))/P(IM,5)
0504         ELSE
0505           IF(I1.EQ.N+1) ZM=V(IM,1)
0506           IF(I1.EQ.N+2) ZM=1.-V(IM,1)
0507           PML=SQRT((V(IM,5)-V(N+1,5)-V(N+2,5))**2-
0508      &    4.*V(N+1,5)*V(N+2,5))
0509           PED=PEM*(0.5*(V(IM,5)-PML+V(I1,5)-V(I2,5))+PML*ZM)/V(IM,5)
0510         ENDIF
0511         IF(MOD(MSTJ(43),2).EQ.1) THEN
0512           PMQTH3=0.5*PARJ(82)
0513           IF(KFLGD2.EQ.22) PMQTH3=0.5*PARJ(83)
0514           PMQ1=(PMTH(1,KFLGD1)**2+PMQTH3**2)/V(I1,5)
0515           PMQ2=(PMTH(1,KFLGD2)**2+PMQTH3**2)/V(I1,5)
0516           ZD=SQRT(MAX(0.,(1.-V(I1,5)/PED**2)*((1.-PMQ1-PMQ2)**2-
0517      &    4.*PMQ1*PMQ2)))
0518           ZH=1.+PMQ1-PMQ2
0519         ELSE
0520           ZD=SQRT(MAX(0.,1.-V(I1,5)/PED**2))
0521           ZH=1.
0522         ENDIF
0523         ZL=0.5*(ZH-ZD)
0524         ZU=0.5*(ZH+ZD)
0525         IF(I1.EQ.N+1.AND.(V(I1,1).LT.ZL.OR.V(I1,1).GT.ZU)) ISL(1)=1
0526         IF(I1.EQ.N+2.AND.(V(I1,1).LT.ZL.OR.V(I1,1).GT.ZU)) ISL(2)=1
0527         IF(KFLDA.EQ.21) V(I1,4)=LOG(ZU*(1.-ZL)/MAX(1E-20,ZL*(1.-ZU)))
0528         IF(KFLDA.NE.21) V(I1,4)=LOG((1.-ZL)/MAX(1E-10,1.-ZU))
0529   320   CONTINUE
0530         IF(ISL(1).EQ.1.AND.ISL(2).EQ.1.AND.ISLM.NE.0) THEN
0531           ISL(3-ISLM)=0
0532           ISLM=3-ISLM
0533         ELSEIF(ISL(1).EQ.1.AND.ISL(2).EQ.1) THEN
0534           ZDR1=MAX(0.,V(N+1,3)/V(N+1,4)-1.)
0535           ZDR2=MAX(0.,V(N+2,3)/V(N+2,4)-1.)
0536           IF(ZDR2.GT.RLU(0)*(ZDR1+ZDR2)) ISL(1)=0
0537           IF(ISL(1).EQ.1) ISL(2)=0
0538           IF(ISL(1).EQ.0) ISLM=1
0539           IF(ISL(2).EQ.0) ISLM=2
0540         ENDIF
0541         IF(ISL(1).EQ.1.OR.ISL(2).EQ.1) GOTO 200
0542       ENDIF
0543       IF(IGM.GT.0.AND.MOD(MSTJ(43),2).EQ.1.AND.(P(N+1,5).GE.
0544      &PMTH(2,KFLD(1)).OR.P(N+2,5).GE.PMTH(2,KFLD(2)))) THEN
0545         PMQ1=V(N+1,5)/V(IM,5)
0546         PMQ2=V(N+2,5)/V(IM,5)
0547         ZD=SQRT(MAX(0.,(1.-V(IM,5)/PEM**2)*((1.-PMQ1-PMQ2)**2-
0548      &  4.*PMQ1*PMQ2)))
0549         ZH=1.+PMQ1-PMQ2
0550         ZL=0.5*(ZH-ZD)
0551         ZU=0.5*(ZH+ZD)
0552         IF(V(IM,1).LT.ZL.OR.V(IM,1).GT.ZU) GOTO 200
0553       ENDIF
0554
0555 C...Accepted branch. Construct four-momentum for initial partons.
0556   330 MAZIP=0
0557       MAZIC=0
0558       IF(NEP.EQ.1) THEN
0559         P(N+1,1)=0.
0560         P(N+1,2)=0.
0561         P(N+1,3)=SQRT(MAX(0.,(P(IPA(1),4)+P(N+1,5))*(P(IPA(1),4)-
0562      &  P(N+1,5))))
0563         P(N+1,4)=P(IPA(1),4)
0564         V(N+1,2)=P(N+1,4)
0565       ELSEIF(IGM.EQ.0.AND.NEP.EQ.2) THEN
0566         PED1=0.5*(V(IM,5)+V(N+1,5)-V(N+2,5))/P(IM,5)
0567         P(N+1,1)=0.
0568         P(N+1,2)=0.
0569         P(N+1,3)=SQRT(MAX(0.,(PED1+P(N+1,5))*(PED1-P(N+1,5))))
0570         P(N+1,4)=PED1
0571         P(N+2,1)=0.
0572         P(N+2,2)=0.
0573         P(N+2,3)=-P(N+1,3)
0574         P(N+2,4)=P(IM,5)-PED1
0575         V(N+1,2)=P(N+1,4)
0576         V(N+2,2)=P(N+2,4)
0577       ELSEIF(NEP.EQ.3) THEN
0578         P(N+1,1)=0.
0579         P(N+1,2)=0.
0580         P(N+1,3)=SQRT(MAX(0.,PA1S))
0581         P(N+2,1)=SQRT(PTS)
0582         P(N+2,2)=0.
0583         P(N+2,3)=0.5*(PA3S-PA2S-PA1S)/P(N+1,3)
0584         P(N+3,1)=-P(N+2,1)
0585         P(N+3,2)=0.
0586         P(N+3,3)=-(P(N+1,3)+P(N+2,3))
0587         V(N+1,2)=P(N+1,4)
0588         V(N+2,2)=P(N+2,4)
0589         V(N+3,2)=P(N+3,4)
0590
0591 C...Construct transverse momentum for ordinary branching in shower.
0592       ELSE
0593         ZM=V(IM,1)
0594         PZM=SQRT(MAX(0.,(PEM+P(IM,5))*(PEM-P(IM,5))))
0595         PMLS=(V(IM,5)-V(N+1,5)-V(N+2,5))**2-4.*V(N+1,5)*V(N+2,5)
0596         IF(PZM.LE.0.) THEN
0597           PTS=0.
0598         ELSEIF(MOD(MSTJ(43),2).EQ.1) THEN
0599           PTS=(PEM**2*(ZM*(1.-ZM)*V(IM,5)-(1.-ZM)*V(N+1,5)-
0600      &    ZM*V(N+2,5))-0.25*PMLS)/PZM**2
0601         ELSE
0602           PTS=PMLS*(ZM*(1.-ZM)*PEM**2/V(IM,5)-0.25)/PZM**2
0603         ENDIF
0604         PT=SQRT(MAX(0.,PTS))
0605
0606 C...Find coefficient of azimuthal asymmetry due to gluon polarization.
0607         HAZIP=0.
0608         IF(MSTJ(49).NE.1.AND.MOD(MSTJ(46),2).EQ.1.AND.K(IM,2).EQ.21.
0609      &  AND.IAU.NE.0) THEN
0610           IF(K(IGM,3).NE.0) MAZIP=1
0611           ZAU=V(IGM,1)
0612           IF(IAU.EQ.IM+1) ZAU=1.-V(IGM,1)
0613           IF(MAZIP.EQ.0) ZAU=0.
0614           IF(K(IGM,2).NE.21) THEN
0615             HAZIP=2.*ZAU/(1.+ZAU**2)
0616           ELSE
0617             HAZIP=(ZAU/(1.-ZAU*(1.-ZAU)))**2
0618           ENDIF
0619           IF(K(N+1,2).NE.21) THEN
0620             HAZIP=HAZIP*(-2.*ZM*(1.-ZM))/(1.-2.*ZM*(1.-ZM))
0621           ELSE
0622             HAZIP=HAZIP*(ZM*(1.-ZM)/(1.-ZM*(1.-ZM)))**2
0623           ENDIF
0624         ENDIF
0625
0626 C...Find coefficient of azimuthal asymmetry due to soft gluon
0627 C...interference.
0628         HAZIC=0.
0629         IF(MSTJ(46).GE.2.AND.(K(N+1,2).EQ.21.OR.K(N+2,2).EQ.21).
0630      &  AND.IAU.NE.0) THEN
0631           IF(K(IGM,3).NE.0) MAZIC=N+1
0632           IF(K(IGM,3).NE.0.AND.K(N+1,2).NE.21) MAZIC=N+2
0633           IF(K(IGM,3).NE.0.AND.K(N+1,2).EQ.21.AND.K(N+2,2).EQ.21.AND.
0634      &    ZM.GT.0.5) MAZIC=N+2
0635           IF(K(IAU,2).EQ.22) MAZIC=0
0636           ZS=ZM
0637           IF(MAZIC.EQ.N+2) ZS=1.-ZM
0638           ZGM=V(IGM,1)
0639           IF(IAU.EQ.IM-1) ZGM=1.-V(IGM,1)
0640           IF(MAZIC.EQ.0) ZGM=1.
0641           HAZIC=(P(IM,5)/P(IGM,5))*SQRT((1.-ZS)*(1.-ZGM)/(ZS*ZGM))
0642           HAZIC=MIN(0.95,HAZIC)
0643         ENDIF
0644       ENDIF
0645
0646 C...Construct kinematics for ordinary branching in shower.
0647   340 IF(NEP.EQ.2.AND.IGM.GT.0) THEN
0648         IF(MOD(MSTJ(43),2).EQ.1) THEN
0649           P(N+1,4)=PEM*V(IM,1)
0650         ELSE
0651           P(N+1,4)=PEM*(0.5*(V(IM,5)-SQRT(PMLS)+V(N+1,5)-V(N+2,5))+
0652      &    SQRT(PMLS)*ZM)/V(IM,5)
0653         ENDIF
0654         PHI=PARU(2)*RLU(0)
0655         P(N+1,1)=PT*COS(PHI)
0656         P(N+1,2)=PT*SIN(PHI)
0657         IF(PZM.GT.0.) THEN
0658           P(N+1,3)=0.5*(V(N+2,5)-V(N+1,5)-V(IM,5)+2.*PEM*P(N+1,4))/PZM
0659         ELSE
0660           P(N+1,3)=0.
0661         ENDIF
0662         P(N+2,1)=-P(N+1,1)
0663         P(N+2,2)=-P(N+1,2)
0664         P(N+2,3)=PZM-P(N+1,3)
0665         P(N+2,4)=PEM-P(N+1,4)
0666         IF(MSTJ(43).LE.2) THEN
0667           V(N+1,2)=(PEM*P(N+1,4)-PZM*P(N+1,3))/P(IM,5)
0668           V(N+2,2)=(PEM*P(N+2,4)-PZM*P(N+2,3))/P(IM,5)
0669         ENDIF
0670       ENDIF
0671
0672 C...Rotate and boost daughters.
0673       IF(IGM.GT.0) THEN
0674         IF(MSTJ(43).LE.2) THEN
0675           BEX=P(IGM,1)/P(IGM,4)
0676           BEY=P(IGM,2)/P(IGM,4)
0677           BEZ=P(IGM,3)/P(IGM,4)
0678           GA=P(IGM,4)/P(IGM,5)
0679           GABEP=GA*(GA*(BEX*P(IM,1)+BEY*P(IM,2)+BEZ*P(IM,3))/(1.+GA)-
0680      &    P(IM,4))
0681         ELSE
0682           BEX=0.
0683           BEY=0.
0684           BEZ=0.
0685           GA=1.
0686           GABEP=0.
0687         ENDIF
0688         THE=ULANGL(P(IM,3)+GABEP*BEZ,SQRT((P(IM,1)+GABEP*BEX)**2+
0689      &  (P(IM,2)+GABEP*BEY)**2))
0690         PHI=ULANGL(P(IM,1)+GABEP*BEX,P(IM,2)+GABEP*BEY)
0691         DO 350 I=N+1,N+2
0692         DP(1)=COS(THE)*COS(PHI)*P(I,1)-SIN(PHI)*P(I,2)+
0693      &  SIN(THE)*COS(PHI)*P(I,3)
0694         DP(2)=COS(THE)*SIN(PHI)*P(I,1)+COS(PHI)*P(I,2)+
0695      &  SIN(THE)*SIN(PHI)*P(I,3)
0696         DP(3)=-SIN(THE)*P(I,1)+COS(THE)*P(I,3)
0697         DP(4)=P(I,4)
0698         DBP=BEX*DP(1)+BEY*DP(2)+BEZ*DP(3)
0699         DGABP=GA*(GA*DBP/(1D0+GA)+DP(4))
0700         P(I,1)=DP(1)+DGABP*BEX
0701         P(I,2)=DP(2)+DGABP*BEY
0702         P(I,3)=DP(3)+DGABP*BEZ
0703   350   P(I,4)=GA*(DP(4)+DBP)
0704       ENDIF
0705
0706 C...Weight with azimuthal distribution, if required.
0707       IF(MAZIP.NE.0.OR.MAZIC.NE.0) THEN
0708         DO 360 J=1,3
0709         DPT(1,J)=P(IM,J)
0710         DPT(2,J)=P(IAU,J)
0711   360   DPT(3,J)=P(N+1,J)
0712         DPMA=DPT(1,1)*DPT(2,1)+DPT(1,2)*DPT(2,2)+DPT(1,3)*DPT(2,3)
0713         DPMD=DPT(1,1)*DPT(3,1)+DPT(1,2)*DPT(3,2)+DPT(1,3)*DPT(3,3)
0714         DPMM=DPT(1,1)**2+DPT(1,2)**2+DPT(1,3)**2
0715         DO 370 J=1,3
0716         DPT(4,J)=DPT(2,J)-DPMA*DPT(1,J)/DPMM
0717   370   DPT(5,J)=DPT(3,J)-DPMD*DPT(1,J)/DPMM
0718         DPT(4,4)=SQRT(DPT(4,1)**2+DPT(4,2)**2+DPT(4,3)**2)
0719         DPT(5,4)=SQRT(DPT(5,1)**2+DPT(5,2)**2+DPT(5,3)**2)
0720         IF(MIN(DPT(4,4),DPT(5,4)).GT.0.1*PARJ(82)) THEN
0721           CAD=(DPT(4,1)*DPT(5,1)+DPT(4,2)*DPT(5,2)+
0722      &    DPT(4,3)*DPT(5,3))/(DPT(4,4)*DPT(5,4))
0723           IF(MAZIP.NE.0) THEN
0724             IF(1.+HAZIP*(2.*CAD**2-1.).LT.RLU(0)*(1.+ABS(HAZIP)))
0725      &      GOTO 340
0726           ENDIF
0727           IF(MAZIC.NE.0) THEN
0728             IF(MAZIC.EQ.N+2) CAD=-CAD
0729             IF((1.-HAZIC)*(1.-HAZIC*CAD)/(1.+HAZIC**2-2.*HAZIC*CAD).
0730      &      LT.RLU(0)) GOTO 340
0731           ENDIF
0732         ENDIF
0733       ENDIF
0734
0735 C...Continue loop over partons that may branch, until none left.
0736       IF(IGM.GE.0) K(IM,1)=14
0737       N=N+NEP
0738       NEP=2
0739       IF(N.GT.MSTU(4)-MSTU(32)-5) THEN
0740         CALL LUERRM(11,'(LUSHOW:) no more memory left in LUJETS')
0741         IF(MSTU(21).GE.1) N=NS
0742         IF(MSTU(21).GE.1) RETURN
0743       ENDIF
0744       GOTO 140
0745
0746 C...Set information on imagined shower initiator.
0747   380 IF(NPA.GE.2) THEN
0748         K(NS+1,1)=11
0749         K(NS+1,2)=94
0750         K(NS+1,3)=IP1
0751         IF(IP2.GT.0.AND.IP2.LT.IP1) K(NS+1,3)=IP2
0752         K(NS+1,4)=NS+2
0753         K(NS+1,5)=NS+1+NPA
0754         IIM=1
0755       ELSE
0756         IIM=0
0757       ENDIF
0758
0759 C...Reconstruct string drawing information.
0760       DO 390 I=NS+1+IIM,N
0761       IF(K(I,1).LE.10.AND.K(I,2).EQ.22) THEN
0762         K(I,1)=1
0763       ELSEIF(K(I,1).LE.10) THEN
0764         K(I,4)=MSTU(5)*(K(I,4)/MSTU(5))
0765         K(I,5)=MSTU(5)*(K(I,5)/MSTU(5))
0766       ELSEIF(K(MOD(K(I,4),MSTU(5))+1,2).NE.22) THEN
0767         ID1=MOD(K(I,4),MSTU(5))
0768         IF(K(I,2).GE.1.AND.K(I,2).LE.8) ID1=MOD(K(I,4),MSTU(5))+1
0769         ID2=2*MOD(K(I,4),MSTU(5))+1-ID1
0770         K(I,4)=MSTU(5)*(K(I,4)/MSTU(5))+ID1
0771         K(I,5)=MSTU(5)*(K(I,5)/MSTU(5))+ID2
0772         K(ID1,4)=K(ID1,4)+MSTU(5)*I
0773         K(ID1,5)=K(ID1,5)+MSTU(5)*ID2
0774         K(ID2,4)=K(ID2,4)+MSTU(5)*ID1
0775         K(ID2,5)=K(ID2,5)+MSTU(5)*I
0776       ELSE
0777         ID1=MOD(K(I,4),MSTU(5))
0778         ID2=ID1+1
0779         K(I,4)=MSTU(5)*(K(I,4)/MSTU(5))+ID1
0780         K(I,5)=MSTU(5)*(K(I,5)/MSTU(5))+ID1
0781         K(ID1,4)=K(ID1,4)+MSTU(5)*I
0782         K(ID1,5)=K(ID1,5)+MSTU(5)*I
0783         K(ID2,4)=0
0784         K(ID2,5)=0
0785       ENDIF
0786   390 CONTINUE
0787
0788 C...Transformation from CM frame.
0789       IF(NPA.GE.2) THEN
0790         BEX=PS(1)/PS(4)
0791         BEY=PS(2)/PS(4)
0792         BEZ=PS(3)/PS(4)
0793         GA=PS(4)/PS(5)
0794         GABEP=GA*(GA*(BEX*P(IPA(1),1)+BEY*P(IPA(1),2)+BEZ*P(IPA(1),3))
0795      &  /(1.+GA)-P(IPA(1),4))
0796       ELSE
0797         BEX=0.
0798         BEY=0.
0799         BEZ=0.
0800         GABEP=0.
0801       ENDIF
0802       THE=ULANGL(P(IPA(1),3)+GABEP*BEZ,SQRT((P(IPA(1),1)
0803      &+GABEP*BEX)**2+(P(IPA(1),2)+GABEP*BEY)**2))
0804       PHI=ULANGL(P(IPA(1),1)+GABEP*BEX,P(IPA(1),2)+GABEP*BEY)
0805       IF(NPA.EQ.3) THEN
0806         CHI=ULANGL(COS(THE)*COS(PHI)*(P(IPA(2),1)+GABEP*BEX)+COS(THE)*
0807      &  SIN(PHI)*(P(IPA(2),2)+GABEP*BEY)-SIN(THE)*(P(IPA(2),3)+GABEP*
0808      &  BEZ),-SIN(PHI)*(P(IPA(2),1)+GABEP*BEX)+COS(PHI)*(P(IPA(2),2)+
0809      &  GABEP*BEY))
0810         CALL LUDBRB(NS+1,N,0.,CHI,0D0,0D0,0D0)
0811       ENDIF
0812       DBEX=DBLE(BEX)
0813       DBEY=DBLE(BEY)
0814       DBEZ=DBLE(BEZ)
0815       CALL LUDBRB(NS+1,N,THE,PHI,DBEX,DBEY,DBEZ)
0816
0817 C...Decay vertex of shower.
0818       DO 400 I=NS+1,N
0819       DO 400 J=1,5
0820   400 V(I,J)=V(IP1,J)
0821
0822 C...Delete trivial shower, else connect initiators.
0823       IF(N.EQ.NS+NPA+IIM) THEN
0824         N=NS
0825       ELSE
0826         DO 410 IP=1,NPA
0827         K(IPA(IP),1)=14
0828         K(IPA(IP),4)=K(IPA(IP),4)+NS+IIM+IP
0829         K(IPA(IP),5)=K(IPA(IP),5)+NS+IIM+IP
0830         K(NS+IIM+IP,3)=IPA(IP)
0831         IF(IIM.EQ.1.AND.MSTU(16).NE.2) K(NS+IIM+IP,3)=NS+1
0832         K(NS+IIM+IP,4)=MSTU(5)*IPA(IP)+K(NS+IIM+IP,4)
0833   410   K(NS+IIM+IP,5)=MSTU(5)*IPA(IP)+K(NS+IIM+IP,5)
0834       ENDIF
0835
0836       RETURN
0837       END