CDECK  ID>, HWHEPA.
*CMZ :-        -26/04/91  14.55.44  by  Federico Carminati
*-- Author :    Bryan Webber and Ian Knowles
C-----------------------------------------------------------------------
      SUBROUTINE HWHEPA
C-----------------------------------------------------------------------
C     (Initially polarised) e+e- --> ffbar (f=quark, mu or tau)
C     If IPROC=107: --> gg, distributed as sum of light quarks.
C     If fermion flavour specified mass effects fully included.
C     EVWGT=sig(e+e- --> ffbar) in nb
C-----------------------------------------------------------------------
      INCLUDE 'HERWIG59.INC'
      DOUBLE PRECISION HWRGEN,HWRUNI,HWUPCM,HWUAEM,Q2NOW,Q2LST,FACTR,
     & VF2,VF,CLF(7),PRAN,PQWT,PMAX,PTHETA,SINTH2,CPHI,SPHI,C2PHI,S2PHI,
     & PPHI,SINTH,PCM,PP(5),EWGT
      INTEGER ID1,ID2,IDF,IQ,IQ1,I
      EXTERNAL HWRGEN,HWRUNI,HWUPCM,HWUAEM
      SAVE Q2LST,FACTR,ID1,ID2,VF2,VF,CLF,EWGT
      DATA Q2LST/0./
      IF (GENEV) THEN
        IF (ID2.EQ.0) THEN
C Choose quark flavour
          PRAN=TQWT*HWRGEN(0)
          PQWT=0.
          DO 10 IQ=1,MAXFL
          PQWT=PQWT+CLQ(1,IQ)
          IF (PQWT.GT.PRAN) GOTO 11
   10     CONTINUE
          IQ=MAXFL
   11     IQ1=MAPQ(IQ)
          DO 20 I=1,7
   20     CLF(I)=CLQ(I,IQ)
        ELSE
          IQ1=ID1
        ENDIF
C Label particles, assign outgoing particle masses
        IDHW(NHEP+1)=200
        IDHEP(NHEP+1)=23
        ISTHEP(NHEP+1)=110
        IF (ID1.EQ.7) THEN
          IDHW(NHEP+2)=13
          IDHW(NHEP+3)=13
          IDHEP(NHEP+2)=21
          IDHEP(NHEP+3)=21
          PHEP(5,NHEP+2)=RMASS(13)
          PHEP(5,NHEP+3)=RMASS(13)
        ELSE
          IDHW(NHEP+2)=IQ1
          IDHW(NHEP+3)=IQ1+6
          IDHEP(NHEP+2)=IDPDG(IQ1)
          IDHEP(NHEP+3)=-IDHEP(NHEP+2)
          PHEP(5,NHEP+2)=RMASS(IQ1)
          PHEP(5,NHEP+3)=RMASS(IQ1)
        ENDIF
        ISTHEP(NHEP+2)=113
        ISTHEP(NHEP+3)=114
        JMOHEP(1,NHEP+1)=1
        IF (JDAHEP(1,1).NE.0) JMOHEP(1,NHEP+1)=JDAHEP(1,1)
        JMOHEP(2,NHEP+1)=2
        IF (JDAHEP(1,2).NE.0) JMOHEP(1,NHEP+1)=JDAHEP(1,2)
        JMOHEP(1,NHEP+2)=NHEP+1
        JMOHEP(2,NHEP+2)=NHEP+3
        JMOHEP(1,NHEP+3)=NHEP+1
        JMOHEP(2,NHEP+3)=NHEP+2
        JDAHEP(1,NHEP+1)=NHEP+2
        JDAHEP(2,NHEP+1)=NHEP+3
        JDAHEP(1,NHEP+2)=0
        JDAHEP(2,NHEP+2)=NHEP+3
        JDAHEP(1,NHEP+3)=0
        JDAHEP(2,NHEP+3)=NHEP+2
C Generate polar and azimuthal angular distributions:
C  CLF(1)*(1+(VF*COSTH)**2)+CLF(2)*(1-VF**2)+CLF(3)*2.*VF*COSTH
C +(VF*SINTH)**2*(CLF(4)*COS(2*PHI-PHI1-PHI2)
C                +CLF(6)*SIN(2*PHI-PHI1-PHI2))
        PMAX=CLF(1)*(1.+VF2)+CLF(2)*(1.-VF2)+ABS(CLF(3))*2.*VF
  30    COSTH=HWRUNI(0,-ONE, ONE)
        PTHETA=CLF(1)*(1.+VF2*COSTH**2)+CLF(2)*(1.-VF2)
     &        +CLF(3)*2.*VF*COSTH
        IF (PTHETA.LT.PMAX*HWRGEN(1)) GOTO 30
        IF (IDHW(1).GT.IDHW(2)) COSTH=-COSTH
        SINTH2=1.-COSTH**2
        IF (TPOL) THEN
          PMAX=PTHETA+VF2*SINTH2*SQRT(CLF(4)**2+CLF(6)**2)
  40      CALL HWRAZM(ONE,CPHI,SPHI)
          C2PHI=2.*CPHI**2-1.
          S2PHI=2.*CPHI*SPHI
          PPHI=PTHETA+(CLF(4)*(C2PHI*COSS+S2PHI*SINS)
     &                +CLF(6)*(S2PHI*COSS-C2PHI*SINS))*VF2*SINTH2
          IF (PPHI.LT.PMAX*HWRGEN(1)) GOTO 40
        ELSE
          CALL HWRAZM(ONE,CPHI,SPHI)
        ENDIF
C Construct final state 4-mommenta
        CALL HWVEQU(5,PHEP(1,3),PHEP(1,NHEP+1))
        PCM=HWUPCM(PHEP(5,NHEP+1),PHEP(5,NHEP+2),PHEP(5,NHEP+3))
C PP is momentum of track NHEP+2 in CoM (track NHEP+1) frame
        SINTH=SQRT(SINTH2)
        PP(5)=PHEP(5,NHEP+2)
        PP(1)=PCM*SINTH*CPHI
        PP(2)=PCM*SINTH*SPHI
        PP(3)=PCM*COSTH
        PP(4)=SQRT(PCM**2+PP(5)**2)
        CALL HWULOB(PHEP(1,NHEP+1),PP(1),PHEP(1,NHEP+2))
        CALL HWVDIF(4,PHEP(1,NHEP+1),PHEP(1,NHEP+2),PHEP(1,NHEP+3))
C Set production vertices
        CALL HWVZRO(4,VHEP(1,NHEP+2))
        CALL HWVEQU(4,VHEP(1,NHEP+2),VHEP(1,NHEP+3))
        NHEP=NHEP+3
      ELSE
        EMSCA=PHEP(5,3)
        Q2NOW=EMSCA**2
        IF (Q2NOW.NE.Q2LST) THEN
C Calculate coefficients for cross-section
          EMSCA=PHEP(5,3)
          Q2LST=Q2NOW
          FACTR=PIFAC*GEV2NB*HWUAEM(Q2NOW)**2/Q2NOW
          ID1=MOD(IPROC,10)
          ID2=MOD(ID1,7)
          IF (ID2.EQ.0) THEN
            CALL HWUEEC(1)
            VF2=1.
            VF=1.
            EWGT=FACTR*FLOAT(NCOLO)*TQWT*4./3.
          ELSE
            IF (IPROC.LT.150) THEN
              IDF=ID1
              FACTR=FACTR*FLOAT(NCOLO)
            ELSE
              ID1=2*ID1+119
              IDF=ID1-110
            ENDIF
            IF (EMSCA.LE.2.*RMASS(ID1)) then
              EWGT=0.
            ELSE
              CALL HWUCFF(11,IDF,Q2NOW,CLF(1))
              VF2=1.-4.*RMASS(ID1)**2/Q2NOW
              VF=SQRT(VF2)
              EWGT=FACTR*VF*(CLF(1)*(1.+VF2/3.)+CLF(2)*(1.-VF2))
            ENDIF
          ENDIF
        ENDIF
        EVWGT=EWGT
      ENDIF
  999 END
