      SUBROUTINE PYRAND
C...Double precision and integer declarations.
      IMPLICIT DOUBLE PRECISION(A-H, O-Z)
      IMPLICIT INTEGER(I-N)
      INTEGER PYK,PYCHGE,PYCOMP
C...Parameter statement to help give large particle numbers.
      PARAMETER (KSUSY1=1000000,KSUSY2=2000000,KTECHN=3000000,
     &KEXCIT=4000000,KDIMEN=5000000)
C...User process initialization and event commonblocks.
      INTEGER MAXPUP
      PARAMETER (MAXPUP=100)
      INTEGER IDBMUP,PDFGUP,PDFSUP,IDWTUP,NPRUP,LPRUP
      DOUBLE PRECISION EBMUP,XSECUP,XERRUP,XMAXUP
      COMMON/HEPRUP/IDBMUP(2),EBMUP(2),PDFGUP(2),PDFSUP(2),
     &IDWTUP,NPRUP,XSECUP(MAXPUP),XERRUP(MAXPUP),XMAXUP(MAXPUP),
     &LPRUP(MAXPUP)
      INTEGER MAXNUP
      PARAMETER (MAXNUP=500)
      INTEGER NUP,IDPRUP,IDUP,ISTUP,MOTHUP,ICOLUP
      DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,PUP,VTIMUP,SPINUP
      COMMON/HEPEUP/NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,IDUP(MAXNUP),
     &ISTUP(MAXNUP),MOTHUP(2,MAXNUP),ICOLUP(2,MAXNUP),PUP(5,MAXNUP),
     &VTIMUP(MAXNUP),SPINUP(MAXNUP)
      SAVE /HEPRUP/,/HEPEUP/
C...Commonblocks.
      COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
      COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
      COMMON/PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
      COMMON/PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200)
      COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
      COMMON/PYINT1/MINT(400),VINT(400)
      COMMON/PYINT2/ISET(500),KFPR(500,2),COEF(500,20),ICOL(40,4,2)
      COMMON/PYINT3/XSFX(2,-40:40),ISIG(1000,3),SIGH(1000)
      COMMON/PYINT4/MWID(500),WIDS(500,5)
      COMMON/PYINT5/NGENPD,NGEN(0:500,3),XSEC(0:500,3)
      COMMON/PYINT7/SIGT(0:6,0:6,0:5)
      COMMON/PYMSSM/IMSS(0:99),RMSS(0:99)
      SAVE /PYDAT1/,/PYDAT2/,/PYDAT3/,/PYSUBS/,/PYPARS/,/PYINT1/,
     &/PYINT2/,/PYINT3/,/PYINT4/,/PYINT5/,/PYINT7/,/PYMSSM/
C...Local arrays.
      DIMENSION XPQ(-25:25),PMM(2),PDIF(4),BHAD(4),PMMN(2)
C...Parameters and data used in elastic/diffractive treatment.
      DATA EPS/0.0808D0/, ALP/0.25D0/, CRES/2D0/, PMRC/1.062D0/,
     &SMP/0.880D0/, BHAD/2.3D0,1.4D0,1.4D0,0.23D0/
C...Initial values, specifically for (first) semihard interaction.
      MINT(10)=0
      MINT(17)=0
      MINT(18)=0
      VINT(97)=1D0
      VINT(143)=1D0
      VINT(144)=1D0
      VINT(157)=0D0
      VINT(158)=0D0
      MFAIL=0
      IF(MSTP(171).EQ.1.AND.MSTP(172).EQ.2) MFAIL=1
      ISUB=0
      ISTSB=0
      LOOP=0
  100 LOOP=LOOP+1
      MINT(51)=0
      MINT(143)=1
C...Start by assuming incoming photon is entering subprocess.
      IF(MINT(11).EQ.22) THEN
         MINT(15)=22
         VINT(307)=VINT(3)**2
      ENDIF
      IF(MINT(12).EQ.22) THEN
         MINT(16)=22
         VINT(308)=VINT(4)**2
      ENDIF
      MINT(103)=MINT(11)
      MINT(104)=MINT(12)
C...Choice of process type - first event of pileup.
      INMULT=0
      IF(MINT(82).EQ.1.AND.ISUB.GE.91.AND.ISUB.LE.96) THEN
      ELSEIF(MINT(82).EQ.1) THEN
C...For gamma-p or gamma-gamma first pick between alternatives.
        IGA=0
        IF(MINT(121).GT.1) CALL PYSAVE(4,IGA)
        MINT(122)=IGA
C...For real gamma + gamma with different nature, flip at random.
        IF(MINT(11).EQ.22.AND.MINT(12).EQ.22.AND.MINT(123).GE.4.AND.
     &  MSTP(14).LE.10.AND.PYR(0).GT.0.5D0) THEN
          MINTSV=MINT(41)
          MINT(41)=MINT(42)
          MINT(42)=MINTSV
          MINTSV=MINT(45)
          MINT(45)=MINT(46)
          MINT(46)=MINTSV
          MINTSV=MINT(107)
          MINT(107)=MINT(108)
          MINT(108)=MINTSV
          IF(MINT(47).EQ.2.OR.MINT(47).EQ.3) MINT(47)=5-MINT(47)
        ENDIF
C...Pick process type, possibly by user process machinery.
C...(If the latter, also event will be picked here.)
        IF(MINT(111).EQ.11.AND.IABS(IDWTUP).EQ.2.AND.LOOP.GE.2) THEN
          CALL UPEVNT
        ELSEIF(MINT(111).EQ.11.AND.IABS(IDWTUP).GE.3) THEN
          CALL UPEVNT
          ISUB=0
  110     ISUB=ISUB+1
          IF((ISET(ISUB).NE.11.OR.KFPR(ISUB,2).NE.IDPRUP).AND.
     &    ISUB.LT.500) GOTO 110
        ELSE
          RSUB=XSEC(0,1)*PYR(0)
          DO 120 I=1,500
            IF(MSUB(I).NE.1) GOTO 120
            ISUB=I
            RSUB=RSUB-XSEC(I,1)
            IF(RSUB.LE.0D0) GOTO 130
  120     CONTINUE
  130     IF(ISUB.EQ.95) ISUB=96
          IF(ISUB.EQ.96) INMULT=1
          IF(ISET(ISUB).EQ.11) THEN
            IDPRUP=KFPR(ISUB,2)
            CALL UPEVNT
          ENDIF
        ENDIF
C...Choice of inclusive process type - pileup events.
      ELSEIF(MINT(82).GE.2.AND.ISUB.EQ.0) THEN
        RSUB=VINT(131)*PYR(0)
        ISUB=96
        IF(RSUB.GT.SIGT(0,0,5)) ISUB=94
        IF(RSUB.GT.SIGT(0,0,5)+SIGT(0,0,4)) ISUB=93
        IF(RSUB.GT.SIGT(0,0,5)+SIGT(0,0,4)+SIGT(0,0,3)) ISUB=92
        IF(RSUB.GT.SIGT(0,0,5)+SIGT(0,0,4)+SIGT(0,0,3)+SIGT(0,0,2))
     &  ISUB=91
        IF(ISUB.EQ.96) INMULT=1
      ENDIF
C...Choice of photon energy and flux factor inside lepton.
      IF(MINT(141).NE.0.OR.MINT(142).NE.0) THEN
        CALL PYGAGA(3,WTGAGA)
        IF(ISUB.GE.131.AND.ISUB.LE.140) THEN
          CKIN(3)=MAX(VINT(285),VINT(154))
          CKIN(1)=2D0*CKIN(3)
        ENDIF
C...When necessary set direct/resolved photon by hand.
      ELSEIF(MINT(15).EQ.22.OR.MINT(16).EQ.22) THEN
        IF(MINT(15).EQ.22.AND.MINT(41).EQ.2) MINT(15)=0
        IF(MINT(16).EQ.22.AND.MINT(42).EQ.2) MINT(16)=0
      ENDIF
C...Restrict direct*resolved processes to pTmin >= Q,
C...to avoid doublecounting  with DIS.
      IF(MSTP(18).EQ.3.AND.ISUB.GE.131.AND.ISUB.LE.136) THEN
        IF(MINT(15).EQ.22) THEN
          CKIN(3)=MAX(VINT(285),VINT(154),ABS(VINT(3)))
        ELSE
          CKIN(3)=MAX(VINT(285),VINT(154),ABS(VINT(4)))
        ENDIF
        CKIN(1)=2D0*CKIN(3)
      ENDIF
C...Set up for multiple interactions.
      IF(INMULT.EQ.1) CALL PYMULT(2)
C...Loopback point for minimum bias in photon physics.
      LOOP2=0
  140 LOOP2=LOOP2+1
      IF(MINT(82).EQ.1) NGEN(0,1)=NGEN(0,1)+MINT(143)
      IF(MINT(82).EQ.1) NGEN(ISUB,1)=NGEN(ISUB,1)+MINT(143)
      IF(ISUB.EQ.96.AND.LOOP2.EQ.1.AND.MINT(82).EQ.1)
     &NGEN(97,1)=NGEN(97,1)+MINT(143)
      MINT(1)=ISUB
      ISTSB=ISET(ISUB)
C...Random choice of flavour for some SUSY processes.
      IF(ISUB.GE.201.AND.ISUB.LE.301) THEN
C...~e_L ~nu_e or ~mu_L ~nu_mu.
        IF(ISUB.EQ.210) THEN
          KFPR(ISUB,1)=KSUSY1+11+2*INT(0.5D0+PYR(0))
          KFPR(ISUB,2)=KFPR(ISUB,1)+1
C...~nu_e ~nu_e(bar) or ~nu_mu ~nu_mu(bar).
        ELSEIF(ISUB.EQ.213) THEN
          KFPR(ISUB,1)=KSUSY1+12+2*INT(0.5D0+PYR(0))
          KFPR(ISUB,2)=KFPR(ISUB,1)
C...~q ~chi/~g; ~q = ~d, ~u, ~s, ~c or ~b.
        ELSEIF(ISUB.GE.246.AND.ISUB.LE.259) THEN
          IF(ISUB.GE.258) THEN
            RKF=4D0
          ELSE
            RKF=5D0
          ENDIF
          IF(MOD(ISUB,2).EQ.0) THEN
            KFPR(ISUB,1)=KSUSY1+1+INT(RKF*PYR(0))
          ELSE
            KFPR(ISUB,1)=KSUSY2+1+INT(RKF*PYR(0))
          ENDIF
C...~q1 ~q2; ~q = ~d, ~u, ~s, or ~c.
        ELSEIF(ISUB.GE.271.AND.ISUB.LE.276) THEN
          IF(ISUB.EQ.271.OR.ISUB.EQ.274) THEN
            KSU1=KSUSY1
            KSU2=KSUSY1
          ELSEIF(ISUB.EQ.272.OR.ISUB.EQ.275) THEN
            KSU1=KSUSY2
            KSU2=KSUSY2
          ELSEIF(PYR(0).LT.0.5D0) THEN
            KSU1=KSUSY1
            KSU2=KSUSY2
          ELSE
            KSU1=KSUSY2
            KSU2=KSUSY1
          ENDIF
          KFPR(ISUB,1)=KSU1+1+INT(4D0*PYR(0))
          KFPR(ISUB,2)=KSU2+1+INT(4D0*PYR(0))
C...~q ~q(bar);  ~q = ~d, ~u, ~s, or ~c.
        ELSEIF(ISUB.EQ.277.OR.ISUB.EQ.279) THEN
          KFPR(ISUB,1)=KSUSY1+1+INT(4D0*PYR(0))
          KFPR(ISUB,2)=KFPR(ISUB,1)
        ELSEIF(ISUB.EQ.278.OR.ISUB.EQ.280) THEN
          KFPR(ISUB,1)=KSUSY2+1+INT(4D0*PYR(0))
          KFPR(ISUB,2)=KFPR(ISUB,1)
C...~q1 ~q2; ~q = ~d, ~u, ~s, or ~c.
        ELSEIF(ISUB.GE.281.AND.ISUB.LE.286) THEN
          IF(ISUB.EQ.281.OR.ISUB.EQ.284) THEN
            KSU1=KSUSY1
            KSU2=KSUSY1
          ELSEIF(ISUB.EQ.282.OR.ISUB.EQ.285) THEN
            KSU1=KSUSY2
            KSU2=KSUSY2
          ELSEIF(PYR(0).LT.0.5D0) THEN
            KSU1=KSUSY1
            KSU2=KSUSY2
          ELSE
            KSU1=KSUSY2
            KSU2=KSUSY1
          ENDIF
          IF(ISUB.EQ.281.OR.ISUB.LE.283) THEN
            RKF=5D0
          ELSE
            RKF=4D0
          ENDIF
          KFPR(ISUB,2)=KSU2+1+INT(RKF*PYR(0))
        ENDIF
      ENDIF
C...Find resonances (explicit or implicit in cross-section).
      MINT(72)=0
      KFR1=0
      IF(ISTSB.EQ.1.OR.ISTSB.EQ.3.OR.ISTSB.EQ.5) THEN
        KFR1=KFPR(ISUB,1)
      ELSEIF(ISUB.EQ.24.OR.ISUB.EQ.25.OR.ISUB.EQ.110.OR.ISUB.EQ.165.OR.
     &  ISUB.EQ.171.OR.ISUB.EQ.176) THEN
        KFR1=23
      ELSEIF(ISUB.EQ.23.OR.ISUB.EQ.26.OR.ISUB.EQ.166.OR.ISUB.EQ.172.OR.
     &  ISUB.EQ.177) THEN
        KFR1=24
      ELSEIF(ISUB.GE.71.AND.ISUB.LE.77) THEN
        KFR1=25
        IF(MSTP(46).EQ.5) THEN
          KFR1=89
          PMAS(89,1)=PARP(45)
          PMAS(89,2)=PARP(45)**3/(96D0*PARU(1)*PARP(47)**2)
        ENDIF
      ELSEIF(ISUB.EQ.194) THEN
        KFR1=KTECHN+113
      ELSEIF(ISUB.EQ.195) THEN
        KFR1=KTECHN+213
      ELSEIF(ISUB.GE.361.AND.ISUB.LE.368) THEN
        KFR1=KTECHN+113
      ELSEIF(ISUB.GE.370.AND.ISUB.LE.377) THEN
        KFR1=KTECHN+213
      ENDIF
      CKMX=CKIN(2)
      IF(CKMX.LE.0D0) CKMX=VINT(1)
      KCR1=PYCOMP(KFR1)
      IF(KFR1.NE.0) THEN
        IF(CKIN(1).GT.PMAS(KCR1,1)+20D0*PMAS(KCR1,2).OR.
     &  CKMX.LT.PMAS(KCR1,1)-20D0*PMAS(KCR1,2)) KFR1=0
      ENDIF
      IF(KFR1.NE.0) THEN
        TAUR1=PMAS(KCR1,1)**2/VINT(2)
        IF(KFR1.EQ.KTECHN+113) THEN
          CALL PYTECM(S1,S2)
          TAUR1=S1/VINT(2)
        ENDIF
        GAMR1=PMAS(KCR1,1)*PMAS(KCR1,2)/VINT(2)
        MINT(72)=1
        MINT(73)=KFR1
        VINT(73)=TAUR1
        VINT(74)=GAMR1
      ENDIF
      IF(ISUB.EQ.141.OR.ISUB.EQ.194.OR.(ISUB.GE.364.AND.ISUB.LE.368))
     $THEN
        KFR2=23
        IF(ISUB.EQ.194) THEN
          KFR2=KTECHN+223
        ELSEIF(ISUB.GE.364.AND.ISUB.LE.368) THEN
          KFR2=KTECHN+223
        ENDIF
        KCR2=PYCOMP(KFR2)
        TAUR2=PMAS(KCR2,1)**2/VINT(2)
        IF(KFR2.EQ.KTECHN+223) THEN
          CALL PYTECM(S1,S2)
          TAUR2=S2/VINT(2)
        ENDIF
        GAMR2=PMAS(KCR2,1)*PMAS(KCR2,2)/VINT(2)
        IF(CKIN(1).GT.PMAS(KCR2,1)+20D0*PMAS(KCR2,2).OR.
     &  CKMX.LT.PMAS(KCR2,1)-20D0*PMAS(KCR2,2)) KFR2=0
        IF(KFR2.NE.0.AND.KFR1.NE.0) THEN
          MINT(72)=2
          MINT(74)=KFR2
          VINT(75)=TAUR2
          VINT(76)=GAMR2
        ELSEIF(KFR2.NE.0) THEN
          KFR1=KFR2
          TAUR1=TAUR2
          GAMR1=GAMR2
          MINT(72)=1
          MINT(73)=KFR1
          VINT(73)=TAUR1
          VINT(74)=GAMR1
        ENDIF
      ENDIF
C...Find product masses and minimum pT of process,
C...optionally with broadening according to a truncated Breit-Wigner.
      VINT(63)=0D0
      VINT(64)=0D0
      MINT(71)=0
      VINT(71)=CKIN(3)
      IF(MINT(82).GE.2) VINT(71)=0D0
      VINT(80)=1D0
      IF(ISTSB.EQ.2.OR.ISTSB.EQ.4) THEN
        NBW=0
        DO 160 I=1,2
          PMMN(I)=0D0
          IF(KFPR(ISUB,I).EQ.0) THEN
          ELSEIF(MSTP(42).LE.0.OR.PMAS(PYCOMP(KFPR(ISUB,I)),2).LT.
     &      PARP(41)) THEN
            VINT(62+I)=PMAS(PYCOMP(KFPR(ISUB,I)),1)**2
          ELSE
            NBW=NBW+1
C...This prevents SUSY/t particles from becoming too light.
            KFLW=KFPR(ISUB,I)
            IF(KFLW/KSUSY1.EQ.1.OR.KFLW/KSUSY1.EQ.2) THEN
              KCW=PYCOMP(KFLW)
              PMMN(I)=PMAS(KCW,1)
              DO 150 IDC=MDCY(KCW,2),MDCY(KCW,2)+MDCY(KCW,3)-1
                IF(MDME(IDC,1).GT.0.AND.BRAT(IDC).GT.1E-4) THEN
                  PMSUM=PMAS(PYCOMP(KFDP(IDC,1)),1)+
     &            PMAS(PYCOMP(KFDP(IDC,2)),1)
                  IF(KFDP(IDC,3).NE.0) PMSUM=PMSUM+
     &            PMAS(PYCOMP(KFDP(IDC,3)),1)
                  PMMN(I)=MIN(PMMN(I),PMSUM)
                ENDIF
  150         CONTINUE
            ELSEIF(KFLW.EQ.6) THEN
              PMMN(I)=PMAS(24,1)+PMAS(5,1)
            ENDIF
          ENDIF
  160   CONTINUE
        IF(NBW.GE.1) THEN
          CKIN41=CKIN(41)
          CKIN43=CKIN(43)
          CKIN(41)=MAX(PMMN(1),CKIN(41))
          CKIN(43)=MAX(PMMN(2),CKIN(43))
          CALL PYOFSH(4,0,KFPR(ISUB,1),KFPR(ISUB,2),0D0,PQM3,PQM4)
          CKIN(41)=CKIN41
          CKIN(43)=CKIN43
          IF(MINT(51).EQ.1) THEN
            IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
            IF(MFAIL.EQ.1) THEN
              MSTI(61)=1
              RETURN
            ENDIF
            GOTO 100
          ENDIF
          VINT(63)=PQM3**2
          VINT(64)=PQM4**2
        ENDIF
        IF(MIN(VINT(63),VINT(64)).LT.CKIN(6)**2) MINT(71)=1
        IF(MINT(71).EQ.1) VINT(71)=MAX(CKIN(3),CKIN(5))
      ENDIF
C...Prepare for additional variable choices in 2 -> 3.
      IF(ISTSB.EQ.5) THEN
        VINT(201)=0D0
        IF(KFPR(ISUB,2).GT.0) VINT(201)=PMAS(PYCOMP(KFPR(ISUB,2)),1)
        VINT(206)=VINT(201)
        VINT(204)=PMAS(23,1)
        IF(ISUB.EQ.124.OR.ISUB.EQ.351) VINT(204)=PMAS(24,1)
        IF(ISUB.EQ.352) VINT(204)=PMAS(PYCOMP(9900024),1)
        IF(ISUB.EQ.121.OR.ISUB.EQ.122.OR.ISUB.EQ.181.OR.ISUB.EQ.182.OR.
     &  ISUB.EQ.186.OR.ISUB.EQ.187) VINT(204)=VINT(201)
        VINT(209)=VINT(204)
      ENDIF
C...Select incoming VDM particle (rho/omega/phi/J/psi).
      IF(ISTSB.NE.0.AND.(MINT(101).GE.2.OR.MINT(102).GE.2).AND.
     &(MINT(123).EQ.2.OR.MINT(123).EQ.3.OR.MINT(123).EQ.7)) THEN
        VRN=PYR(0)*SIGT(0,0,5)
        IF(MINT(101).LE.1) THEN
          I1MN=0
          I1MX=0
        ELSE
          I1MN=1
          I1MX=MINT(101)
        ENDIF
        IF(MINT(102).LE.1) THEN
          I2MN=0
          I2MX=0
        ELSE
          I2MN=1
          I2MX=MINT(102)
        ENDIF
        DO 180 I1=I1MN,I1MX
          KFV1=110*I1+3
          DO 170 I2=I2MN,I2MX
            KFV2=110*I2+3
            VRN=VRN-SIGT(I1,I2,5)
            IF(VRN.LE.0D0) GOTO 190
  170     CONTINUE
  180   CONTINUE
  190   IF(MINT(101).GE.2) MINT(103)=KFV1
        IF(MINT(102).GE.2) MINT(104)=KFV2
      ENDIF
      IF(ISTSB.EQ.0) THEN
C...Elastic scattering or single or double diffractive scattering.
C...Select incoming particle (rho/omega/phi/J/psi for VDM) and mass.
        MINT(103)=MINT(11)
        MINT(104)=MINT(12)
        PMM(1)=VINT(3)
        PMM(2)=VINT(4)
        IF(MINT(101).GE.2.OR.MINT(102).GE.2) THEN
          JJ=ISUB-90
          VRN=PYR(0)*SIGT(0,0,JJ)
          IF(MINT(101).LE.1) THEN
            I1MN=0
            I1MX=0
          ELSE
            I1MN=1
            I1MX=MINT(101)
          ENDIF
          IF(MINT(102).LE.1) THEN
            I2MN=0
            I2MX=0
          ELSE
            I2MN=1
            I2MX=MINT(102)
          ENDIF
          DO 210 I1=I1MN,I1MX
            KFV1=110*I1+3
            DO 200 I2=I2MN,I2MX
              KFV2=110*I2+3
              VRN=VRN-SIGT(I1,I2,JJ)
              IF(VRN.LE.0D0) GOTO 220
  200       CONTINUE
  210     CONTINUE
  220     IF(MINT(101).GE.2) THEN
            MINT(103)=KFV1
            PMM(1)=PYMASS(KFV1)
          ENDIF
          IF(MINT(102).GE.2) THEN
            MINT(104)=KFV2
            PMM(2)=PYMASS(KFV2)
          ENDIF
        ENDIF
        VINT(67)=PMM(1)
        VINT(68)=PMM(2)
C...Select mass for GVMD states (rejecting previous assignment).
        Q0S=4D0*PARP(15)**2
        Q1S=4D0*VINT(154)**2
        LOOP3=0
  230   LOOP3=LOOP3+1
        DO 240 JT=1,2
          IF(MINT(106+JT).EQ.3) THEN
            PS=VINT(2+JT)**2
            PMM(JT)=(Q0S+PS)*(Q1S+PS)/
     &      (Q0S+PYR(0)*(Q1S-Q0S)+PS)-PS
            IF(MINT(102+JT).GE.333) PMM(JT)=PMM(JT)-
     &      PMAS(PYCOMP(113),1)+PMAS(PYCOMP(MINT(102+JT)),1)
          ENDIF
  240   CONTINUE
        IF(PMM(1)+PMM(2)+PARP(104).GE.VINT(1)) THEN
          IF(LOOP3.LT.100.AND.(MINT(107).EQ.3.OR.MINT(108).EQ.3))
     &    GOTO 230
          GOTO 100
        ENDIF
C...Side/sides of diffractive system.
        MINT(17)=0
        MINT(18)=0
        IF(ISUB.EQ.92.OR.ISUB.EQ.94) MINT(17)=1
        IF(ISUB.EQ.93.OR.ISUB.EQ.94) MINT(18)=1
C...Find masses of particles and minimal masses of diffractive states.
        DO 250 JT=1,2
          PDIF(JT)=PMM(JT)
          VINT(68+JT)=PDIF(JT)
          IF(MINT(16+JT).EQ.1) PDIF(JT)=PDIF(JT)+PARP(102)
  250   CONTINUE
        SH=VINT(2)
        SQM1=PMM(1)**2
        SQM2=PMM(2)**2
        SQM3=PDIF(1)**2
        SQM4=PDIF(2)**2
        SMRES1=(PMM(1)+PMRC)**2
        SMRES2=(PMM(2)+PMRC)**2
C...Find elastic slope and lower limit diffractive slope.
        IHA=MAX(2,IABS(MINT(103))/110)
        IF(IHA.GE.5) IHA=1
        IHB=MAX(2,IABS(MINT(104))/110)
        IF(IHB.GE.5) IHB=1
        IF(ISUB.EQ.91) THEN
          BMN=2D0*BHAD(IHA)+2D0*BHAD(IHB)+4D0*SH**EPS-4.2D0
        ELSEIF(ISUB.EQ.92) THEN
          BMN=MAX(2D0,2D0*BHAD(IHB))
        ELSEIF(ISUB.EQ.93) THEN
          BMN=MAX(2D0,2D0*BHAD(IHA))
        ELSEIF(ISUB.EQ.94) THEN
          BMN=2D0*ALP*4D0
        ENDIF
C...Determine maximum possible t range and coefficient of generation.
        SQLA12=(SH-SQM1-SQM2)**2-4D0*SQM1*SQM2
        SQLA34=(SH-SQM3-SQM4)**2-4D0*SQM3*SQM4
        THA=SH-(SQM1+SQM2+SQM3+SQM4)+(SQM1-SQM2)*(SQM3-SQM4)/SH
        THB=SQRT(MAX(0D0,SQLA12))*SQRT(MAX(0D0,SQLA34))/SH
        THC=(SQM3-SQM1)*(SQM4-SQM2)+(SQM1+SQM4-SQM2-SQM3)*
     &  (SQM1*SQM4-SQM2*SQM3)/SH
        THL=-0.5D0*(THA+THB)
        THU=THC/THL
        THRND=EXP(MAX(-50D0,BMN*(THL-THU)))-1D0
C...Select diffractive mass/masses according to dm^2/m^2.
        LOOP3=0
  260   LOOP3=LOOP3+1
        DO 270 JT=1,2
          IF(MINT(16+JT).EQ.0) THEN
            PDIF(2+JT)=PDIF(JT)
          ELSE
            PMMIN=PDIF(JT)
            PMMAX=MAX(VINT(2+JT),VINT(1)-PDIF(3-JT))
            PDIF(2+JT)=PMMIN*(PMMAX/PMMIN)**PYR(0)
          ENDIF
  270   CONTINUE
        SQM3=PDIF(3)**2
        SQM4=PDIF(4)**2
C..Additional mass factors, including resonance enhancement.
        IF(PDIF(3)+PDIF(4).GE.VINT(1)) THEN
          IF(LOOP3.LT.100) GOTO 260
          GOTO 100
        ENDIF
        IF(ISUB.EQ.92) THEN
          FSD=(1D0-SQM3/SH)*(1D0+CRES*SMRES1/(SMRES1+SQM3))
          IF(FSD.LT.PYR(0)*(1D0+CRES)) GOTO 260
        ELSEIF(ISUB.EQ.93) THEN
          FSD=(1D0-SQM4/SH)*(1D0+CRES*SMRES2/(SMRES2+SQM4))
          IF(FSD.LT.PYR(0)*(1D0+CRES)) GOTO 260
        ELSEIF(ISUB.EQ.94) THEN
          FDD=(1D0-(PDIF(3)+PDIF(4))**2/SH)*(SH*SMP/
     &    (SH*SMP+SQM3*SQM4))*(1D0+CRES*SMRES1/(SMRES1+SQM3))*
     &    (1D0+CRES*SMRES2/(SMRES2+SQM4))
          IF(FDD.LT.PYR(0)*(1D0+CRES)**2) GOTO 260
        ENDIF
C...Select t according to exp(Bmn*t) and correct to right slope.
        TH=THU+LOG(1D0+THRND*PYR(0))/BMN
        IF(ISUB.GE.92) THEN
          IF(ISUB.EQ.92) THEN
            BADD=2D0*ALP*LOG(SH/SQM3)
            IF(BHAD(IHB).LT.1D0) BADD=MAX(0D0,BADD+2D0*BHAD(IHB)-2D0)
          ELSEIF(ISUB.EQ.93) THEN
            BADD=2D0*ALP*LOG(SH/SQM4)
            IF(BHAD(IHA).LT.1D0) BADD=MAX(0D0,BADD+2D0*BHAD(IHA)-2D0)
          ELSEIF(ISUB.EQ.94) THEN
            BADD=2D0*ALP*(LOG(EXP(4D0)+SH/(ALP*SQM3*SQM4))-4D0)
          ENDIF
          IF(EXP(MAX(-50D0,BADD*(TH-THU))).LT.PYR(0)) GOTO 260
        ENDIF
C...Check whether m^2 and t choices are consistent.
        SQLA34=(SH-SQM3-SQM4)**2-4D0*SQM3*SQM4
        THA=SH-(SQM1+SQM2+SQM3+SQM4)+(SQM1-SQM2)*(SQM3-SQM4)/SH
        THB=SQRT(MAX(0D0,SQLA12))*SQRT(MAX(0D0,SQLA34))/SH
        IF(THB.LE.1D-8) GOTO 260
        THC=(SQM3-SQM1)*(SQM4-SQM2)+(SQM1+SQM4-SQM2-SQM3)*
     &  (SQM1*SQM4-SQM2*SQM3)/SH
        THLM=-0.5D0*(THA+THB)
        THUM=THC/THLM
        IF(TH.LT.THLM.OR.TH.GT.THUM) GOTO 260
C...Information to output.
        VINT(21)=1D0
        VINT(22)=0D0
        VINT(23)=MIN(1D0,MAX(-1D0,(THA+2D0*TH)/THB))
        VINT(45)=TH
        VINT(59)=2D0*SQRT(MAX(0D0,-(THC+THA*TH+TH**2)))/THB
        VINT(63)=PDIF(3)**2
        VINT(64)=PDIF(4)**2
        VINT(283)=PMM(1)**2/4D0
        VINT(284)=PMM(2)**2/4D0
C...Note: in the following, by In is meant the integral over the
C...quantity multiplying coefficient cn.
C...Choose tau according to h1(tau)/tau, where
C...h1(tau) = c1 + I1/I2*c2*1/tau + I1/I3*c3*1/(tau+tau_R) +
C...I1/I4*c4*tau/((s*tau-m^2)^2+(m*Gamma)^2) +
C...I1/I5*c5*1/(tau+tau_R') +
C...I1/I6*c6*tau/((s*tau-m'^2)^2+(m'*Gamma')^2) +
C...I1/I7*c7*tau/(1.-tau), and
C...c1 + c2 + c3 + c4 + c5 + c6 + c7 = 1.
      ELSEIF(ISTSB.GE.1.AND.ISTSB.LE.5) THEN
        CALL PYKLIM(1)
        IF(MINT(51).NE.0) THEN
          IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
          IF(MFAIL.EQ.1) THEN
            MSTI(61)=1
            RETURN
          ENDIF
          GOTO 100
        ENDIF
        RTAU=PYR(0)
        MTAU=1
        IF(RTAU.GT.COEF(ISUB,1)) MTAU=2
        IF(RTAU.GT.COEF(ISUB,1)+COEF(ISUB,2)) MTAU=3
        IF(RTAU.GT.COEF(ISUB,1)+COEF(ISUB,2)+COEF(ISUB,3)) MTAU=4
        IF(RTAU.GT.COEF(ISUB,1)+COEF(ISUB,2)+COEF(ISUB,3)+COEF(ISUB,4))
     &  MTAU=5
        IF(RTAU.GT.COEF(ISUB,1)+COEF(ISUB,2)+COEF(ISUB,3)+COEF(ISUB,4)+
     &  COEF(ISUB,5)) MTAU=6
        IF(RTAU.GT.COEF(ISUB,1)+COEF(ISUB,2)+COEF(ISUB,3)+COEF(ISUB,4)+
     &  COEF(ISUB,5)+COEF(ISUB,6)) MTAU=7
        CALL PYKMAP(1,MTAU,PYR(0))
C...2 -> 3, 4 processes:
C...Choose tau' according to h4(tau,tau')/tau', where
C...h4(tau,tau') = c1 + I1/I2*c2*(1 - tau/tau')^3/tau' +
C...I1/I3*c3*1/(1 - tau'), and c1 + c2 + c3 = 1.
        IF(ISTSB.GE.3.AND.ISTSB.LE.5) THEN
          CALL PYKLIM(4)
          IF(MINT(51).NE.0) THEN
            IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
            IF(MFAIL.EQ.1) THEN
              MSTI(61)=1
              RETURN
            ENDIF
            GOTO 100
          ENDIF
          RTAUP=PYR(0)
          MTAUP=1
          IF(RTAUP.GT.COEF(ISUB,18)) MTAUP=2
          IF(RTAUP.GT.COEF(ISUB,18)+COEF(ISUB,19)) MTAUP=3
          CALL PYKMAP(4,MTAUP,PYR(0))
        ENDIF
C...Choose y* according to h2(y*), where
C...h2(y*) = I0/I1*c1*(y*-y*min) + I0/I2*c2*(y*max-y*) +
C...I0/I3*c3*1/cosh(y*) + I0/I4*c4*1/(1-exp(y*-y*max)) +
C...I0/I5*c5*1/(1-exp(-y*-y*min)), I0 = y*max-y*min,
C...and c1 + c2 + c3 + c4 + c5 = 1.
        CALL PYKLIM(2)
        IF(MINT(51).NE.0) THEN
          IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
          IF(MFAIL.EQ.1) THEN
            MSTI(61)=1
            RETURN
          ENDIF
          GOTO 100
        ENDIF
        RYST=PYR(0)
        MYST=1
        IF(RYST.GT.COEF(ISUB,8)) MYST=2
        IF(RYST.GT.COEF(ISUB,8)+COEF(ISUB,9)) MYST=3
        IF(RYST.GT.COEF(ISUB,8)+COEF(ISUB,9)+COEF(ISUB,10)) MYST=4
        IF(RYST.GT.COEF(ISUB,8)+COEF(ISUB,9)+COEF(ISUB,10)+
     &  COEF(ISUB,11)) MYST=5
        CALL PYKMAP(2,MYST,PYR(0))
C...2 -> 2 processes:
C...Choose cos(theta-hat) (cth) according to h3(cth), where
C...h3(cth) = c0 + I0/I1*c1*1/(A - cth) + I0/I2*c2*1/(A + cth) +
C...I0/I3*c3*1/(A - cth)^2 + I0/I4*c4*1/(A + cth)^2,
C...A = 1 + 2*(m3*m4/sh)^2 (= 1 for massless products),
C...and c0 + c1 + c2 + c3 + c4 = 1.
        CALL PYKLIM(3)
        IF(MINT(51).NE.0) THEN
          IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
          IF(MFAIL.EQ.1) THEN
            MSTI(61)=1
            RETURN
          ENDIF
          GOTO 100
        ENDIF
        IF(ISTSB.EQ.2.OR.ISTSB.EQ.4) THEN
          RCTH=PYR(0)
          MCTH=1
          IF(RCTH.GT.COEF(ISUB,13)) MCTH=2
          IF(RCTH.GT.COEF(ISUB,13)+COEF(ISUB,14)) MCTH=3
          IF(RCTH.GT.COEF(ISUB,13)+COEF(ISUB,14)+COEF(ISUB,15)) MCTH=4
          IF(RCTH.GT.COEF(ISUB,13)+COEF(ISUB,14)+COEF(ISUB,15)+
     &    COEF(ISUB,16)) MCTH=5
          CALL PYKMAP(3,MCTH,PYR(0))
        ENDIF
C...2 -> 3 : select pT1, phi1, pT2, phi2, y3 for 3 outgoing.
        IF(ISTSB.EQ.5) THEN
          CALL PYKMAP(5,0,0D0)
          IF(MINT(51).NE.0) THEN
            IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
            IF(MFAIL.EQ.1) THEN
              MSTI(61)=1
              RETURN
            ENDIF
            GOTO 100
          ENDIF
        ENDIF
C...DIS as f + gamma* -> f process: set dummy values.
      ELSEIF(ISTSB.EQ.8) THEN
        VINT(21)=0.9D0
        VINT(22)=0D0
        VINT(23)=0D0
        VINT(47)=0D0
        VINT(48)=0D0
C...Low-pT or multiple interactions (first semihard interaction).
      ELSEIF(ISTSB.EQ.9) THEN
        CALL PYMULT(3)
        ISUB=MINT(1)
C...Study user-defined process: kinematics plus weight.
      ELSEIF(ISTSB.EQ.11) THEN
        IF(IDWTUP.GT.0.AND.XWGTUP.LT.0D0) CALL
     &  PYERRM(26,'(PYRAND:) Negative XWGTUP for user process')
        MSTI(51)=0
        IF(NUP.LE.0) THEN
          MINT(51)=2
          MSTI(51)=1
          IF(MINT(82).EQ.1) THEN
            NGEN(0,1)=NGEN(0,1)-1
            NGEN(ISUB,1)=NGEN(ISUB,1)-1
          ENDIF
          IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
          RETURN
        ENDIF
C...Extract cross section event weight.
        IF(IABS(IDWTUP).EQ.1.OR.IABS(IDWTUP).EQ.4) THEN
          SIGS=1D-9*XWGTUP
        ELSE
          SIGS=1D-9*XSECUP(KFPR(ISUB,1))
        ENDIF
        IF(IABS(IDWTUP).GE.1.AND.IABS(IDWTUP).LE.3) THEN
          VINT(97)=SIGN(1D0,XWGTUP)
        ELSE
          VINT(97)=1D-9*XWGTUP
        ENDIF
C...Construct 'trivial' kinematical variables needed.
        KFL1=IDUP(1)
        KFL2=IDUP(2)
        VINT(41)=PUP(4,1)/EBMUP(1)
        VINT(42)=PUP(4,2)/EBMUP(2)
        VINT(21)=VINT(41)*VINT(42)
        VINT(22)=0.5D0*LOG(VINT(41)/VINT(42))
        VINT(44)=VINT(21)*VINT(2)
        VINT(43)=SQRT(MAX(0D0,VINT(44)))
        VINT(55)=SCALUP
        IF(SCALUP.LE.0D0) VINT(55)=VINT(43)
        VINT(56)=VINT(55)**2
        VINT(57)=AQEDUP
        VINT(58)=AQCDUP
C...Construct other kinematical variables needed (approximately).
        VINT(23)=0D0
        VINT(26)=VINT(21)
        VINT(45)=-0.5D0*VINT(44)
        VINT(46)=-0.5D0*VINT(44)
        VINT(49)=VINT(43)
        VINT(50)=VINT(44)
        VINT(51)=VINT(55)
        VINT(52)=VINT(56)
        VINT(53)=VINT(55)
        VINT(54)=VINT(56)
        VINT(25)=0D0
        VINT(48)=0D0
        IF(ISTUP(1).NE.-1.OR.ISTUP(2).NE.-1) CALL PYERRM(26,
     &  '(PYRAND:) unacceptable ISTUP code for incoming particles')
        DO 280 IUP=3,NUP
          IF(ISTUP(IUP).LT.1.OR.ISTUP(IUP).GT.3) CALL PYERRM(26,
     &    '(PYRAND:) unacceptable ISTUP code for particles')
          IF(ISTUP(IUP).EQ.1) VINT(25)=VINT(25)+2D0*(PUP(5,IUP)**2+
     &    PUP(1,IUP)**2+PUP(2,IUP)**2)/VINT(2)
          IF(ISTUP(IUP).EQ.1) VINT(48)=VINT(48)+0.5D0*(PUP(1,IUP)**2+
     &    PUP(2,IUP)**2)
  280   CONTINUE
        VINT(47)=SQRT(VINT(48))
      ENDIF
C...Choose azimuthal angle.
      VINT(24)=0D0
      IF(ISTSB.NE.11) VINT(24)=PARU(2)*PYR(0)
C...Check against user cuts on kinematics at parton level.
      MINT(51)=0
      IF((ISUB.LE.90.OR.ISUB.GT.100).AND.ISTSB.LE.10) CALL PYKLIM(0)
      IF(MINT(51).NE.0) THEN
        IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
        IF(MFAIL.EQ.1) THEN
          MSTI(61)=1
          RETURN
        ENDIF
        GOTO 100
      ENDIF
      IF(MINT(82).EQ.1.AND.MSTP(141).GE.1.AND.ISTSB.LE.10) THEN
        MCUT=0
        IF(MSUB(91)+MSUB(92)+MSUB(93)+MSUB(94)+MSUB(95).EQ.0)
     &  CALL PYKCUT(MCUT)
        IF(MCUT.NE.0) THEN
          IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
          IF(MFAIL.EQ.1) THEN
            MSTI(61)=1
            RETURN
          ENDIF
          GOTO 100
        ENDIF
      ENDIF
C...Calculate differential cross-section for different subprocesses.
      IF(ISTSB.LE.10) CALL PYSIGH(NCHN,SIGS)
      SIGSOR=SIGS
      SIGLPT=SIGT(0,0,5)*VINT(315)*VINT(316)
C...Multiply cross section by lepton -> photon flux factor.
      IF(MINT(141).NE.0.OR.MINT(142).NE.0) THEN
        SIGS=WTGAGA*SIGS
        DO 290 ICHN=1,NCHN
          SIGH(ICHN)=WTGAGA*SIGH(ICHN)
  290   CONTINUE
        SIGLPT=WTGAGA*SIGLPT
      ENDIF
C...Multiply cross-section by user-defined weights.
      IF(MSTP(173).EQ.1) THEN
        SIGS=PARP(173)*SIGS
        DO 300 ICHN=1,NCHN
          SIGH(ICHN)=PARP(173)*SIGH(ICHN)
  300   CONTINUE
        SIGLPT=PARP(173)*SIGLPT
      ENDIF
      WTXS=1D0
      SIGSWT=SIGS
      VINT(99)=1D0
      VINT(100)=1D0
      IF(MINT(82).EQ.1.AND.MSTP(142).GE.1) THEN
        IF(ISUB.NE.96.AND.MSUB(91)+MSUB(92)+MSUB(93)+MSUB(94)+
     &  MSUB(95).EQ.0) CALL PYEVWT(WTXS)
        SIGSWT=WTXS*SIGS
        VINT(99)=WTXS
        IF(MSTP(142).EQ.1) VINT(100)=1D0/WTXS
      ENDIF
C...Calculations for Monte Carlo estimate of all cross-sections.
      IF(MINT(82).EQ.1.AND.ISUB.LE.90.OR.ISUB.GE.96) THEN
        IF(MSTP(142).LE.1) THEN
          XSEC(ISUB,2)=XSEC(ISUB,2)+SIGS
        ELSE
          XSEC(ISUB,2)=XSEC(ISUB,2)+SIGSWT
        ENDIF
      ELSEIF(MINT(82).EQ.1) THEN
        XSEC(ISUB,2)=XSEC(ISUB,2)+SIGS
      ENDIF
      IF((ISUB.EQ.95.OR.ISUB.EQ.96).AND.LOOP2.EQ.1.AND.
     &MINT(82).EQ.1) XSEC(97,2)=XSEC(97,2)+SIGLPT
C...Multiple interactions: store results of cross-section calculation.
      IF(MINT(50).EQ.1.AND.MSTP(82).GE.3) THEN
        VINT(153)=SIGSOR
        CALL PYMULT(4)
      ENDIF
C...Ratio of actual to maximum cross section.
      IF(ISTSB.NE.11) THEN
        VIOL=SIGSWT/XSEC(ISUB,1)
        IF(ISUB.EQ.96.AND.MSTP(173).EQ.1) VIOL=VIOL/PARP(174)
      ELSEIF(IDWTUP.EQ.1.OR.IDWTUP.EQ.2) THEN
        VIOL=XWGTUP/XMAXUP(KFPR(ISUB,1))
      ELSEIF(IDWTUP.EQ.-1.OR.IDWTUP.EQ.-2) THEN
        VIOL=ABS(XWGTUP)/ABS(XMAXUP(KFPR(ISUB,1)))
      ELSE
        VIOL=1D0
      ENDIF
C...Check that weight not negative.
      IF(MSTP(123).LE.0) THEN
        IF(VIOL.LT.-1D-3) THEN
          WRITE(MSTU(11),5000) VIOL,NGEN(0,3)+1
          IF(MSTP(122).GE.1) WRITE(MSTU(11),5100) ISUB,VINT(21),
     &    VINT(22),VINT(23),VINT(26)
          STOP
        ENDIF
      ELSE
        IF(VIOL.LT.MIN(-1D-3,VINT(109))) THEN
          VINT(109)=VIOL
          WRITE(MSTU(11),5200) VIOL,NGEN(0,3)+1
          IF(MSTP(122).GE.1) WRITE(MSTU(11),5100) ISUB,VINT(21),
     &    VINT(22),VINT(23),VINT(26)
        ENDIF
      ENDIF
C...Weighting using estimate of maximum of differential cross-section.
      IF(MFAIL.EQ.0.AND.ISUB.NE.95.AND.ISUB.NE.96) THEN
        IF(VIOL.LT.PYR(0)) THEN
          IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
          IF(ISUB.GE.91.AND.ISUB.LE.94) ISUB=0
          GOTO 100
        ENDIF
      ELSEIF(MFAIL.EQ.0) THEN
        RATND=SIGLPT/XSEC(95,1)
        VIOL=VIOL/RATND
        IF(LOOP2.EQ.1.AND.RATND.LT.PYR(0)) THEN
          IF(VIOL.GT.PYR(0).AND.MINT(82).EQ.1.AND.MSUB(95).EQ.1.AND.
     &    (ISUB.LE.90.OR.ISUB.GE.95)) NGEN(95,1)=NGEN(95,1)+MINT(143)
          IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
          ISUB=0
          GOTO 100
        ENDIF
        IF(VIOL.LT.PYR(0)) THEN
          GOTO 140
        ENDIF
      ELSEIF(ISUB.NE.95.AND.ISUB.NE.96) THEN
        IF(VIOL.LT.PYR(0)) THEN
          MSTI(61)=1
          IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
          RETURN
        ENDIF
      ELSE
        RATND=SIGLPT/XSEC(95,1)
        IF(LOOP.EQ.1.AND.RATND.LT.PYR(0)) THEN
          MSTI(61)=1
          IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
          RETURN
        ENDIF
        VIOL=VIOL/RATND
        IF(VIOL.LT.PYR(0)) THEN
          IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
          GOTO 100
        ENDIF
      ENDIF
C...Check for possible violation of estimated maximum of differential
C...cross-section used in weighting.
      IF(MSTP(123).LE.0) THEN
        IF(VIOL.GT.1D0) THEN
          WRITE(MSTU(11),5300) VIOL,NGEN(0,3)+1
          IF(MSTP(122).GE.2) WRITE(MSTU(11),5100) ISUB,VINT(21),
     &    VINT(22),VINT(23),VINT(26)
          STOP
        ENDIF
      ELSEIF(MSTP(123).EQ.1) THEN
        IF(VIOL.GT.VINT(108)) THEN
          VINT(108)=VIOL
          IF(VIOL.GT.1.0001D0) THEN
            MINT(10)=1
            WRITE(MSTU(11),5400) VIOL,NGEN(0,3)+1
            IF(MSTP(122).GE.2) WRITE(MSTU(11),5100) ISUB,VINT(21),
     &      VINT(22),VINT(23),VINT(26)
          ENDIF
        ENDIF
      ELSEIF(VIOL.GT.VINT(108)) THEN
        VINT(108)=VIOL
        IF(VIOL.GT.1D0) THEN
          MINT(10)=1
          WRITE(MSTU(11),5400) VIOL,NGEN(0,3)+1
          IF(ISTSB.EQ.11.AND.(IABS(IDWTUP).EQ.1.OR.IABS(IDWTUP).EQ.2))
     &    THEN
            XMAXUP(KFPR(ISUB,1))=VIOL*XMAXUP(KFPR(ISUB,1))
            IF(KFPR(ISUB,1).LE.9) THEN
              WRITE(MSTU(11),5800) KFPR(ISUB,1),XMAXUP(KFPR(ISUB,1))
            ELSEIF(KFPR(ISUB,1).LE.99) THEN
              WRITE(MSTU(11),5900) KFPR(ISUB,1),XMAXUP(KFPR(ISUB,1))
            ELSE
              WRITE(MSTU(11),6000) KFPR(ISUB,1),XMAXUP(KFPR(ISUB,1))
            ENDIF
          ENDIF
          IF(ISTSB.NE.11.OR.IABS(IDWTUP).EQ.1) THEN
            XDIF=XSEC(ISUB,1)*(VIOL-1D0)
            XSEC(ISUB,1)=XSEC(ISUB,1)+XDIF
            IF(MSUB(ISUB).EQ.1.AND.(ISUB.LE.90.OR.ISUB.GT.96))
     &      XSEC(0,1)=XSEC(0,1)+XDIF
            IF(MSTP(122).GE.2) WRITE(MSTU(11),5100) ISUB,VINT(21),
     &      VINT(22),VINT(23),VINT(26)
            IF(ISUB.LE.9) THEN
              WRITE(MSTU(11),5500) ISUB,XSEC(ISUB,1)
            ELSEIF(ISUB.LE.99) THEN
              WRITE(MSTU(11),5600) ISUB,XSEC(ISUB,1)
            ELSE
              WRITE(MSTU(11),5700) ISUB,XSEC(ISUB,1)
            ENDIF
          ENDIF
          VINT(108)=1D0
        ENDIF
      ENDIF
C...Multiple interactions: choose impact parameter.
      VINT(148)=1D0
      IF(MINT(50).EQ.1.AND.(ISUB.LE.90.OR.ISUB.GE.96).AND.
     &MSTP(82).GE.3) THEN
        CALL PYMULT(5)
        IF(VINT(150).LT.PYR(0)) THEN
          IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
          IF(MFAIL.EQ.1) THEN
            MSTI(61)=1
            RETURN
          ENDIF
          GOTO 100
        ENDIF
      ENDIF
      IF(MINT(82).EQ.1) NGEN(0,2)=NGEN(0,2)+1
      IF(MINT(82).EQ.1.AND.MSUB(95).EQ.1) THEN
        IF(ISUB.LE.90.OR.ISUB.GE.95) NGEN(95,1)=NGEN(95,1)+MINT(143)
        IF(ISUB.LE.90.OR.ISUB.GE.96) NGEN(96,2)=NGEN(96,2)+1
      ENDIF
      IF(ISUB.LE.90.OR.ISUB.GE.96) MINT(31)=MINT(31)+1
C...Choose flavour of reacting partons (and subprocess).
      IF(ISTSB.GE.11) GOTO 320
      RSIGS=SIGS*PYR(0)
      QT2=VINT(48)
      RQQBAR=PARP(87)*(1D0-(QT2/(QT2+(PARP(88)*PARP(82)*
     &(VINT(1)/PARP(89))**PARP(90))**2))**2)
      IF(ISUB.NE.95.AND.(ISUB.NE.96.OR.MSTP(82).LE.1.OR.
     &PYR(0).GT.RQQBAR)) THEN
        DO 310 ICHN=1,NCHN
          KFL1=ISIG(ICHN,1)
          KFL2=ISIG(ICHN,2)
          MINT(2)=ISIG(ICHN,3)
          RSIGS=RSIGS-SIGH(ICHN)
          IF(RSIGS.LE.0D0) GOTO 320
  310   CONTINUE
C...Multiple interactions: choose qqbar preferentially at small pT.
      ELSEIF(ISUB.EQ.96) THEN
        MINT(105)=MINT(103)
        MINT(109)=MINT(107)
        CALL PYSPLI(MINT(11),21,KFL1,KFLDUM)
        MINT(105)=MINT(104)
        MINT(109)=MINT(108)
        CALL PYSPLI(MINT(12),21,KFL2,KFLDUM)
        MINT(1)=11
        MINT(2)=1
        IF(KFL1.EQ.KFL2.AND.PYR(0).LT.0.5D0) MINT(2)=2
C...Low-pT: choose string drawing configuration.
      ELSE
        KFL1=21
        KFL2=21
        RSIGS=6D0*PYR(0)
        MINT(2)=1
        IF(RSIGS.GT.1D0) MINT(2)=2
        IF(RSIGS.GT.2D0) MINT(2)=3
      ENDIF
C...Reassign QCD process. Partons before initial state radiation.
  320 IF(MINT(2).GT.10) THEN
        MINT(1)=MINT(2)/10
        MINT(2)=MOD(MINT(2),10)
      ENDIF
      IF(MINT(82).EQ.1.AND.MSTP(111).GE.0) NGEN(MINT(1),2)=
     &NGEN(MINT(1),2)+1
      MINT(15)=KFL1
      MINT(16)=KFL2
      MINT(13)=MINT(15)
      MINT(14)=MINT(16)
      VINT(141)=VINT(41)
      VINT(142)=VINT(42)
      VINT(151)=0D0
      VINT(152)=0D0
C...Calculate x value of photon for parton inside photon inside e.
      DO 350 JT=1,2
        MINT(18+JT)=0
        VINT(154+JT)=0D0
        MSPLI=0
        IF(JT.EQ.1.AND.MINT(43).LE.2) MSPLI=1
        IF(JT.EQ.2.AND.MOD(MINT(43),2).EQ.1) MSPLI=1
        IF(IABS(MINT(14+JT)).LE.8.OR.MINT(14+JT).EQ.21) MSPLI=MSPLI+1
        IF(MSPLI.EQ.2) THEN
          KFLH=MINT(14+JT)
          XHRD=VINT(140+JT)
          Q2HRD=VINT(54)
          MINT(105)=MINT(102+JT)
          MINT(109)=MINT(106+JT)
          VINT(120)=VINT(2+JT)
          IF(MSTP(57).LE.1) THEN
            CALL PYPDFU(22,XHRD,Q2HRD,XPQ)
          ELSE
            CALL PYPDFL(22,XHRD,Q2HRD,XPQ)
          ENDIF
          WTMX=4D0*XPQ(KFLH)
          IF(MSTP(13).EQ.2) THEN
            Q2PMS=Q2HRD/PMAS(11,1)**2
            WTMX=WTMX*LOG(MAX(2D0,Q2PMS*(1D0-XHRD)/XHRD**2))
          ENDIF
  330     XE=XHRD**PYR(0)
          XG=MIN(1D0-1D-10,XHRD/XE)
          IF(MSTP(57).LE.1) THEN
            CALL PYPDFU(22,XG,Q2HRD,XPQ)
          ELSE
            CALL PYPDFL(22,XG,Q2HRD,XPQ)
          ENDIF
          WT=(1D0+(1D0-XE)**2)*XPQ(KFLH)
          IF(MSTP(13).EQ.2) WT=WT*LOG(MAX(2D0,Q2PMS*(1D0-XE)/XE**2))
          IF(WT.LT.PYR(0)*WTMX) GOTO 330
          MINT(18+JT)=1
          VINT(154+JT)=XE
          DO 340 KFLS=-25,25
            XSFX(JT,KFLS)=XPQ(KFLS)
  340     CONTINUE
        ENDIF
  350 CONTINUE
C...Pick scale where photon is resolved.
      Q0S=PARP(15)**2
      Q1S=VINT(154)**2
      VINT(283)=0D0
      IF(MINT(107).EQ.3) THEN
        IF(MSTP(66).EQ.1) THEN
          VINT(283)=Q0S*(VINT(54)/Q0S)**PYR(0)
        ELSEIF(MSTP(66).EQ.2) THEN
          PS=VINT(3)**2
          Q2EFF=VINT(54)*((Q0S+PS)/(VINT(54)+PS))*
     &    EXP(PS*(VINT(54)-Q0S)/((VINT(54)+PS)*(Q0S+PS)))
          Q2INT=SQRT(Q0S*Q2EFF)
          VINT(283)=Q2INT*(VINT(54)/Q2INT)**PYR(0)
        ELSEIF(MSTP(66).EQ.3) THEN
          VINT(283)=Q0S*(Q1S/Q0S)**PYR(0)
        ELSEIF(MSTP(66).GE.4) THEN
          PS=0.25D0*VINT(3)**2
          VINT(283)=(Q0S+PS)*(Q1S+PS)/
     &    (Q0S+PYR(0)*(Q1S-Q0S)+PS)-PS
        ENDIF
      ENDIF
      VINT(284)=0D0
      IF(MINT(108).EQ.3) THEN
        IF(MSTP(66).EQ.1) THEN
          VINT(284)=Q0S*(VINT(54)/Q0S)**PYR(0)
        ELSEIF(MSTP(66).EQ.2) THEN
          PS=VINT(4)**2
          Q2EFF=VINT(54)*((Q0S+PS)/(VINT(54)+PS))*
     &    EXP(PS*(VINT(54)-Q0S)/((VINT(54)+PS)*(Q0S+PS)))
          Q2INT=SQRT(Q0S*Q2EFF)
          VINT(284)=Q2INT*(VINT(54)/Q2INT)**PYR(0)
        ELSEIF(MSTP(66).EQ.3) THEN
          VINT(284)=Q0S*(Q1S/Q0S)**PYR(0)
        ELSEIF(MSTP(66).GE.4) THEN
          PS=0.25D0*VINT(4)**2
          VINT(284)=(Q0S+PS)*(Q1S+PS)/
     &    (Q0S+PYR(0)*(Q1S-Q0S)+PS)-PS
        ENDIF
      ENDIF
      IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
C...Format statements for differential cross-section maximum violations.
 5000 FORMAT(/1X,'Error: negative cross-section fraction',1P,D11.3,1X,
     &'in event',1X,I7,'D0'/1X,'Execution stopped!')
 5100 FORMAT(1X,'ISUB = ',I3,'; Point of violation:'/1X,'tau =',1P,
     &D11.3,', y* =',D11.3,', cthe = ',0P,F11.7,', tau'' =',1P,D11.3)
 5200 FORMAT(/1X,'Warning: negative cross-section fraction',1P,D11.3,1X,
     &'in event',1X,I7)
 5300 FORMAT(/1X,'Error: maximum violated by',1P,D11.3,1X,
     &'in event',1X,I7,'D0'/1X,'Execution stopped!')
 5400 FORMAT(/1X,'Advisory warning: maximum violated by',1P,D11.3,1X,
     &'in event',1X,I7)
 5500 FORMAT(1X,'XSEC(',I1,',1) increased to',1P,D11.3)
 5600 FORMAT(1X,'XSEC(',I2,',1) increased to',1P,D11.3)
 5700 FORMAT(1X,'XSEC(',I3,',1) increased to',1P,D11.3)
 5800 FORMAT(1X,'XMAXUP(',I1,') increased to',1P,D11.3)
 5900 FORMAT(1X,'XMAXUP(',I2,') increased to',1P,D11.3)
 6000 FORMAT(1X,'XMAXUP(',I3,') increased to',1P,D11.3)
      RETURN
      END
