PROGRAM HWIGPR
   

C---COMMON BLOCKS ARE INCLUDED AS FILE HERWIG65.INC
      INCLUDE 'HERWIG65.INC'
      INTEGER N 
      EXTERNAL HWUDAT
C---MAX NUMBER OF EVENTS THIS RUN
      MAXEV= 1000
C---BEAM PARTICLES
      PART1='P'
      PART2='P'
      open (unit = 1, file = "output.dat")

  
C---BEAM MOMENTA 
      PBEAM1=7000.
      PBEAM2=7000.
C---PROCESS
      IPROC=1353
C---INITIALISE OTHER COMMON BLOCKS
      CALL HWIGIN
C---USER CAN RESET PARAMETERS AT THIS POINT, OTHERWISE DEFAULT
C   VALUES IN HWIGIN WILL BE USED.
      NOWGT = .TRUE.
      PRVTX = .FALSE.
      MAXPR = 5
C      PTMIN=100.
      PTMIN = 5
C---COMPUTE PARAMETER-DEPENDENT CONSTANTS
      CALL HWUINC
C---CALL HWUSTA TO MAKE ANY PARTICLE STABLE
      CALL HWUSTA('PI0     ')
C---USER'S INITIAL CALCULATIONS
      CALL HWABEG
C---INITIALISE ELEMENTARY PROCESS
      CALL HWEINI


C---LOOP OVER EVENTS
      DO 100 N=1,MAXEV
C---INITIALISE EVENT
      CALL HWUINE
C---GENERATE HARD SUBPROCESS
      CALL HWEPRO
C---GENERATE PARTON CASCADES
      CALL HWBGEN
C---FINISH EVENT
      CALL HWUFNE
C---USER'S EVENT ANALYSIS

       CALL HWRAP
  100 CONTINUE
C---TERMINATE ELEMENTARY PROCESS
      CALL HWEFIN
C---USER'S TERMINAL CALCULATIONS
      CALL HWAEND
      close(1)

      STOP
      END

      SUBROUTINE HWABEG
C     USER'S ROUTINE FOR INITIALIZATION
C----------------------------------------------------------------------
      END
C----------------------------------------------------------------------
      SUBROUTINE HWAEND
C     USER'S ROUTINE FOR TERMINAL CALCULATIONS, HISTOGRAM OUTPUT, ETC
C----------------------------------------------------------------------
      END
C----------------------------------------------------------------------



      SUBROUTINE HWRAP
C     USER'S ROUTINE TO CALCULATE THE RAPIDITY DISTRIBUTION OF THE Z
C     ZLOW    : MINIMUM VALUE OF Z FOR CUT
C     PTLCUT  : CUT ON LEPTON MOMENTUM
C     ETACUT  : CUT ON LEPTON RAPIDITY
C----------------------------------------------------------------------
      INCLUDE 'HERWIG65.INC'
      INTEGER IST, I
      LOGICAL ZCUT, PTTEST, PLUS, MINUS
      DOUBLE PRECISION ZLOW, PX, PY, PT, PTLP,PTLM, ETAL, 
     &  PTLCUT, MUPLUS, MUMINUS, E, PZ, ETA,
     &  ETAFNL, ETACUT 

      
      ETACUT = 50
      ZLOW = 40
      PTLCUT = 5
      IF (IERROR.NE.0) RETURN

      ZCUT = .FALSE.
      PLUS = .FALSE.
      MINUS =.FALSE. 
         DO 30 I=1, NHEP
            IST = ISTHEP(I)
C MAKE CUT ON Z
            IF((IST .EQ. 120) .AND. (IDHEP(I) .EQ. 23)) THEN
              IF(PHEP(5,I) .GT. ZLOW)  THEN
                E = PHEP(4,I)
                PZ= PHEP(3,I)
                ETAFNL =  0.5*LOG((E+PZ)/(E-PZ))
                ZCUT = .TRUE.
              ENDIF
            ENDIF
C MAKES CUT ON PT and Y FOR MU-MINUS 
            IF ((IST .EQ. 190) .AND.(IDHEP(I) .EQ. 13) ) THEN
               PX = PHEP(1,I)
               PY = PHEP(2,I)
               E = PHEP(4,I)
               PZ = PHEP(3,I)
               ETA = 0.5*LOG((E+PZ)/(E-PZ))
               PTLM = SQRT(PX*PX + PY*PY)
               IF((PTLM .GT. PTLCUT) .AND. (ABS(ETA) .LT.ETACUT )) THEN
                 PLUS = .TRUE.
               ENDIF
            ENDIF
C MAKES CUT ON PT and Y FOR MU-PLUS
            IF ((IST .EQ. 190) .AND.(IDHEP(I) .EQ. -13) ) THEN
               PX = PHEP(1,I)
               PY = PHEP(2,I)
               E = PHEP(4,I)
               PZ = PHEP(3,I)
               ETA = 0.5*LOG((E+PZ)/(E-PZ))
               PTLP = SQRT(PX*PX + PY*PY)
               IF((PTLP .GT. PTLCUT) .AND. (ABS(ETA) .LT.ETACUT)) THEN
                  MINUS = .TRUE.
               ENDIF
            ENDIF
   30    CONTINUE
         IF(ZCUT) THEN
C ONYLY WRITE EVENTS THAT PASS LEPTON TEST
            IF(PLUS .AND. MINUS) THEN
              WRITE(1,*), ETAFNL
            ENDIF
          ENDIF
      END
C--------------------------------------------------------------------------