pythiaeRHIC/pythia/pysspa.f

0001
0002 C*********************************************************************
0003
0004 C...PYSSPA
0005 C...Generates spacelike parton showers.
0006
0007       SUBROUTINE PYSSPA(IPU1,IPU2)
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/PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200)
0018       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0019       COMMON/PYINT1/MINT(400),VINT(400)
0020       COMMON/PYINT2/ISET(500),KFPR(500,2),COEF(500,20),ICOL(40,4,2)
0021       COMMON/PYINT3/XSFX(2,-40:40),ISIG(1000,3),SIGH(1000)
0022       SAVE /PYJETS/,/PYDAT1/,/PYDAT2/,/PYSUBS/,/PYPARS/,/PYINT1/,
0023      &/PYINT2/,/PYINT3/
0024 C...Local arrays and data.
0025       DIMENSION KFLS(4),IS(2),XS(2),ZS(2),Q2S(2),TEVCSV(2),TEVESV(2),
0026      &XFS(2,-25:25),XFA(-25:25),XFB(-25:25),XFN(-25:25),WTAPC(-25:25),
0027      &WTAPE(-25:25),WTSF(-25:25),THE2(2),ALAM(2),DQ2(3),DPC(3),DPD(4),
0028      &DPB(4),ROBO(5),MORE(2),KFBEAM(2),Q2MNCS(2),KCFI(2),NFIS(2),
0029      &THEFIS(2,2),ISFI(2),DPHI(2),MCESV(2)
0030       DATA IS/2*0/
0031
0032 C...Read out basic information; set global Q^2 scale.
0033       IPUS1=IPU1
0034       IPUS2=IPU2
0035       ISUB=MINT(1)
0036       Q2MX=VINT(56)
0037       VINT2R=VINT(2)*VINT(143)*VINT(144)
0038       IF(ISET(ISUB).EQ.2.OR.ISET(ISUB).EQ.9.OR.ISET(ISUB).EQ.11) Q2MX=
0039      &MIN(VINT2R,PARP(67)*VINT(56))
0040       FCQ2MX=1D0
0041
0042 C...Define which processes ME corrections have been implemented for.
0043       MECOR=0
0044       IF(MSTP(68).EQ.1.OR.MSTP(68).EQ.3) THEN
0045         IF(ISUB.EQ.1.OR.ISUB.EQ.2.OR.ISUB.EQ.141.OR.ISUB.EQ.142.OR.
0046      &  ISUB.EQ.144) MECOR=1
0047         IF(ISUB.EQ.102.OR.ISUB.EQ.152.OR.ISUB.EQ.157) MECOR=2
0048         IF(ISUB.EQ.3.OR.ISUB.EQ.151.OR.ISUB.EQ.156) MECOR=3
0049       ENDIF
0050
0051 C...Initialize QCD evolution and check phase space.
0052       Q2MNC=PARP(62)**2
0053       Q2MNCS(1)=Q2MNC
0054       Q2MNCS(2)=Q2MNC
0055       IF(MINT(107).EQ.2.AND.MSTP(66).EQ.2) THEN
0056         Q0S=PARP(15)**2
0057         PS=VINT(3)**2
0058         Q2EFF=VINT(54)*((Q0S+PS)/(VINT(54)+PS))*
0059      &  EXP(PS*(VINT(54)-Q0S)/((VINT(54)+PS)*(Q0S+PS)))
0060         Q2INT=SQRT(Q0S*Q2EFF)
0061         Q2MNCS(1)=MAX(Q2MNC,Q2INT)
0062       ELSEIF(MINT(107).EQ.3.AND.MSTP(66).GE.1) THEN
0063         Q2MNCS(1)=MAX(Q2MNC,VINT(283))
0064       ENDIF
0065       IF(MINT(108).EQ.2.AND.MSTP(66).EQ.2) THEN
0066         Q0S=PARP(15)**2
0067         PS=VINT(4)**2
0068         Q2EFF=VINT(54)*((Q0S+PS)/(VINT(54)+PS))*
0069      &  EXP(PS*(VINT(54)-Q0S)/((VINT(54)+PS)*(Q0S+PS)))
0070         Q2INT=SQRT(Q0S*Q2EFF)
0071         Q2MNCS(2)=MAX(Q2MNC,Q2INT)
0072       ELSEIF(MINT(108).EQ.3.AND.MSTP(66).GE.1) THEN
0073         Q2MNCS(2)=MAX(Q2MNC,VINT(284))
0074       ENDIF
0075       MCEV=0
0076       ALAMS=PARU(112)
0077       PARU(112)=PARP(61)
0078       FQ2C=1D0
0079       TCMX=0D0
0080       IF(MINT(47).GE.2.AND.(MINT(47).LT.5.OR.MSTP(12).GE.1)) THEN
0081         MCEV=1
0082         IF(MSTP(64).EQ.1) FQ2C=PARP(63)
0083         IF(MSTP(64).EQ.2) FQ2C=PARP(64)
0084         TCMX=LOG(FQ2C*Q2MX/PARP(61)**2)
0085         IF(Q2MX.LT.MAX(Q2MNC,2D0*PARP(61)**2).OR.TCMX.LT.0.2D0)
0086      &  MCEV=0
0087       ENDIF
0088
0089 C...Initialize QED evolution and check phase space.
0090       MEEV=0
0091       XEE=1D-10
0092       SPME=PMAS(11,1)**2
0093       IF(IABS(MINT(11)).EQ.13.OR.IABS(MINT(12)).EQ.13)
0094      &SPME=PMAS(13,1)**2
0095       IF(IABS(MINT(11)).EQ.15.OR.IABS(MINT(12)).EQ.15)
0096      &SPME=PMAS(15,1)**2
0097       Q2MNE=MAX(PARP(68)**2,2D0*SPME)
0098       TEMX=0D0
0099       FWTE=10D0
0100       IF(MINT(45).EQ.3.OR.MINT(46).EQ.3) THEN
0101         MEEV=1
0102         TEMX=LOG(Q2MX/SPME)
0103         IF(Q2MX.LE.Q2MNE.OR.TEMX.LT.0.2D0) MEEV=0
0104       ENDIF
0105       IF(MSTP(61).GE.2.AND.MCEV.EQ.1.AND.MEEV.EQ.0) THEN
0106         MEEV=2
0107         TEMX=TCMX
0108         FWTE=1D0
0109       ENDIF
0110       IF(MCEV.EQ.0.AND.MEEV.EQ.0) RETURN
0111
0112 C...Loopback point in case of failure to reconstruct kinematics.
0113       NS=N
0114       LOOP=0
0115       MNT352=MINT(352)
0116       MNT353=MINT(353)
0117       VNT352=VINT(352)
0118       VNT353=VINT(353)
0119   100 LOOP=LOOP+1
0120       IF(LOOP.GT.100) THEN
0121         MINT(51)=1
0122         RETURN
0123       ENDIF
0124       N=NS
0125       MINT(352)=MNT352
0126       MINT(353)=MNT353
0127       VINT(352)=VNT352
0128       VINT(353)=VNT353
0129
0130 C...Initial values: flavours, momenta, virtualities.
0131       DO 120 JT=1,2
0132         MORE(JT)=1
0133         KFBEAM(JT)=MINT(10+JT)
0134         IF(MINT(18+JT).EQ.1)KFBEAM(JT)=22
0135         KFLS(JT)=MINT(14+JT)
0136         KFLS(JT+2)=KFLS(JT)
0137         XS(JT)=VINT(40+JT)
0138         IF(MINT(18+JT).EQ.1) XS(JT)=VINT(40+JT)/VINT(154+JT)
0139         IF(MINT(31).GE.2) XS(JT)=XS(JT)/VINT(142+JT)
0140         ZS(JT)=1D0
0141         Q2S(JT)=FCQ2MX*Q2MX
0142         DQ2(JT)=0D0
0143         TEVCSV(JT)=TCMX
0144         ALAM(JT)=PARP(61)
0145         THE2(JT)=1D0
0146         TEVESV(JT)=TEMX
0147         MCESV(JT)=0
0148 C...Calculate initial parton distribution weights.
0149         MINT(105)=MINT(102+JT)
0150         MINT(109)=MINT(106+JT)
0151         VINT(120)=VINT(2+JT)
0152         IF(XS(JT).LT.1D0-XEE) THEN
0153           IF(MINT(31).GE.2) MINT(30)=JT
0154           IF(MSTP(57).LE.1) THEN
0155             CALL PYPDFU(KFBEAM(JT),XS(JT),Q2S(JT),XFB)
0156           ELSE
0157             CALL PYPDFL(KFBEAM(JT),XS(JT),Q2S(JT),XFB)
0158           ENDIF
0159         ENDIF
0160         DO 110 KFL=-25,25
0161           XFS(JT,KFL)=XFB(KFL)
0162   110   CONTINUE
0163 C...Special kinematics check for c/b quarks (that g -> c cbar or
0164 C...b bbar kinematically possible).
0165       KFLCB=IABS(KFLS(JT))
0166       IF(KFBEAM(JT).NE.22.AND.(KFLCB.EQ.4.OR.KFLCB.EQ.5)) THEN
0167         IF(XS(JT).GT.0.9D0*Q2S(JT)/(PMAS(KFLCB,1)**2+Q2S(JT))) THEN
0168           MINT(51)=1
0169           RETURN
0170         ENDIF
0171       ENDIF
0172   120 CONTINUE
0173       DSH=VINT(44)
0174       IF(ISET(ISUB).GE.3.AND.ISET(ISUB).LE.5) DSH=VINT(26)*VINT(2)
0175
0176 C...Find if interference with final state partons.
0177       MFIS=0
0178       IF(MSTP(67).GE.1.AND.MSTP(67).LE.3) MFIS=MSTP(67)
0179       IF(MFIS.NE.0) THEN
0180         DO 140 I=1,2
0181           KCFI(I)=0
0182           KCA=PYCOMP(IABS(KFLS(I)))
0183           IF(KCA.NE.0) KCFI(I)=KCHG(KCA,2)*ISIGN(1,KFLS(I))
0184           NFIS(I)=0
0185           IF(KCFI(I).NE.0) THEN
0186             IF(I.EQ.1) IPFS=IPUS1
0187             IF(I.EQ.2) IPFS=IPUS2
0188             DO 130 J=1,2
0189               ICSI=MOD(K(IPFS,3+J),MSTU(5))
0190               IF(ICSI.GT.0.AND.ICSI.NE.IPUS1.AND.ICSI.NE.IPUS2.AND.
0191      &        (KCFI(I).EQ.(-1)**(J+1).OR.KCFI(I).EQ.2)) THEN
0192                 NFIS(I)=NFIS(I)+1
0193                 THEFIS(I,NFIS(I))=PYANGL(P(ICSI,3),SQRT(P(ICSI,1)**2+
0194      &          P(ICSI,2)**2))
0195                 IF(I.EQ.2) THEFIS(I,NFIS(I))=PARU(1)-THEFIS(I,NFIS(I))
0196               ENDIF
0197   130       CONTINUE
0198           ENDIF
0199   140   CONTINUE
0200         IF(NFIS(1)+NFIS(2).EQ.0) MFIS=0
0201       ENDIF
0202
0203 C...Pick up leg with highest virtuality.
0204       JTOLD=1
0205   150 N=N+1
0206       JT=1
0207       IF(N.GT.NS+1.AND.Q2S(2).GT.Q2S(1)) JT=2
0208       IF(N.EQ.NS+2.AND.JT.EQ.JTOLD) JT=3-JT
0209       IF(MORE(JT).EQ.0) JT=3-JT
0210       JTOLD=JT
0211       KFLB=KFLS(JT)
0212       XB=XS(JT)
0213       DO 160 KFL=-25,25
0214         XFB(KFL)=XFS(JT,KFL)
0215   160 CONTINUE
0216       DSHR=2D0*SQRT(DSH)
0217       DSHZ=DSH/ZS(JT)
0218
0219 C...Check if allowed to branch.
0220       MCEV=0
0221       IF(IABS(KFLB).LE.10.OR.KFLB.EQ.21) THEN
0222         MCEV=1
0223         XEC=MAX(PARP(65)*DSHR/VINT2R,XB*(1D0/(1D0-PARP(66))-1D0))
0224         IF(XB.GE.1D0-2D0*XEC) MCEV=0
0225       ENDIF
0226       MEEV=0
0227       IF(MINT(44+JT).EQ.3) THEN
0228         MEEV=1
0229         IF(XB.GE.1D0-2D0*XEE) MEEV=0
0230         IF((IABS(KFLB).LE.10.OR.KFLB.EQ.21).AND.XB.GE.1D0-2D0*XEC)
0231      &  MEEV=0
0232 C***Currently kill QED shower for resolved photoproduction.
0233         IF(MINT(18+JT).EQ.1) MEEV=0
0234 C***Currently kill shower for W inside electron.
0235         IF(IABS(KFLB).EQ.24) THEN
0236           MCEV=0
0237           MEEV=0
0238         ENDIF
0239       ENDIF
0240       IF(MSTP(61).GE.2.AND.MCEV.EQ.1.AND.MEEV.EQ.0.AND.IABS(KFLB).LE.10)
0241      &MEEV=2
0242       IF(MCEV.EQ.0.AND.MEEV.EQ.0) THEN
0243         Q2B=0D0
0244         GOTO 260
0245       ENDIF
0246
0247 C...Maximum Q2 with or without Q2 ordering. Effective Lambda and n_f.
0248       Q2B=Q2S(JT)
0249       TEVCB=TEVCSV(JT)
0250       TEVEB=TEVESV(JT)
0251       IF(MSTP(62).LE.1) THEN
0252         IF(ZS(JT).GT.0.99999D0) THEN
0253           Q2B=Q2S(JT)
0254         ELSE
0255           Q2B=0.5D0*(1D0/ZS(JT)+1D0)*Q2S(JT)+0.5D0*(1D0/ZS(JT)-1D0)*
0256      &    (Q2S(3-JT)-DSH+SQRT((DSH+Q2S(1)+Q2S(2))**2+
0257      &    8D0*Q2S(1)*Q2S(2)*ZS(JT)/(1D0-ZS(JT))))
0258         ENDIF
0259         IF(MCEV.EQ.1) TEVCB=LOG(FQ2C*Q2B/ALAM(JT)**2)
0260         IF(MEEV.EQ.1) TEVEB=LOG(Q2B/SPME)
0261       ENDIF
0262       IF(MCEV.EQ.1) THEN
0263         ALSDUM=PYALPS(FQ2C*Q2B)
0264         TEVCB=TEVCB+2D0*LOG(ALAM(JT)/PARU(117))
0265         ALAM(JT)=PARU(117)
0266         B0=(33D0-2D0*MSTU(118))/6D0
0267       ENDIF
0268       IF(MEEV.EQ.2) TEVEB=TEVCB
0269       TEVCBS=TEVCB
0270       TEVEBS=TEVEB
0271
0272 C...Select side for interference with final state partons.
0273       IF(MFIS.GE.1.AND.N.LE.NS+2) THEN
0274         IFI=N-NS
0275         ISFI(IFI)=0
0276         IF(IABS(KCFI(IFI)).EQ.1.AND.NFIS(IFI).EQ.1) THEN
0277           ISFI(IFI)=1
0278         ELSEIF(KCFI(IFI).EQ.2.AND.NFIS(IFI).EQ.1) THEN
0279           IF(PYR(0).GT.0.5D0) ISFI(IFI)=1
0280         ELSEIF(KCFI(IFI).EQ.2.AND.NFIS(IFI).EQ.2) THEN
0281           ISFI(IFI)=1
0282           IF(PYR(0).GT.0.5D0) ISFI(IFI)=2
0283         ENDIF
0284       ENDIF
0285
0286 C...Calculate preweighting factor for ME-corrected processes.
0287       IF(MECOR.GE.1) CALL PYMEMX(MECOR,WTFF,WTGF,WTFG,WTGG)
0288
0289 C...Calculate Altarelli-Parisi weights.
0290       DO 170 KFL=-25,25
0291         WTAPC(KFL)=0D0
0292         WTAPE(KFL)=0D0
0293         WTSF(KFL)=0D0
0294   170 CONTINUE
0295 C...q -> q (g or gamma emission), g -> q.
0296       IF(IABS(KFLB).LE.10) THEN
0297         WTAPC(KFLB)=(8D0/3D0)*LOG((1D0-XEC-XB)*(XB+XEC)/(XEC*(1D0-XEC)))
0298         WTAPC(21)=0.5D0*(XB/(XB+XEC)-XB/(1D0-XEC))
0299         EQ2=1D0/9D0
0300         IF(MOD(IABS(KFLB),2).EQ.0) EQ2=4D0*EQ2
0301         IF(MEEV.EQ.2) WTAPE(KFLB)=2.*EQ2*LOG((1D0-XEC-XB)*(XB+XEC)/
0302      &  (XEC*(1D0-XEC)))
0303         IF(MECOR.GE.1.AND.(N.EQ.NS+1.OR.N.EQ.NS+2)) THEN
0304           WTAPC(KFLB)=WTFF*WTAPC(KFLB)
0305           WTAPC(21)=WTGF*WTAPC(21)
0306           WTAPE(KFLB)=WTFF*WTAPE(KFLB)
0307         ENDIF
0308 C...f -> f, gamma -> f.
0309       ELSEIF(IABS(KFLB).LE.20) THEN
0310         WTAPF1=LOG((1D0-XEE-XB)*(XB+XEE)/(XEE*(1D0-XEE)))
0311         WTAPF2=LOG((1D0-XEE-XB)*(1D0-XEE)/(XEE*(XB+XEE)))
0312         WTAPE(KFLB)=2D0*(WTAPF1+WTAPF2)
0313         IF(MSTP(12).GE.1) WTAPE(22)=XB/(XB+XEE)-XB/(1D0-XEE)
0314         IF(MECOR.GE.1.AND.(N.EQ.NS+1.OR.N.EQ.NS+2)) THEN
0315           WTAPE(KFLB)=WTFF*WTAPE(KFLB)
0316           WTAPE(22)=WTGF*WTAPE(22)
0317         ENDIF
0318 C...f -> g, g -> g.
0319       ELSEIF(KFLB.EQ.21) THEN
0320         WTAPQ=(16D0/3D0)*(SQRT((1D0-XEC)/XB)-SQRT((XB+XEC)/XB))
0321         DO 180 KFL=1,MSTP(58)
0322           WTAPC(KFL)=WTAPQ
0323           WTAPC(-KFL)=WTAPQ
0324   180   CONTINUE
0325         WTAPC(21)=6D0*LOG((1D0-XEC-XB)/XEC)
0326         IF(MECOR.GE.1.AND.(N.EQ.NS+1.OR.N.EQ.NS+2)) THEN
0327           DO 190 KFL=1,MSTP(58)
0328             WTAPC(KFL)=WTFG*WTAPC(KFL)
0329             WTAPC(-KFL)=WTFG*WTAPC(-KFL)
0330   190     CONTINUE
0331           WTAPC(21)=WTGG*WTAPC(21)
0332         ENDIF
0333 C...f -> gamma, W+, W-.
0334       ELSEIF(KFLB.EQ.22) THEN
0335         WTAPF=LOG((1D0-XEE-XB)*(1D0-XEE)/(XEE*(XB+XEE)))/XB
0336         WTAPE(11)=WTAPF
0337         WTAPE(-11)=WTAPF
0338         IF(MECOR.GE.1.AND.(N.EQ.NS+1.OR.N.EQ.NS+2)) THEN
0339           WTAPE(11)=WTFG*WTAPE(11)
0340           WTAPE(-11)=WTFG*WTAPE(-11)
0341         ENDIF
0342       ELSEIF(KFLB.EQ.24) THEN
0343         WTAPE(-11)=1D0/(4D0*PARU(102))*LOG((1D0-XEE-XB)*(1D0-XEE)/
0344      &  (XEE*(XB+XEE)))/XB
0345       ELSEIF(KFLB.EQ.-24) THEN
0346         WTAPE(11)=1D0/(4D0*PARU(102))*LOG((1D0-XEE-XB)*(1D0-XEE)/
0347      &  (XEE*(XB+XEE)))/XB
0348       ENDIF
0349
0350 C...Calculate parton distribution weights and sum.
0351       NTRY=0
0352   200 NTRY=NTRY+1
0353       IF(NTRY.GT.500) THEN
0354         MINT(51)=1
0355         RETURN
0356       ENDIF
0357       WTSUMC=0D0
0358       WTSUME=0D0
0359       XFBO=MAX(1D-10,XFB(KFLB))
0360       DO 210 KFL=-25,25
0361         WTSF(KFL)=XFB(KFL)/XFBO
0362         WTSUMC=WTSUMC+WTAPC(KFL)*WTSF(KFL)
0363         WTSUME=WTSUME+WTAPE(KFL)*WTSF(KFL)
0364   210 CONTINUE
0365       WTSUMC=MAX(0.0001D0,WTSUMC)
0366       WTSUME=MAX(0.0001D0/FWTE,WTSUME)
0367
0368 C...Choose new t: fix alpha_s, alpha_s(Q^2), alpha_s(k_T^2).
0369       NTRY2=0
0370   220 NTRY2=NTRY2+1
0371       IF(NTRY2.GT.500) THEN
0372         MINT(51)=1
0373         RETURN
0374       ENDIF
0375       IF(MCEV.EQ.1) THEN
0376         IF(MSTP(64).LE.0) THEN
0377           TEVCB=TEVCB+LOG(PYR(0))*PARU(2)/(PARU(111)*WTSUMC)
0378         ELSEIF(MSTP(64).EQ.1) THEN
0379           TEVCB=TEVCB*EXP(MAX(-50D0,LOG(PYR(0))*B0/WTSUMC))
0380         ELSE
0381           TEVCB=TEVCB*EXP(MAX(-50D0,LOG(PYR(0))*B0/(5D0*WTSUMC)))
0382         ENDIF
0383       ENDIF
0384       IF(MEEV.EQ.1) THEN
0385         TEVEB=TEVEB*EXP(MAX(-50D0,LOG(PYR(0))*PARU(2)/
0386      &  (PARU(101)*FWTE*WTSUME*TEMX)))
0387       ELSEIF(MEEV.EQ.2) THEN
0388         TEVEB=TEVEB+LOG(PYR(0))*PARU(2)/(PARU(101)*WTSUME)
0389       ENDIF
0390
0391 C...Translate t into Q2 scale; choose between QCD and QED evolution.
0392   230 IF(MCEV.EQ.1) Q2CB=ALAM(JT)**2*EXP(MAX(-50D0,TEVCB))/FQ2C
0393       IF(MEEV.EQ.1) Q2EB=SPME*EXP(MAX(-50D0,TEVEB))
0394       IF(MEEV.EQ.2) Q2EB=ALAM(JT)**2*EXP(MAX(-50D0,TEVEB))/FQ2C
0395 C...Ensure that Q2 is above threshold for charm/bottom.
0396       KFLCB=IABS(KFLB)
0397       IF(KFBEAM(JT).NE.22.AND.(KFLCB.EQ.4.OR.KFLCB.EQ.5).AND.
0398      &MCEV.EQ.1) THEN
0399         IF(Q2CB.LT.PMAS(KFLCB,1)**2) THEN
0400           Q2CB=1.1D0*PMAS(KFLCB,1)**2
0401           TEVCB=LOG(FQ2C*Q2B/ALAM(JT)**2)
0402           FCQ2MX=MIN(2D0,1.05D0*FCQ2MX)
0403         ENDIF
0404       ENDIF
0405       IF(KFBEAM(JT).NE.22.AND.(KFLCB.EQ.4.OR.KFLCB.EQ.5).AND.
0406      &MEEV.EQ.2) THEN
0407         IF(Q2EB.LT.PMAS(KFLCB,1)**2) MEEV=0
0408       ENDIF
0409       MCE=0
0410       IF(MCEV.EQ.0.AND.MEEV.EQ.0) THEN
0411       ELSEIF(MCEV.EQ.1.AND.MEEV.EQ.0) THEN
0412         IF(Q2CB.GT.Q2MNCS(JT)) MCE=1
0413       ELSEIF(MCEV.EQ.0.AND.MEEV.EQ.1) THEN
0414         IF(Q2EB.GT.Q2MNE) MCE=2
0415       ELSEIF(MCEV.EQ.0.AND.MEEV.EQ.2) THEN
0416         IF(Q2EB.GT.Q2MNCS(JT)) MCE=2
0417       ELSEIF(MCEV.EQ.1.AND.MEEV.EQ.2) THEN
0418         IF(Q2CB.GT.Q2EB.AND.Q2CB.GT.Q2MNCS(JT)) MCE=1
0419         IF(Q2EB.GT.Q2CB.AND.Q2EB.GT.Q2MNCS(JT)) MCE=2
0420       ELSEIF(Q2MNCS(JT).GT.Q2MNE) THEN
0421         MCE=1
0422         IF(Q2EB.GT.Q2CB.OR.Q2CB.LE.Q2MNCS(JT)) MCE=2
0423         IF(MCE.EQ.2.AND.Q2EB.LE.Q2MNE) MCE=0
0424       ELSE
0425         MCE=2
0426         IF(Q2CB.GT.Q2EB.OR.Q2EB.LE.Q2MNE) MCE=1
0427         IF(MCE.EQ.1.AND.Q2CB.LE.Q2MNCS(JT)) MCE=0
0428       ENDIF
0429
0430 C...Evolution possibly ended. Update t values.
0431       IF(MCE.EQ.0) THEN
0432         Q2B=0D0
0433         GOTO 260
0434       ELSEIF(MCE.EQ.1) THEN
0435         Q2B=Q2CB
0436         Q2REF=FQ2C*Q2B
0437         IF(MEEV.EQ.1) TEVEB=LOG(Q2B/SPME)
0438         IF(MEEV.EQ.2) TEVEB=LOG(FQ2C*Q2B/ALAM(JT)**2)
0439       ELSE
0440         Q2B=Q2EB
0441         Q2REF=Q2B
0442         IF(MCEV.EQ.1) TEVCB=LOG(FQ2C*Q2B/ALAM(JT)**2)
0443       ENDIF
0444
0445 C...Select flavour for branching parton.
0446       IF(MCE.EQ.1) WTRAN=PYR(0)*WTSUMC
0447       IF(MCE.EQ.2) WTRAN=PYR(0)*WTSUME
0448       KFLA=-25
0449   240 KFLA=KFLA+1
0450       IF(MCE.EQ.1) WTRAN=WTRAN-WTAPC(KFLA)*WTSF(KFLA)
0451       IF(MCE.EQ.2) WTRAN=WTRAN-WTAPE(KFLA)*WTSF(KFLA)
0452       IF(KFLA.LE.24.AND.WTRAN.GT.0D0) GOTO 240
0453       IF(KFLA.EQ.25) THEN
0454         Q2B=0D0
0455         GOTO 260
0456       ENDIF
0457
0458 C...Choose z value and corrective weight.
0459       WTZ=0D0
0460 C...q -> q + g or q -> q + gamma.
0461       IF(IABS(KFLA).LE.10.AND.IABS(KFLB).LE.10) THEN
0462         Z=1D0-((1D0-XB-XEC)/(1D0-XEC))*
0463      &  (XEC*(1D0-XEC)/((XB+XEC)*(1D0-XB-XEC)))**PYR(0)
0464         WTZ=0.5D0*(1D0+Z**2)
0465 C...q -> g + q.
0466       ELSEIF(IABS(KFLA).LE.10.AND.KFLB.EQ.21) THEN
0467         Z=XB/(SQRT(XB+XEC)+PYR(0)*(SQRT(1D0-XEC)-SQRT(XB+XEC)))**2
0468         WTZ=0.5D0*(1D0+(1D0-Z)**2)*SQRT(Z)
0469 C...f -> f + gamma.
0470       ELSEIF(IABS(KFLA).LE.20.AND.IABS(KFLB).LE.20) THEN
0471         IF(WTAPF1.GT.PYR(0)*(WTAPF1+WTAPF2)) THEN
0472           Z=1D0-((1D0-XB-XEE)/(1D0-XEE))*
0473      &    (XEE*(1D0-XEE)/((XB+XEE)*(1D0-XB-XEE)))**PYR(0)
0474         ELSE
0475           Z=XB+XB*(XEE/(1D0-XEE))*
0476      &    ((1D0-XB-XEE)*(1D0-XEE)/(XEE*(XB+XEE)))**PYR(0)
0477         ENDIF
0478         WTZ=0.5D0*(1D0+Z**2)*(Z-XB)/(1D0-XB)
0479 C...f -> gamma + f.
0480       ELSEIF(IABS(KFLA).LE.20.AND.KFLB.EQ.22) THEN
0481         Z=XB+XB*(XEE/(1D0-XEE))*
0482      &  ((1D0-XB-XEE)*(1D0-XEE)/(XEE*(XB+XEE)))**PYR(0)
0483         WTZ=0.5D0*(1D0+(1D0-Z)**2)*XB*(Z-XB)/Z
0484 C...f -> W+- + f.
0485       ELSEIF(IABS(KFLA).LE.20.AND.IABS(KFLB).EQ.24) THEN
0486         Z=XB+XB*(XEE/(1D0-XEE))*
0487      &  ((1D0-XB-XEE)*(1D0-XEE)/(XEE*(XB+XEE)))**PYR(0)
0488         WTZ=0.5D0*(1D0+(1D0-Z)**2)*(XB*(Z-XB)/Z)*
0489      &  (Q2B/(Q2B+PMAS(24,1)**2))
0490 C...g -> q + qbar.
0491       ELSEIF(KFLA.EQ.21.AND.IABS(KFLB).LE.10) THEN
0492         Z=XB/(1D0-XEC)+PYR(0)*(XB/(XB+XEC)-XB/(1D0-XEC))
0493         WTZ=1D0-2D0*Z*(1D0-Z)
0494 C...g -> g + g.
0495       ELSEIF(KFLA.EQ.21.AND.KFLB.EQ.21) THEN
0496         Z=1D0/(1D0+((1D0-XEC-XB)/XB)*(XEC/(1D0-XEC-XB))**PYR(0))
0497         WTZ=(1D0-Z*(1D0-Z))**2
0498 C...gamma -> f + fbar.
0499       ELSEIF(KFLA.EQ.22.AND.IABS(KFLB).LE.20) THEN
0500         Z=XB/(1D0-XEE)+PYR(0)*(XB/(XB+XEE)-XB/(1D0-XEE))
0501         WTZ=1D0-2D0*Z*(1D0-Z)
0502       ENDIF
0503       IF(MCE.EQ.2.AND.MEEV.EQ.1) WTZ=(WTZ/FWTE)*(TEVEB/TEMX)
0504
0505 C...Option with resummation of soft gluon emission as effective z shift.
0506       IF(MCE.EQ.1) THEN
0507         IF(MSTP(65).GE.1) THEN
0508           RSOFT=6D0
0509           IF(KFLB.NE.21) RSOFT=8D0/3D0
0510           Z=Z*(TEVCB/TEVCSV(JT))**(RSOFT*XEC/((XB+XEC)*B0))
0511           IF(Z.LE.XB) GOTO 220
0512         ENDIF
0513
0514 C...Option with alpha_s(k_T^2): demand k_T^2 > cutoff, reweight.
0515         IF(MSTP(64).GE.2) THEN
0516           IF((1D0-Z)*Q2B.LT.Q2MNCS(JT)) GOTO 220
0517           ALPRAT=TEVCB/(TEVCB+LOG(1D0-Z))
0518           IF(ALPRAT.LT.5D0*PYR(0)) GOTO 220
0519           IF(ALPRAT.GT.5D0) WTZ=WTZ*ALPRAT/5D0
0520         ENDIF
0521       ENDIF
0522
0523 C...Remove kinematically impossible branchings.
0524       UHAT=Q2B-DSH*(1D0-Z)/Z
0525       IF(MSTP(68).GE.0.AND.UHAT.GT.0D0) GOTO 220
0526
0527 C...Select phi angle of branching at random.
0528       PHIBR=PARU(2)*PYR(0)
0529
0530 C...Matrix-element corrections for some processes.
0531       IF(MECOR.GE.1.AND.(N.EQ.NS+1.OR.N.EQ.NS+2)) THEN
0532         IF(IABS(KFLA).LE.20.AND.IABS(KFLB).LE.20) THEN
0533           CALL PYMEWT(MECOR,1,Q2B,Z,PHIBR,WTME)
0534           WTZ=WTZ*WTME/WTFF
0535         ELSEIF((KFLA.EQ.21.OR.KFLA.EQ.22).AND.IABS(KFLB).LE.20) THEN
0536           CALL PYMEWT(MECOR,2,Q2B,Z,PHIBR,WTME)
0537           WTZ=WTZ*WTME/WTGF
0538         ELSEIF(IABS(KFLA).LE.20.AND.(KFLB.EQ.21.OR.KFLB.EQ.22)) THEN
0539           CALL PYMEWT(MECOR,3,Q2B,Z,PHIBR,WTME)
0540           WTZ=WTZ*WTME/WTFG
0541         ELSEIF(KFLA.EQ.21.AND.KFLB.EQ.21) THEN
0542           CALL PYMEWT(MECOR,4,Q2B,Z,PHIBR,WTME)
0543           WTZ=WTZ*WTME/WTGG
0544         ENDIF
0545       ENDIF
0546
0547 C...Impose angular constraint in first branching from interference
0548 C...with final state partons.
0549       IF(MCE.EQ.1) THEN
0550         IF(MFIS.GE.1.AND.N.LE.NS+2.AND.NTRY2.LT.200) THEN
0551           THE2D=(4D0*Q2B)/(DSH*(1D0-Z))
0552           IF(N.EQ.NS+1.AND.ISFI(1).GE.1) THEN
0553             IF(THE2D.GT.THEFIS(1,ISFI(1))**2) GOTO 220
0554           ELSEIF(N.EQ.NS+2.AND.ISFI(2).GE.1) THEN
0555             IF(THE2D.GT.THEFIS(2,ISFI(2))**2) GOTO 220
0556           ENDIF
0557         ENDIF
0558
0559 C...Option with angular ordering requirement.
0560         IF(MSTP(62).GE.3.AND.NTRY2.LT.200) THEN
0561           THE2T=(4D0*Z**2*Q2B)/(4D0*Z**2*Q2B+(1D0-Z)*XB**2*VINT2R)
0562           IF(THE2T.GT.THE2(JT)) GOTO 220
0563         ENDIF
0564       ENDIF
0565
0566 C...Weighting with new parton distributions.
0567       MINT(105)=MINT(102+JT)
0568       MINT(109)=MINT(106+JT)
0569       VINT(120)=VINT(2+JT)
0570       IF(MINT(31).GE.2) MINT(30)=JT
0571       IF(MSTP(57).LE.1) THEN
0572         CALL PYPDFU(KFBEAM(JT),XB,Q2REF,XFN)
0573       ELSE
0574         CALL PYPDFL(KFBEAM(JT),XB,Q2REF,XFN)
0575       ENDIF
0576       XFBN=XFN(KFLB)
0577       IF(XFBN.LT.1D-20) THEN
0578         IF(KFLA.EQ.KFLB) THEN
0579           TEVCB=TEVCBS
0580           TEVEB=TEVEBS
0581           WTAPC(KFLB)=0D0
0582           WTAPE(KFLB)=0D0
0583           GOTO 200
0584         ELSEIF(MCE.EQ.1.AND.TEVCBS-TEVCB.GT.0.2D0) THEN
0585           TEVCB=0.5D0*(TEVCBS+TEVCB)
0586           GOTO 230
0587         ELSEIF(MCE.EQ.2.AND.TEVEBS-TEVEB.GT.0.2D0) THEN
0588           TEVEB=0.5D0*(TEVEBS+TEVEB)
0589           GOTO 230
0590         ELSE
0591           XFBN=1D-10
0592           XFN(KFLB)=XFBN
0593         ENDIF
0594       ENDIF
0595       DO 250 KFL=-25,25
0596         XFB(KFL)=XFN(KFL)
0597   250 CONTINUE
0598       XA=XB/Z
0599       IF(MINT(31).GE.2) MINT(30)=JT
0600       IF(MSTP(57).LE.1) THEN
0601         CALL PYPDFU(KFBEAM(JT),XA,Q2REF,XFA)
0602       ELSE
0603         CALL PYPDFL(KFBEAM(JT),XA,Q2REF,XFA)
0604       ENDIF
0605       XFAN=XFA(KFLA)
0606       IF(XFAN.LT.1D-20) GOTO 200
0607       WTSFA=WTSF(KFLA)
0608       IF(WTZ*XFAN/XFBN.LT.PYR(0)*WTSFA) GOTO 200
0609
0610 C...Define two hard scatterers in their CM-frame.
0611   260 IF(N.EQ.NS+2) THEN
0612         DQ2(JT)=Q2B
0613         DPLCM=SQRT((DSH+DQ2(1)+DQ2(2))**2-4D0*DQ2(1)*DQ2(2))/DSHR
0614         DO 280 JR=1,2
0615           I=NS+JR
0616           IF(JR.EQ.1) IPO=IPUS1
0617           IF(JR.EQ.2) IPO=IPUS2
0618           DO 270 J=1,5
0619             K(I,J)=0
0620             P(I,J)=0D0
0621             V(I,J)=0D0
0622   270     CONTINUE
0623           K(I,1)=14
0624           K(I,2)=KFLS(JR+2)
0625           K(I,4)=IPO
0626           K(I,5)=IPO
0627           P(I,3)=DPLCM*(-1)**(JR+1)
0628           P(I,4)=(DSH+DQ2(3-JR)-DQ2(JR))/DSHR
0629           P(I,5)=-SQRT(DQ2(JR))
0630           K(IPO,1)=14
0631           K(IPO,3)=I
0632           K(IPO,4)=MOD(K(IPO,4),MSTU(5))+MSTU(5)*I
0633           K(IPO,5)=MOD(K(IPO,5),MSTU(5))+MSTU(5)*I
0634   280   CONTINUE
0635
0636 C...Find maximum allowed mass of timelike parton.
0637       ELSEIF(N.GT.NS+2) THEN
0638         JR=3-JT
0639         DQ2(3)=Q2B
0640         DPC(1)=P(IS(1),4)
0641         DPC(2)=P(IS(2),4)
0642         DPC(3)=0.5D0*(ABS(P(IS(1),3))+ABS(P(IS(2),3)))
0643         DPD(1)=DSH+DQ2(JR)+DQ2(JT)
0644         DPD(2)=DSHZ+DQ2(JR)+DQ2(3)
0645         DPD(3)=SQRT(DPD(1)**2-4D0*DQ2(JR)*DQ2(JT))
0646         DPD(4)=SQRT(DPD(2)**2-4D0*DQ2(JR)*DQ2(3))
0647         IKIN=0
0648         IF(Q2S(JR).GE.0.25D0*Q2MNC.AND.DPD(1)-DPD(3).GE.
0649      &  1D-10*DPD(1)) IKIN=1
0650         IF(IKIN.EQ.0) DMSMA=(DQ2(JT)/ZS(JT)-DQ2(3))*
0651      &  (DSH/(DSH+DQ2(JT))-DSH/(DSHZ+DQ2(3)))
0652         IF(IKIN.EQ.1) DMSMA=(DPD(1)*DPD(2)-DPD(3)*DPD(4))/
0653      &  (2D0*DQ2(JR))-DQ2(JT)-DQ2(3)
0654
0655 C...Generate timelike parton shower (if required).
0656         IT=N
0657         DO 290 J=1,5
0658           K(IT,J)=0
0659           P(IT,J)=0D0
0660           V(IT,J)=0D0
0661   290   CONTINUE
0662 C...f -> f + g (gamma).
0663         IF(IABS(KFLB).LE.20.AND.IABS(KFLS(JT+2)).LE.20) THEN
0664           K(IT,2)=21
0665           IF(MCESV(JT).EQ.2.OR.IABS(KFLB).GE.11) K(IT,2)=22
0666 C...f -> g (gamma, W+-) + f.
0667         ELSEIF(IABS(KFLB).LE.20.AND.IABS(KFLS(JT+2)).GT.20) THEN
0668           K(IT,2)=KFLB
0669           IF(KFLS(JT+2).EQ.24) THEN
0670             K(IT,2)=-12
0671           ELSEIF(KFLS(JT+2).EQ.-24) THEN
0672             K(IT,2)=12
0673           ENDIF
0674 C...g (gamma) -> f + fbar, g + g.
0675         ELSE
0676           K(IT,2)=-KFLS(JT+2)
0677           IF(KFLS(JT+2).GT.20) K(IT,2)=KFLS(JT+2)
0678         ENDIF
0679         K(IT,1)=3
0680         IF((IABS(K(IT,2)).GE.11.AND.IABS(K(IT,2)).LE.18).OR.
0681      &  IABS(K(IT,2)).EQ.22) K(IT,1)=1
0682         P(IT,5)=PYMASS(K(IT,2))
0683         IF(DMSMA.LE.P(IT,5)**2) GOTO 100
0684         IF(MSTP(63).GE.1.AND.MCESV(JT).EQ.1) THEN
0685           MSTJ48=MSTJ(48)
0686           PARJ85=PARJ(85)
0687           P(IT,4)=(DSHZ-DSH-P(IT,5)**2)/DSHR
0688           P(IT,3)=SQRT(P(IT,4)**2-P(IT,5)**2)
0689           IF(MSTP(63).EQ.1) THEN
0690             Q2TIM=DMSMA
0691           ELSEIF(MSTP(63).EQ.2) THEN
0692             Q2TIM=MIN(DMSMA,PARP(71)*Q2S(JT))
0693           ELSE
0694             Q2TIM=DMSMA
0695             MSTJ(48)=1
0696             IF(IKIN.EQ.0) DPT2=DMSMA*(DSHZ+DQ2(3))/(DSH+DQ2(JT))
0697             IF(IKIN.EQ.1) DPT2=DMSMA*(0.5D0*DPD(1)*DPD(2)+0.5D0*DPD(3)*
0698      &      DPD(4)-DQ2(JR)*(DQ2(JT)+DQ2(3)))/(4D0*DSH*DPC(3)**2)
0699             PARJ(85)=SQRT(MAX(0D0,DPT2))*
0700      &      (1D0/P(IT,4)+1D0/P(IS(JT),4))
0701           ENDIF
0702           CALL PYSHOW(IT,0,SQRT(Q2TIM))
0703           MSTJ(48)=MSTJ48
0704           PARJ(85)=PARJ85
0705           IF(N.GE.IT+1) P(IT,5)=P(IT+1,5)
0706         ENDIF
0707
0708 C...Reconstruct kinematics of branching: timelike parton shower.
0709         DMS=P(IT,5)**2
0710         IF(IKIN.EQ.0) DPT2=(DMSMA-DMS)*(DSHZ+DQ2(3))/(DSH+DQ2(JT))
0711         IF(IKIN.EQ.1) DPT2=(DMSMA-DMS)*(0.5D0*DPD(1)*DPD(2)+
0712      &  0.5D0*DPD(3)*DPD(4)-DQ2(JR)*(DQ2(JT)+DQ2(3)+DMS))/
0713      &  (4D0*DSH*DPC(3)**2)
0714         IF(DPT2.LT.0D0) GOTO 100
0715         DPB(1)=(0.5D0*DPD(2)-DPC(JR)*(DSHZ+DQ2(JR)-DQ2(JT)-DMS)/
0716      &  DSHR)/DPC(3)-DPC(3)
0717         P(IT,1)=SQRT(DPT2)
0718         P(IT,3)=DPB(1)*(-1)**(JT+1)
0719         P(IT,4)=SQRT(DPT2+DPB(1)**2+DMS)
0720         IF(N.GE.IT+1) THEN
0721           DPB(1)=SQRT(DPB(1)**2+DPT2)
0722           DPB(2)=SQRT(DPB(1)**2+DMS)
0723           DPB(3)=P(IT+1,3)
0724           DPB(4)=SQRT(DPB(3)**2+DMS)
0725           DBEZ=(DPB(4)*DPB(1)-DPB(3)*DPB(2))/(DPB(4)*DPB(2)-DPB(3)*
0726      &    DPB(1))
0727           CALL PYROBO(IT+1,N,0D0,0D0,0D0,0D0,DBEZ)
0728           THE=PYANGL(P(IT,3),P(IT,1))
0729           CALL PYROBO(IT+1,N,THE,0D0,0D0,0D0,0D0)
0730         ENDIF
0731
0732 C...Reconstruct kinematics of branching: spacelike parton.
0733         DO 300 J=1,5
0734           K(N+1,J)=0
0735           P(N+1,J)=0D0
0736           V(N+1,J)=0D0
0737   300   CONTINUE
0738         K(N+1,1)=14
0739         K(N+1,2)=KFLB
0740         P(N+1,1)=P(IT,1)
0741         P(N+1,3)=P(IT,3)+P(IS(JT),3)
0742         P(N+1,4)=P(IT,4)+P(IS(JT),4)
0743         P(N+1,5)=-SQRT(DQ2(3))
0744
0745 C...Define colour flow of branching.
0746         K(IS(JT),3)=N+1
0747         K(IT,3)=N+1
0748         IM1=N+1
0749         IM2=N+1
0750 C...f -> f + gamma (Z, W).
0751         IF(IABS(K(IT,2)).GE.22) THEN
0752           K(IT,1)=1
0753           ID1=IS(JT)
0754           ID2=IS(JT)
0755 C...f -> gamma (Z, W) + f.
0756         ELSEIF(IABS(K(IS(JT),2)).GE.22) THEN
0757           ID1=IT
0758           ID2=IT
0759 C...gamma -> q + qbar, g + g.
0760         ELSEIF(K(N+1,2).EQ.22) THEN
0761           ID1=IS(JT)
0762           ID2=IT
0763           IM1=ID2
0764           IM2=ID1
0765 C...q -> q + g.
0766         ELSEIF(K(N+1,2).GT.0.AND.K(N+1,2).NE.21.AND.K(IT,2).EQ.21) THEN
0767           ID1=IT
0768           ID2=IS(JT)
0769 C...q -> g + q.
0770         ELSEIF(K(N+1,2).GT.0.AND.K(N+1,2).NE.21) THEN
0771           ID1=IS(JT)
0772           ID2=IT
0773 C...qbar -> qbar + g.
0774         ELSEIF(K(N+1,2).LT.0.AND.K(IT,2).EQ.21) THEN
0775           ID1=IS(JT)
0776           ID2=IT
0777 C...qbar -> g + qbar.
0778         ELSEIF(K(N+1,2).LT.0) THEN
0779           ID1=IT
0780           ID2=IS(JT)
0781 C...g -> g + g; g -> q + qbar.
0782         ELSEIF((K(IT,2).EQ.21.AND.PYR(0).GT.0.5D0).OR.K(IT,2).LT.0) THEN
0783           ID1=IS(JT)
0784           ID2=IT
0785         ELSE
0786           ID1=IT
0787           ID2=IS(JT)
0788         ENDIF
0789         IF(IM1.EQ.N+1) K(IM1,4)=K(IM1,4)+ID1
0790         IF(IM2.EQ.N+1) K(IM2,5)=K(IM2,5)+ID2
0791         K(ID1,4)=K(ID1,4)+MSTU(5)*IM1
0792         K(ID2,5)=K(ID2,5)+MSTU(5)*IM2
0793         IF(ID1.NE.ID2) THEN
0794           K(ID1,5)=K(ID1,5)+MSTU(5)*ID2
0795           K(ID2,4)=K(ID2,4)+MSTU(5)*ID1
0796         ENDIF
0797         N=N+1
0798         IF(K(IT,1).EQ.1) THEN
0799           K(IT,4)=0
0800           K(IT,5)=0
0801         ENDIF
0802
0803 C...Boost to new CM-frame.
0804         DBSVX=(P(N,1)+P(IS(JR),1))/(P(N,4)+P(IS(JR),4))
0805         DBSVZ=(P(N,3)+P(IS(JR),3))/(P(N,4)+P(IS(JR),4))
0806         IF(DBSVX**2+DBSVZ**2.GE.1D0) GOTO 100
0807         CALL PYROBO(NS+1,N,0D0,0D0,-DBSVX,0D0,-DBSVZ)
0808         IR=N+(JT-1)*(IS(1)-N)
0809         CALL PYROBO(NS+1,N,-PYANGL(P(IR,3),P(IR,1)),DPHI(JT),
0810      &  0D0,0D0,0D0)
0811
0812 C...Global statistics.
0813         MINT(352)=MINT(352)+1
0814         VINT(352)=VINT(352)+SQRT(P(IT,1)**2+P(IT,2)**2)
0815         IF (MINT(352).EQ.1) VINT(357)=SQRT(P(IT,1)**2+P(IT,2)**2)
0816       ENDIF
0817
0818 C...Update kinematics variables.
0819       IS(JT)=N
0820       DQ2(JT)=Q2B
0821       IF(MSTP(62).GE.3.AND.NTRY2.LT.200) THE2(JT)=THE2T
0822       DSH=DSHZ
0823
0824 C...Save quantities; loop back.
0825       Q2S(JT)=Q2B
0826       DPHI(JT)=PHIBR
0827       MCESV(JT)=MCE
0828       IF((MCEV.EQ.1.AND.Q2B.GE.0.25D0*Q2MNC).OR.
0829      &(MEEV.EQ.1.AND.Q2B.GE.Q2MNE)) THEN
0830         KFLS(JT+2)=KFLS(JT)
0831         KFLS(JT)=KFLA
0832         XS(JT)=XA
0833         ZS(JT)=Z
0834         DO 310 KFL=-25,25
0835           XFS(JT,KFL)=XFA(KFL)
0836   310   CONTINUE
0837         TEVCSV(JT)=TEVCB
0838         TEVESV(JT)=TEVEB
0839       ELSE
0840         MORE(JT)=0
0841         IF(JT.EQ.1) IPU1=N
0842         IF(JT.EQ.2) IPU2=N
0843       ENDIF
0844       IF(N.GT.MSTU(4)-MSTU(32)-10) THEN
0845         CALL PYERRM(11,'(PYSSPA:) no more memory left in PYJETS')
0846         IF(MSTU(21).GE.1) N=NS
0847         IF(MSTU(21).GE.1) RETURN
0848       ENDIF
0849       IF(MORE(1).EQ.1.OR.MORE(2).EQ.1) GOTO 150
0850
0851 C...Boost hard scattering partons to frame of shower initiators.
0852       DO 320 J=1,3
0853         ROBO(J+2)=(P(NS+1,J)+P(NS+2,J))/(P(NS+1,4)+P(NS+2,4))
0854   320 CONTINUE
0855       K(N+2,1)=1
0856       DO 330 J=1,5
0857         P(N+2,J)=P(NS+1,J)
0858   330 CONTINUE
0859       CALL PYROBO(N+2,N+2,0D0,0D0,-ROBO(3),-ROBO(4),-ROBO(5))
0860       ROBO(2)=PYANGL(P(N+2,1),P(N+2,2))
0861       ROBO(1)=PYANGL(P(N+2,3),SQRT(P(N+2,1)**2+P(N+2,2)**2))
0862       IMIN=MINT(83)+5
0863       IF(MINT(31).GE.2) IMIN=MIN(IPUS1,IPUS2)
0864       CALL PYROBO(IMIN,NS,0D0,-ROBO(2),0D0,0D0,0D0)
0865       CALL PYROBO(IMIN,NS,ROBO(1),ROBO(2),ROBO(3),ROBO(4),ROBO(5))
0866
0867 C...Store user information. Reset Lambda value.
0868       IF(MINT(31).LE.1) THEN
0869         K(IPU1,3)=MINT(83)+3
0870         K(IPU2,3)=MINT(83)+4
0871       ELSE
0872         K(IPU1,3)=MINT(83)+1
0873         K(IPU2,3)=MINT(83)+2
0874       ENDIF
0875       DO 340 JT=1,2
0876         MINT(12+JT)=KFLS(JT)
0877         VINT(140+JT)=XS(JT)
0878         IF(MINT(18+JT).EQ.1) VINT(140+JT)=VINT(154+JT)*XS(JT)
0879         IF(MINT(31).GE.2) VINT(140+JT)=VINT(140+JT)*VINT(142+JT)
0880   340 CONTINUE
0881       PARU(112)=ALAMS
0882
0883       RETURN
0884       END