*
* $Id: hwhew2.F,v 1.1.1.1 1996/03/08 17:02:14 mclareni Exp $
*
* $Log: hwhew2.F,v $
* Revision 1.1.1.1  1996/03/08 17:02:14  mclareni
* Herwig58
*
*
*CMZ :          29/08/94  11.51.47  by  Unknown
*-- Author :
CDECK  ID>, HWHEW2.
*CMZ :-        -26/04/91  13.22.25  by  Federico Carminati
*-- Author :    Zoltan Kunszt, modified by Bryan Webber
C-----------------------------------------------------------------------
      SUBROUTINE HWHEW2(NPART,PPCM,H,CH,D)
C PCM SHOULD BE DEFINED SUCH THAT ALL 4-MOMENTA ARE OUTGOING.
C CONVENTION FOR PCM AND P IS THAT DIRECTION 1 =BEAM, COMPONENT
C 4 = ENERGY AND COMPONENT 2 AND 3 ARE TRANSVERSE COMPONENTS.
C THUS INCOMING MOMENTA SHOULD CORRESPOND TO OUTGOING MOMENTA
C OF NEGATIVE ENERGY.
C PCM IS FILLED BY PHASE SPACE MONTE CARLO.
C I1-I7 HERE REFER TO HOW PCM INDEXING IS MAPPED TO OUR STANDARD
C 1-6=GLUON,GLUON,Q,QBAR,QP,QPBAR ORDERING `
      IMPLICIT NONE
      DOUBLE PRECISION ZERO, ONE, TWO, THREE, HALF
      PARAMETER (ZERO=0.D0, ONE=1.D0, TWO=2.D0,
     &           THREE=3.D0, HALF=0.5D0)
      INTEGER J,L,IJ,II,JJ,I,NPART,IP1,IPP1
      COMPLEX PT5,ZT,Z1,ZI,ZP,ZQ,ZD,ZPS,ZQS,
     & ZDPM,ZDMP,H(7,7),CH(7,7),D(7,7)
      DOUBLE PRECISION PPCM(5,7),P(5,7),WRN(7),EPS,Q1,Q2,QP,QM,
     & P1,P2,PP,PM,DMP,DPM,PT,QT,PTI,QTI
      EPS=0.0000001
      ZI=CMPLX(0.,1.)
      Z1=CMPLX(1.,0.)
C FOLLOWING DO LOOP IS TO CONVERT TO OUR STANDARD INDEXING
      DO 1 L=1,NPART
      DO 1 IJ=1,4
1     P(IJ,L)=PPCM(IJ,L)
      DO 2 II=1,7
      WRN(II)=ONE
      IF(P(4,II).LT.ZERO) WRN(II)=-ONE
      DO 2 JJ=1,4
      P(JJ,II)=WRN(II)*P(JJ,II)
    2 CONTINUE
C THE ABOVE CHECKS FOR MOMENTA WITH NEGATIVE ENERGY,INNER PRODUCTS
C ARE EXPRESSED DIFFERENTLY FOR DIFFERENT CASES
      DO 11 I=1,NPART-1
      IP1=I+1
      DO 11 J=IP1,NPART
      Q1=P(4,I)+P(1,I)
      QP=0.0
      IF(Q1.GT.EPS)QP=SQRT(Q1)
      Q2=P(4,I)-P(1,I)
      QM=0.0
      IF(Q2.GT.EPS)QM=SQRT(Q2)
      P1=P(4,J)+P(1,J)
      PP=0.
      IF(P1.GT.EPS)PP=SQRT(P1)
      P2=P(4,J)-P(1,J)
      PM=0.
      IF(P2.GT.EPS)PM=SQRT(P2)
      DMP=PM*QP
      ZDMP=CMPLX(DMP,ZERO)
      DPM=PP*QM
      ZDPM=CMPLX(DPM,ZERO)
C NOTE THAT IN OUR INNER PRODUCT NOTATION WE ARE COMPUTING <P,Q>
      PT=SQRT(P(2,J)**2+P(3,J)**2)
      QT=SQRT(P(2,I)**2+P(3,I)**2)
      IF(PT.GT.EPS) GO TO 99
      ZP=Z1
      GO TO 98
   99 PTI=ONE/PT
      ZP=CMPLX(PTI*P(2,J),PTI*P(3,J))
   98 ZPS=CONJG(ZP)
      IF(QT.GT.EPS) GO TO 89
      ZQ=Z1
      GO TO 88
   89 QTI=ONE/QT
      ZQ=CMPLX(QTI*P(2,I),QTI*P(3,I))
   88 ZQS=CONJG(ZQ)
      ZT=Z1
      IF(WRN(I).LT.0) ZT=ZT*ZI
      IF(WRN(J).LT.0) ZT=ZT*ZI
      H(J,I)=(ZDMP*ZP-ZDPM*ZQ)*ZT
      CH(J,I)=(ZDMP*ZPS-ZDPM*ZQS)*ZT
      ZD=H(J,I)*CH(J,I)
      PT5=CMPLX(.5,0.)
      D(J,I)=PT5*ZD
   11 CONTINUE
      DO 60 I=1,NPART-1
      IPP1=I+1
      DO 60 J=IPP1,NPART
      H(I,J)=-H(J,I)
      CH(I,J)=-CH(J,I)
   60 D(I,J)=D(J,I)
      RETURN
      END
