PYTHIA/pythia/pyrand.F

   1
   2 C*********************************************************************
   3
   4       SUBROUTINE PYRAND
   5
   6 C...Generates quantities characterizing the high-pT scattering at the
   7 C...parton level according to the matrix elements. Chooses incoming,
   8 C...reacting partons, their momentum fractions and one of the possible
   9 C...subprocesses.
  10       COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
  11       COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
  12       COMMON/PYSUBS/MSEL,MSUB(200),KFIN(2,-40:40),CKIN(200)
  13       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
  14       COMMON/PYINT1/MINT(400),VINT(400)
  15       COMMON/PYINT2/ISET(200),KFPR(200,2),COEF(200,20),ICOL(40,4,2)
  16       COMMON/PYINT3/XSFX(2,-40:40),ISIG(1000,3),SIGH(1000)
  17       COMMON/PYINT4/WIDP(21:40,0:40),WIDE(21:40,0:40),WIDS(21:40,3)
  18       COMMON/PYINT5/NGEN(0:200,3),XSEC(0:200,3)
  19       COMMON/PYINT7/SIGT(0:6,0:6,0:5)
  20       COMMON/PYINT9/DXSEC(0:200)
  21       DOUBLE PRECISION DXSEC
  22       COMMON/PYUPPR/NUP,KUP(20,7),PUP(20,5),NFUP,IFUP(10,2),Q2UP(0:10)
  23       SAVE /LUDAT1/,/LUDAT2/
  24       SAVE /PYSUBS/,/PYPARS/,/PYINT1/,/PYINT2/,/PYINT3/,/PYINT4/,
  25      &/PYINT5/,/PYINT7/,/PYINT9/,/PYUPPR/
  26       DIMENSION XPQ(-25:25),PMM(2),PDIF(4),BHAD(4)
  27
  28 C...Parameters and data used in elastic/diffractive treatment.
  29       DATA EPS/0.0808/, ALP/0.25/, CRES/2./, PMRC/1.062/, SMP/0.880/
  30       DATA BHAD/2.3,1.4,1.4,0.23/
  31
  32 C...Initial values, specifically for (first) semihard interaction.
  33       MINT(10)=0
  34       MINT(17)=0
  35       MINT(18)=0
  36       VINT(143)=1.
  37       VINT(144)=1.
  38       MFAIL=0
  39       IF(MSTP(171).EQ.1.AND.MSTP(172).EQ.2) MFAIL=1
  40       ISUB=0
  41       LOOP=0
  42   100 LOOP=LOOP+1
  43       MINT(51)=0
  44
  45 C...Choice of process type - first event of pileup.
  46       IF(MINT(82).EQ.1.AND.(ISUB.LE.90.OR.ISUB.GT.96)) THEN
  47
  48 C...For gamma-p or gamma-gamma first pick between alternatives.
  49         IF(MINT(121).GT.1) CALL PYSAVE(4,IGA)
  50         MINT(122)=IGA
  51
  52 C...For gamma + gamma with different nature, flip at random.
  53         IF(MINT(11).EQ.22.AND.MINT(12).EQ.22.AND.MINT(123).GE.4.AND.
  54      &  RLU(0).GT.0.5) THEN
  55           MINTSV=MINT(41)
  56           MINT(41)=MINT(42)
  57           MINT(42)=MINTSV
  58           MINTSV=MINT(45)
  59           MINT(45)=MINT(46)
  60           MINT(46)=MINTSV
  61           MINTSV=MINT(107)
  62           MINT(107)=MINT(108)
  63           MINT(108)=MINTSV
  64           IF(MINT(47).EQ.2.OR.MINT(47).EQ.3) MINT(47)=5-MINT(47)
  65         ENDIF
  66
  67 C...Pick process type.
  68         RSUB=XSEC(0,1)*RLU(0)
  69         DO 110 I=1,200
  70         IF(MSUB(I).NE.1) GOTO 110
  71         ISUB=I
  72         RSUB=RSUB-XSEC(I,1)
  73         IF(RSUB.LE.0.) GOTO 120
  74   110   CONTINUE
  75   120   IF(ISUB.EQ.95) ISUB=96
  76         IF(ISUB.EQ.96) CALL PYMULT(2)
  77
  78 C...Choice of inclusive process type - pileup events.
  79       ELSEIF(MINT(82).GE.2.AND.ISUB.EQ.0) THEN
  80         RSUB=VINT(131)*RLU(0)
  81         ISUB=96
  82         IF(RSUB.GT.SIGT(0,0,5)) ISUB=94
  83         IF(RSUB.GT.SIGT(0,0,5)+SIGT(0,0,4)) ISUB=93
  84         IF(RSUB.GT.SIGT(0,0,5)+SIGT(0,0,4)+SIGT(0,0,3)) ISUB=92
  85         IF(RSUB.GT.SIGT(0,0,5)+SIGT(0,0,4)+SIGT(0,0,3)+SIGT(0,0,2))
  86      &  ISUB=91
  87         IF(ISUB.EQ.96) CALL PYMULT(2)
  88       ENDIF
  89       IF(MINT(82).EQ.1) NGEN(0,1)=NGEN(0,1)+1
  90       IF(MINT(82).EQ.1) NGEN(ISUB,1)=NGEN(ISUB,1)+1
  91       IF(ISUB.EQ.96.AND.LOOP.EQ.1.AND.MINT(82).EQ.1)
  92      &NGEN(97,1)=NGEN(97,1)+1
  93       MINT(1)=ISUB
  94       ISTSB=ISET(ISUB)
  95
  96 C...Find resonances (explicit or implicit in cross-section).
  97       MINT(72)=0
  98       KFR1=0
  99       IF(ISTSB.EQ.1.OR.ISTSB.EQ.3.OR.ISTSB.EQ.5) THEN
 100         KFR1=KFPR(ISUB,1)
 101       ELSEIF(ISUB.EQ.24.OR.ISUB.EQ.25.OR.ISUB.EQ.110.OR.ISUB.EQ.165.OR.
 102      &ISUB.EQ.171.OR.ISUB.EQ.176) THEN
 103         KFR1=23
 104       ELSEIF(ISUB.EQ.23.OR.ISUB.EQ.26.OR.ISUB.EQ.166.OR.ISUB.EQ.172.OR.
 105      &ISUB.EQ.177) THEN
 106         KFR1=24
 107       ELSEIF(ISUB.GE.71.AND.ISUB.LE.77) THEN
 108         KFR1=25
 109         IF(MSTP(46).EQ.5) THEN
 110           KFR1=30
 111           PMAS(30,1)=PARP(45)
 112           PMAS(30,2)=PARP(45)**3/(96.*PARU(1)*PARP(47)**2)
 113         ENDIF
 114       ENDIF
 115       CKMX=CKIN(2)
 116       IF(CKMX.LE.0.) CKMX=VINT(1)
 117       IF(KFR1.NE.0) THEN
 118         IF(CKIN(1).GT.PMAS(KFR1,1)+20.*PMAS(KFR1,2).OR.
 119      &  CKMX.LT.PMAS(KFR1,1)-20.*PMAS(KFR1,2)) KFR1=0
 120       ENDIF
 121       IF(KFR1.NE.0) THEN
 122         TAUR1=PMAS(KFR1,1)**2/VINT(2)
 123         GAMR1=PMAS(KFR1,1)*PMAS(KFR1,2)/VINT(2)
 124         MINT(72)=1
 125         MINT(73)=KFR1
 126         VINT(73)=TAUR1
 127         VINT(74)=GAMR1
 128       ENDIF
 129       IF(ISUB.EQ.141) THEN
 130         KFR2=23
 131         TAUR2=PMAS(KFR2,1)**2/VINT(2)
 132         GAMR2=PMAS(KFR2,1)*PMAS(KFR2,2)/VINT(2)
 133         IF(CKIN(1).GT.PMAS(KFR2,1)+20.*PMAS(KFR2,2).OR.
 134      &  CKMX.LT.PMAS(KFR2,1)-20.*PMAS(KFR2,2)) KFR2=0
 135         IF(KFR2.NE.0.AND.KFR1.NE.0) THEN
 136           MINT(72)=2
 137           MINT(74)=KFR2
 138           VINT(75)=TAUR2
 139           VINT(76)=GAMR2
 140         ELSEIF(KFR2.NE.0) THEN
 141           KFR1=KFR2
 142           TAUR1=TAUR2
 143           GAMR1=GAMR2
 144           MINT(72)=1
 145           MINT(73)=KFR1
 146           VINT(73)=TAUR1
 147           VINT(74)=GAMR1
 148         ENDIF
 149       ENDIF
 150
 151 C...Find product masses and minimum pT of process,
 152 C...optionally with broadening according to a truncated Breit-Wigner.
 153       VINT(63)=0.
 154       VINT(64)=0.
 155       MINT(71)=0
 156       VINT(71)=CKIN(3)
 157       IF(MINT(82).GE.2) VINT(71)=0.
 158       VINT(80)=1.
 159       IF(ISTSB.EQ.2.OR.ISTSB.EQ.4) THEN
 160         NBW=0
 161         DO 130 I=1,2
 162         IF(KFPR(ISUB,I).EQ.0) THEN
 163         ELSEIF(MSTP(42).LE.0.OR.PMAS(LUCOMP(KFPR(ISUB,I)),2).LT.
 164      &  PARP(41)) THEN
 165           VINT(62+I)=PMAS(LUCOMP(KFPR(ISUB,I)),1)**2
 166         ELSE
 167           NBW=NBW+1
 168         ENDIF
 169   130   CONTINUE
 170         IF(NBW.GE.1) THEN
 171           CALL PYOFSH(4,0,KFPR(ISUB,1),KFPR(ISUB,2),0.,PQM3,PQM4)
 172           IF(MINT(51).EQ.1) THEN
 173             IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
 174             IF(MFAIL.EQ.1) THEN
 175               MSTI(61)=1
 176               RETURN
 177             ENDIF
 178             GOTO 100
 179           ENDIF
 180           VINT(63)=PQM3**2
 181           VINT(64)=PQM4**2
 182         ENDIF
 183         IF(MIN(VINT(63),VINT(64)).LT.CKIN(6)**2) MINT(71)=1
 184         IF(MINT(71).EQ.1) VINT(71)=MAX(CKIN(3),CKIN(5))
 185       ELSEIF(ISTSB.EQ.6) THEN
 186         CALL PYOFSH(6,0,KFPR(ISUB,1),KFPR(ISUB,2),0.,PQM3,PQM4)
 187         IF(MINT(51).EQ.1) THEN
 188           IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
 189           IF(MFAIL.EQ.1) THEN
 190             MSTI(61)=1
 191             RETURN
 192           ENDIF
 193           GOTO 100
 194         ENDIF
 195         VINT(63)=PQM3**2
 196         VINT(64)=PQM4**2
 197       ENDIF
 198
 199 C...Prepare for additional variable choices in 2 -> 3.
 200       IF(ISTSB.EQ.5) THEN
 201         VINT(201)=0.
 202         IF(KFPR(ISUB,2).GT.0) VINT(201)=PMAS(KFPR(ISUB,2),1)
 203         VINT(206)=VINT(201)
 204         VINT(204)=PMAS(23,1)
 205         IF(ISUB.EQ.124) VINT(204)=PMAS(24,1)
 206         IF(ISUB.EQ.121.OR.ISUB.EQ.122.OR.ISUB.EQ.181.OR.ISUB.EQ.182.OR.
 207      &  ISUB.EQ.186.OR.ISUB.EQ.187) VINT(204)=VINT(201)
 208         VINT(209)=VINT(204)
 209       ENDIF
 210
 211 C...Select incoming VDM particle (rho/omega/phi/J/psi).
 212       IF(ISTSB.NE.0.AND.(MINT(101).GE.2.OR.MINT(102).GE.2).AND.
 213      &(MINT(123).EQ.2.OR.MINT(123).EQ.5.OR.MINT(123).EQ.7)) THEN
 214         VRN=RLU(0)*SIGT(0,0,5)
 215         IF(MINT(101).LE.1) THEN
 216           I1MN=0
 217           I1MX=0
 218         ELSE
 219           I1MN=1
 220           I1MX=MINT(101)
 221         ENDIF
 222         IF(MINT(102).LE.1) THEN
 223           I2MN=0
 224           I2MX=0
 225         ELSE
 226           I2MN=1
 227           I2MX=MINT(102)
 228         ENDIF
 229         DO 150 I1=I1MN,I1MX
 230         KFV1=110*I1+3
 231         DO 140 I2=I2MN,I2MX
 232         KFV2=110*I2+3
 233         VRN=VRN-SIGT(I1,I2,5)
 234         IF(VRN.LE.0.) GOTO 160
 235   140   CONTINUE
 236   150   CONTINUE
 237   160   IF(MINT(101).GE.2) MINT(103)=KFV1
 238         IF(MINT(102).GE.2) MINT(104)=KFV2
 239       ENDIF
 240
 241       IF(ISTSB.EQ.0) THEN
 242 C...Elastic scattering or single or double diffractive scattering.
 243
 244 C...Select incoming particle (rho/omega/phi/J/psi for VDM) and mass.
 245         MINT(103)=MINT(11)
 246         MINT(104)=MINT(12)
 247         PMM(1)=VINT(3)
 248         PMM(2)=VINT(4)
 249         IF(MINT(101).GE.2.OR.MINT(102).GE.2) THEN
 250           JJ=ISUB-90
 251           VRN=RLU(0)*SIGT(0,0,JJ)
 252           IF(MINT(101).LE.1) THEN
 253             I1MN=0
 254             I1MX=0
 255           ELSE
 256             I1MN=1
 257             I1MX=MINT(101)
 258           ENDIF
 259           IF(MINT(102).LE.1) THEN
 260             I2MN=0
 261             I2MX=0
 262           ELSE
 263             I2MN=1
 264             I2MX=MINT(102)
 265           ENDIF
 266           DO 180 I1=I1MN,I1MX
 267           KFV1=110*I1+3
 268           DO 170 I2=I2MN,I2MX
 269           KFV2=110*I2+3
 270           VRN=VRN-SIGT(I1,I2,JJ)
 271           IF(VRN.LE.0.) GOTO 190
 272   170     CONTINUE
 273   180     CONTINUE
 274   190     IF(MINT(101).GE.2) THEN
 275             MINT(103)=KFV1
 276             PMM(1)=ULMASS(KFV1)
 277           ENDIF
 278           IF(MINT(102).GE.2) THEN
 279             MINT(104)=KFV2
 280             PMM(2)=ULMASS(KFV2)
 281           ENDIF
 282         ENDIF
 283
 284 C...Side/sides of diffractive system.
 285         MINT(17)=0
 286         MINT(18)=0
 287         IF(ISUB.EQ.92.OR.ISUB.EQ.94) MINT(17)=1
 288         IF(ISUB.EQ.93.OR.ISUB.EQ.94) MINT(18)=1
 289
 290 C...Find masses of particles and minimal masses of diffractive states.
 291         DO 200 JT=1,2
 292         PDIF(JT)=PMM(JT)
 293         VINT(66+JT)=PDIF(JT)
 294         IF(MINT(16+JT).EQ.1) PDIF(JT)=PDIF(JT)+PARP(102)
 295   200   CONTINUE
 296         SH=VINT(2)
 297         SQM1=PMM(1)**2
 298         SQM2=PMM(2)**2
 299         SQM3=PDIF(1)**2
 300         SQM4=PDIF(2)**2
 301         SMRES1=(PMM(1)+PMRC)**2
 302         SMRES2=(PMM(2)+PMRC)**2
 303
 304 C...Find elastic slope and lower limit diffractive slope.
 305         IHA=MAX(2,IABS(MINT(103))/110)
 306         IF(IHA.GE.5) IHA=1
 307         IHB=MAX(2,IABS(MINT(104))/110)
 308         IF(IHB.GE.5) IHB=1
 309         IF(ISUB.EQ.91) THEN
 310           BMN=2.*BHAD(IHA)+2.*BHAD(IHB)+4.*SH**EPS-4.2
 311         ELSEIF(ISUB.EQ.92) THEN
 312           BMN=MAX(2.,2.*BHAD(IHB))
 313         ELSEIF(ISUB.EQ.93) THEN
 314           BMN=MAX(2.,2.*BHAD(IHA))
 315         ELSEIF(ISUB.EQ.94) THEN
 316           BMN=2.*ALP*4.
 317         ENDIF
 318
 319 C...Determine maximum possible t range and coefficient of generation.
 320         SQLA12=(SH-SQM1-SQM2)**2-4.*SQM1*SQM2
 321         SQLA34=(SH-SQM3-SQM4)**2-4.*SQM3*SQM4
 322         THA=SH-(SQM1+SQM2+SQM3+SQM4)+(SQM1-SQM2)*(SQM3-SQM4)/SH
 323         THB=SQRT(MAX(0.,SQLA12))*SQRT(MAX(0.,SQLA34))/SH
 324         THC=(SQM3-SQM1)*(SQM4-SQM2)+(SQM1+SQM4-SQM2-SQM3)*
 325      &  (SQM1*SQM4-SQM2*SQM3)/SH
 326         THL=-0.5*(THA+THB)
 327         THU=THC/THL
 328         THRND=EXP(MAX(-50.,BMN*(THL-THU)))-1.
 329
 330 C...Select diffractive mass/masses according to dm^2/m^2.
 331   210   DO 220 JT=1,2
 332         IF(MINT(16+JT).EQ.0) THEN
 333           PDIF(2+JT)=PDIF(JT)
 334         ELSE
 335           PMMIN=PDIF(JT)
 336           PMMAX=MAX(VINT(2+JT),VINT(1)-PDIF(3-JT))
 337           PDIF(2+JT)=PMMIN*(PMMAX/PMMIN)**RLU(0)
 338         ENDIF
 339   220   CONTINUE
 340         SQM3=PDIF(3)**2
 341         SQM4=PDIF(4)**2
 342
 343 C..Additional mass factors, including resonance enhancement.
 344         IF(PDIF(3)+PDIF(4).GE.VINT(1)) GOTO 210
 345         IF(ISUB.EQ.92) THEN
 346           FSD=(1.-SQM3/SH)*(1.+CRES*SMRES1/(SMRES1+SQM3))
 347           IF(FSD.LT.RLU(0)*(1.+CRES)) GOTO 210
 348         ELSEIF(ISUB.EQ.93) THEN
 349           FSD=(1.-SQM4/SH)*(1.+CRES*SMRES2/(SMRES2+SQM4))
 350           IF(FSD.LT.RLU(0)*(1.+CRES)) GOTO 210
 351         ELSEIF(ISUB.EQ.94) THEN
 352           FDD=(1.-(PDIF(3)+PDIF(4))**2/SH)*(SH*SMP/(SH*SMP+SQM3*SQM4))*
 353      &    (1.+CRES*SMRES1/(SMRES1+SQM3))*(1.+CRES*SMRES2/(SMRES2+SQM4))
 354           IF(FDD.LT.RLU(0)*(1.+CRES)**2) GOTO 210
 355         ENDIF
 356
 357 C...Select t according to exp(Bmn*t) and correct to right slope.
 358         TH=THU+LOG(1.+THRND*RLU(0))/BMN
 359         IF(ISUB.GE.92) THEN
 360           IF(ISUB.EQ.92) THEN
 361             BADD=2.*ALP*LOG(SH/SQM3)
 362             IF(BHAD(IHB).LT.1.) BADD=MAX(0.,BADD+2.*BHAD(IHB)-2.)
 363           ELSEIF(ISUB.EQ.93) THEN
 364             BADD=2.*ALP*LOG(SH/SQM4)
 365             IF(BHAD(IHA).LT.1.) BADD=MAX(0.,BADD+2.*BHAD(IHA)-2.)
 366           ELSEIF(ISUB.EQ.94) THEN
 367             BADD=2.*ALP*(LOG(EXP(4.)+SH/(ALP*SQM3*SQM4))-4.)
 368           ENDIF
 369           IF(EXP(MAX(-50.,BADD*(TH-THU))).LT.RLU(0)) GOTO 210
 370         ENDIF
 371
 372 C...Check whether m^2 and t choices are consistent.
 373         SQLA34=(SH-SQM3-SQM4)**2-4.*SQM3*SQM4
 374         THA=SH-(SQM1+SQM2+SQM3+SQM4)+(SQM1-SQM2)*(SQM3-SQM4)/SH
 375         THB=SQRT(MAX(0.,SQLA12))*SQRT(MAX(0.,SQLA34))/SH
 376         IF(THB.LE.1E-8) GOTO 210
 377         THC=(SQM3-SQM1)*(SQM4-SQM2)+(SQM1+SQM4-SQM2-SQM3)*
 378      &  (SQM1*SQM4-SQM2*SQM3)/SH
 379         THLM=-0.5*(THA+THB)
 380         THUM=THC/THLM
 381         IF(TH.LT.THLM.OR.TH.GT.THUM) GOTO 210
 382
 383 C...Information to output.
 384         VINT(21)=1.
 385         VINT(22)=0.
 386         VINT(23)=MIN(1.,MAX(-1.,(THA+2.*TH)/THB))
 387         VINT(45)=TH
 388         VINT(59)=2.*SQRT(MAX(0.,-(THC+THA*TH+TH**2)))/THB
 389         VINT(63)=PDIF(3)**2
 390         VINT(64)=PDIF(4)**2
 391
 392 C...Note: in the following, by In is meant the integral over the
 393 C...quantity multiplying coefficient cn.
 394 C...Choose tau according to h1(tau)/tau, where
 395 C...h1(tau) = c1 + I1/I2*c2*1/tau + I1/I3*c3*1/(tau+tau_R) +
 396 C...I1/I4*c4*tau/((s*tau-m^2)^2+(m*Gamma)^2) +
 397 C...I1/I5*c5*1/(tau+tau_R') +
 398 C...I1/I6*c6*tau/((s*tau-m'^2)^2+(m'*Gamma')^2) +
 399 C...I1/I7*c7*tau/(1.-tau), and
 400 C...c1 + c2 + c3 + c4 + c5 + c6 + c7 = 1.
 401       ELSEIF(ISTSB.GE.1.AND.ISTSB.LE.6) THEN
 402         CALL PYKLIM(1)
 403         IF(MINT(51).NE.0) THEN
 404           IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
 405           IF(MFAIL.EQ.1) THEN
 406             MSTI(61)=1
 407             RETURN
 408           ENDIF
 409           GOTO 100
 410         ENDIF
 411         RTAU=RLU(0)
 412         MTAU=1
 413         IF(RTAU.GT.COEF(ISUB,1)) MTAU=2
 414         IF(RTAU.GT.COEF(ISUB,1)+COEF(ISUB,2)) MTAU=3
 415         IF(RTAU.GT.COEF(ISUB,1)+COEF(ISUB,2)+COEF(ISUB,3)) MTAU=4
 416         IF(RTAU.GT.COEF(ISUB,1)+COEF(ISUB,2)+COEF(ISUB,3)+COEF(ISUB,4))
 417      &  MTAU=5
 418         IF(RTAU.GT.COEF(ISUB,1)+COEF(ISUB,2)+COEF(ISUB,3)+COEF(ISUB,4)+
 419      &  COEF(ISUB,5)) MTAU=6
 420         IF(RTAU.GT.COEF(ISUB,1)+COEF(ISUB,2)+COEF(ISUB,3)+COEF(ISUB,4)+
 421      &  COEF(ISUB,5)+COEF(ISUB,6)) MTAU=7
 422         CALL PYKMAP(1,MTAU,RLU(0))
 423
 424 C...2 -> 3, 4 processes:
 425 C...Choose tau' according to h4(tau,tau')/tau', where
 426 C...h4(tau,tau') = c1 + I1/I2*c2*(1 - tau/tau')^3/tau' +
 427 C...I1/I3*c3*1/(1 - tau'), and c1 + c2 + c3 = 1.
 428         IF(ISTSB.GE.3.AND.ISTSB.LE.5) THEN
 429           CALL PYKLIM(4)
 430           IF(MINT(51).NE.0) THEN
 431             IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
 432             IF(MFAIL.EQ.1) THEN
 433               MSTI(61)=1
 434               RETURN
 435             ENDIF
 436             GOTO 100
 437           ENDIF
 438           RTAUP=RLU(0)
 439           MTAUP=1
 440           IF(RTAUP.GT.COEF(ISUB,18)) MTAUP=2
 441           IF(RTAUP.GT.COEF(ISUB,18)+COEF(ISUB,19)) MTAUP=3
 442           CALL PYKMAP(4,MTAUP,RLU(0))
 443         ENDIF
 444
 445 C...Choose y* according to h2(y*), where
 446 C...h2(y*) = I0/I1*c1*(y*-y*min) + I0/I2*c2*(y*max-y*) +
 447 C...I0/I3*c3*1/cosh(y*) + I0/I4*c4*1/(1-exp(y*-y*max)) +
 448 C...I0/I5*c5*1/(1-exp(-y*-y*min)), I0 = y*max-y*min,
 449 C...and c1 + c2 + c3 + c4 + c5 = 1.
 450         CALL PYKLIM(2)
 451         IF(MINT(51).NE.0) THEN
 452           IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
 453           IF(MFAIL.EQ.1) THEN
 454             MSTI(61)=1
 455             RETURN
 456           ENDIF
 457           GOTO 100
 458         ENDIF
 459         RYST=RLU(0)
 460         MYST=1
 461         IF(RYST.GT.COEF(ISUB,8)) MYST=2
 462         IF(RYST.GT.COEF(ISUB,8)+COEF(ISUB,9)) MYST=3
 463         IF(RYST.GT.COEF(ISUB,8)+COEF(ISUB,9)+COEF(ISUB,10)) MYST=4
 464         IF(RYST.GT.COEF(ISUB,8)+COEF(ISUB,9)+COEF(ISUB,10)+
 465      &  COEF(ISUB,11)) MYST=5
 466         CALL PYKMAP(2,MYST,RLU(0))
 467
 468 C...2 -> 2 processes:
 469 C...Choose cos(theta-hat) (cth) according to h3(cth), where
 470 C...h3(cth) = c0 + I0/I1*c1*1/(A - cth) + I0/I2*c2*1/(A + cth) +
 471 C...I0/I3*c3*1/(A - cth)^2 + I0/I4*c4*1/(A + cth)^2,
 472 C...A = 1 + 2*(m3*m4/sh)^2 (= 1 for massless products),
 473 C...and c0 + c1 + c2 + c3 + c4 = 1.
 474         CALL PYKLIM(3)
 475         IF(MINT(51).NE.0) THEN
 476           IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
 477           IF(MFAIL.EQ.1) THEN
 478             MSTI(61)=1
 479             RETURN
 480           ENDIF
 481           GOTO 100
 482         ENDIF
 483         IF(ISTSB.EQ.2.OR.ISTSB.EQ.4.OR.ISTSB.EQ.6) THEN
 484           RCTH=RLU(0)
 485           MCTH=1
 486           IF(RCTH.GT.COEF(ISUB,13)) MCTH=2
 487           IF(RCTH.GT.COEF(ISUB,13)+COEF(ISUB,14)) MCTH=3
 488           IF(RCTH.GT.COEF(ISUB,13)+COEF(ISUB,14)+COEF(ISUB,15)) MCTH=4
 489           IF(RCTH.GT.COEF(ISUB,13)+COEF(ISUB,14)+COEF(ISUB,15)+
 490      &    COEF(ISUB,16)) MCTH=5
 491           CALL PYKMAP(3,MCTH,RLU(0))
 492         ENDIF
 493
 494 C...2 -> 3 : select pT1, phi1, pT2, phi2, y3 for 3 outgoing.
 495         IF(ISTSB.EQ.5) THEN
 496           CALL PYKMAP(5,0,0.)
 497           IF(MINT(51).NE.0) THEN
 498             IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
 499             IF(MFAIL.EQ.1) THEN
 500               MSTI(61)=1
 501               RETURN
 502             ENDIF
 503             GOTO 100
 504           ENDIF
 505         ENDIF
 506
 507 C...Low-pT or multiple interactions (first semihard interaction).
 508       ELSEIF(ISTSB.EQ.9) THEN
 509         CALL PYMULT(3)
 510         ISUB=MINT(1)
 511
 512 C...Generate user-defined process: kinematics plus weight.
 513       ELSEIF(ISTSB.EQ.11) THEN
 514         MSTI(51)=0
 515         CALL PYUPEV(ISUB,SIGS)
 516         IF(NUP.LE.0) THEN
 517           MINT(51)=2
 518           MSTI(51)=1
 519           IF(MINT(82).EQ.1) THEN
 520             NGEN(0,1)=NGEN(0,1)-1
 521             NGEN(0,2)=NGEN(0,2)-1
 522             NGEN(ISUB,1)=NGEN(ISUB,1)-1
 523           ENDIF
 524           IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
 525           RETURN
 526         ENDIF
 527
 528 C...Construct 'trivial' kinematical variables needed.
 529         KFL1=KUP(1,2)
 530         KFL2=KUP(2,2)
 531         VINT(41)=2.*PUP(1,4)/VINT(1)
 532         VINT(42)=2.*PUP(2,4)/VINT(1)
 533         VINT(21)=VINT(41)*VINT(42)
 534         VINT(22)=0.5*LOG(VINT(41)/VINT(42))
 535         VINT(44)=VINT(21)*VINT(2)
 536         VINT(43)=SQRT(MAX(0.,VINT(44)))
 537         VINT(56)=Q2UP(0)
 538         VINT(55)=SQRT(MAX(0.,VINT(56)))
 539
 540 C...Construct other kinematical variables needed (approximately).
 541         VINT(23)=0.
 542         VINT(26)=VINT(21)
 543         VINT(45)=-0.5*VINT(44)
 544         VINT(46)=-0.5*VINT(44)
 545         VINT(49)=VINT(43)
 546         VINT(50)=VINT(44)
 547         VINT(51)=VINT(55)
 548         VINT(52)=VINT(56)
 549         VINT(53)=VINT(55)
 550         VINT(54)=VINT(56)
 551         VINT(25)=0.
 552         VINT(48)=0.
 553         DO 230 IUP=3,NUP
 554         IF(KUP(IUP,1).EQ.1) VINT(25)=VINT(25)+2.*(PUP(IUP,5)**2+
 555      &  PUP(IUP,1)**2+PUP(IUP,2)**2)/VINT(1)
 556         IF(KUP(IUP,1).EQ.1) VINT(48)=VINT(48)+0.5*(PUP(IUP,1)**2+
 557      &  PUP(IUP,2)**2)
 558   230   CONTINUE
 559         VINT(47)=SQRT(VINT(48))
 560
 561 C...Calculate structure function weights.
 562         IF(MINT(47).GE.2) THEN
 563           DO 250 I=3-MIN(2,MINT(45)),MIN(2,MINT(46))
 564           MINT(105)=MINT(102+I)
 565           MINT(109)=MINT(106+I)
 566           IF(MSTP(57).LE.1) THEN
 567             CALL PYSTFU(MINT(10+I),VINT(40+I),Q2UP(0),XPQ)
 568           ELSE
 569             CALL PYSTFL(MINT(10+I),VINT(40+I),Q2UP(0),XPQ)
 570           ENDIF
 571           DO 240 KFL=-25,25
 572           XSFX(I,KFL)=XPQ(KFL)
 573   240     CONTINUE
 574   250     CONTINUE
 575         ENDIF
 576       ENDIF
 577
 578 C...Choose azimuthal angle.
 579       VINT(24)=PARU(2)*RLU(0)
 580
 581 C...Check against user cuts on kinematics at parton level.
 582       MINT(51)=0
 583       IF((ISUB.LE.90.OR.ISUB.GT.100).AND.ISTSB.LE.10) CALL PYKLIM(0)
 584       IF(MINT(51).NE.0) THEN
 585           IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
 586           IF(MFAIL.EQ.1) THEN
 587             MSTI(61)=1
 588             RETURN
 589           ENDIF
 590           GOTO 100
 591         ENDIF
 592       IF(MINT(82).EQ.1.AND.MSTP(141).GE.1.AND.ISTSB.LE.10) THEN
 593         MCUT=0
 594         IF(MSUB(91)+MSUB(92)+MSUB(93)+MSUB(94)+MSUB(95).EQ.0)
 595      &  CALL PYKCUT(MCUT)
 596         IF(MCUT.NE.0) THEN
 597           IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
 598           IF(MFAIL.EQ.1) THEN
 599             MSTI(61)=1
 600             RETURN
 601           ENDIF
 602           GOTO 100
 603         ENDIF
 604       ENDIF
 605
 606 C...Calculate differential cross-section for different subprocesses.
 607       IF(ISTSB.LE.10) CALL PYSIGH(NCHN,SIGS)
 608       SIGSOR=SIGS
 609       SIGLPT=SIGT(0,0,5)
 610
 611 C...Multiply cross-section by user-defined weights.
 612       IF(MSTP(173).EQ.1) THEN
 613         SIGS=PARP(173)*SIGS
 614         DO 260 ICHN=1,NCHN
 615         SIGH(ICHN)=PARP(173)*SIGH(ICHN)
 616   260   CONTINUE
 617         SIGLPT=PARP(173)*SIGLPT
 618       ENDIF
 619       WTXS=1.
 620       SIGSWT=SIGS
 621       VINT(99)=1.
 622       VINT(100)=1.
 623       IF(MINT(82).EQ.1.AND.MSTP(142).GE.1) THEN
 624         IF(ISUB.NE.96.AND.MSUB(91)+MSUB(92)+MSUB(93)+MSUB(94)+
 625      &  MSUB(95).EQ.0) CALL PYEVWT(WTXS)
 626         SIGSWT=WTXS*SIGS
 627         VINT(99)=WTXS
 628         IF(MSTP(142).EQ.1) VINT(100)=1./WTXS
 629       ENDIF
 630
 631 C...Calculations for Monte Carlo estimate of all cross-sections.
 632       IF(MINT(82).EQ.1.AND.ISUB.LE.90.OR.ISUB.GE.96) THEN
 633         IF(MSTP(142).LE.1) THEN
 634           XSEC(ISUB,2)=XSEC(ISUB,2)+SIGS
 635           DXSEC(ISUB)=DXSEC(ISUB)+SIGS
 636         ELSE
 637           XSEC(ISUB,2)=XSEC(ISUB,2)+SIGSWT
 638           DXSEC(ISUB)=DXSEC(ISUB)+SIGSWT
 639         ENDIF
 640       ELSEIF(MINT(82).EQ.1) THEN
 641         XSEC(ISUB,2)=XSEC(ISUB,2)+SIGS
 642         DXSEC(ISUB)=DXSEC(ISUB)+SIGS
 643       ENDIF
 644       IF((ISUB.EQ.95.OR.ISUB.EQ.96).AND.LOOP.EQ.1.AND.MINT(82).EQ.1)
 645      &THEN
 646         XSEC(97,2)=XSEC(97,2)+SIGLPT
 647         DXSEC(97)=DXSEC(97)+SIGLPT
 648       ENDIF
 649
 650 C...Multiple interactions: store results of cross-section calculation.
 651       IF(MINT(50).EQ.1.AND.MSTP(82).GE.3) THEN
 652         VINT(153)=SIGSOR
 653         CALL PYMULT(4)
 654       ENDIF
 655
 656 C...Check that weight not negative.
 657       VIOL=SIGSWT/XSEC(ISUB,1)
 658       IF(ISUB.EQ.96.AND.MSTP(173).EQ.1) VIOL=VIOL/PARP(174)
 659       IF(MSTP(123).LE.0) THEN
 660         IF(VIOL.LT.-1E-3) THEN
 661           WRITE(MSTU(11),5000) VIOL,NGEN(0,3)+1
 662           WRITE(MSTU(11),5100) ISUB,VINT(21),VINT(22),VINT(23),VINT(26)
 663           STOP
 664         ENDIF
 665       ELSE
 666         IF(VIOL.LT.MIN(-1E-3,VINT(109))) THEN
 667           VINT(109)=VIOL
 668           WRITE(MSTU(11),5200) VIOL,NGEN(0,3)+1
 669           WRITE(MSTU(11),5100) ISUB,VINT(21),VINT(22),VINT(23),VINT(26)
 670         ENDIF
 671       ENDIF
 672
 673 C...Weighting using estimate of maximum of differential cross-section.
 674       IF(MFAIL.EQ.0) THEN
 675         IF(VIOL.LT.RLU(0)) THEN
 676           IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
 677           GOTO 100
 678         ENDIF
 679       ELSEIF(ISUB.NE.95.AND.ISUB.NE.96) THEN
 680         IF(VIOL.LT.RLU(0)) THEN
 681           MSTI(61)=1
 682           IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
 683           RETURN
 684         ENDIF
 685       ELSE
 686         RATND=SIGLPT/XSEC(95,1)
 687         IF(LOOP.EQ.1.AND.RATND.LT.RLU(0)) THEN
 688           MSTI(61)=1
 689           IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
 690           RETURN
 691         ENDIF
 692         VIOL=VIOL/RATND
 693         IF(VIOL.LT.RLU(0)) THEN
 694           IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
 695           GOTO 100
 696         ENDIF
 697       ENDIF
 698
 699 C...Check for possible violation of estimated maximum of differential
 700 C...cross-section used in weighting.
 701       IF(MSTP(123).LE.0) THEN
 702         IF(VIOL.GT.1.) THEN
 703           WRITE(MSTU(11),5300) VIOL,NGEN(0,3)+1
 704           WRITE(MSTU(11),5100) ISUB,VINT(21),VINT(22),VINT(23),VINT(26)
 705           STOP
 706         ENDIF
 707       ELSEIF(MSTP(123).EQ.1) THEN
 708         IF(VIOL.GT.VINT(108)) THEN
 709           VINT(108)=VIOL
 710           IF(VIOL.GT.1.) THEN
 711             MINT(10)=1
 712             WRITE(MSTU(11),5400) VIOL,NGEN(0,3)+1
 713             WRITE(MSTU(11),5100) ISUB,VINT(21),VINT(22),VINT(23),
 714      &      VINT(26)
 715           ENDIF
 716         ENDIF
 717       ELSEIF(VIOL.GT.VINT(108)) THEN
 718         VINT(108)=VIOL
 719         IF(VIOL.GT.1.) THEN
 720           MINT(10)=1
 721           XDIF=XSEC(ISUB,1)*(VIOL-1.)
 722           XSEC(ISUB,1)=XSEC(ISUB,1)+XDIF
 723           IF(MSUB(ISUB).EQ.1.AND.(ISUB.LE.90.OR.ISUB.GT.96))
 724      &    XSEC(0,1)=XSEC(0,1)+XDIF
 725           WRITE(MSTU(11),5400) VIOL,NGEN(0,3)+1
 726           WRITE(MSTU(11),5100) ISUB,VINT(21),VINT(22),VINT(23),VINT(26)
 727           IF(ISUB.LE.9) THEN
 728             WRITE(MSTU(11),5500) ISUB,XSEC(ISUB,1)
 729           ELSEIF(ISUB.LE.99) THEN
 730             WRITE(MSTU(11),5600) ISUB,XSEC(ISUB,1)
 731           ELSE
 732             WRITE(MSTU(11),5700) ISUB,XSEC(ISUB,1)
 733           ENDIF
 734           VINT(108)=1.
 735         ENDIF
 736       ENDIF
 737
 738 C...Multiple interactions: choose impact parameter.
 739       VINT(148)=1.
 740       IF(MINT(50).EQ.1.AND.(ISUB.LE.90.OR.ISUB.GE.96).AND.MSTP(82).GE.3)
 741      &THEN
 742         CALL PYMULT(5)
 743         IF(VINT(150).LT.RLU(0)) THEN
 744           IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
 745           IF(MFAIL.EQ.1) THEN
 746             MSTI(61)=1
 747             RETURN
 748           ENDIF
 749           GOTO 100
 750         ENDIF
 751       ENDIF
 752       IF(MINT(82).EQ.1) NGEN(0,2)=NGEN(0,2)+1
 753       IF(MINT(82).EQ.1.AND.MSUB(95).EQ.1) THEN
 754         IF(ISUB.LE.90.OR.ISUB.GE.95) NGEN(95,1)=NGEN(95,1)+1
 755         IF(ISUB.LE.90.OR.ISUB.GE.96) NGEN(96,2)=NGEN(96,2)+1
 756       ENDIF
 757       IF(ISUB.LE.90.OR.ISUB.GE.96) MINT(31)=MINT(31)+1
 758
 759 C...Choose flavour of reacting partons (and subprocess).
 760       IF(ISTSB.GE.11) GOTO 280
 761       RSIGS=SIGS*RLU(0)
 762       QT2=VINT(48)
 763       RQQBAR=PARP(87)*(1.-(QT2/(QT2+(PARP(88)*PARP(82))**2))**2)
 764       IF(ISUB.NE.95.AND.(ISUB.NE.96.OR.MSTP(82).LE.1.OR.
 765      &RLU(0).GT.RQQBAR)) THEN
 766         DO 270 ICHN=1,NCHN
 767         KFL1=ISIG(ICHN,1)
 768         KFL2=ISIG(ICHN,2)
 769         MINT(2)=ISIG(ICHN,3)
 770         RSIGS=RSIGS-SIGH(ICHN)
 771         IF(RSIGS.LE.0.) GOTO 280
 772   270   CONTINUE
 773
 774 C...Multiple interactions: choose qq~ preferentially at small pT.
 775       ELSEIF(ISUB.EQ.96) THEN
 776         MINT(105)=MINT(103)
 777         MINT(109)=MINT(107)
 778         CALL PYSPLI(MINT(11),21,KFL1,KFLDUM)
 779         MINT(105)=MINT(104)
 780         MINT(109)=MINT(108)
 781         CALL PYSPLI(MINT(12),21,KFL2,KFLDUM)
 782         MINT(1)=11
 783         MINT(2)=1
 784         IF(KFL1.EQ.KFL2.AND.RLU(0).LT.0.5) MINT(2)=2
 785
 786 C...Low-pT: choose string drawing configuration.
 787       ELSE
 788         KFL1=21
 789         KFL2=21
 790         RSIGS=6.*RLU(0)
 791         MINT(2)=1
 792         IF(RSIGS.GT.1.) MINT(2)=2
 793         IF(RSIGS.GT.2.) MINT(2)=3
 794       ENDIF
 795
 796 C...Reassign QCD process. Partons before initial state radiation.
 797   280 IF(MINT(2).GT.10) THEN
 798         MINT(1)=MINT(2)/10
 799         MINT(2)=MOD(MINT(2),10)
 800       ENDIF
 801       IF(MINT(82).EQ.1.AND.MSTP(111).GE.0) NGEN(MINT(1),2)=
 802      &NGEN(MINT(1),2)+1
 803       MINT(15)=KFL1
 804       MINT(16)=KFL2
 805       MINT(13)=MINT(15)
 806       MINT(14)=MINT(16)
 807       VINT(141)=VINT(41)
 808       VINT(142)=VINT(42)
 809       VINT(151)=0.
 810       VINT(152)=0.
 811
 812 C...Calculate x value of photon for parton inside photon inside e.
 813       DO 310 JT=1,2
 814       MINT(18+JT)=0
 815       VINT(154+JT)=0.
 816       MSPLI=0
 817       IF(JT.EQ.1.AND.MINT(43).LE.2) MSPLI=1
 818       IF(JT.EQ.2.AND.MOD(MINT(43),2).EQ.1) MSPLI=1
 819       IF(IABS(MINT(14+JT)).LE.8.OR.MINT(14+JT).EQ.21) MSPLI=MSPLI+1
 820       IF(MSPLI.EQ.2) THEN
 821         KFLH=MINT(14+JT)
 822         XHRD=VINT(140+JT)
 823         Q2HRD=VINT(54)
 824         MINT(105)=MINT(102+JT)
 825         MINT(109)=MINT(106+JT)
 826         IF(MSTP(57).LE.1) THEN
 827           CALL PYSTFU(22,XHRD,Q2HRD,XPQ)
 828         ELSE
 829           CALL PYSTFL(22,XHRD,Q2HRD,XPQ)
 830         ENDIF
 831         WTMX=4.*XPQ(KFLH)
 832         IF(MSTP(13).EQ.2) THEN
 833           Q2PMS=Q2HRD/PMAS(11,1)**2
 834           WTMX=WTMX*LOG(MAX(2.,Q2PMS*(1.-XHRD)/XHRD**2))
 835         ENDIF
 836   290   XE=XHRD**RLU(0)
 837         XG=MIN(0.999999,XHRD/XE)
 838         IF(MSTP(57).LE.1) THEN
 839           CALL PYSTFU(22,XG,Q2HRD,XPQ)
 840         ELSE
 841           CALL PYSTFL(22,XG,Q2HRD,XPQ)
 842         ENDIF
 843         WT=(1.+(1.-XE)**2)*XPQ(KFLH)
 844         IF(MSTP(13).EQ.2) WT=WT*LOG(MAX(2.,Q2PMS*(1.-XE)/XE**2))
 845         IF(WT.LT.RLU(0)*WTMX) GOTO 290
 846         MINT(18+JT)=1
 847         VINT(154+JT)=XE
 848         DO 300 KFLS=-25,25
 849         XSFX(JT,KFLS)=XPQ(KFLS)
 850   300   CONTINUE
 851       ENDIF
 852   310 CONTINUE
 853
 854 C...Pick scale where photon is resolved.
 855       IF(MINT(107).EQ.3) VINT(283)=PARP(15)**2*
 856      &(VINT(54)/PARP(15)**2)**RLU(0)
 857       IF(MINT(108).EQ.3) VINT(284)=PARP(15)**2*
 858      &(VINT(54)/PARP(15)**2)**RLU(0)
 859       IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
 860
 861 C...Format statements for differential cross-section maximum violations.
 862  5000 FORMAT(1X,'Error: negative cross-section fraction',1P,E11.3,1X,
 863      &'in event',1X,I7,'.'/1X,'Execution stopped!')
 864  5100 FORMAT(1X,'ISUB = ',I3,'; Point of violation:'/1X,'tau =',1P,
 865      &E11.3,', y* =',E11.3,', cthe = ',0P,F11.7,', tau'' =',1P,E11.3)
 866  5200 FORMAT(1X,'Warning: negative cross-section fraction',1P,E11.3,1X,
 867      &'in event',1X,I7)
 868  5300 FORMAT(1X,'Error: maximum violated by',1P,E11.3,1X,
 869      &'in event',1X,I7,'.'/1X,'Execution stopped!')
 870  5400 FORMAT(1X,'Warning: maximum violated by',1P,E11.3,1X,
 871      &'in event',1X,I7)
 872  5500 FORMAT(1X,'XSEC(',I1,',1) increased to',1P,E11.3)
 873  5600 FORMAT(1X,'XSEC(',I2,',1) increased to',1P,E11.3)
 874  5700 FORMAT(1X,'XSEC(',I3,',1) increased to',1P,E11.3)
 875
 876       RETURN
 877       END