PROGRAM MATCHEX

C...Double precision and integer declarations
      IMPLICIT DOUBLE PRECISION(A-H, O-Z)
      IMPLICIT INTEGER(I-N)
      INTEGER PYK,PYCHGE,PYCOMP
C...PYTHIA Commonblocks used
      COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
      COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
      COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
      COMMON/PYINT2/ISET(500),KFPR(500,2),COEF(500,20),ICOL(40,4,2)
      COMMON/PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200)
C...User process initialization commonblock
      PARAMETER (MAXPUP=100)
      INTEGER IDBMUP,PDFGUP,PDFSUP,IDWTUP,NPRUP,LPRUP
      COMMON/HEPRUP/IDBMUP(2),EBMUP(2),PDFGUP(2),PDFSUP(2),
     &IDWTUP,NPRUP,XSECUP(MAXPUP),XERRUP(MAXPUP),XMAXUP(MAXPUP),
     &LPRUP(MAXPUP)

C...Select internal processes to run
      MSEL=0
c   gb->tH+
      MSUB(161)=1
c   gg->tbH+
      MSUB(401)=1
c   qqbar->tbH+
      MSUB(402)=1

C...Set parameters
c   tanb
      PARU(141) = 30d0
c   H+ mass
      PMAS(37,1) = 250d0

C...Set factorization scale
      MSTP(32)=12
      MSTP(39)=8
c   factorization scale
      PARP(193) = 0.25d0*((PYMASS(37)+PYMASS(6))/2d0)**2
c   renormalization scale
      PARP(194) = PYMASS(37)**2

C...Initial state radiation onshell
      MSTP(63)=0

C...BW shape of H+ and top masses in double counting term
      MSTP(190)=1

C   Set number of tries to find highest xsect with 2->3 variables
      MSTP(129)=1000

C...CM Energy (LHC)
      ECM=1.4d4

C...Set incoming beams to protons
      IDBMUP(1)=2212
      IDBMUP(2)=2212
C...Set beam energies
      EBMUP(1)=0.5*ECM
      EBMUP(2)=0.5*ECM

C...Initialize PYTHIA
      CALL PYINIT('USER',' ',' ',0d0)

C...Start PYTHIA loop (generate events)
      NPOS=0
      NNEG=0
      NITER=1000
      DO 20 J=1,NITER

C     Generate an event
        CALL PYEVNT

C     Determine negative or positive event
         IF(MSTI(1).EQ.161.OR.MSTI(1).EQ.401.OR.MSTI(1).EQ.402) THEN
           NPOS=NPOS+1
         ELSE IF(KFPR(MSTI(1),2).EQ.LPRUP(1))THEN
c        See the PYTHIA 6.2 manual, p.222
           NNEG=NNEG+1
         END IF

 20   CONTINUE

      CALL PYSTAT(1)

      WRITE(*,*) 'Number of pos events:',NPOS
      WRITE(*,*) 'Number of neg events:',NNEG

      END