*
* $Id: hwdbos.F,v 1.1.1.1 1996/03/08 17:02:11 mclareni Exp $
*
* $Log: hwdbos.F,v $
* Revision 1.1.1.1  1996/03/08 17:02:11  mclareni
* Herwig58
*
*
*CMZ :          29/08/94  11.51.46  by  Unknown
*-- Author :
CDECK  ID>, HWDBOS.
*CMZ :-        -26/04/91  14.55.44  by  Federico Carminati
*-- Author :    Mike Seymour
C-----------------------------------------------------------------------
      SUBROUTINE HWDBOS(IBOS)
C     DECAY GAUGE BOSONS (ALREADY FOUND BY HWDHAD)
C     USES SPIN DENSITY MATRIX IN RHOHEP (1ST CMPT=>-VE,2=>LONG,3=>+VE)
C     IF BOSON CAME FROM HIGGS DECAY, GIVE BOTH THE SAME HELICITY (EPR)
C     IF BOSON CAME FROM W+1JET, GIVE IT THE CORRECT DECAY CORRELATIONS
C-----------------------------------------------------------------------
#include "herwig58/herwig58.inc"
      INTEGER HWRINT,IBOS,IPAIR,ICMF,IOPT,IHEL,IMOTH,I,IQRK,IANT
      DOUBLE PRECISION R(3,3),HWRGEN,HWRUNI,CV,CA,BR,PCM,HWUPCM,PBOS(5),
     &  PMAX,PROB,HWULDO
      LOGICAL QUARKS
      IF (IDHW(IBOS).LT.198.OR.IDHW(IBOS).GT.200)
     &  CALL HWWARN('HWDBOS',101,*999)
      QUARKS=.FALSE.
C---SEE IF IT IS PART OF A PAIR
      IMOTH=JMOHEP(1,IBOS)
      IPAIR=JMOHEP(2,IBOS)
      ICMF=JMOHEP(1,IBOS)
      IF (IDHW(ICMF).EQ.IDHW(IBOS).AND.ISTHEP(ICMF)/10.EQ.12)
     &  ICMF=JMOHEP(1,ICMF)
      IOPT=0
      IF (IPAIR.NE.0) THEN
        IF (JMOHEP(2,IPAIR).NE.IBOS.OR.
     &    IDHW(IPAIR).LT.198.OR.IDHW(IPAIR).GT.200) IPAIR=0
      ENDIF
      IF (IPAIR.GT.0) IOPT=1
C---SELECT DECAY PRODUCTS
   10 CALL HWDBOZ(IDHW(IBOS),IDN(1),IDN(2),CV,CA,BR,IOPT)
C---W + 1JET DECAYS ARE NOW HANDLED HERE !
      IF (IPRO.EQ.21) THEN
        IQRK=IDHW(JMOHEP(1,ICMF))
        IANT=IDHW(JMOHEP(2,ICMF))
        IF (IQRK.EQ.13 .AND. IANT.LE.6) THEN
          IQRK=JMOHEP(2,ICMF)
          IANT=JDAHEP(2,ICMF)
        ELSEIF (IQRK.EQ.13) THEN
          IQRK=JDAHEP(2,ICMF)
          IANT=JMOHEP(2,ICMF)
        ELSEIF (IANT.EQ.13 .AND. IQRK.LE.6) THEN
          IQRK=JMOHEP(1,ICMF)
          IANT=JDAHEP(2,ICMF)
        ELSEIF (IANT.EQ.13) THEN
          IQRK=JDAHEP(2,ICMF)
          IANT=JMOHEP(1,ICMF)
        ELSEIF (IQRK.GT.IANT) THEN
          IQRK=JMOHEP(2,ICMF)
          IANT=JMOHEP(1,ICMF)
        ELSE
          IQRK=JMOHEP(1,ICMF)
          IANT=JMOHEP(2,ICMF)
        ENDIF
        PHEP(5,NHEP+1)=RMASS(IDN(1))
        PHEP(5,NHEP+2)=RMASS(IDN(2))
        PCM=HWUPCM(PHEP(5,IBOS),PHEP(5,NHEP+1),PHEP(5,NHEP+2))
        IF (PCM.LT.0) CALL HWWARN('HWDBOS',103,*999)
        PMAX=HWULDO(PHEP(1,IANT),PHEP(1,IBOS))**2+
     &       HWULDO(PHEP(1,IQRK),PHEP(1,IBOS))**2
 1      CALL HWDTWO(PHEP(1,IBOS),PHEP(1,NHEP+1),PHEP(1,NHEP+2),
     &              PCM,TWO,.TRUE.)
        PROB=HWULDO(PHEP(1,IANT),PHEP(1,NHEP+1))**2+
     &       HWULDO(PHEP(1,IQRK),PHEP(1,NHEP+2))**2
        IF (PROB.GT.PMAX.OR.PROB.LT.0) CALL HWWARN('HWDBOS',104,*999)
        IF (PMAX*HWRGEN(0).GT.PROB) GOTO 1
      ELSE
C---SELECT HELICITY, UNLESS IT IS THE SECOND OF A HIGGS DECAY (EPR)
      IF (IPAIR.NE.IBOS .OR. IDHW(ICMF).NE.201) THEN
      IF (RHOHEP(1,IBOS)+RHOHEP(2,IBOS)+RHOHEP(3,IBOS).LE.0) THEN
C---COPY PARENT HELICITY IF IT WAS A GAUGE BOSON
        IF (IDHW(IMOTH).GE.198.AND.IDHW(IMOTH).LE.200) THEN
          CALL HWVEQU(3,RHOHEP(1,IMOTH),RHOHEP(1,IBOS))
          IF (RHOHEP(1,IBOS)+RHOHEP(2,IBOS)+RHOHEP(3,IBOS).GT.0)
     &    GO TO 20
        ENDIF
        CALL HWWARN('HWDBOS',1,*999)
        RHOHEP(1,IBOS)=1.
        RHOHEP(2,IBOS)=1.
        RHOHEP(3,IBOS)=1.
      ENDIF
 20   IHEL=HWRINT(1,3)
      IF (HWRGEN(0).GT.RHOHEP(IHEL,IBOS)) GOTO 20
      ENDIF
C---SELECT DIRECTION OF FERMION
 30   COSTH=HWRUNI(0,-ONE,ONE)
      IF (IHEL.EQ.1 .AND. 0.25*(1+COSTH)**2.LT.HWRGEN(0)) GOTO 30
      IF (IHEL.EQ.2 .AND.      (1-COSTH**2).LT.HWRGEN(0)) GOTO 30
      IF (IHEL.EQ.3 .AND. 0.25*(1-COSTH)**2.LT.HWRGEN(0)) GOTO 30
C---GENERATE DECAY RELATIVE TO Z-AXIS
      PHEP(5,NHEP+1)=RMASS(IDN(1))
      PHEP(5,NHEP+2)=RMASS(IDN(2))
      PCM=HWUPCM(PHEP(5,IBOS),PHEP(5,NHEP+1),PHEP(5,NHEP+2))
      IF (PCM.LT.0) CALL HWWARN('HWDBOS',102,*999)
      CALL HWRAZM(PCM*SQRT(1-COSTH**2),PHEP(1,NHEP+1),PHEP(2,NHEP+1))
      PHEP(3,NHEP+1)=PCM*COSTH
      PHEP(4,NHEP+1)=SQRT(PHEP(5,NHEP+1)**2+PCM**2)
C---ROTATE SO THAT Z-AXIS BECOMES BOSON'S DIRECTION IN ORIGINAL CM FRAME
      CALL HWULOF(PHEP(1,ICMF),PHEP(1,IBOS),PBOS)
      CALL HWUROT(PBOS, ONE,ZERO,R)
      CALL HWUROB(R,PHEP(1,NHEP+1),PHEP(1,NHEP+1))
C---BOOST BACK TO LAB
      CALL HWULOB(PHEP(1,IBOS),PHEP(1,NHEP+1),PHEP(1,NHEP+1))
      CALL HWVDIF(4,PHEP(1,IBOS),PHEP(1,NHEP+1),PHEP(1,NHEP+2))
      ENDIF
C---STATUS, IDs AND POINTERS
      ISTHEP(IBOS)=195
      DO 50 I=1,2
        ISTHEP(NHEP+I)=193
        IDHW(NHEP+I)=IDN(I)
        IDHEP(NHEP+I)=IDPDG(IDN(I))
        JDAHEP(I,IBOS)=NHEP+I
        JMOHEP(1,NHEP+I)=IBOS
        JMOHEP(2,NHEP+I)=JMOHEP(1,IBOS)
 50   CONTINUE
      NHEP=NHEP+2
      IF (IDN(1).LE.12) THEN
        ISTHEP(NHEP-1)=113
        ISTHEP(NHEP)=114
        JMOHEP(2,NHEP)=NHEP-1
        JDAHEP(2,NHEP)=NHEP-1
        JMOHEP(2,NHEP-1)=NHEP
        JDAHEP(2,NHEP-1)=NHEP
        QUARKS=.TRUE.
      ENDIF
C---IF FIRST OF A PAIR, DO SECOND DECAY
      IF (IPAIR.NE.0 .AND. IPAIR.NE.IBOS) THEN
        IBOS=IPAIR
        GOTO 10
      ENDIF
C---IF QUARK DECAY, HADRONIZE
      IF (QUARKS) THEN
        EMSCA=PHEP(5,IBOS)
        CALL HWBGEN
        CALL HWDHQK
        CALL HWCFOR
        CALL HWCDEC
      ENDIF
 999  END