C
      SUBROUTINE HSUSER(ICALL,X,Y,Q2)
      IMPLICIT DOUBLE PRECISION (A-H,M,O-Z)
      COMMON /HSOPTN/ INT2(5),INT3(15),ISAM2(5),ISAM3(15),
     *                IOPLOT,IPRINT,ICUT
      COMMON /HSNUME/ SIGTOT,SIGTRR,SIGG(20),SIGGRR(20),NEVENT,NEVE(20)
      COMMON /HSKNST/ PI,ALPHA,ALP1PI,ALP2PI,ALP4PI,E,GF,SXNORM,SX1NRM
      COMMON /HSELAB/ SP,EELE,PELE,EPRO,PPRO
      COMMON /HSPARM/ POLARI,LLEPT,LQUA
      COMMON /HSCUTS/ XMIN,XMAX,Q2MIN,Q2MAX,YMIN,YMAX,WMIN,GMIN
      COMMON /HSIKP/  S,T,U,SS,TS,US,DKP,DKPS,DKQ,DKQS
      COMMON /HSKPXY/ XX,YY
      COMMON /HSCHNN/ ICHNN
      PARAMETER (NMXHEP=2000)
      COMMON /HEPEVT/ NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
     &                JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),
     &                PHEP(5,NMXHEP),VHKK(4,NMXHEP)
C...HBook declarations
      INTEGER NWPAWC,CHANNELS,LIMITS,CHANNELCOUNTER
      REAL*4 HMEMORY
      REAL*4 HI,CHANNELSUM
      REAL*8 HINT,HBINW
      PARAMETER (NWPAWC=100000)
      PARAMETER (CHANNELS=50,LIMITS=CHANNELS+1)
      COMMON/PAWC/HMEMORY(NWPAWC)
      REAL*4 XBINS(LIMITS)
      INTEGER NID(1:200)
      DATA NID/200*0/
      REAL*8 BINW(1:200)
      DATA BINW/200*0./
      COMMON /HBOOKC/NID,BINW

      GOTO (100,200,300) ICALL

C...Initialization call 
 100  CONTINUE
C...INITIALIZE HBOOK
      CALL HLIMIT(NWPAWC)
C...Initialize limits for histograms and cuts
      XMINH=1D-5
      XMAXH=1D0
      YMINH=0D0
      YMAXH=1D0
      Q2MINH=0.1D0
      Q2MAXH=1D4
      EGMINH=0D0
      EGMAXH=35D0
      THETGN=0D0
      THETGX=5D0
      CTHGMX=DCOS(PI/180D0*THETGN)
      CTHGMN=DCOS(PI/180D0*THETGX)
      XMINL=DLOG10(XMINH)
      XMAXL=DLOG10(XMAXH)
      Q2MINL=DLOG10(Q2MINH)
      Q2MAXL=DLOG10(Q2MAXH)
C...Initialize and counters
      EGMIN=1D0
      THGMND=3D0
      THGMXD=177D0
      THGMIN=THGMND/180D0*PI
      THGMAX=THGMXD/180D0*PI
      PTMIN=0.5D0
      NTOTAL=0
      NINCUT=0
      NREJEC=0
C...Book histograms
      CALL HBOOK1(11,'<d[S]>/<d log x>',
     . CHANNELS,SNGL(XMINL),SNGL(XMAXL),0.)
      BINW(11)=HBINW(11)
      CALL HBOOK1(12,'<d[S]>/<d y>',
     . CHANNELS,SNGL(YMINH),SNGL(YMAXH),0.)
      BINW(12)=HBINW(12)
      CALL HBOOK1(13,'<d[S]>/<d log Q2>',
     . CHANNELS,SNGL(Q2MINL),SNGL(Q2MAXL),0.)
      BINW(13)=HBINW(13)
      CALL HBOOK1(21,'<d[S]>/<d Egamma>',
     . CHANNELS,SNGL(EGMINH),SNGL(EGMAXH),0.)
      BINW(21)=HBINW(21)
      CALL HBOOK1(22,'<d[S]>/<d cos theta gamma>',
     . CHANNELS,SNGL(CTHGMN),SNGL(CTHGMX),0.)
      BINW(22)=HBINW(22)
C...Fill header record of common HEPEVT
C     CALL HSHEAD
      NEVMOD=10000
      WRITE(*,*) ' START EVENT GENERATION '
      RETURN

C...Event scoring
 200  CONTINUE
      DSIGHB=SIGTOT/DFLOAT(NEVENT)
      DSIGH1=DSIGHB/BINW(11)
      DSIGH2=DSIGHB/BINW(12)
      DSIGH3=DSIGHB/BINW(13)
      DSIGH4=DSIGHB/BINW(21)
      DSIGH5=DSIGHB/BINW(22)
      CALL HFF1(11,NID(11),SNGL(DLOG10(XX)),SNGL(DSIGH1))
      CALL HFF1(12,NID(12),SNGL(YY),SNGL(DSIGH2))
      CALL HFF1(13,NID(13),SNGL(DLOG10(Q2)),SNGL(DSIGH3))
      ETG=PHEP(4,3)
      THG=0D0
      PTG=DABS(PHEP(1,3)**2+PHEP(2,3)**2)
      PTG=DSQRT(PTG)
      IF (PTG.LE.DABS(PHEP(3,3)).AND.PHEP(3,3).NE.0D0)
     &THG=DATAN(PTG/PHEP(3,3))
      CTHG=DCOS(THG)
      CALL HFF1(21,NID(21),SNGL(ETG),SNGL(DSIGH4))
      CALL HFF1(22,NID(22),SNGL(CTHG),SNGL(DSIGH4))
C...Print event kinematics
      IF(MOD(NEVHEP,NEVMOD).EQ.0) THEN
        WRITE(6,211) NEVHEP,NHEP,X,Y,Q2
        DO 220 I=1,NHEP
        WRITE(6,212) IDHEP(I),PHEP(1,I),PHEP(2,I),PHEP(3,I),PHEP(4,I)
 220    CONTINUE
        WRITE(6,213)  S,T,U,  SS,TS,US,  DKP,DKPS,DKQ,DKQS
 211    FORMAT(/,' NEVHEP = ',I6,/,' NHEP = ',I3,/,
     &  '  X = ',E12.5,'   Y = ',E12.5,/,
     &  ' Q2 = ',E12.5)
 212    FORMAT(I6,4D16.9)
 213    FORMAT(3D16.9/3D16.9/4D16.9)
      ENDIF
C...Count radiative events inside cuts
      NTOTAL=NTOTAL+1
      IF (ICHNN.GT.5) THEN
        EGAMMA=PHEP(4,3)
        PTGAMM=DSQRT(PHEP(1,3)**2+PHEP(2,3)**2)
        THGAMM=DASIN(PTGAMM/EGAMMA)
        IF (THGAMM.GT.THGMIN.AND.THGAMM.LT.THGMAX.AND.
     &      EGAMMA.GT.EGMIN.AND.PTGAMM.GT.PTMIN) NINCUT=NINCUT+1
      ENDIF
      RETURN

C...Final Call
 300  CONTINUE
C...Save histograms
      CALL HRPUT(0,'HERACLES.PAW','N')
C...Terminating record on HEPEVT
C     CALL HSTERM
C...Print final statistics of accepted events
      SIGICT=SIGTOT*DFLOAT(NINCUT)/DFLOAT(NTOTAL)
      WRITE(6,3001) NTOTAL,SIGTOT,EGMIN,THGMND,THGMXD,NINCUT,SIGICT
 3001 FORMAT(/,' Number of events   = ',I8,/
     &        ,' total cross section   = ',E15.5,//
     &        ,' Radiative events with E > ',F8.2,' GeV',/
     &        ,' and ',F8.2,' < theta < ',F8.2,/ 
     &        ,' events inside cuts = ',I8,/
     &        ,' cross section in cuts = ',E15.5)
      RETURN
      END
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C---FILL HEADER RECORD OF COMMON HEPEVT
C
      SUBROUTINE HSHEAD
      IMPLICIT DOUBLE PRECISION (A-H,M,O-Z)
C---------------------------------------------------------------------
      COMMON /HSOPTN/ INT2(5),INT3(15),ISAM2(5),ISAM3(15),
     *                IOPLOT,IPRINT,ICUT
      COMMON /HSUNTS/ LUNTES,LUNDAT,LUNIN,LUNOUT,LUNRND,LUNPD6,LUNPD7
      COMMON /VGASIO/ NINP,NOUTP
      COMMON /HSVGLP/ NPOVEG,NUMINT,NPHYP
      COMMON /HSRDIO/ ISDINP,ISDOUT
      COMMON /HSPARL/ LPAR(20),LPARIN(12),IPART
      COMMON /HSNUME/ SIGTOT,SIGTRR,SIGG(20),SIGGRR(20),NEVENT,NEVE(20)
      COMMON /HSELAB/ SP,EELE,PELE,EPRO,PPRO
      COMMON /HSPARM/ POLARI,LLEPT,LQUA
      COMMON /HSGSW/  SW,CW,SW2,CW2
     *              ,MW,MZ,MH,ME,MMY,MTAU,MU,MD,MS,MC,MB,MT
     *              ,MW2,MZ2,MH2,ME2,MMY2,MTAU2,MU2,MD2,MS2,MC2,MB2,MT2
      COMMON /PYSTFUC/ PYSTOP,PYSLAM,NPYMOD,NPYMAX,NPYMIN
      REAL*4           PYSTOP,PYSLAM
      COMMON /HSCUTS/ XMIN,XMAX,Q2MIN,Q2MAX,YMIN,YMAX,WMIN,GMIN
      COMMON /HSIRCT/ DELEPS,DELTA,EGMIN,IOPEGM
      PARAMETER (NMXHEP=2000)
      COMMON /HEPEVT/ NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
     &                JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),
     &                PHEP(5,NMXHEP),VHKK(4,NMXHEP)
      DIMENSION NOW(8)

      NEVHEP=-1
C---RUN AND GENERATOR'S IDENTIFIER
      ISTHEP(1)=1
      IDHEP(1)=0
C---DATE AND TIME
C      CALL DATIM(NOW)
       DO 11 INW=1,8
   11  NOW(INW)=0
      JMOHEP(1,1)=NOW(6)+NOW(7)*100+(NOW(8)-1900)*10000
      JMOHEP(2,1)=NOW(3)+NOW(4)*100+NOW(5)*10000
C---VERSION NUMBER AND DATE OF LAST CHANGE
      JDAHEP(1,1)=0406
      JDAHEP(2,1)=961127

C---REAL PARAMETERS FROM INPUT AND HSPRLG
      PHEP(1,2)=EELE
      PHEP(1,3)=POLARI
      PHEP(1,4)=EPRO
      PHEP(1,5)=XMIN
      PHEP(1,6)=XMAX
      PHEP(1,7)=YMIN
      PHEP(1,8)=YMAX
      PHEP(1,9)=Q2MIN
      PHEP(1,10)=WMIN
      PHEP(1,11)=EGMIN
      PHEP(1,12)=MW
      PHEP(1,13)=MZ
      PHEP(1,14)=MH
      PHEP(1,15)=MT

C---OPTION FLAGS
      ISTHEP(2)=LLEPT
      ISTHEP(3)=ICUT
      DO 1 I=1,12
    1 ISTHEP(I+3)=LPARIN(I)
      ISTHEP(16)=IPART
      ISTHEP(17)=NPYMIN
      ISTHEP(18)=NPYMAX
      ISTHEP(19)=LUNIN
      ISTHEP(20)=LUNOUT
      ISTHEP(21)=LUNTES
      ISTHEP(22)=LUNRND
      ISTHEP(23)=LUNDAT
      ISTHEP(24)=NINP
      ISTHEP(25)=NOUTP
      DO 2 I=1,5
    2 ISTHEP(25+I)=INT2(I)
      DO 3 I=1,15
    3 ISTHEP(30+I)=INT3(I)
      ISTHEP(46)=NPOVEG
      ISTHEP(47)=NUMINT
      ISTHEP(48)=NPHYP
      DO 4 I=1,5
    4 ISTHEP(48+I)=ISAM2(I)
      DO 5 I=1,15
    5 ISTHEP(53+I)=ISAM3(I)
      ISTHEP(69)=IOPLOT
      ISTHEP(70)=IPRINT
      ISTHEP(71)=ISDINP
      ISTHEP(72)=ISDOUT
      ISTHEP(73)=NEVENT
      NHEP=73

      RETURN
      END
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C---FILL TERMINATING RECORD OF COMMON HEPEVT
C
      SUBROUTINE HSTERM
      IMPLICIT DOUBLE PRECISION (A-H,M,O-Z)
C---------------------------------------------------------------------
      COMMON /HSNUME/ SIGTOT,SIGTRR,SIGG(20),SIGGRR(20),NEVENT,NEVE(20)
      COMMON /HSGSW/  SW,CW,SW2,CW2
     *              ,MW,MZ,MH,ME,MMY,MTAU,MU,MD,MS,MC,MB,MT
     *              ,MW2,MZ2,MH2,ME2,MMY2,MTAU2,MU2,MD2,MS2,MC2,MB2,MT2
      COMMON /PYSTFUC/ PYSTOP,PYSLAM,NPYMOD,NPYMAX,NPYMIN
      REAL*4           PYSTOP,PYSLAM
      COMMON /HSCUTS/ XMIN,XMAX,Q2MIN,Q2MAX,YMIN,YMAX,WMIN,GMIN
      PARAMETER (NMXHEP=2000)
      COMMON /HEPEVT/ NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
     &                JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),
     &                PHEP(5,NMXHEP),VHKK(4,NMXHEP)

      NEVHEP=-2

C---CROSS SECTIONS
      PHEP(1,1)=SIGTOT/1D3
      PHEP(1,2)=SIGTRR/1D3
      DO 1 I=1,20
      PHEP(1,2+I)=SIGG(I)/1D3
    1 PHEP(1,22+I)=SIGGRR(I)/1D3
C---MASSES, SW2
      PHEP(1,43)=SW2
      PHEP(1,44)=MW
C---EVENT NUMBERS
      ISTHEP(1)=NEVENT
      DO 2 I=1,20
    2 ISTHEP(1+I)=NEVE(I)
      NHEP=44

      RETURN
      END

************************************************************************

C FUNCTION RETURNS THE BIN WIDTH OF A HISTOGRAM
      FUNCTION HBINW(ID)
      IMPLICIT NONE
      REAL*8 HBINW
      INTEGER ID
      CHARACTER*80 CHTITL
      INTEGER NX,NY,NWT,LOC
      REAL*4 XMI,XMA,YMI,YMA
      CALL HGIVE(ID,CHTITL,NX,XMI,XMA,NY,YMI,YMA,NWT,LOC)
      IF (NY.EQ.0) THEN
         HBINW=(DBLE(XMA)-DBLE(XMI))/DBLE(NX)
      ELSE
         HBINW=(DBLE(XMA)-DBLE(XMI))*(DBLE(YMA)-DBLE(YMI))
     .        /(DBLE(NX)*DBLE(NY))
      ENDIF
      RETURN
      END