PYTHIA/pythia/pysspa.F

   1
   2 C*********************************************************************
   3
   4       SUBROUTINE PYSSPA(IPU1,IPU2)
   5
   6 C...Generates spacelike parton showers.
   7       IMPLICIT DOUBLE PRECISION(D)
   8       COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5)
   9       COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
  10       COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
  11       COMMON/PYSUBS/MSEL,MSUB(200),KFIN(2,-40:40),CKIN(200)
  12       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
  13       COMMON/PYINT1/MINT(400),VINT(400)
  14       COMMON/PYINT2/ISET(200),KFPR(200,2),COEF(200,20),ICOL(40,4,2)
  15       COMMON/PYINT3/XSFX(2,-40:40),ISIG(1000,3),SIGH(1000)
  16       SAVE /LUJETS/,/LUDAT1/,/LUDAT2/
  17       SAVE /PYSUBS/,/PYPARS/,/PYINT1/,/PYINT2/,/PYINT3/
  18       DIMENSION KFLS(4),IS(2),XS(2),ZS(2),Q2S(2),TEVCSV(2),TEVESV(2),
  19      &XFS(2,-25:25),XFA(-25:25),XFB(-25:25),XFN(-25:25),WTAPC(-25:25),
  20      &WTAPE(-25:25),WTSF(-25:25),THE2(2),ALAM(2),DQ2(3),DPC(3),DPD(4),
  21      &DPB(4),ROBO(5),MORE(2),KFBEAM(2),Q2MNCS(2),KCFI(2),NFIS(2),
  22      &THEFIS(2,2),ISFI(2)
  23       DATA IS/2*0/
  24
  25 C...Read out basic information; set global Q^2 scale.
  26       IPUS1=IPU1
  27       IPUS2=IPU2
  28       ISUB=MINT(1)
  29       Q2MX=VINT(56)
  30       IF(ISET(ISUB).EQ.2) Q2MX=PARP(67)*VINT(56)
  31
  32 C...Initialize QCD evolution and check phase space.
  33       Q2MNC=PARP(62)**2
  34       Q2MNCS(1)=Q2MNC
  35       IF(MSTP(66).EQ.1.AND.MINT(107).EQ.3)
  36      &Q2MNCS(1)=MAX(Q2MNC,VINT(283))
  37       Q2MNCS(2)=Q2MNC
  38       IF(MSTP(66).EQ.1.AND.MINT(108).EQ.3)
  39      &Q2MNCS(2)=MAX(Q2MNC,VINT(284))
  40       MCEV=0
  41       XEC0=2.*PARP(65)/VINT(1)
  42       ALAMS=PARU(112)
  43       PARU(112)=PARP(61)
  44       FQ2C=1.
  45       TCMX=0.
  46       IF(MINT(47).GE.2.AND.(MINT(47).NE.5.OR.MSTP(12).GE.1)) THEN
  47         MCEV=1
  48         IF(MSTP(64).EQ.1) FQ2C=PARP(63)
  49         IF(MSTP(64).EQ.2) FQ2C=PARP(64)
  50         TCMX=LOG(FQ2C*Q2MX/PARP(61)**2)
  51         IF(Q2MX.LT.MAX(Q2MNC,2.*PARP(61)**2).OR.TCMX.LT.0.2)
  52      &  MCEV=0
  53       ENDIF
  54
  55 C...Initialize QED evolution and check phase space.
  56       Q2MNE=PARP(68)**2
  57       MEEV=0
  58       XEE=1E-6
  59       SPME=PMAS(11,1)**2
  60       TEMX=0.
  61       FWTE=10.
  62       IF(MINT(45).EQ.3.OR.MINT(46).EQ.3) THEN
  63         MEEV=1
  64         TEMX=LOG(Q2MX/SPME)
  65         IF(Q2MX.LE.Q2MNE.OR.TEMX.LT.0.2) MEEV=0
  66       ENDIF
  67       IF(MCEV.EQ.0.AND.MEEV.EQ.0) RETURN
  68
  69 C...Initial values: flavours, momenta, virtualities.
  70       NS=N
  71   100 N=NS
  72       DO 120 JT=1,2
  73       MORE(JT)=1
  74       KFBEAM(JT)=MINT(10+JT)
  75       IF(MINT(18+JT).EQ.1)KFBEAM(JT)=22
  76       KFLS(JT)=MINT(14+JT)
  77       KFLS(JT+2)=KFLS(JT)
  78       XS(JT)=VINT(40+JT)
  79       IF(MINT(18+JT).EQ.1) XS(JT)=VINT(40+JT)/VINT(154+JT)
  80       ZS(JT)=1.
  81       Q2S(JT)=Q2MX
  82       TEVCSV(JT)=TCMX
  83       ALAM(JT)=PARP(61)
  84       THE2(JT)=100.
  85       TEVESV(JT)=TEMX
  86       DO 110 KFL=-25,25
  87       XFS(JT,KFL)=XSFX(JT,KFL)
  88   110 CONTINUE
  89   120 CONTINUE
  90       DSH=VINT(44)
  91       IF(ISET(ISUB).GE.3.AND.ISET(ISUB).LE.5) DSH=VINT(26)*VINT(2)
  92
  93 C...Find if interference with final state partons.
  94       MFIS=0
  95       IF(MSTP(67).GE.1.AND.MSTP(67).LE.3) MFIS=MSTP(67)
  96       IF(MFIS.NE.0) THEN
  97         DO 140 I=1,2
  98         KCFI(I)=0
  99         KCA=LUCOMP(IABS(KFLS(I)))
 100         IF(KCA.NE.0) KCFI(I)=KCHG(KCA,2)*ISIGN(1,KFLS(I))
 101         NFIS(I)=0
 102         IF(KCFI(I).NE.0) THEN
 103           IF(I.EQ.1) IPFS=IPUS1
 104           IF(I.EQ.2) IPFS=IPUS2
 105           DO 130 J=1,2
 106           ICSI=MOD(K(IPFS,3+J),MSTU(5))
 107           IF(ICSI.GT.0.AND.ICSI.NE.IPUS1.AND.ICSI.NE.IPUS2.AND.
 108      &    (KCFI(I).EQ.(-1)**(J+1).OR.KCFI(I).EQ.2)) THEN
 109             NFIS(I)=NFIS(I)+1
 110             THEFIS(I,NFIS(I))=ULANGL(P(ICSI,3),SQRT(P(ICSI,1)**2+
 111      &      P(ICSI,2)**2))
 112             IF(I.EQ.2) THEFIS(I,NFIS(I))=PARU(1)-THEFIS(I,NFIS(I))
 113           ENDIF
 114   130     CONTINUE
 115         ENDIF
 116   140   CONTINUE
 117         IF(NFIS(1)+NFIS(2).EQ.0) MFIS=0
 118       ENDIF
 119
 120 C...Pick up leg with highest virtuality.
 121   150 N=N+1
 122       JT=1
 123       IF(N.GT.NS+1.AND.Q2S(2).GT.Q2S(1)) JT=2
 124       IF(MORE(JT).EQ.0) JT=3-JT
 125       KFLB=KFLS(JT)
 126       XB=XS(JT)
 127       DO 160 KFL=-25,25
 128       XFB(KFL)=XFS(JT,KFL)
 129   160 CONTINUE
 130       DSHR=2D0*SQRT(DSH)
 131       DSHZ=DSH/DBLE(ZS(JT))
 132
 133 C...Check if allowed to branch.
 134       MCEV=0
 135       IF(IABS(KFLB).LE.10.OR.KFLB.EQ.21) THEN
 136         MCEV=1
 137         XEC=MAX(XEC0,XB*(1./(1.-PARP(66))-1.))
 138         IF(XB.GE.1.-2.*XEC) MCEV=0
 139       ENDIF
 140       MEEV=0
 141       IF(MINT(44+JT).EQ.3) THEN
 142         MEEV=1
 143         IF(XB.GE.1.-2.*XEE) MEEV=0
 144         IF((IABS(KFLB).LE.10.OR.KFLB.EQ.21).AND.XB.GE.1.-2.*XEC) MEEV=0
 145 C***Currently kill QED shower for resolved photoproduction.
 146         IF(MINT(18+JT).EQ.1) MEEV=0
 147 C***Currently kill shower for W inside electron.
 148         IF(IABS(KFLB).EQ.24) THEN
 149           MCEV=0
 150           MEEV=0
 151         ENDIF
 152       ENDIF
 153       IF(MCEV.EQ.0.AND.MEEV.EQ.0) THEN
 154         Q2B=0.
 155         GOTO 250
 156       ENDIF
 157
 158 C...Maximum Q2 with or without Q2 ordering. Effective Lambda and n_f.
 159       Q2B=Q2S(JT)
 160       TEVCB=TEVCSV(JT)
 161       TEVEB=TEVESV(JT)
 162       IF(MSTP(62).LE.1) THEN
 163         IF(ZS(JT).GT.0.99999) THEN
 164           Q2B=Q2S(JT)
 165         ELSE
 166           Q2B=0.5*(1./ZS(JT)+1.)*Q2S(JT)+0.5*(1./ZS(JT)-1.)*(Q2S(3-JT)-
 167      &    SNGL(DSH)+SQRT((SNGL(DSH)+Q2S(1)+Q2S(2))**2+8.*Q2S(1)*Q2S(2)*
 168      &    ZS(JT)/(1.-ZS(JT))))
 169         ENDIF
 170         IF(MCEV.EQ.1) TEVCB=LOG(FQ2C*Q2B/ALAM(JT)**2)
 171         IF(MEEV.EQ.1) TEVEB=LOG(Q2B/SPME)
 172       ENDIF
 173       IF(MCEV.EQ.1) THEN
 174         ALSDUM=ULALPS(FQ2C*Q2B)
 175         TEVCB=TEVCB+2.*LOG(ALAM(JT)/PARU(117))
 176         ALAM(JT)=PARU(117)
 177         B0=(33.-2.*MSTU(118))/6.
 178       ENDIF
 179       TEVCBS=TEVCB
 180       TEVEBS=TEVEB
 181
 182 C...Select side for interference with final state partons.
 183       IF(MFIS.GE.1.AND.N.LE.NS+2) THEN
 184         IFI=N-NS
 185         ISFI(IFI)=0
 186         IF(IABS(KCFI(IFI)).EQ.1.AND.NFIS(IFI).EQ.1) THEN
 187           ISFI(IFI)=1
 188         ELSEIF(KCFI(IFI).EQ.2.AND.NFIS(IFI).EQ.1) THEN
 189           IF(RLU(0).GT.0.5) ISFI(IFI)=1
 190         ELSEIF(KCFI(IFI).EQ.2.AND.NFIS(IFI).EQ.2) THEN
 191           ISFI(IFI)=1
 192           IF(RLU(0).GT.0.5) ISFI(IFI)=2
 193         ENDIF
 194       ENDIF
 195
 196 C...Calculate Altarelli-Parisi weights.
 197       DO 170 KFL=-25,25
 198       WTAPC(KFL)=0.
 199       WTAPE(KFL)=0.
 200       WTSF(KFL)=0.
 201   170 CONTINUE
 202 C...q -> q, g -> q.
 203       IF(IABS(KFLB).LE.10) THEN
 204         WTAPC(KFLB)=(8./3.)*LOG((1.-XEC-XB)*(XB+XEC)/(XEC*(1.-XEC)))
 205         WTAPC(21)=0.5*(XB/(XB+XEC)-XB/(1.-XEC))
 206 C...f -> f, gamma -> f.
 207       ELSEIF(IABS(KFLB).LE.20) THEN
 208         WTAPF1=LOG((1.-XEE-XB)*(XB+XEE)/(XEE*(1.-XEE)))
 209         WTAPF2=LOG((1.-XEE-XB)*(1.-XEE)/(XEE*(XB+XEE)))
 210         WTAPE(KFLB)=2.*(WTAPF1+WTAPF2)
 211         IF(MSTP(12).GE.1) WTAPE(22)=XB/(XB+XEE)-XB/(1.-XEE)
 212 C...f -> g, g -> g.
 213       ELSEIF(KFLB.EQ.21) THEN
 214         WTAPQ=(16./3.)*(SQRT((1.-XEC)/XB)-SQRT((XB+XEC)/XB))
 215         DO 180 KFL=1,MSTP(58)
 216         WTAPC(KFL)=WTAPQ
 217         WTAPC(-KFL)=WTAPQ
 218   180   CONTINUE
 219         WTAPC(21)=6.*LOG((1.-XEC-XB)/XEC)
 220 C...f -> gamma, W+, W-.
 221       ELSEIF(KFLB.EQ.22) THEN
 222         WTAPF=LOG((1.-XEE-XB)*(1.-XEE)/(XEE*(XB+XEE)))/XB
 223         WTAPE(11)=WTAPF
 224         WTAPE(-11)=WTAPF
 225       ELSEIF(KFLB.EQ.24) THEN
 226         WTAPE(-11)=1./(4.*PARU(102))*LOG((1.-XEE-XB)*(1.-XEE)/
 227      &  (XEE*(XB+XEE)))/XB
 228       ELSEIF(KFLB.EQ.-24) THEN
 229         WTAPE(11)=1./(4.*PARU(102))*LOG((1.-XEE-XB)*(1.-XEE)/
 230      &  (XEE*(XB+XEE)))/XB
 231       ENDIF
 232
 233 C...Calculate structure function weights and sum.
 234       NTRY=0
 235   190 NTRY=NTRY+1
 236       IF(NTRY.GT.500) THEN
 237         MINT(51)=1
 238         RETURN
 239       ENDIF
 240       WTSUMC=0.
 241       WTSUME=0.
 242       XFBO=MAX(1E-10,XFB(KFLB))
 243       DO 200 KFL=-25,25
 244       WTSF(KFL)=XFB(KFL)/XFBO
 245       WTSUMC=WTSUMC+WTAPC(KFL)*WTSF(KFL)
 246       WTSUME=WTSUME+WTAPE(KFL)*WTSF(KFL)
 247   200 CONTINUE
 248       WTSUMC=MAX(0.0001,WTSUMC)
 249       WTSUME=MAX(0.0001/FWTE,WTSUME)
 250
 251 C...Choose new t: fix alpha_s, alpha_s(Q^2), alpha_s(k_T^2).
 252       NTRY2=0
 253   210 NTRY2=NTRY2+1
 254       IF(NTRY2.GT.500) THEN
 255         MINT(51)=1
 256         RETURN
 257       ENDIF
 258       IF(MCEV.EQ.1) THEN
 259         IF(MSTP(64).LE.0) THEN
 260           TEVCB=TEVCB+LOG(RLU(0))*PARU(2)/(PARU(111)*WTSUMC)
 261         ELSEIF(MSTP(64).EQ.1) THEN
 262           TEVCB=TEVCB*EXP(MAX(-50.,LOG(RLU(0))*B0/WTSUMC))
 263         ELSE
 264           TEVCB=TEVCB*EXP(MAX(-50.,LOG(RLU(0))*B0/(5.*WTSUMC)))
 265         ENDIF
 266       ENDIF
 267       IF(MEEV.EQ.1) THEN
 268         TEVEB=TEVEB*EXP(MAX(-50.,LOG(RLU(0))*PARU(2)/
 269      &  (PARU(101)*FWTE*WTSUME*TEMX)))
 270       ENDIF
 271
 272 C...Translate t into Q2 scale; choose between QCD and QED evolution.
 273   220 IF(MCEV.EQ.1) Q2CB=ALAM(JT)**2*EXP(MAX(-50.,TEVCB))/FQ2C
 274       IF(MEEV.EQ.1) Q2EB=SPME*EXP(MAX(-50.,TEVEB))
 275       MCE=0
 276       IF(MCEV.EQ.0.AND.MEEV.EQ.0) THEN
 277       ELSEIF(MCEV.EQ.1.AND.MEEV.EQ.0) THEN
 278         IF(Q2CB.GT.Q2MNCS(JT)) MCE=1
 279       ELSEIF(MCEV.EQ.0.AND.MEEV.EQ.1) THEN
 280         IF(Q2EB.GT.Q2MNE) MCE=2
 281       ELSEIF(Q2MNCS(JT).GT.Q2MNE) THEN
 282         MCE=1
 283         IF(Q2EB.GT.Q2CB.OR.Q2CB.LE.Q2MNCS(JT)) MCE=2
 284         IF(MCE.EQ.2.AND.Q2EB.LE.Q2MNE) MCE=0
 285       ELSE
 286         MCE=2
 287         IF(Q2CB.GT.Q2EB.OR.Q2EB.LE.Q2MNE) MCE=1
 288         IF(MCE.EQ.1.AND.Q2CB.LE.Q2MNCS(JT)) MCE=0
 289       ENDIF
 290
 291 C...Evolution possibly ended. Update t values.
 292       IF(MCE.EQ.0) THEN
 293         Q2B=0.
 294         GOTO 250
 295       ELSEIF(MCE.EQ.1) THEN
 296         Q2B=Q2CB
 297         Q2REF=FQ2C*Q2B
 298         IF(MEEV.EQ.1) TEVEB=LOG(Q2B/SPME)
 299       ELSE
 300         Q2B=Q2EB
 301         Q2REF=Q2B
 302         IF(MCEV.EQ.1) TEVCB=LOG(FQ2C*Q2B/ALAM(JT)**2)
 303       ENDIF
 304
 305 C...Select flavour for branching parton.
 306       IF(MCE.EQ.1) WTRAN=RLU(0)*WTSUMC
 307       IF(MCE.EQ.2) WTRAN=RLU(0)*WTSUME
 308       KFLA=-25
 309   230 KFLA=KFLA+1
 310       IF(MCE.EQ.1) WTRAN=WTRAN-WTAPC(KFLA)*WTSF(KFLA)
 311       IF(MCE.EQ.2) WTRAN=WTRAN-WTAPE(KFLA)*WTSF(KFLA)
 312       IF(KFLA.LE.24.AND.WTRAN.GT.0.) GOTO 230
 313       IF(KFLA.EQ.25) THEN
 314         Q2B=0.
 315         GOTO 250
 316       ENDIF
 317
 318 C...Choose z value and corrective weight.
 319       WTZ=0.
 320 C...q -> q + g.
 321       IF(IABS(KFLA).LE.10.AND.IABS(KFLB).LE.10) THEN
 322         Z=1.-((1.-XB-XEC)/(1.-XEC))*
 323      &  (XEC*(1.-XEC)/((XB+XEC)*(1.-XB-XEC)))**RLU(0)
 324         WTZ=0.5*(1.+Z**2)
 325 C...q -> g + q.
 326       ELSEIF(IABS(KFLA).LE.10.AND.KFLB.EQ.21) THEN
 327         Z=XB/(SQRT(XB+XEC)+RLU(0)*(SQRT(1.-XEC)-SQRT(XB+XEC)))**2
 328         WTZ=0.5*(1.+(1.-Z)**2)*SQRT(Z)
 329 C...f -> f + gamma.
 330       ELSEIF(IABS(KFLA).LE.20.AND.IABS(KFLB).LE.20) THEN
 331         IF(WTAPF1.GT.RLU(0)*(WTAPF1+WTAPF2)) THEN
 332           Z=1.-((1.-XB-XEE)/(1.-XEE))*
 333      &    (XEE*(1.-XEE)/((XB+XEE)*(1.-XB-XEE)))**RLU(0)
 334         ELSE
 335           Z=XB+XB*(XEE/(1.-XEE))*
 336      &    ((1.-XB-XEE)*(1.-XEE)/(XEE*(XB+XEE)))**RLU(0)
 337         ENDIF
 338         WTZ=0.5*(1.+Z**2)*(Z-XB)/(1.-XB)
 339 C...f -> gamma + f.
 340       ELSEIF(IABS(KFLA).LE.20.AND.KFLB.EQ.22) THEN
 341         Z=XB+XB*(XEE/(1.-XEE))*
 342      &  ((1.-XB-XEE)*(1.-XEE)/(XEE*(XB+XEE)))**RLU(0)
 343         WTZ=0.5*(1.+(1.-Z)**2)*XB*(Z-XB)/Z
 344 C...f -> W+- + f'.
 345       ELSEIF(IABS(KFLA).LE.20.AND.IABS(KFLB).EQ.24) THEN
 346         Z=XB+XB*(XEE/(1.-XEE))*
 347      &  ((1.-XB-XEE)*(1.-XEE)/(XEE*(XB+XEE)))**RLU(0)
 348         WTZ=0.5*(1.+(1.-Z)**2)*(XB*(Z-XB)/Z)*(Q2B/(Q2B+PMAS(24,1)**2))
 349 C...g -> q + q~.
 350       ELSEIF(KFLA.EQ.21.AND.IABS(KFLB).LE.10) THEN
 351         Z=XB/(1.-XEC)+RLU(0)*(XB/(XB+XEC)-XB/(1.-XEC))
 352         WTZ=1.-2.*Z*(1.-Z)
 353 C...g -> g + g.
 354       ELSEIF(KFLA.EQ.21.AND.KFLB.EQ.21) THEN
 355         Z=1./(1.+((1.-XEC-XB)/XB)*(XEC/(1.-XEC-XB))**RLU(0))
 356         WTZ=(1.-Z*(1.-Z))**2
 357 C...gamma -> f + f~.
 358       ELSEIF(KFLA.EQ.22.AND.IABS(KFLB).LE.20) THEN
 359         Z=XB/(1.-XEE)+RLU(0)*(XB/(XB+XEE)-XB/(1.-XEE))
 360         WTZ=1.-2.*Z*(1.-Z)
 361       ENDIF
 362       IF(MCE.EQ.2) WTZ=(WTZ/FWTE)*(TEVEB/TEMX)
 363
 364 C...Option with resummation of soft gluon emission as effective z shift.
 365       IF(MCE.EQ.1) THEN
 366         IF(MSTP(65).GE.1) THEN
 367           RSOFT=6.
 368           IF(KFLB.NE.21) RSOFT=8./3.
 369           Z=Z*(TEVCB/TEVCSV(JT))**(RSOFT*XEC/((XB+XEC)*B0))
 370           IF(Z.LE.XB) GOTO 210
 371         ENDIF
 372
 373 C...Option with alpha_s(k_T^2): demand k_T^2 > cutoff, reweight.
 374         IF(MSTP(64).GE.2) THEN
 375           IF((1.-Z)*Q2B.LT.Q2MNCS(JT)) GOTO 210
 376           ALPRAT=TEVCB/(TEVCB+LOG(1.-Z))
 377           IF(ALPRAT.LT.5.*RLU(0)) GOTO 210
 378           IF(ALPRAT.GT.5.) WTZ=WTZ*ALPRAT/5.
 379         ENDIF
 380
 381 C...Impose angular constraint in first branching from interference
 382 C...with final state partons.
 383         IF(MFIS.GE.1.AND.N.LE.NS+2.AND.NTRY2.LT.200) THEN
 384           THE2D=(4.*Q2B)/(DSH*(1.-Z))
 385           IF(N.EQ.NS+1.AND.ISFI(1).GE.1) THEN
 386             IF(THE2D.GT.THEFIS(1,ISFI(1))**2) GOTO 210
 387           ELSEIF(N.EQ.NS+2.AND.ISFI(2).GE.1) THEN
 388             IF(THE2D.GT.THEFIS(2,ISFI(2))**2) GOTO 210
 389           ENDIF
 390         ENDIF
 391
 392 C...Option with angular ordering requirement.
 393         IF(MSTP(62).GE.3.AND.NTRY2.LT.200) THEN
 394           THE2T=(4.*Z**2*Q2B)/(VINT(2)*(1.-Z)*XB**2)
 395           IF(THE2T.GT.THE2(JT)) GOTO 210
 396         ENDIF
 397       ENDIF
 398
 399 C...Weighting with new structure functions.
 400       MINT(105)=MINT(102+JT)
 401       MINT(109)=MINT(106+JT)
 402       IF(MSTP(57).LE.1) THEN
 403         CALL PYSTFU(KFBEAM(JT),XB,Q2REF,XFN)
 404       ELSE
 405         CALL PYSTFL(KFBEAM(JT),XB,Q2REF,XFN)
 406       ENDIF
 407       XFBN=XFN(KFLB)
 408       IF(XFBN.LT.1E-20) THEN
 409         IF(KFLA.EQ.KFLB) THEN
 410           TEVCB=TEVCBS
 411           TEVEB=TEVEBS
 412           WTAPC(KFLB)=0.
 413           WTAPE(KFLB)=0.
 414           GOTO 190
 415         ELSEIF(MCE.EQ.1.AND.TEVCBS-TEVCB.GT.0.2) THEN
 416           TEVCB=0.5*(TEVCBS+TEVCB)
 417           GOTO 220
 418         ELSEIF(MCE.EQ.2.AND.TEVEBS-TEVEB.GT.0.2) THEN
 419           TEVEB=0.5*(TEVEBS+TEVEB)
 420           GOTO 220
 421         ELSE
 422           XFBN=1E-10
 423           XFN(KFLB)=XFBN
 424         ENDIF
 425       ENDIF
 426       DO 240 KFL=-25,25
 427       XFB(KFL)=XFN(KFL)
 428   240 CONTINUE
 429       XA=XB/Z
 430       IF(MSTP(57).LE.1) THEN
 431         CALL PYSTFU(KFBEAM(JT),XA,Q2REF,XFA)
 432       ELSE
 433         CALL PYSTFL(KFBEAM(JT),XA,Q2REF,XFA)
 434       ENDIF
 435       XFAN=XFA(KFLA)
 436       IF(XFAN.LT.1E-20) GOTO 190
 437       WTSFA=WTSF(KFLA)
 438       IF(WTZ*XFAN/XFBN.LT.RLU(0)*WTSFA) GOTO 190
 439
 440 C...Define two hard scatterers in their CM-frame.
 441   250 IF(N.EQ.NS+2) THEN
 442         DQ2(JT)=Q2B
 443         DPLCM=SQRT((DSH+DQ2(1)+DQ2(2))**2-4D0*DQ2(1)*DQ2(2))/DSHR
 444         DO 270 JR=1,2
 445         I=NS+JR
 446         IF(JR.EQ.1) IPO=IPUS1
 447         IF(JR.EQ.2) IPO=IPUS2
 448         DO 260 J=1,5
 449         K(I,J)=0
 450         P(I,J)=0.
 451         V(I,J)=0.
 452   260   CONTINUE
 453         K(I,1)=14
 454         K(I,2)=KFLS(JR+2)
 455         K(I,4)=IPO
 456         K(I,5)=IPO
 457         P(I,3)=DPLCM*(-1)**(JR+1)
 458         P(I,4)=(DSH+DQ2(3-JR)-DQ2(JR))/DSHR
 459         P(I,5)=-SQRT(SNGL(DQ2(JR)))
 460         K(IPO,1)=14
 461         K(IPO,3)=I
 462         K(IPO,4)=MOD(K(IPO,4),MSTU(5))+MSTU(5)*I
 463         K(IPO,5)=MOD(K(IPO,5),MSTU(5))+MSTU(5)*I
 464   270   CONTINUE
 465
 466 C...Find maximum allowed mass of timelike parton.
 467       ELSEIF(N.GT.NS+2) THEN
 468         JR=3-JT
 469         DQ2(3)=Q2B
 470         DPC(1)=P(IS(1),4)
 471         DPC(2)=P(IS(2),4)
 472         DPC(3)=0.5*(ABS(P(IS(1),3))+ABS(P(IS(2),3)))
 473         DPD(1)=DSH+DQ2(JR)+DQ2(JT)
 474         DPD(2)=DSHZ+DQ2(JR)+DQ2(3)
 475         DPD(3)=SQRT(DPD(1)**2-4D0*DQ2(JR)*DQ2(JT))
 476         DPD(4)=SQRT(DPD(2)**2-4D0*DQ2(JR)*DQ2(3))
 477         IKIN=0
 478         IF(Q2S(JR).GE.0.25*Q2MNC.AND.DPD(1)-DPD(3).GE.
 479      &  1D-10*DPD(1)) IKIN=1
 480         IF(IKIN.EQ.0) DMSMA=(DQ2(JT)/DBLE(ZS(JT))-DQ2(3))*(DSH/
 481      &  (DSH+DQ2(JT))-DSH/(DSHZ+DQ2(3)))
 482         IF(IKIN.EQ.1) DMSMA=(DPD(1)*DPD(2)-DPD(3)*DPD(4))/(2.*
 483      &  DQ2(JR))-DQ2(JT)-DQ2(3)
 484
 485 C...Generate timelike parton shower (if required).
 486         IT=N
 487         DO 280 J=1,5
 488         K(IT,J)=0
 489         P(IT,J)=0.
 490         V(IT,J)=0.
 491   280   CONTINUE
 492         K(IT,1)=3
 493 C...f -> f + g (gamma).
 494         IF(IABS(KFLB).LE.20.AND.IABS(KFLS(JT+2)).LE.20) THEN
 495           K(IT,2)=21
 496           IF(IABS(KFLB).GE.11) K(IT,2)=22
 497 C...f -> g (gamma, W+-) + f.
 498         ELSEIF(IABS(KFLB).LE.20.AND.IABS(KFLS(JT+2)).GT.20) THEN
 499           K(IT,2)=KFLB
 500           IF(KFLS(JT+2).EQ.24) THEN
 501             K(IT,2)=-12
 502           ELSEIF(KFLS(JT+2).EQ.-24) THEN
 503             K(IT,2)=12
 504           ENDIF
 505 C...g (gamma) -> f + f~, g + g.
 506         ELSE
 507           K(IT,2)=-KFLS(JT+2)
 508           IF(KFLS(JT+2).GT.20) K(IT,2)=KFLS(JT+2)
 509         ENDIF
 510         P(IT,5)=ULMASS(K(IT,2))
 511         IF(SNGL(DMSMA).LE.P(IT,5)**2) GOTO 100
 512         IF(MSTP(63).GE.1.AND.MCE.EQ.1) THEN
 513           MSTJ48=MSTJ(48)
 514           PARJ85=PARJ(85)
 515           P(IT,4)=(DSHZ-DSH-P(IT,5)**2)/DSHR
 516           P(IT,3)=SQRT(P(IT,4)**2-P(IT,5)**2)
 517           IF(MSTP(63).EQ.1) THEN
 518             Q2TIM=DMSMA
 519           ELSEIF(MSTP(63).EQ.2) THEN
 520             Q2TIM=MIN(SNGL(DMSMA),PARP(71)*Q2S(JT))
 521           ELSE
 522             Q2TIM=DMSMA
 523             MSTJ(48)=1
 524             IF(IKIN.EQ.0) DPT2=DMSMA*(DSHZ+DQ2(3))/(DSH+DQ2(JT))
 525             IF(IKIN.EQ.1) DPT2=DMSMA*(0.5*DPD(1)*DPD(2)+0.5*DPD(3)*
 526      &      DPD(4)-DQ2(JR)*(DQ2(JT)+DQ2(3)))/(4.*DSH*DPC(3)**2)
 527             PARJ(85)=SQRT(MAX(0.,SNGL(DPT2)))*
 528      &      (1./P(IT,4)+1./P(IS(JT),4))
 529           ENDIF
 530           CALL LUSHOW(IT,0,SQRT(Q2TIM))
 531           MSTJ(48)=MSTJ48
 532           PARJ(85)=PARJ85
 533           IF(N.GE.IT+1) P(IT,5)=P(IT+1,5)
 534         ENDIF
 535
 536 C...Reconstruct kinematics of branching: timelike parton shower.
 537         DMS=P(IT,5)**2
 538         IF(IKIN.EQ.0) DPT2=(DMSMA-DMS)*(DSHZ+DQ2(3))/(DSH+DQ2(JT))
 539         IF(IKIN.EQ.1) DPT2=(DMSMA-DMS)*(0.5*DPD(1)*DPD(2)+0.5*DPD(3)*
 540      &  DPD(4)-DQ2(JR)*(DQ2(JT)+DQ2(3)+DMS))/(4.*DSH*DPC(3)**2)
 541         IF(DPT2.LT.0.) GOTO 100
 542         DPB(1)=(0.5*DPD(2)-DPC(JR)*(DSHZ+DQ2(JR)-DQ2(JT)-DMS)/
 543      &  DSHR)/DPC(3)-DPC(3)
 544         P(IT,1)=SQRT(SNGL(DPT2))
 545         P(IT,3)=DPB(1)*(-1)**(JT+1)
 546         P(IT,4)=SQRT(DPT2+DPB(1)**2+DMS)
 547         IF(N.GE.IT+1) THEN
 548           DPB(1)=SQRT(DPB(1)**2+DPT2)
 549           DPB(2)=SQRT(DPB(1)**2+DMS)
 550           DPB(3)=P(IT+1,3)
 551           DPB(4)=SQRT(DPB(3)**2+DMS)
 552           DBEZ=(DPB(4)*DPB(1)-DPB(3)*DPB(2))/(DPB(4)*DPB(2)-DPB(3)*
 553      &    DPB(1))
 554           CALL LUDBRB(IT+1,N,0.,0.,0D0,0D0,DBEZ)
 555           THE=ULANGL(P(IT,3),P(IT,1))
 556           CALL LUDBRB(IT+1,N,THE,0.,0D0,0D0,0D0)
 557         ENDIF
 558
 559 C...Reconstruct kinematics of branching: spacelike parton.
 560         DO 290 J=1,5
 561         K(N+1,J)=0
 562         P(N+1,J)=0.
 563         V(N+1,J)=0.
 564   290   CONTINUE
 565         K(N+1,1)=14
 566         K(N+1,2)=KFLB
 567         P(N+1,1)=P(IT,1)
 568         P(N+1,3)=P(IT,3)+P(IS(JT),3)
 569         P(N+1,4)=P(IT,4)+P(IS(JT),4)
 570         P(N+1,5)=-SQRT(SNGL(DQ2(3)))
 571
 572 C...Define colour flow of branching.
 573         K(IS(JT),3)=N+1
 574         K(IT,3)=N+1
 575         IM1=N+1
 576         IM2=N+1
 577 C...f -> f + gamma (Z, W).
 578         IF(IABS(K(IT,2)).GE.22) THEN
 579           K(IT,1)=1
 580           ID1=IS(JT)
 581           ID2=IS(JT)
 582 C...f -> gamma (Z, W) + f.
 583         ELSEIF(IABS(K(IS(JT),2)).GE.22) THEN
 584           ID1=IT
 585           ID2=IT
 586 C...gamma -> q + q~, g + g.
 587         ELSEIF(K(N+1,2).EQ.22) THEN
 588           ID1=IS(JT)
 589           ID2=IT
 590           IM1=ID2
 591           IM2=ID1
 592 C...q -> q + g.
 593         ELSEIF(K(N+1,2).GT.0.AND.K(N+1,2).NE.21.AND.K(IT,2).EQ.21) THEN
 594           ID1=IT
 595           ID2=IS(JT)
 596 C...q -> g + q.
 597         ELSEIF(K(N+1,2).GT.0.AND.K(N+1,2).NE.21) THEN
 598           ID1=IS(JT)
 599           ID2=IT
 600 C...q~ -> q~ + g.
 601         ELSEIF(K(N+1,2).LT.0.AND.K(IT,2).EQ.21) THEN
 602           ID1=IS(JT)
 603           ID2=IT
 604 C...q~ -> g + q~.
 605         ELSEIF(K(N+1,2).LT.0) THEN
 606           ID1=IT
 607           ID2=IS(JT)
 608 C...g -> g + g; g -> q + q~.
 609         ELSEIF((K(IT,2).EQ.21.AND.RLU(0).GT.0.5).OR.K(IT,2).LT.0) THEN
 610           ID1=IS(JT)
 611           ID2=IT
 612         ELSE
 613           ID1=IT
 614           ID2=IS(JT)
 615         ENDIF
 616         IF(IM1.EQ.N+1) K(IM1,4)=K(IM1,4)+ID1
 617         IF(IM2.EQ.N+1) K(IM2,5)=K(IM2,5)+ID2
 618         K(ID1,4)=K(ID1,4)+MSTU(5)*IM1
 619         K(ID2,5)=K(ID2,5)+MSTU(5)*IM2
 620         IF(ID1.NE.ID2) THEN
 621           K(ID1,5)=K(ID1,5)+MSTU(5)*ID2
 622           K(ID2,4)=K(ID2,4)+MSTU(5)*ID1
 623         ENDIF
 624         N=N+1
 625
 626 C...Boost to new CM-frame.
 627         DBSVX=DBLE((P(N,1)+P(IS(JR),1))/(P(N,4)+P(IS(JR),4)))
 628         DBSVZ=DBLE((P(N,3)+P(IS(JR),3))/(P(N,4)+P(IS(JR),4)))
 629         IF(DBSVX**2+DBSVZ**2.GE.1D0) GOTO 100
 630         CALL LUDBRB(NS+1,N,0.,0.,-DBSVX,0D0,-DBSVZ)
 631         IR=N+(JT-1)*(IS(1)-N)
 632         CALL LUDBRB(NS+1,N,-ULANGL(P(IR,3),P(IR,1)),PARU(2)*RLU(0),
 633      &  0D0,0D0,0D0)
 634       ENDIF
 635
 636 C...Update kinematics variables.
 637       IS(JT)=N
 638       DQ2(JT)=Q2B
 639       IF(MSTP(62).GE.3) THE2(JT)=THE2T
 640       DSH=DSHZ
 641
 642 C...Save quantities; loop back.
 643       Q2S(JT)=Q2B
 644       IF((MCEV.EQ.1.AND.Q2B.GE.0.25*Q2MNC).OR.
 645      &(MEEV.EQ.1.AND.Q2B.GE.Q2MNE)) THEN
 646         KFLS(JT+2)=KFLS(JT)
 647         KFLS(JT)=KFLA
 648         XS(JT)=XA
 649         ZS(JT)=Z
 650         DO 300 KFL=-25,25
 651         XFS(JT,KFL)=XFA(KFL)
 652   300   CONTINUE
 653         TEVCSV(JT)=TEVCB
 654         TEVESV(JT)=TEVEB
 655       ELSE
 656         MORE(JT)=0
 657         IF(JT.EQ.1) IPU1=N
 658         IF(JT.EQ.2) IPU2=N
 659       ENDIF
 660       IF(N.GT.MSTU(4)-MSTU(32)-10) THEN
 661         CALL LUERRM(11,'(PYSSPA:) no more memory left in LUJETS')
 662         IF(MSTU(21).GE.1) N=NS
 663         IF(MSTU(21).GE.1) RETURN
 664       ENDIF
 665       IF(MORE(1).EQ.1.OR.MORE(2).EQ.1) GOTO 150
 666
 667 C...Boost hard scattering partons to frame of shower initiators.
 668       DO 310 J=1,3
 669       ROBO(J+2)=(P(NS+1,J)+P(NS+2,J))/(P(NS+1,4)+P(NS+2,4))
 670   310 CONTINUE
 671       K(N+2,1)=1
 672       DO 320 J=1,5
 673       P(N+2,J)=P(NS+1,J)
 674   320 CONTINUE
 675       ROBOT=ROBO(3)**2+ROBO(4)**2+ROBO(5)**2
 676       IF(ROBOT.GE.0.999999) THEN
 677         ROBOT=1.00001*SQRT(ROBOT)
 678         ROBO(3)=ROBO(3)/ROBOT
 679         ROBO(4)=ROBO(4)/ROBOT
 680         ROBO(5)=ROBO(5)/ROBOT
 681       ENDIF
 682       CALL LUDBRB(N+2,N+2,0.,0.,-DBLE(ROBO(3)),-DBLE(ROBO(4)),
 683      &-DBLE(ROBO(5)))
 684       ROBO(2)=ULANGL(P(N+2,1),P(N+2,2))
 685       ROBO(1)=ULANGL(P(N+2,3),SQRT(P(N+2,1)**2+P(N+2,2)**2))
 686       CALL LUDBRB(MINT(83)+5,NS,ROBO(1),ROBO(2),DBLE(ROBO(3)),
 687      &DBLE(ROBO(4)),DBLE(ROBO(5)))
 688
 689 C...Store user information. Reset Lambda value.
 690       K(IPU1,3)=MINT(83)+3
 691       K(IPU2,3)=MINT(83)+4
 692       DO 330 JT=1,2
 693       MINT(12+JT)=KFLS(JT)
 694       VINT(140+JT)=XS(JT)
 695       IF(MINT(18+JT).EQ.1) VINT(140+JT)=VINT(154+JT)*XS(JT)
 696   330 CONTINUE
 697       PARU(112)=ALAMS
 698
 699       RETURN
 700       END