CDECK  ID>, HWHEGG.
*CMZ :-        -19/03/92  10.13.56  by  Mike Seymour
*-- Author :    Mike Seymour
C-----------------------------------------------------------------------
      SUBROUTINE HWHEGG
C----------------------------------------------------------------------
C     HARD PROCESS: EE --> EEGAMGAM --> EEFFBAR/WW
C     MEAN EVENT WEIGHT = CROSS-SECTION IN NB
C     AFTER CUTS ON PT AND MASS OF CENTRE-OF-MASS SYSTEM
C     AND COS(THETA) IN CENTRE-OF-MASS SYSTEM
C     AND TIMES BRANCHING FRACTION IF WW
C-----------------------------------------------------------------------
      INCLUDE 'HERWIG59.INC'
      DOUBLE PRECISION HWRGEN,HWULDO,EMSQ,BETA,S,T,U,TMIN,TMAX,TRAT,
     & DSDT,PROB,X,Z(2),ZMIN,ZMAX,PCMIN,PCMAX,PCFAC,PLOGMI,PLOGMA,PTCMF,
     & Q,PC,BLOG,EMCMIN,EMCMAX,EMLMIN,EMLMAX,WGT(6),RWGT,CV,CA,BR,QT(2),
     & QX(2),QY(2),PX,PY,ROOTS,DOT,A,B,C,SHAT,PCF(2),PCM(2),PCMAC,ZZ(2),
     & COLFAC
      INTEGER I,IGAM,ID,IDL,ID1,ID2,IHEP,JHEP,NADD,NTRY,NQ,JGAM
      LOGICAL HWRLOG
      EXTERNAL HWRGEN,HWULDO,HWRLOG
      SAVE S,BETA,X,ID,NQ,WGT,EMLMIN,EMLMAX,PCFAC,PLOGMA,PLOGMI,SHAT,
     &  PCF,PCM,Z,PCMAC,NADD
      IF (IERROR.NE.0) RETURN
C---INITIALIZE LOCAL COPIES OF EMMIN,EMMAX
      IF (FSTWGT) THEN
        EMLMIN=EMMIN
        EMLMAX=EMMAX
      ENDIF
      IF (.NOT.GENEV) THEN
C---CHOOSE Z1,Z2 AND CALCULATE SUB-PROCESS CROSS-SECTION
        EVWGT=0
C-----FIND FINAL STATE PARTICLES
        IHPRO=MOD(IPROC,100)
        IF (IHPRO.EQ.0) THEN
          ID=1
          NQ=6
          COLFAC=FLOAT(NCOLO)
          NADD=6
        ELSEIF (IHPRO.LE.6) THEN
          ID=IHPRO
          NQ=1
          COLFAC=FLOAT(NCOLO)
          NADD=6
          Q=QFCH(ID)
        ELSEIF (IHPRO.LE.9) THEN
          ID=119+2*(IHPRO-6)
          NQ=1
          COLFAC=1.
          NADD=6
          Q=QFCH(ID-110)
        ELSEIF (IHPRO.LE.10) THEN
          ID=198
          NQ=1
          NADD=1
        ELSE
          CALL HWWARN('HWHEGG',200,*999)
        ENDIF
C-----SPLIT ELECTRONS TO PHOTONS
        NHEP=3
        GAMWT=1
        S=2*HWULDO(PHEP(1,1),PHEP(1,2))
        ROOTS=SQRT(S)
        EMCMIN=MAX(EMLMIN,MAX(2*RMASS(ID),PTMIN))
        EMCMAX=MIN(EMLMAX,ROOTS)
        IF (EMCMIN.GT.EMCMAX) RETURN
        ZMIN=EMCMIN**2/S
        ZMAX=1-PHEP(5,1)/PHEP(4,1)
        IF (ZMIN.GT.ZMAX) RETURN
        CALL HWEGAM(1,ZMIN,ZMAX,.TRUE.)
        Z(1)=PHEP(4,NHEP-1)/PHEP(4,1)
        ZMIN=EMCMIN**2/(Z(1)*S)
        ZMAX=MIN(EMCMAX**2/(Z(1)*S), ONE-PHEP(5,2)/PHEP(4,2))
        IF (ZMIN.GT.ZMAX) RETURN
        CALL HWEGAM(2,ZMIN,ZMAX,.TRUE.)
        Z(2)=PHEP(4,NHEP-1)/PHEP(4,2)
        EMSCA=PHEP(5,3)
        SHAT=EMSCA**2
C-----REMOVE LOG TERMS FROM WEIGHT, CALCULATE NEW ONES FROM PT LIMITS
        GAMWT=GAMWT/(0.5*LOG((1-Z(1))*S/(Z(1)*PHEP(5,1)**2))
     &              *0.5*LOG((1-Z(2))*Z(1)*S/(Z(2)*PHEP(5,2)**2)))
        PCF(1)=Z(1)*PHEP(5,1)
        PCF(2)=Z(2)*PHEP(5,2)
        PCFAC=SQRT(PCF(1)*PCF(2))
        PCM(1)=(1-Z(1))*PHEP(4,1)
        PCM(2)=(1-Z(2))*PHEP(4,2)
        PCMAC=SQRT(PCM(1)*PCM(2))
        PCMIN=MAX(PTMIN,MAX(PCF(1),PCF(2)))
        PCMAX=MIN( MIN(PTMAX,PHEP(5,3)) , MIN(PCM(1),PCM(2)) )
        IF (PCMIN.GT.PCMAX) RETURN
        PLOGMI=(LOG(PCMIN/PCFAC))**2
        PLOGMA=(LOG(PCMAX/PCFAC))**2
        GAMWT=GAMWT*(PLOGMA-PLOGMI)
C-----CALCULATE CROSS-SECTION
        DO 10 IDL=1,NQ
          WGT(IDL)=EVWGT
          IF (IHPRO.EQ.0) THEN
            ID=IDL
            Q=QFCH(ID)
          ENDIF
          EMSQ=RMASS(ID)**2
          X=4*EMSQ/SHAT
          IF (X.GT.1) GOTO 10
          BETA=SQRT(1-X)
          BLOG=LOG((1+BETA*CTMAX)/(1-BETA*CTMAX))/BETA
          IF (IHPRO.LE.9) THEN
            EVWGT=EVWGT+GEV2NB*4*PIFAC*COLFAC*Q**4*ALPHEM**2*BETA
     &           /SHAT * GAMWT * ( (1+X-0.5*X**2)*BLOG
     &                     - CTMAX*(1+X**2/(CTMAX**2*(X-1)+1)) )
            WGT(IDL)=EVWGT
          ELSE
            CALL HWDBOZ(198,ID1,ID2,CV,CA,BR,1)
            CALL HWDBOZ(199,ID1,ID2,CV,CA,BR,1)
            EVWGT=EVWGT + GEV2NB*6*PIFAC*ALPHEM**2*BETA/SHAT*BR
     &        * GAMWT * (-(  X-0.5*X**2)*BLOG
     &                     + CTMAX*(1+(X**2+16/3.)/(CTMAX**2*(X-1)+1)) )
          ENDIF
 10     CONTINUE
C-----GAMWT MUST BE RESET TO ONE, SINCE IT IS REAPPLIED LATER!
        GAMWT=ONE
      ELSE
C---GENERATE EVENT
C-----CHOOSE PT OF THE CMF
        PTCMF=PCFAC*EXP(SQRT(HWRGEN(0)*(PLOGMA-PLOGMI)+PLOGMI))
C-----CHOOSE WHICH PHOTON USUALLY HAS SMALLER PT
        NTRY=0
 20     IGAM=1
        IF (LOG(PCM(1)/PCF(1)).LT.HWRGEN(1)*2*LOG(PCMAC/PCFAC)) IGAM=2
        JGAM=3-IGAM
C-----CHOOSE ITS PT
 30     NTRY=NTRY+1
        IF (NTRY.GT.NBTRY) CALL HWWARN('HWHEGG',100,*999)
        QT(IGAM)=(PCM(IGAM)/PCF(IGAM))**HWRGEN(2)
        PROB=(QT(IGAM)**2/(QT(IGAM)**2+1))**2
        QT(IGAM)=QT(IGAM)*PCF(IGAM)
        IF (HWRLOG(1-PROB)) GOTO 30
C-----CHOOSE ITS DIRECTION
        CALL HWRAZM(QT(IGAM),QX(IGAM),QY(IGAM))
C-----CALCULATE THE OTHER PHOTON'S PT
        QX(JGAM)=PTCMF-QX(IGAM)
        QY(JGAM)=     -QY(IGAM)
        QT(JGAM)=SQRT(QX(JGAM)**2+QY(JGAM)**2)
        IF (QT(JGAM).LT.PCF(JGAM).OR.QT(JGAM).GT.PCM(JGAM)) GOTO 20
C-----APPLY A RANDOM ROTATION AROUND THE BEAM AXIS
        CALL HWRAZM(ONE,PX,PY)
        IF (PX.EQ.0) PX=1D-20
        QX(1)=(QX(1)*PX   -QY(1)*PY)
        QY(1)=(QY(1)      +QX(1)*PY)/PX
        QX(2)=(QX(2)*PX   -QY(2)*PY)
        QY(2)=(QY(2)      +QX(2)*PY)/PX
C-----RECONSTRUCT MOMENTA
        IF (QT(IGAM).GT.QT(JGAM)) THEN
          IGAM=3-IGAM
          JGAM=3-JGAM
        ENDIF
        DOT=-Z(JGAM)*S+SHAT+2*(QX(1)*QX(2)+QY(1)*QY(2))
C-------SOLVE QUADRATIC IN Z(IGAM) TO FIND ELECTRON ENERGIES
        A=S*(S*Z(JGAM)+QT(JGAM)**2)
        B=S*DOT*(1+Z(JGAM))
        C=DOT**2+S*QT(IGAM)**2*(1-Z(JGAM))**2-4*QT(IGAM)**2*QT(JGAM)**2
        IF (B**2.LT.4*A*C) GOTO 20
        ZZ(IGAM)=(-B+SQRT(B**2-4*A*C))/(2*A)
        IF (ZZ(IGAM).LT.0 .OR. ZZ(IGAM).GT.1-Z(IGAM)) GOTO 20
        ZZ(JGAM)=1-Z(JGAM)
C-------REJECT AGAINST PHOTON DISTRIBUTION FUNCTION
        PROB=((1+ZZ(IGAM)**2)/(1-ZZ(IGAM)))/((1+(1-Z(IGAM))**2)/Z(IGAM))
     &      *((1+ZZ(JGAM)**2)/(1-ZZ(JGAM)))/((1+(1-Z(JGAM))**2)/Z(JGAM))
        IF (HWRLOG(1-PROB)) GOTO 20
C-------RECONSTRUCT ALL OTHER VARIABLES
        DO 40 I=1,2
          IGAM=2*I+3
          PHEP(1,IGAM)=QX(I)
          PHEP(2,IGAM)=QY(I)
          PHEP(4,IGAM)=ZZ(I)*PHEP(4,I)
          PHEP(5,IGAM)=RMASS(IDHW(IGAM))
C---------IF MOMENTUM CANNOT BE CONSERVED TRY AGAIN
          IF (PHEP(4,IGAM)**2-PHEP(5,IGAM)**2-QT(I)**2 .LT. 0) GOTO 20
          PHEP(3,IGAM)=SIGN(SQRT(PHEP(4,IGAM)**2-PHEP(5,IGAM)**2-
     &      QT(I)**2),PHEP(3,IGAM))
          CALL HWVDIF(4,PHEP(1,I),PHEP(1,IGAM),PHEP(1,IGAM-1))
          CALL HWUMAS(PHEP(1,IGAM-1))
 40     CONTINUE
C-----TIDY UP EVENT RECORD
        NHEP=NHEP+1
        IDHW(NHEP)=IDHW(3)
        IDHEP(NHEP)=IDHEP(3)
        ISTHEP(NHEP)=110
        CALL HWVSUM(4,PHEP(1,4),PHEP(1,6),PHEP(1,NHEP))
        CALL HWVSUM(4,PHEP(1,1),PHEP(1,2),PHEP(1,3))
        CALL HWUMAS(PHEP(1,NHEP))
        CALL HWUMAS(PHEP(1,3))
        JMOHEP(1,NHEP)=4
        JMOHEP(2,NHEP)=6
        JMOHEP(1,3)=0
        JMOHEP(2,3)=0
C-----CHOOSE FINAL STATE QUARK
        IF (IHPRO.EQ.0) THEN
          RWGT=HWRGEN(2)*EVWGT
          ID=1
          DO 50 IDL=1,NQ
            IF (RWGT.GT.WGT(IDL)) ID=IDL+1
 50       CONTINUE
          EMSQ=RMASS(ID)**2
          X=4*EMSQ/SHAT
          BETA=SQRT(1-X)
        ENDIF
C-----CHOOSE T (WHERE T = MANDELSTAM_T - EMSQ)
        TMIN=-SHAT/2
        TMAX=-SHAT/2*(1-BETA*CTMAX)
        TRAT=TMAX/TMIN
        NTRY=0
        IF (IHPRO.LE.9) THEN
C-------FOR FFBAR, CHOOSE T ACCORDING TO -SHAT/T
 60       NTRY=NTRY+1
          IF (NTRY.GT.NBTRY) CALL HWWARN('HWHEGG',101,*999)
          T=TRAT**HWRGEN(3)*TMIN
          U=-T-SHAT
C-------REWEIGHT TO CORRECT DISTRIBUTION
          DSDT=(T*U-2*EMSQ*(T+2*EMSQ))/T**2
     &        +( 2*EMSQ*(SHAT-4*EMSQ))/(T*U)
     &        +(T*U-2*EMSQ*(U+2*EMSQ))/U**2
          PROB=-DSDT*T/SHAT / (1 + 2*X - 2*X**2)
          IF (HWRLOG(1-PROB)) GOTO 60
        ELSE
C-------FOR WW, CHOOSE T ACCORDING TO (SHAT/T)**2
 70       NTRY=NTRY+1
          IF (NTRY.GT.NBTRY) CALL HWWARN('HWHEGG',102,*999)
          T=TMAX/(1-(1-TRAT)*HWRGEN(4))
          U=-T-SHAT
C-------REWEIGHT TO CORRECT DISTRIBUTION
          DSDT=( 3*(T*U)**2 - SHAT*T*U*(4*SHAT+6*EMSQ)
     &      + SHAT**2*(2*SHAT**2+6*EMSQ**2) ) / (T*U)**2
          PROB=DSDT*(T/SHAT)**2 / (4.75 - 1.5*X + 1.5*X**2)
          IF (HWRLOG(1-PROB)) GOTO 70
        ENDIF
C-----SYMMETRIZE IN T,U
        IF (HWRLOG(HALF)) T=U
C-----FILL EVENT RECORD
        COSTH=(1+2*T/SHAT)/BETA
        PC=0.5*BETA*PHEP(5,NHEP)
        PHEP(5,NHEP+1)=RMASS(ID)
        PHEP(5,NHEP+2)=RMASS(ID)
        CALL HWDTWO(PHEP(1,NHEP),PHEP(1,NHEP+1),PHEP(1,NHEP+2),
     &              PC,COSTH,.TRUE.)
        DO 80 I=1,2
          IHEP=NHEP+I
          JHEP=NHEP+3-I
          ISTHEP(IHEP)=190
          IF (IHPRO.LE.6) ISTHEP(IHEP)=112+I
          IDHW(IHEP)=ID+NADD*(I-1)
          IDHEP(IHEP)=IDPDG(IDHW(IHEP))
          JDAHEP(I,NHEP)=IHEP
          JMOHEP(1,IHEP)=NHEP
          JMOHEP(2,IHEP)=JHEP
          JDAHEP(2,IHEP)=JHEP
          IF (IHPRO.EQ.10) THEN
            RHOHEP(1,IHEP)=0.3333
            RHOHEP(2,IHEP)=0.3333
            RHOHEP(3,IHEP)=0.3333
          ENDIF
 80     CONTINUE
        NHEP=NHEP+2
      ENDIF
 999  END
