*
* $Id: hwegam.F,v 1.1.1.1 1996/03/08 17:02:12 mclareni Exp $
*
* $Log: hwegam.F,v $
* Revision 1.1.1.1  1996/03/08 17:02:12  mclareni
* Herwig58
*
*
*CMZ :          29/08/94  11.51.46  by  Unknown
*-- Author :
CDECK  ID>, HWEGAM.
*CMZ :-        -26/04/91  11.11.55  by  Bryan Webber
*-- Author :    Bryan Webber & Luca Stanco
C-----------------------------------------------------------------------
      SUBROUTINE HWEGAM(IHEP,GAMWT,ZMI,ZMA,WWA)
C     GENERATES A PHOTON IN WEIZSACKER-WILLIAMS (WWA=.TRUE.) OR
C     ELSE EQUIVALENT PHOTON APPROX FROM
C     INCOMING E+, E-, MU+ OR MU-
C-----------------------------------------------------------------------
C     MODIFIED 5/8/92 BY MHS TO RESET MOMENTA IF SECOND BEAM IS SPLIT
C       TO LEPTON-PHOTON IN LEPTON-LEPTON COLLISIONS (IE IF IPRO.GE.90)
C     MODIFIED 12/10/93 BY BRW TO INCLUDE EXACT KINEMATICS IF .NOT.WWA
C-----------------------------------------------------------------------
#include "herwig58/herwig58.inc"
      INTEGER IHEP,IHADIS,HQ,NTRY
      LOGICAL WWA
      DOUBLE PRECISION HWRGEN,EGMIN,ZMIN,ZMAX,ZGAM,GAMWT,SS,
     & ZMI,ZMA,PPL,PMI,QT2,Q2,QQMIN,QQMAX,HWRUNI,S0
      DATA EGMIN/5.D0/
      IF (IERROR.NE.0)  RETURN
      IF (IHEP.LT.1.OR.IHEP.GT.2) CALL HWWARN('HWEGAM',500,*999)
      SS=PHEP(5,3)
      IF (IHEP.EQ.1) THEN
        IHADIS=2
      ELSE
        IHADIS=1
        IF (JDAHEP(1,IHADIS).NE.0) IHADIS=JDAHEP(1,IHADIS)
      ENDIF
C---DEFINE LIMITS FOR GAMMA MOMENTUM FRACTION
      IF (ZMI.LE.0D0 .OR. ZMA.GT.1D0) THEN
        IF (IPRO.EQ.13.OR.IPRO.EQ.14) THEN
          S0 = EMMIN**2
        ELSEIF(IPRO.EQ.15.OR.IPRO.EQ.18.OR.IPRO.EQ.22.OR.
     &         IPRO.EQ.50.OR.IPRO.EQ.55)THEN
          S0 = 4.D0*PTMIN**2
        ELSEIF (IPRO.EQ.17.OR.IPRO.EQ.51) THEN
          HQ = MOD(IPROC,100)
          S0 = 4.D0*(PTMIN**2+RMASS(HQ)**2)
        ELSEIF (IPRO.EQ.16.OR.IPRO.EQ.19.OR.IPRO.EQ.95) THEN
          S0 = MAX(2*RMASS(1),RMASS(201)-GAMMAX*GAMH)**2
        ELSEIF (IPRO.EQ.23) THEN
          S0 = MAX(2*RMASS(1),RMASS(201)-GAMMAX*GAMH)**2
          S0 = (PTMIN+SQRT(PTMIN**2+S0))**2
        ELSEIF (IPRO.EQ.20) THEN
          S0 = RMASS(201)**2
        ELSEIF (IPRO.EQ.21) THEN
          S0 = (PTMIN+SQRT(PTMIN**2+RMASS(198)**2))**2
        ELSEIF (IPRO.EQ.52) THEN
          HQ = MOD(IPROC,100)
          S0 = (PTMIN+SQRT(PTMIN**2+RMASS(HQ)**2))**2
        ELSEIF (IPRO.EQ.90) THEN
          S0 = Q2MIN
        ELSEIF (IPRO.EQ.91.OR.IPRO.EQ.92) THEN
          S0 = Q2MIN+4.D0*PTMIN**2
          HQ = MOD(IPROC,100)
          IF (HQ.GT.0) S0 = S0+4.D0*RMASS(HQ)**2
          IF (IPRO.EQ.91) S0 = MAX(S0,EMMIN**2)
        ELSE
          S0 = 0
        ENDIF
        IF (S0.GT.0) THEN
          S0 = (SQRT(S0)+ABS(PHEP(5,IHADIS)))**2-PHEP(5,IHADIS)**2
          ZMIN = S0 / (SS**2 - PHEP(5,IHEP)**2 - PHEP(5,IHADIS)**2)
          ZMAX = ONE
        ELSE
C---UNKNOWN PROCESS: USE ENERGY CUTOFF, AND WARN USER
          IF (FSTWGT) CALL HWWARN('HWEGAM',1,*999)
          ZMIN = EGMIN / PHEP(4,IHEP)
          ZMAX = ONE
        ENDIF
      ELSE
        ZMIN=ZMI
        ZMAX=ZMA
      ENDIF
C---APPLY USER DEFINED CUTS YBMIN,YBMAX AND INDIRECT LIMITS ON Z
      IF (.NOT.WWA) THEN
        ZMIN=MAX(ZMIN,YBMIN,SQRT(Q2WWMN)/ABS(PHEP(3,IHEP)))
        ZMAX=MIN(ZMAX,YBMAX)
        IF (ZMIN.GT.ZMAX) THEN
          GAMWT=ZERO
          RETURN
        ENDIF
      ENDIF
C---GENERATE GAMMA MOMENTUM FRACTION
      NTRY=0
   10 NTRY=NTRY+1
      IF (NTRY.GT.NETRY) CALL HWWARN('HWEGAM',51,*999)
      ZGAM=(ZMIN/ZMAX)**HWRGEN(1)*ZMAX
      IF (ONE+(ONE-ZGAM)**2.LT.TWO*HWRGEN(2)) GO TO 10
      IF (WWA) THEN
        GAMWT = GAMWT * .5*ALPHEM/PIFAC *
     +       LOG((ONE-ZGAM)/ZGAM*(SS/PHEP(5,IHEP))**2) *
     +       (TWO*LOG(ZMAX/ZMIN)-TWO*(ZMAX-ZMIN)+HALF*(ZMAX**2-ZMIN**2))
      ELSE
C---Q2WWMN AND Q2WWMX ARE USER-DEFINED LIMITS IN THE Q**2 INTEGRATION
        QQMAX=MIN(Q2WWMX,(ZGAM*PHEP(3,IHEP))**2)
        QQMIN=MAX(Q2WWMN,(PHEP(5,IHEP)*ZGAM)**2/(1.-ZGAM))
        IF (QQMIN.GT.QQMAX) CALL HWWARN('HWEGAM',50,*10)
        Q2=EXP(HWRUNI(0,LOG(QQMIN),LOG(QQMAX)))
        GAMWT = GAMWT * .5*ALPHEM/PIFAC * LOG(QQMAX/QQMIN) *
     +       (TWO*LOG(ZMAX/ZMIN)-TWO*(ZMAX-ZMIN)+HALF*(ZMAX**2-ZMIN**2))
      ENDIF
      IF (GAMWT.LT.ZERO) GAMWT=ZERO
C---FILL PHOTON
      NHEP=NHEP+1
      IDHW(NHEP)=59
      ISTHEP(NHEP)=3
      IDHEP(NHEP)=22
      JMOHEP(1,NHEP)=IHEP
      JMOHEP(2,NHEP)=0
      JDAHEP(1,NHEP)=0
      JDAHEP(2,NHEP)=0
      JDAHEP(1,IHEP)=NHEP
      IF (WWA) THEN
C---FOR COLLINEAR KINEMATICS, ZGAM IS THE ENERGY FRACTION
        PHEP(4,NHEP)=PHEP(4,IHEP)*ZGAM
        PHEP(3,NHEP)=PHEP(3,IHEP)-SIGN(SQRT(
     &     (PHEP(4,IHEP)-PHEP(4,NHEP))**2-PHEP(5,IHEP)**2),PHEP(3,IHEP))
        PHEP(2,NHEP)=0
        PHEP(1,NHEP)=0
        CALL HWUMAS(PHEP(1,NHEP))
      ELSE
C---FOR EXACT KINEMATICS, ZGAM IS TAKEN TO BE FRACTION OF (E+PZ)
        PPL=ZGAM*(ABS(PHEP(3,IHEP))+PHEP(4,IHEP))
        QT2=(ONE-ZGAM)*Q2-(ZGAM*PHEP(5,IHEP))**2
        PMI=(QT2-Q2)/PPL
        PHEP(5,NHEP)=-SQRT(Q2)
        PHEP(4,NHEP)=(PPL+PMI)/TWO
        PHEP(3,NHEP)=SIGN((PPL-PMI)/TWO,PHEP(3,IHEP))
        CALL HWRAZM(SQRT(QT2),PHEP(1,NHEP),PHEP(2,NHEP))
      ENDIF
C---UPDATE OVERALL CM FRAME
      JMOHEP(IHEP,3)=NHEP
      CALL HWVDIF(4,PHEP(1,3),PHEP(1,IHEP),PHEP(1,3))
      CALL HWVSUM(4,PHEP(1,NHEP),PHEP(1,3),PHEP(1,3))
      CALL HWUMAS(PHEP(1,3))
C---FILL OUTGOING LEPTON
      NHEP=NHEP+1
      IDHW(NHEP)=IDHW(IHEP)
      ISTHEP(NHEP)=1
      IDHEP(NHEP)=IDHEP(IHEP)
      JMOHEP(1,NHEP)=IHEP
      JMOHEP(2,NHEP)=0
      JDAHEP(1,NHEP)=0
      JDAHEP(2,NHEP)=0
      JDAHEP(2,IHEP)=NHEP
      CALL HWVDIF(4,PHEP(1,IHEP),PHEP(1,NHEP-1),PHEP(1,NHEP))
      PHEP(5,NHEP)=PHEP(5,IHEP)
 999  END