*
* $Id: hwccut.F,v 1.1.1.1 1996/03/08 17:02:10 mclareni Exp $
*
* $Log: hwccut.F,v $
* Revision 1.1.1.1  1996/03/08 17:02:10  mclareni
* Herwig58
*
*
*CMZ :          29/08/94  11.51.46  by  Unknown
*-- Author :
CDECK  ID>, HWCCUT.
*CMZ :-        -26/04/91  14.29.39  by  Federico Carminati
*-- Author :    Bryan Webber
C------------------------------------------------------------------------
      SUBROUTINE HWCCUT(JHEP,KHEP,PCL,BTCLUS,SPLIT)
C     CUTS INTO TWO A CLUSTER OF MOMENTUM
C     PCL FORMED FROM PARTONS JHEP AND KHEP
C------------------------------------------------------------------------
#include "herwig58/herwig58.inc"
      DOUBLE PRECISION HWREXQ,HWRGEN,HWUPCM,HWVDOT,EMC,QM1,QM2,EMX,EMY,
     & QM3,PCX,PCY,PXY,RCM,PCL(5),AX(5),PA(5),PB(5),PC(5)
      INTEGER HWRINT,JHEP,KHEP,LHEP,MHEP,ID1,ID2,ID3,NTRY,J
      LOGICAL BTCLUS,SPLIT
      IF (IERROR.NE.0) RETURN
      EMC=PCL(5)
      ID1=IDHW(JHEP)
      ID2=IDHW(KHEP)
      QM1=RMASS(ID1)
      QM2=RMASS(ID2)
      SPLIT=.TRUE.
      NTRY=0
   10 NTRY=NTRY+1
      IF (NTRY.GT.100) THEN
C---SPLITTING FAILED 100 TIMES - ASSUME NO PHASE SPACE
        SPLIT=.FALSE.
        RETURN
      ENDIF
      IF (BTCLUS) THEN
        SPLIT=.FALSE.
        ID3=HWRINT(1,2)
        QM3=RMASS(ID3)
C---SPLIT BEAM AND TARGET CLUSTERS AS SOFT CLUSTERS
        IF (EMC.LE.QM1+QM2+2*QM3) GOTO 10
        EMX=QM1+QM3+HWREXQ(BTCLM,EMC-QM1-QM2-2*QM3)
        EMY=QM2+QM3+HWREXQ(BTCLM,EMC-QM1-QM2-2*QM3)
        IF (EMX+EMY.GE.EMC) GOTO 10
        PCX=HWUPCM(EMX,QM1,QM3)
        PCY=HWUPCM(EMY,QM2,QM3)
      ELSE
C---CHOOSE FRAGMENT MASSES
        PXY=EMC-QM1-QM2
   12   EMX=QM1+PXY*HWRGEN(0)**PSPLT
        EMY=QM2+PXY*HWRGEN(1)**PSPLT
        IF (EMX+EMY.GE.EMC) GO TO 12
C---U,D,S PAIR PRODUCTION WITH WEIGHTS QWT
   15   ID3=HWRINT(1,3)
        IF (QWT(ID3).LT.HWRGEN(3)) GO TO 15
        QM3=RMASS(ID3)
        PCX=HWUPCM(EMX,QM1,QM3)
        IF (PCX.LT.0.) GO TO 10
        PCY=HWUPCM(EMY,QM2,QM3)
        IF (PCY.LT.0.) GO TO 10
      ENDIF
      PXY=HWUPCM(EMC,EMX,EMY)
C---BOOST ANTIQUARK TO C.M. FRAME TO FIND AXIS
      CALL HWULOF(PCL,PHEP(1,KHEP),AX)
      RCM=1./SQRT(HWVDOT(3,AX(1),AX(1)))
      CALL HWVSCA(3,RCM,AX,AX)
C---CONSTRUCT NEW C.M. MOMENTA (COLLINEAR)
      CALL HWVSCA(3,PXY,AX,PC)
      PC(4)=SQRT(PXY**2+EMY**2)
      PC(5)=EMY
      CALL HWVSCA(3,PCY,AX,PA)
      PA(4)=SQRT(PCY**2+QM2**2)
      PA(5)=QM2
      CALL HWULOB(PC,PA,PB)
      CALL HWVDIF(4,PC,PB,PA)
      PA(5)=QM3
      LHEP=NHEP+1
      MHEP=NHEP+2
      CALL HWULOB(PCL,PB,PHEP(1,KHEP))
      CALL HWULOB(PCL,PA,PHEP(1,MHEP))
      CALL HWVSCA(3,-ONE,PC,PC)
      PC(4)=EMC-PC(4)
      PC(5)=EMX
      CALL HWVSCA(3,PCX,AX,PA)
      PA(4)=SQRT(PCX**2+QM3**2)
      CALL HWULOB(PC,PA,PB)
      CALL HWULOB(PCL,PB,PHEP(1,LHEP))
      DO 45 J=1,4
   45 PHEP(J,JHEP)=PCL(J)-PHEP(J,KHEP)-PHEP(J,LHEP)-PHEP(J,MHEP)
      PHEP(5,JHEP)=QM1
      IDHW(LHEP)=ID3+6
      IDHW(MHEP)=ID3
      IDHEP(MHEP)= IDPDG(ID3)
      IDHEP(LHEP)=-IDPDG(ID3)
      ISTHEP(LHEP)=151
      ISTHEP(MHEP)=151
      JMOHEP(2,JHEP)=LHEP
      JDAHEP(2,KHEP)=MHEP
      JMOHEP(1,LHEP)=JMOHEP(1,KHEP)
      JMOHEP(2,LHEP)=MHEP
      JDAHEP(1,LHEP)=0
      JDAHEP(2,LHEP)=JHEP
      JMOHEP(1,MHEP)=JMOHEP(1,JHEP)
      JMOHEP(2,MHEP)=KHEP
      JDAHEP(1,MHEP)=0
      JDAHEP(2,MHEP)=LHEP
      NHEP=NHEP+2
  999 END