C Driver, to be modified by the user (keeping the read statements).
C Works only with Les Houches interface (so enter negative process numbers)
      PROGRAM HWIGPR
      INCLUDE 'HERWIG65.INC'
      INTEGER N,NSTEP,I,JPR0,JPR
C QQIN IS THE EVENT FILE
      CHARACTER*50 QQIN
      COMMON/VVJIN/QQIN
      REAL*8 TMPLAM
      INTEGER IPDF
      REAL*8 flav1,flav2,HWRGEN
      EXTERNAL HWRGEN
c
      WRITE(*,*)'Enter filename for events'
      READ(*,*)QQIN
      WRITE(*,*)'Enter maximum number of events to generate'
      WRITE(*,*)'MUST coincide with the number of events on tape'
      READ(*,*)MAXEV
      WRITE(*,*)'Enter 0 to use Herwig default PDFs'
      WRITE(*,*)'      1 to use PDFs from library'
      READ(*,*)IPDF
C OUTPUT THE RESULTS AFTER ANY NSTEP EVENTS
      NSTEP=20000
C---BEAM PARTICLES
      WRITE(*,*)'Enter colliding particles (PART1, PART2)'
      READ(*,*)PART1,PART2
C---BEAM MOMENTA
      WRITE(*,*)'Enter beam momenta (PBEAM1, PBEAM2)'
      READ(*,*)PBEAM1,PBEAM2
C---PROCESS
      WRITE(*,*)'Enter process number (IPROC) (IPROC<0 for Les Houches)'
      READ(*,*)IPROC
C---INITIALISE OTHER COMMON BLOCKS
      CALL HWIGIN
C---USER CAN RESET PARAMETERS AT
C   THIS POINT, OTHERWISE DEFAULT
C   VALUES IN HWIGIN WILL BE USED.
C
C LHSOFT=.FALSE. SWITCHES UNDERLYING EVENT OFF
      LHSOFT=IPROC.GT.-10000
      WRITE(*,*)
      IF (LHSOFT) THEN
         WRITE(*,*)'Underlying event WILL be generated'
      ELSE
         WRITE(*,*)'Underlying event WILL NOT be generated'
      ENDIF
      WRITE(*,*)
C
C Select W/Z boson decay modes
      MODBOS(1)=4
C      MODBOS(2)=3
C
      IF(IPDF.EQ.1)THEN
         DO I=1,2
            WRITE(*,*)'   Incoming particle # ',I
            WRITE(*,*)'Enter PDF group name (AUTPDF)'
            READ(*,*)AUTPDF(I)
            WRITE(*,*)'Enter PDF set number (MODPDF)'
            READ(*,*)MODPDF(I)
         ENDDO
      ENDIF
      WRITE(*,*)'Enter Lambda_QCD, <0 for Herwig default'
      READ(*,*)TMPLAM
      IF(TMPLAM.GE.0.D0)QCDLAM=TMPLAM
C---CHECK PROCESS CODE
      JPR0=MOD(ABS(IPROC),10000)
      JPR=JPR0/100
      IF (JPR.EQ.13.OR.JPR.EQ.14) THEN
         IF(JPR0.EQ.1350.OR.JPR0.EQ.1396.OR.JPR0.EQ.1399)THEN
           WRITE(*,*)'Enter M_GAMMA*(INF), M_GAMMA*(SUP)'
           READ(*,*)EMMIN,EMMAX
         ELSEIF(JPR0.EQ.1397)THEN
           WRITE(*,*)'Enter Z0 mass, width, and GammaX'
           READ(*,*)RMASS(200),GAMZ,GAMMAX
         ELSEIF(JPR0.EQ.1497)THEN
           WRITE(*,*)'Enter W+ mass, width, and GammaX'
           READ(*,*)RMASS(198),GAMW,GAMMAX
           RMASS(199)=RMASS(198)
         ELSEIF(JPR0.EQ.1498)THEN
           WRITE(*,*)'Enter W- mass, width, and GammaX'
           READ(*,*)RMASS(199),GAMW,GAMMAX
           RMASS(198)=RMASS(199)
         ELSEIF(JPR0.EQ.1499)THEN
           WRITE(*,*)'Enter W mass, width, and GammaX'
           READ(*,*)RMASS(199),GAMW,GAMMAX
           RMASS(198)=RMASS(199)
         ENDIF
      ELSEIF (JPR.EQ.28) THEN
         WRITE(*,*)'Enter W+, W-, Z0 masses'
         READ(*,*)RMASS(198),RMASS(199),RMASS(200)
      ELSEIF (JPR.EQ.16) THEN
         WRITE(*,*)'Enter Higgs boson and top masses'
         READ(*,*)RMASS(201),RMASS(6)
      ELSEIF (JPR.EQ.17) THEN
         IF(ABS(IPROC).EQ.1705.OR.ABS(IPROC).EQ.11705)THEN
           WRITE(*,*)'Enter bottom mass'
           READ(*,*)RMASS(5)
         ELSEIF(ABS(IPROC).EQ.1706.OR.ABS(IPROC).EQ.11706)THEN
           WRITE(*,*)'Enter top mass'
           READ(*,*)RMASS(6)
         ENDIF
      ELSE
         WRITE(*,*)'   Bad process code IPROC =',IPROC
         CALL HWWARN('HWIGPR',502,*999)
      ENDIF
      WRITE(*,*)'Enter quark (d,u,s,c,b) and gluon masses'
      READ(*,*)RMASS(1),RMASS(2),RMASS(3),
     #         RMASS(4),RMASS(5),RMASS(13)
      DO I=1,5
         RMASS(I+6)=RMASS(I)
      ENDDO
C NO SOFT ME CORRECTIONS (ALREADY INCLUDED IN MC@NLO)
      IF(IPROC.LT.0)THEN
        SOFTME=.FALSE.
      ELSE
        SOFTME=.TRUE.
      ENDIF
C HERWIG RUNS IN WEIGHTED MODE
      NOWGT=.FALSE.
C NEGATIVE WEIGHTS ALLOWED
      NEGWTS=.TRUE.
      MAXPR=2
      MAXER=MAXEV/100
      LRSUD=0
      LWSUD=77
C IN THE CASE HERWIG PDFS ARE USED, ADOPT MRST
      NSTRU=8
      PRVTX=.FALSE.
      PTMIN=0.5
      NRN(1)=1973774260
      NRN(2)=1099242306
C THE FOLLOWING SHOULD BE IRRELEVANT (ONLY USED IN UNWEIGHTED MODE)
      WGTMAX=1.000001D0
      AVABW=1.000001D0
C FOR TOP PRODUCTION (HARMLESS ELSEWHERE)
      RLTIM(6)=1.D-23
      RLTIM(12)=1.D-23
C---B FRAGMENTATION PARAMETERS (FOR B PRODUCTION ONLY)
      IF(ABS(IPROC).EQ.1705.OR.ABS(IPROC).EQ.11705)THEN
        PSPLT(2)=0.5
      ENDIF
C---COMPUTE PARAMETER-DEPENDENT CONSTANTS
      CALL HWUINC
C---CALL HWUSTA TO MAKE ANY PARTICLE STABLE
      CALL HWUSTA('PI0     ')
C---USER'S INITIAL CALCULATIONS
      write(6,*) 'CALLING INITIAL HWABEG'
      CALL HWABEG
C---INITIALISE ELEMENTARY PROCESS
      CALL HWEINI
C---LOOP OVER EVENTS
      DO 100 N=1,MAXEV

c---DSW. Set decay modes on event by event basis :
         flav1=HWRGEN(0)
         flav2=HWRGEN(1)
         if (flav1.GE.(2.0/3.0)) then
            MODBOS(1)=4
         else if (flav1.GE.(1.0/3.0)) then
            MODBOS(1)=3
         else 
            MODBOS(1)=2
         endif
         if (flav2.GE.(2.0/3.0)) then
            MODBOS(2)=4
         else if (flav2.GE.(1.0/3.0)) then
            MODBOS(2)=3
         else 
            MODBOS(2)=2
         endif
c        write(6,*) 'MODBOS(1),(2) = ',MODBOS(1),' , ',MODBOS(2)

C---INITIALISE EVENT
         CALL HWUINE
C---GENERATE HARD SUBPROCESS
         CALL HWEPRO
C---GENERATE PARTON CASCADES
         CALL HWBGEN
C---DO HEAVY OBJECT DECAYS
         CALL HWDHOB
C---DO CLUSTER FORMATION
         CALL HWCFOR
C---DO CLUSTER DECAYS
         CALL HWCDEC
C---DO UNSTABLE PARTICLE DECAYS
         CALL HWDHAD
C---DO HEAVY FLAVOUR HADRON DECAYS
         CALL HWDHVY
C---ADD SOFT UNDERLYING EVENT IF NEEDED
         CALL HWMEVT
C---FINISH EVENT
         CALL HWUFNE
C---USER'S EVENT ANALYSIS
         CALL HWANAL
         IF(MOD(NEVHEP,NSTEP).EQ.0) THEN
            WRITE(*,*)'# of events processed=',NEVHEP
            CALL HWAEND
         ENDIF
  100 CONTINUE
C---TERMINATE ELEMENTARY PROCESS
      CALL HWEFIN
C---USER'S TERMINAL CALCULATIONS
      WRITE(*,*)'# of events processed=',NEVHEP
      CALL HWAEND
 999  STOP
      END