PROGRAM MATCHEX

C...Double precision and integer declarations
      IMPLICIT DOUBLE PRECISION(A-H, O-Z)
      IMPLICIT INTEGER(I-N)
C...PYTHIA Commonblocks used
      COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
      COMMON/PYINT2/ISET(500),KFPR(500,2),COEF(500,20),ICOL(40,4,2)
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
      CALL PYGIVE('MSEL=0')
c   gb->tH+
      CALL PYGIVE('MSUB(161)=1')
c   gg->tbH+
      CALL PYGIVE('MSUB(401)=1')
c   qqbar->tbH+
      CALL PYGIVE('MSUB(402)=1')

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

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

C...Turn off multiple interactions (to run faster)
      CALL PYGIVE('MSTP(81)=0')

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