pyinpr.f

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pyinpr.f

 
C*********************************************************************
 
C...PYINPR
C...Selects partonic subprocesses to be included in the simulation.
 
      SUBROUTINE pyinpr
 
C...Double precision and integer declarations.
      IMPLICIT DOUBLE PRECISION(a-h, o-z)
      IMPLICIT INTEGER(i-n)
      INTEGER pyk,pychge,pycomp
 
C...User process initialization commonblock.
      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)
      SAVE /heprup/
 
C...Commonblocks and character variables.
      common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
      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/pyint6/proc(0:500)
      CHARACTER proc*28
      SAVE /pydat1/,/pydat3/,/pysubs/,/pypars/,/pyint1/,/pyint2/,
     &/pyint6/
      CHARACTER chipr*10
 
C...Reset processes to be included.
      IF(msel.NE.0) THEN
        DO 100 i=1,500
          msub(i)=0
  100   CONTINUE
      ENDIF
 
C...Set running pTmin scale.
      IF(mstp(82).LE.1) THEN
        ptmrun=parp(81)*(vint(1)/parp(89))**parp(90)
      ELSE
        ptmrun=parp(82)*(vint(1)/parp(89))**parp(90)
      ENDIF
 
C...Begin by assuming incoming photon to enter subprocess.
      IF(mint(11).EQ.22) mint(15)=22
      IF(mint(12).EQ.22) mint(16)=22
 
C...For e-gamma with MSTP(14)=10 allow mixture of VMD and anomalous.
      IF(mint(121).EQ.2.AND.mstp(14).EQ.10) THEN
        msub(10)=1
        mint(123)=mint(122)+1
 
C...For gamma-p or gamma-gamma with MSTP(14) = 10, 20, 25 or 30
C...allow mixture.
C...Here also set a few parameters otherwise normally not touched.
      ELSEIF(mint(121).GT.1) THEN
 
C...Parton distributions dampened at small Q2; go to low energies,
C...alpha_s <1; no minimum pT cut-off a priori.
        IF(mstp(18).EQ.2) THEN
          mstp(57)=3
          parp(2)=2d0
          paru(115)=1d0
          ckin(5)=0.2d0
          ckin(6)=0.2d0
        ENDIF
 
C...Define pT cut-off parameters and whether run involves low-pT.
        ptmvmd=ptmrun
        vint(154)=ptmvmd
        ptmdir=ptmvmd
        IF(mstp(18).EQ.2) ptmdir=parp(15)
        ptmano=ptmvmd
        IF(mstp(15).EQ.5) ptmano=0.60d0+
     &  0.125d0*log(1d0+0.10d0*vint(1))**2
        iptl=1
        IF(vint(285).GT.max(ptmvmd,ptmdir,ptmano)) iptl=0
        IF(msel.EQ.2) iptl=1
 
C...Set up for p/gamma * gamma; real or virtual photons.
        IF(mint(121).EQ.3.OR.mint(121).EQ.6.OR.(mint(121).EQ.4.AND.
     &  mstp(14).EQ.30)) THEN
 
C...Set up for p/VMD * VMD.
        IF(mint(122).EQ.1) THEN
          mint(123)=2
          msub(11)=1
          msub(12)=1
          msub(13)=1
          msub(28)=1
          msub(53)=1
          msub(68)=1
          IF(iptl.EQ.1) msub(95)=1
          IF(msel.EQ.2) THEN
            msub(91)=1
            msub(92)=1
            msub(93)=1
            msub(94)=1
          ENDIF
          IF(iptl.EQ.1) ckin(3)=0d0
 
C...Set up for p/VMD * direct gamma.
        ELSEIF(mint(122).EQ.2) THEN
          mint(123)=0
          IF(mint(121).EQ.6) mint(123)=5
          msub(131)=1
          msub(132)=1
          msub(135)=1
          msub(136)=1
          IF(iptl.EQ.1) ckin(3)=ptmdir
 
C...Set up for p/VMD * anomalous gamma.
        ELSEIF(mint(122).EQ.3) THEN
          mint(123)=3
          IF(mint(121).EQ.6) mint(123)=7
          msub(11)=1
          msub(12)=1
          msub(13)=1
          msub(28)=1
          msub(53)=1
          msub(68)=1
          IF(iptl.EQ.1) msub(95)=1
          IF(msel.EQ.2) THEN
            msub(91)=1
            msub(92)=1
            msub(93)=1
            msub(94)=1
          ENDIF
          IF(iptl.EQ.1) ckin(3)=0d0
 
C...Set up for DIS * p.
        ELSEIF(mint(122).EQ.4.AND.(iabs(mint(11)).GT.100.OR.
     &  iabs(mint(12)).GT.100)) THEN
          mint(123)=8
          IF(iptl.EQ.1) msub(99)=1
 
C...Set up for direct * direct gamma (switch off leptons).
        ELSEIF(mint(122).EQ.4) THEN
          mint(123)=0
          msub(137)=1
          msub(138)=1
          msub(139)=1
          msub(140)=1
          DO 110 ii=mdcy(22,2),mdcy(22,2)+mdcy(22,3)-1
            IF(iabs(kfdp(ii,1)).GE.10) mdme(ii,1)=min(0,mdme(ii,1))
  110     CONTINUE
          IF(iptl.EQ.1) ckin(3)=ptmdir
 
C...Set up for direct * anomalous gamma.
        ELSEIF(mint(122).EQ.5) THEN
          mint(123)=6
          msub(131)=1
          msub(132)=1
          msub(135)=1
          msub(136)=1
          IF(iptl.EQ.1) ckin(3)=ptmano
 
C...Set up for anomalous * anomalous gamma.
        ELSEIF(mint(122).EQ.6) THEN
          mint(123)=3
          msub(11)=1
          msub(12)=1
          msub(13)=1
          msub(28)=1
          msub(53)=1
          msub(68)=1
          IF(iptl.EQ.1) msub(95)=1
          IF(msel.EQ.2) THEN
            msub(91)=1
            msub(92)=1
            msub(93)=1
            msub(94)=1
          ENDIF
          IF(iptl.EQ.1) ckin(3)=0d0
        ENDIF
 
C...Set up for gamma* * gamma*; virtual photons = dir, VMD, anom.
        ELSEIF(mint(121).EQ.9.OR.mint(121).EQ.13) THEN
 
C...Set up for direct * direct gamma (switch off leptons).
        IF(mint(122).EQ.1) THEN
          mint(123)=0
          msub(137)=1
          msub(138)=1
          msub(139)=1
          msub(140)=1
          DO 120 ii=mdcy(22,2),mdcy(22,2)+mdcy(22,3)-1
            IF(iabs(kfdp(ii,1)).GE.10) mdme(ii,1)=min(0,mdme(ii,1))
  120     CONTINUE
          IF(iptl.EQ.1) ckin(3)=ptmdir
 
C...Set up for direct * VMD and VMD * direct gamma.
        ELSEIF(mint(122).EQ.2.OR.mint(122).EQ.4) THEN
          mint(123)=5
          msub(131)=1
          msub(132)=1
          msub(135)=1
          msub(136)=1
          IF(iptl.EQ.1) ckin(3)=ptmdir
 
C...Set up for direct * anomalous and anomalous * direct gamma.
        ELSEIF(mint(122).EQ.3.OR.mint(122).EQ.7) THEN
          mint(123)=6
          msub(131)=1
          msub(132)=1
          msub(135)=1
          msub(136)=1
          IF(iptl.EQ.1) ckin(3)=ptmano
 
C...Set up for VMD*VMD.
        ELSEIF(mint(122).EQ.5) THEN
          mint(123)=2
          msub(11)=1
          msub(12)=1
          msub(13)=1
          msub(28)=1
          msub(53)=1
          msub(68)=1
          IF(iptl.EQ.1) msub(95)=1
          IF(msel.EQ.2) THEN
            msub(91)=1
            msub(92)=1
            msub(93)=1
            msub(94)=1
          ENDIF
          IF(iptl.EQ.1) ckin(3)=0d0
 
C...Set up for VMD * anomalous and anomalous * VMD gamma.
        ELSEIF(mint(122).EQ.6.OR.mint(122).EQ.8) THEN
          mint(123)=7
          msub(11)=1
          msub(12)=1
          msub(13)=1
          msub(28)=1
          msub(53)=1
          msub(68)=1
          IF(iptl.EQ.1) msub(95)=1
          IF(msel.EQ.2) THEN
            msub(91)=1
            msub(92)=1
            msub(93)=1
            msub(94)=1
          ENDIF
          IF(iptl.EQ.1) ckin(3)=0d0
 
C...Set up for anomalous * anomalous gamma.
        ELSEIF(mint(122).EQ.9) THEN
          mint(123)=3
          msub(11)=1
          msub(12)=1
          msub(13)=1
          msub(28)=1
          msub(53)=1
          msub(68)=1
          IF(iptl.EQ.1) msub(95)=1
          IF(msel.EQ.2) THEN
            msub(91)=1
            msub(92)=1
            msub(93)=1
            msub(94)=1
          ENDIF
          IF(iptl.EQ.1) ckin(3)=0d0
 
C...Set up for DIS * VMD and VMD * DIS gamma.
        ELSEIF(mint(122).EQ.10.OR.mint(122).EQ.12) THEN
          mint(123)=8
          IF(iptl.EQ.1) msub(99)=1
 
C...Set up for DIS * anomalous and anomalous * DIS gamma.
        ELSEIF(mint(122).EQ.11.OR.mint(122).EQ.13) THEN
          mint(123)=9
          IF(iptl.EQ.1) msub(99)=1
        ENDIF
 
C...Set up for gamma* * p; virtual photons = dir, res.
        ELSEIF(mint(121).EQ.2) THEN
 
C...Set up for direct * p.
        IF(mint(122).EQ.1) THEN
          mint(123)=0
          msub(131)=1
          msub(132)=1
          msub(135)=1
          msub(136)=1
          IF(iptl.EQ.1) ckin(3)=ptmdir
 
C...Set up for resolved * p.
        ELSEIF(mint(122).EQ.2) THEN
          mint(123)=1
          msub(11)=1
          msub(12)=1
          msub(13)=1
          msub(28)=1
          msub(53)=1
          msub(68)=1
          IF(iptl.EQ.1) msub(95)=1
          IF(msel.EQ.2) THEN
            msub(91)=1
            msub(92)=1
            msub(93)=1
            msub(94)=1
          ENDIF
          IF(iptl.EQ.1) ckin(3)=0d0
        ENDIF
 
C...Set up for gamma* * gamma*; virtual photons = dir, res.
        ELSEIF(mint(121).EQ.4) THEN
 
C...Set up for direct * direct gamma (switch off leptons).
        IF(mint(122).EQ.1) THEN
          mint(123)=0
          msub(137)=1
          msub(138)=1
          msub(139)=1
          msub(140)=1
          DO 130 ii=mdcy(22,2),mdcy(22,2)+mdcy(22,3)-1
            IF(iabs(kfdp(ii,1)).GE.10) mdme(ii,1)=min(0,mdme(ii,1))
  130     CONTINUE
          IF(iptl.EQ.1) ckin(3)=ptmdir
 
C...Set up for direct * resolved and resolved * direct gamma.
        ELSEIF(mint(122).EQ.2.OR.mint(122).EQ.3) THEN
          mint(123)=5
          msub(131)=1
          msub(132)=1
          msub(135)=1
          msub(136)=1
          IF(iptl.EQ.1) ckin(3)=ptmdir
 
C...Set up for resolved * resolved gamma.
        ELSEIF(mint(122).EQ.4) THEN
          mint(123)=2
          msub(11)=1
          msub(12)=1
          msub(13)=1
          msub(28)=1
          msub(53)=1
          msub(68)=1
          IF(iptl.EQ.1) msub(95)=1
          IF(msel.EQ.2) THEN
            msub(91)=1
            msub(92)=1
            msub(93)=1
            msub(94)=1
          ENDIF
          IF(iptl.EQ.1) ckin(3)=0d0
        ENDIF
 
C...End of special set up for gamma-p and gamma-gamma.
        ENDIF
        ckin(1)=2d0*ckin(3)
      ENDIF
 
C...Flavour information for individual beams.
      DO 140 i=1,2
        mint(40+i)=1
        IF(mint(123).GE.1.AND.mint(10+i).EQ.22) mint(40+i)=2
        IF(iabs(mint(10+i)).GT.100) mint(40+i)=2
        mint(44+i)=mint(40+i)
        IF(mstp(11).GE.1.AND.(iabs(mint(10+i)).EQ.11.OR.
     &  iabs(mint(10+i)).EQ.13.OR.iabs(mint(10+i)).EQ.15)) mint(44+i)=3
  140 CONTINUE
 
C...If two real gammas, whereof one direct, pick the first.
C...For two virtual photons, keep requested order.
      IF(mint(11).EQ.22.AND.mint(12).EQ.22) THEN
        IF(mstp(14).LE.10.AND.mint(123).GE.4.AND.mint(123).LE.6) THEN
          mint(41)=1
          mint(45)=1
        ELSEIF(mstp(14).EQ.12.OR.mstp(14).EQ.13.OR.mstp(14).EQ.22.OR.
     &  mstp(14).EQ.26.OR.mstp(14).EQ.27) THEN
          mint(41)=1
          mint(45)=1
        ELSEIF(mstp(14).EQ.14.OR.mstp(14).EQ.17.OR.mstp(14).EQ.23.OR.
     &  mstp(14).EQ.28.OR.mstp(14).EQ.29) THEN
          mint(42)=1
          mint(46)=1
        ELSEIF((mstp(14).EQ.20.OR.mstp(14).EQ.30).AND.(mint(122).EQ.2
     &  .OR.mint(122).EQ.3.OR.mint(122).EQ.10.OR.mint(122).EQ.11)) THEN
          mint(41)=1
          mint(45)=1
        ELSEIF((mstp(14).EQ.20.OR.mstp(14).EQ.30).AND.(mint(122).EQ.4
     &  .OR.mint(122).EQ.7.OR.mint(122).EQ.12.OR.mint(122).EQ.13)) THEN
          mint(42)=1
          mint(46)=1
        ELSEIF(mstp(14).EQ.25.AND.mint(122).EQ.2) THEN
          mint(41)=1
          mint(45)=1
        ELSEIF(mstp(14).EQ.25.AND.mint(122).EQ.3) THEN
          mint(42)=1
          mint(46)=1
        ENDIF
      ELSEIF(mint(11).EQ.22.OR.mint(12).EQ.22) THEN
        IF(mstp(14).EQ.26.OR.mstp(14).EQ.28.OR.mint(122).EQ.4) THEN
          IF(mint(11).EQ.22) THEN
            mint(41)=1
            mint(45)=1
          ELSE
            mint(42)=1
            mint(46)=1
          ENDIF
        ENDIF
        IF(mint(123).GE.4.AND.mint(123).LE.7) CALL pyerrm(26,
     &  '(PYINPR:) unallowed MSTP(14) code for single photon')
      ENDIF
 
C...Flavour information on combination of incoming particles.
      mint(43)=2*mint(41)+mint(42)-2
      mint(44)=mint(43)
      IF(mint(123).LE.0) THEN
        IF(mint(11).EQ.22) mint(43)=mint(43)+2
        IF(mint(12).EQ.22) mint(43)=mint(43)+1
      ELSEIF(mint(123).LE.3) THEN
        IF(mint(11).EQ.22) mint(44)=mint(44)-2
        IF(mint(12).EQ.22) mint(44)=mint(44)-1
      ELSEIF(mint(11).EQ.22.AND.mint(12).EQ.22) THEN
        mint(43)=4
        mint(44)=1
      ENDIF
      mint(47)=2*min(2,mint(45))+min(2,mint(46))-2
      IF(min(mint(45),mint(46)).EQ.3) mint(47)=5
      IF(mint(45).EQ.1.AND.mint(46).EQ.3) mint(47)=6
      IF(mint(45).EQ.3.AND.mint(46).EQ.1) mint(47)=7
      mint(50)=0
      IF(mint(41).EQ.2.AND.mint(42).EQ.2.AND.mint(111).NE.12) mint(50)=1
      mint(107)=0
      mint(108)=0
      IF(mint(121).EQ.9.OR.mint(121).EQ.13) THEN
        IF((mint(122).GE.4.AND.mint(122).LE.6).OR.mint(122).EQ.12)
     &  mint(107)=2
        IF((mint(122).GE.7.AND.mint(122).LE.9).OR.mint(122).EQ.13)
     &  mint(107)=3
        IF(mint(122).EQ.10.OR.mint(122).EQ.11) mint(107)=4
        IF(mint(122).EQ.2.OR.mint(122).EQ.5.OR.mint(122).EQ.8.OR.
     &  mint(122).EQ.10) mint(108)=2
        IF(mint(122).EQ.3.OR.mint(122).EQ.6.OR.mint(122).EQ.9.OR.
     &  mint(122).EQ.11) mint(108)=3
        IF(mint(122).EQ.12.OR.mint(122).EQ.13) mint(108)=4
      ELSEIF(mint(121).EQ.4.AND.mstp(14).EQ.25) THEN
        IF(mint(122).GE.3) mint(107)=1
        IF(mint(122).EQ.2.OR.mint(122).EQ.4) mint(108)=1
      ELSEIF(mint(121).EQ.2) THEN
        IF(mint(122).EQ.2.AND.mint(11).EQ.22) mint(107)=1
        IF(mint(122).EQ.2.AND.mint(12).EQ.22) mint(108)=1
      ELSE
        IF(mint(11).EQ.22) THEN
          mint(107)=mint(123)
          IF(mint(123).GE.4) mint(107)=0
          IF(mint(123).EQ.7) mint(107)=2
          IF(mstp(14).EQ.26.OR.mstp(14).EQ.27) mint(107)=4
          IF(mstp(14).EQ.28) mint(107)=2
          IF(mstp(14).EQ.29) mint(107)=3
          IF(mstp(14).EQ.30.AND.mint(121).EQ.4.AND.mint(122).EQ.4)
     &    mint(107)=4
        ENDIF
        IF(mint(12).EQ.22) THEN
          mint(108)=mint(123)
          IF(mint(123).GE.4) mint(108)=mint(123)-3
          IF(mint(123).EQ.7) mint(108)=3
          IF(mstp(14).EQ.26) mint(108)=2
          IF(mstp(14).EQ.27) mint(108)=3
          IF(mstp(14).EQ.28.OR.mstp(14).EQ.29) mint(108)=4
          IF(mstp(14).EQ.30.AND.mint(121).EQ.4.AND.mint(122).EQ.4)
     &    mint(108)=4
        ENDIF
        IF(mint(11).EQ.22.AND.mint(12).EQ.22.AND.(mstp(14).EQ.14.OR.
     &  mstp(14).EQ.17.OR.mstp(14).EQ.18.OR.mstp(14).EQ.23)) THEN
          minttp=mint(107)
          mint(107)=mint(108)
          mint(108)=minttp
        ENDIF
      ENDIF
      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
 
C...Select default processes according to incoming beams
C...(already done for gamma-p and gamma-gamma with
C...MSTP(14) = 10, 20, 25 or 30).
      IF(mint(121).GT.1) THEN
      ELSEIF(msel.EQ.1.OR.msel.EQ.2) THEN
 
        IF(mint(43).EQ.1) THEN
C...Lepton + lepton -> gamma/Z0 or W.
          IF(mint(11)+mint(12).EQ.0) msub(1)=1
          IF(mint(11)+mint(12).NE.0) msub(2)=1
 
        ELSEIF(mint(43).LE.3.AND.mint(123).EQ.0.AND.
     &    (mint(11).EQ.22.OR.mint(12).EQ.22)) THEN
C...Unresolved photon + lepton: Compton scattering.
          msub(133)=1
          msub(134)=1
 
        ELSEIF((mint(123).EQ.8.OR.mint(123).EQ.9).AND.(mint(11).EQ.22
     &  .OR.mint(12).EQ.22)) THEN
C...DIS as pure gamma* + f -> f process.
          msub(99)=1
 
        ELSEIF(mint(43).LE.3) THEN
C...Lepton + hadron: deep inelastic scattering.
          msub(10)=1
 
        ELSEIF(mint(123).EQ.0.AND.mint(11).EQ.22.AND.
     &    mint(12).EQ.22) THEN
C...Two unresolved photons: fermion pair production,
C...exclude lepton pairs.
          DO 150 isub=137,140
            msub(isub)=1
  150     CONTINUE
          DO 160 ii=mdcy(22,2),mdcy(22,2)+mdcy(22,3)-1
            IF(iabs(kfdp(ii,1)).GE.10) mdme(ii,1)=min(0,mdme(ii,1))
  160     CONTINUE
          ptmdir=ptmrun
          IF(mstp(18).EQ.2) ptmdir=parp(15)
          IF(ckin(3).LT.ptmrun.OR.msel.EQ.2) ckin(3)=ptmdir
          ckin(1)=max(ckin(1),2d0*ckin(3))
 
        ELSEIF((mint(123).EQ.0.AND.(mint(11).EQ.22.OR.mint(12).EQ.22))
     &    .OR.(mint(123).GE.4.AND.mint(123).LE.6.AND.mint(11).EQ.22.AND.
     &    mint(12).EQ.22)) THEN
C...Unresolved photon + hadron: photon-parton scattering.
          DO 170 isub=131,136
            msub(isub)=1
  170     CONTINUE
 
        ELSEIF(msel.EQ.1) THEN
C...High-pT QCD processes:
          msub(11)=1
          msub(12)=1
          msub(13)=1
          msub(28)=1
          msub(53)=1
          msub(68)=1
          ptmn=ptmrun
          vint(154)=ptmn
          IF(ckin(3).LT.ptmn) msub(95)=1
          IF(msub(95).EQ.1.AND.mint(50).EQ.0) msub(95)=0
 
        ELSE
C...All QCD processes:
          msub(11)=1
          msub(12)=1
          msub(13)=1
          msub(28)=1
          msub(53)=1
          msub(68)=1
          msub(91)=1
          msub(92)=1
          msub(93)=1
          msub(94)=1
          msub(95)=1
        ENDIF
 
      ELSEIF(msel.GE.4.AND.msel.LE.8) THEN
C...Heavy quark production.
        msub(81)=1
        msub(82)=1
        msub(84)=1
        DO 180 j=1,min(8,mdcy(21,3))
          mdme(mdcy(21,2)+j-1,1)=0
  180   CONTINUE
        mdme(mdcy(21,2)+msel-1,1)=1
        msub(85)=1
        DO 190 j=1,min(12,mdcy(22,3))
          mdme(mdcy(22,2)+j-1,1)=0
  190   CONTINUE
        mdme(mdcy(22,2)+msel-1,1)=1
 
      ELSEIF(msel.EQ.10) THEN
C...Prompt photon production:
        msub(14)=1
        msub(18)=1
        msub(29)=1
 
      ELSEIF(msel.EQ.11) THEN
C...Z0/gamma* production:
        msub(1)=1
 
      ELSEIF(msel.EQ.12) THEN
C...W+/- production:
        msub(2)=1
 
      ELSEIF(msel.EQ.13) THEN
C...Z0 + jet:
        msub(15)=1
        msub(30)=1
 
      ELSEIF(msel.EQ.14) THEN
C...W+/- + jet:
        msub(16)=1
        msub(31)=1
 
      ELSEIF(msel.EQ.15) THEN
C...Z0 & W+/- pair production:
        msub(19)=1
        msub(20)=1
        msub(22)=1
        msub(23)=1
        msub(25)=1
 
      ELSEIF(msel.EQ.16) THEN
C...h0 production:
        msub(3)=1
        msub(102)=1
        msub(103)=1
        msub(123)=1
        msub(124)=1
 
      ELSEIF(msel.EQ.17) THEN
C...h0 & Z0 or W+/- pair production:
        msub(24)=1
        msub(26)=1
 
      ELSEIF(msel.EQ.18) THEN
C...h0 production; interesting processes in e+e-.
        msub(24)=1
        msub(103)=1
        msub(123)=1
        msub(124)=1
 
      ELSEIF(msel.EQ.19) THEN
C...h0, H0 and A0 production; interesting processes in e+e-.
        msub(24)=1
        msub(103)=1
        msub(123)=1
        msub(124)=1
        msub(153)=1
        msub(171)=1
        msub(173)=1
        msub(174)=1
        msub(158)=1
        msub(176)=1
        msub(178)=1
        msub(179)=1
 
      ELSEIF(msel.EQ.21) THEN
C...Z'0 production:
        msub(141)=1
 
      ELSEIF(msel.EQ.22) THEN
C...W'+/- production:
        msub(142)=1
 
      ELSEIF(msel.EQ.23) THEN
C...H+/- production:
        msub(143)=1
 
      ELSEIF(msel.EQ.24) THEN
C...R production:
        msub(144)=1
 
      ELSEIF(msel.EQ.25) THEN
C...LQ (leptoquark) production.
        msub(145)=1
        msub(162)=1
        msub(163)=1
        msub(164)=1
 
      ELSEIF(msel.GE.35.AND.msel.LE.38) THEN
C...Production of one heavy quark (W exchange):
        msub(83)=1
        DO 200 j=1,min(8,mdcy(21,3))
          mdme(mdcy(21,2)+j-1,1)=0
  200   CONTINUE
        mdme(mdcy(21,2)+msel-31,1)=1
 
CMRENNA++Define SUSY alternatives.
      ELSEIF(msel.EQ.39) THEN
C...Turn on all SUSY processes.
        IF(mint(43).EQ.4) THEN
C...Hadron-hadron processes.
          DO 210 i=201,301
            IF(iset(i).GE.0) msub(i)=1
  210     CONTINUE
        ELSEIF(mint(43).EQ.1) THEN
C...Lepton-lepton processes: QED production of squarks.
          DO 220 i=201,214
            msub(i)=1
  220     CONTINUE
          msub(210)=0
          msub(211)=0
          msub(212)=0
          DO 230 i=216,228
            msub(i)=1
  230     CONTINUE
          DO 240 i=261,263
            msub(i)=1
  240     CONTINUE
          msub(277)=1
          msub(278)=1
        ENDIF
 
      ELSEIF(msel.EQ.40) THEN
C...Gluinos and squarks.
        IF(mint(43).EQ.4) THEN
          msub(243)=1
          msub(244)=1
          msub(258)=1
          msub(259)=1
          msub(261)=1
          msub(262)=1
          msub(264)=1
          msub(265)=1
          DO 250 i=271,296
            msub(i)=1
  250     CONTINUE
        ELSEIF(mint(43).EQ.1) THEN
          msub(277)=1
          msub(278)=1
        ENDIF
 
      ELSEIF(msel.EQ.41) THEN
C...Stop production.
        msub(261)=1
        msub(262)=1
        msub(263)=1
        IF(mint(43).EQ.4) THEN
          msub(264)=1
          msub(265)=1
        ENDIF
 
      ELSEIF(msel.EQ.42) THEN
C...Slepton production.
        DO 260 i=201,214
          msub(i)=1
  260   CONTINUE
        IF(mint(43).NE.4) THEN
          msub(210)=0
          msub(211)=0
          msub(212)=0
        ENDIF
 
      ELSEIF(msel.EQ.43) THEN
C...Neutralino/Chargino + Gluino/Squark.
        IF(mint(43).EQ.4) THEN
          DO 270 i=237,242
            msub(i)=1
  270     CONTINUE
          DO 280 i=246,254
            msub(i)=1
  280     CONTINUE
          msub(256)=1
        ENDIF
 
      ELSEIF(msel.EQ.44) THEN
C...Neutralino/Chargino pair production.
        IF(mint(43).EQ.4) THEN
          DO 290 i=216,236
            msub(i)=1
  290     CONTINUE
        ELSEIF(mint(43).EQ.1) THEN
          DO 300 i=216,228
            msub(i)=1
  300     CONTINUE
        ENDIF
 
      ELSEIF(msel.EQ.45) THEN
C...Sbottom production.
        msub(287)=1
        msub(288)=1
        IF(mint(43).EQ.4) THEN
          DO 310 i=281,296
            msub(i)=1
  310     CONTINUE
        ENDIF
 
      ELSEIF(msel.EQ.50) THEN
C...Pair production of technipions and gauge bosons.
        DO 320 i=361,368
          msub(i)=1
  320   CONTINUE
        IF(mint(43).EQ.4) THEN
          DO 330 i=370,377
            msub(i)=1
  330     CONTINUE
        ENDIF
 
      ELSEIF(msel.EQ.51) THEN
C...QCD 2 -> 2 processes with compositeness/technicolor modifications.
        DO 340 i=381,386
          msub(i)=1
  340   CONTINUE
 
      ELSEIF(msel.EQ.61) THEN
C...Charmonium production in colour octet model, with recoiling parton.
        DO 342 i=421,439
          msub(i)=1
 342   CONTINUE
 
      ELSEIF(msel.EQ.62) THEN
C...Bottomonium production in colour octet model, with recoiling parton.
        DO 344 i=461,479
          msub(i)=1
 344   CONTINUE
 
      ELSEIF(msel.EQ.63) THEN
C...Charmonium and bottomonium production in colour octet model.
        DO 346 i=421,439
          msub(i)=1
          msub(i+40)=1
 346   CONTINUE
      ENDIF
 
C...Find heaviest new quark flavour allowed in processes 81-84.
      kflqm=1
      DO 350 i=1,min(8,mdcy(21,3))
        idc=i+mdcy(21,2)-1
        IF(mdme(idc,1).LE.0) goto 350
        kflqm=i
  350 CONTINUE
      IF(mstp(7).GE.1.AND.mstp(7).LE.8.AND.(msel.LE.3.OR.msel.GE.9))
     &kflqm=mstp(7)
      mint(55)=kflqm
      kfpr(81,1)=kflqm
      kfpr(81,2)=kflqm
      kfpr(82,1)=kflqm
      kfpr(82,2)=kflqm
      kfpr(83,1)=kflqm
      kfpr(84,1)=kflqm
      kfpr(84,2)=kflqm
 
C...Find heaviest new fermion flavour allowed in process 85.
      kflfm=1
      DO 360 i=1,min(12,mdcy(22,3))
        idc=i+mdcy(22,2)-1
        IF(mdme(idc,1).LE.0) goto 360
        kflfm=kfdp(idc,1)
  360 CONTINUE
      IF(((mstp(7).GE.1.AND.mstp(7).LE.8).OR.(mstp(7).GE.11.AND.
     &mstp(7).LE.18)).AND.(msel.LE.3.OR.msel.GE.9)) kflfm=mstp(7)
      mint(56)=kflfm
      kfpr(85,1)=kflfm
      kfpr(85,2)=kflfm
 
C...Import relevant information on external user processes.
      IF(mint(111).GE.11) THEN
        ipypr=0
        DO 390 iup=1,nprup
C...Find next empty PYTHIA process number slot and enable it.
  370     ipypr=ipypr+1
          IF(ipypr.GT.500) CALL pyerrm(26,
     &    '(PYINPR.) no more empty slots for user processes')
          IF(iset(ipypr).GE.0.AND.iset(ipypr).LE.9) goto 370
          IF(ipypr.GE.91.AND.ipypr.LE.100) goto 370
          iset(ipypr)=11
C...Overwrite KFPR with references back to process number and ID.
          kfpr(ipypr,1)=iup
          kfpr(ipypr,2)=lprup(iup)
C...Process title.
          WRITE(chipr,'(I10)') lprup(iup)
          ichin=1
          DO 380 ich=1,9
            IF(chipr(ich:ich).EQ.' ') ichin=ich+1
  380     CONTINUE
          proc(ipypr)='User process '//chipr(ichin:10)//' '
C...Switch on process.
          msub(ipypr)=1
  390   CONTINUE
      ENDIF
 
      RETURN
      END