CDECK  ID>, HWDHIG.
*CMZ :-        -24/04/92  14.23.44  by  Mike Seymour
*-- Author :    Mike Seymour
C-----------------------------------------------------------------------
      SUBROUTINE HWDHIG(GAMINP)
C-----------------------------------------------------------------------
C     HIGGS DECAY ROUTINE
C     A) FOR GAMinp=0 FIND AND DECAY HIGGS
C     B) FOR GAMinp>0 CALCULATE TOTAL HIGGS WIDTH
C                     FOR EMH=GAMINP. STORE RESULT IN GAMINP.
C-----------------------------------------------------------------------
      INCLUDE 'HERWIG59.INC'
      DOUBLE PRECISION HWDHGF,HWRGEN,HWRUNI,HWUSQR,HWUPCM,GAMINP,EMH,
     & EMF,COLFAC,ENF,K1,K0,BET0,BET1,GAM0,GAM1,SCLOG,CFAC,XF,EM,GAMLIM,
     & GAM,XW,EMW,XZ,EMZ,YW,YZ,EMI,TAUT,TAUW,WIDHIG,VECDEC,EMB,GAMB,
     & TMIN,TMAX1,EM1,TMAX2,EM2,X1,X2,PROB,PCM,SUMR,SUMI,TAUTR,TAUTI,
     & TAUWR,TAUWI,GFACTR
      INTEGER HWRINT,IHIG,I,IFERM,NLOOK,I1,I2,IPART,IMODE,IDEC,MMAX
      LOGICAL HWRLOG
      EXTERNAL HWDHGF,HWRGEN,HWRUNI,HWUSQR,HWUPCM,HWRINT,HWRLOG
      SAVE GAM,EM,VECDEC
      PARAMETER (NLOOK=100)
      DIMENSION VECDEC(2,0:NLOOK)
      EQUIVALENCE (EMW,RMASS(198)),(EMZ,RMASS(200))
      DATA GAMLIM,GAM,EM/10D0,2*0D0/
C---IF DECAY, FIND HIGGS (HWDHAD WILL HAVE GIVEN IT STATUS=1)
      IF (GAMINP.EQ.0) THEN
        IHIG=0
        DO 10 I=1,NHEP
 10       IF (IHIG.EQ.0.AND.IDHW(I).EQ.201.AND.ISTHEP(I).EQ.1) IHIG=I
        IF (IHIG.EQ.0) CALL HWWARN('HWDHIG',101,*999)
        EMH=PHEP(5,IHIG)
        IF (EMH.LE.0) CALL HWWARN('HWDHIG',102,*999)
        EMSCA=EMH
      ELSE
        EMH=GAMINP
        IF (EMH.LE.0) THEN
          GAMINP=0
          RETURN
        ENDIF
      ENDIF
C---CALCULATE BRANCHING FRACTIONS
C---FERMIONS
C---NLL CORRECTION TO QUARK DECAY RATE (HHG eq 2.6-9)
      ENF=0
      DO 1 I=1,6
 1      IF (2*RMASS(I).LT.EMH) ENF=ENF+1
      K1=5/PIFAC**2
      K0=3/(4*PIFAC**2)
      BET0=(11*CAFAC-2*ENF)/3
      BET1=(34*CAFAC**2-(10*CAFAC+6*CFFAC)*ENF)/3
      GAM0=-8
      GAM1=-404./3+40*ENF/9
      SCLOG=LOG(EMH**2/QCDLAM**2)
      CFAC=1 + ( K1/K0 - 2*GAM0 + GAM0*BET1/BET0**2*LOG(SCLOG)
     &       +   (GAM0*BET1-GAM1*BET0)/BET0**2) / (BET0*SCLOG)
      DO 100 IFERM=1,9
        IF (IFERM.LE.6) THEN
          EMF=RMASS(IFERM)
          XF=(EMF/EMH)**2
          COLFAC=FLOAT(NCOLO)
          IF (EMF.GT.QCDLAM)
     &      EMF=EMF*(LOG(EMH/QCDLAM)/LOG(EMF/QCDLAM))**(GAM0/(2*BET0))
        ELSE
          EMF=RMASS(107+IFERM*2)
          XF=(EMF/EMH)**2
          COLFAC=1
          CFAC=1
        ENDIF
        IF (XF.LT.0.25) THEN
        GFACTR=ALPHEM/(8.*SWEIN*EMW**2)
          BRHIG(IFERM)=COLFAC*GFACTR*EMH*EMF**2 * (1-4*XF)**1.5 * CFAC
        ELSE
          BRHIG(IFERM)=0
        ENDIF
 100  CONTINUE
C---W*W*/Z*Z*
      IF (ABS(EM-EMH).GE.GAMLIM*GAM) THEN
C---OFF EDGE OF LOOK-UP TABLE
        XW=(EMW/EMH)**2
        XZ=(EMZ/EMH)**2
        YW=EMW*GAMW/EMH**2
        YZ=EMZ*GAMZ/EMH**2
        BRHIG(10)=.50*GFACTR * EMH**3 * HWDHGF(XW,YW)
        BRHIG(11)=.25*GFACTR * EMH**3 * HWDHGF(XZ,YZ)
      ELSE
C---LOOK IT UP
        EMI=((EMH-EM)/(GAM*GAMLIM)+1)*NLOOK/2.0
        I1=INT(EMI)
        I2=INT(EMI+1)
        BRHIG(10)=.50*GFACTR * EMH**3 * ( VECDEC(1,I1)*(I2-EMI) +
     &                                    VECDEC(1,I2)*(EMI-I1) )
        BRHIG(11)=.25*GFACTR * EMH**3 * ( VECDEC(2,I1)*(I2-EMI) +
     &                                    VECDEC(2,I2)*(EMI-I1) )
      ENDIF
C---GAMMAGAMMA
      TAUT=(2*RMASS(6)/EMH)**2
      TAUW=(2*EMW/EMH)**2
      CALL HWDHGC(TAUT,TAUTR,TAUTI)
      CALL HWDHGC(TAUW,TAUWR,TAUWI)
      SUMR=4./3*(  - 2*TAUT*( 1 + (1-TAUT)*TAUTR ) ) * ENHANC(6)
     &         +(2 + 3*TAUW*( 1 + (2-TAUW)*TAUWR ) ) * ENHANC(10)
      SUMI=4./3*(  - 2*TAUT*(     (1-TAUT)*TAUTI ) ) * ENHANC(6)
     &         +(    3*TAUW*(     (2-TAUW)*TAUWI ) ) * ENHANC(10)
      BRHIG(12)=GFACTR*.03125*(ALPHEM/PIFAC)**2
     &         *EMH**3 * (SUMR**2 + SUMI**2)
      WIDHIG=0
      DO 200 IPART=1, 12
        IF (IPART.LT.12) BRHIG(IPART)=BRHIG(IPART)*ENHANC(IPART)**2
 200    WIDHIG=WIDHIG+BRHIG(IPART)
      IF (WIDHIG.EQ.0) CALL HWWARN('HWDHIG',103,*999)
      DO 300 IPART=1, 12
 300    BRHIG(IPART)=BRHIG(IPART)/WIDHIG
      IF (EM.NE.RMASS(201)) THEN
C---SET UP W*W*/Z*Z* LOOKUP TABLES
        EM=EMH
        GAM=WIDHIG
        GAMLIM=MAX(GAMLIM,GAMMAX)
        DO 400 I=0,NLOOK
          EMH=(I*2.0/NLOOK-1)*GAM*GAMLIM+EM
          XW=(EMW/EMH)**2
          XZ=(EMZ/EMH)**2
          YW=EMW*GAMW/EMH**2
          YZ=EMZ*GAMZ/EMH**2
          VECDEC(1,I)=HWDHGF(XW,YW)
          VECDEC(2,I)=HWDHGF(XZ,YZ)
 400    CONTINUE
        EMH=EM
      ENDIF
      IF (GAMINP.GT.0) THEN
        GAMINP=WIDHIG
        RETURN
      ENDIF
C---SEE IF USER SPECIFIED A DECAY MODE
      IMODE=MOD(IPROC,100)
C---IF NOT, CHOOSE ONE
      IF (IMODE.LT.1.OR.IMODE.GT.12) THEN
        MMAX=12
        IF (IMODE.LT.1) MMAX=6
 500    IMODE=HWRINT(1,MMAX)
        IF (BRHIG(IMODE).LT.HWRGEN(0)) GOTO 500
      ENDIF
C---SEE IF SPECIFIED DECAY IS POSSIBLE
      IF (BRHIG(IMODE).EQ.0) CALL HWWARN('HWDHIG',104,*999)
      IF (IMODE.LE.6) THEN
        IDEC=IMODE
      ELSEIF (IMODE.LE.9) THEN
        IDEC=107+IMODE*2
      ELSEIF (IMODE.EQ.10) THEN
        IDEC=198
      ELSEIF (IMODE.EQ.11) THEN
        IDEC=200
      ELSEIF (IMODE.EQ.12) THEN
        IDEC=59
      ENDIF
C---STATUS, IDs AND POINTERS
      ISTHEP(IHIG)=195
      DO 600 I=1,2
        ISTHEP(NHEP+I)=193
        IDHW(NHEP+I)=IDEC
        IDHEP(NHEP+I)=IDPDG(IDEC)
        JDAHEP(I,IHIG)=NHEP+I
        JMOHEP(1,NHEP+I)=IHIG
        JMOHEP(2,NHEP+I)=NHEP+(3-I)
        JDAHEP(2,NHEP+I)=NHEP+(3-I)
        PHEP(5,NHEP+I)=RMASS(IDEC)
        IDEC=IDEC+6
        IF (IDEC.EQ.204) IDEC=199
        IF (IDEC.EQ.206) IDEC=200
        IF (IDEC.EQ. 65) IDEC= 59
 600  CONTINUE
C---ALLOW W/Z TO BE OFF-SHELL
      IF (IMODE.EQ.10.OR.IMODE.EQ.11) THEN
        IF (IMODE.EQ.10) THEN
          EMB=EMW
          GAMB=GAMW
        ELSE
          EMB=EMZ
          GAMB=GAMZ
        ENDIF
C---STANDARD MASS DISTRIBUTION
 700    TMIN=ATAN(-EMB/GAMB)
        TMAX1=ATAN((EMH**2/EMB-EMB)/GAMB)
        EM1=HWUSQR(EMB*(GAMB*TAN(HWRUNI(0,TMIN,TMAX1))+EMB))
        TMAX2=ATAN(((EMH-EM1)**2/EMB-EMB)/GAMB)
        EM2=HWUSQR(EMB*(GAMB*TAN(HWRUNI(0,TMIN,TMAX2))+EMB))
        X1=(EM1/EMH)**2
        X2=(EM2/EMH)**2
C---CORRECT MASS DISTRIBUTION
        PROB=HWUSQR(1+X1**2+X2**2-2*X1-2*X2-2*X1*X2)
     &        * ((X1+X2-1)**2 + 8*X1*X2)
        IF (.NOT.HWRLOG(PROB)) GOTO 700
C---CALCULATE SPIN DENSITY MATRIX
        RHOHEP(1,NHEP+1)=4*X1*X2      / (8*X1*X2 + (X1+X2-1)**2)
        RHOHEP(2,NHEP+1)=(X1+X2-1)**2 / (8*X1*X2 + (X1+X2-1)**2)
        RHOHEP(3,NHEP+1)=RHOHEP(1,NHEP+1)
C---SYMMETRIZE DISTRIBUTIONS IN PARTICLES 1,2
        IF (HWRLOG(HALF)) THEN
          PHEP(5,NHEP+1)=EM1
          PHEP(5,NHEP+2)=EM2
        ELSE
          PHEP(5,NHEP+1)=EM2
          PHEP(5,NHEP+2)=EM1
        ENDIF
      ENDIF
C---DO DECAY
      PCM=HWUPCM(EMH,PHEP(5,NHEP+1),PHEP(5,NHEP+2))
      IF (PCM.LT.0) CALL HWWARN('HWDHIG',105,*999)
      CALL HWDTWO(PHEP(1,IHIG),PHEP(1,NHEP+1),PHEP(1,NHEP+2),
     &            PCM,TWO,.TRUE.)
      NHEP=NHEP+2
C---IF QUARK DECAY, HADRONIZE
      IF (IMODE.LE.6) THEN
        ISTHEP(NHEP-1)=113
        ISTHEP(NHEP)=114
        CALL HWBGEN
        CALL HWDHQK
        CALL HWCFOR
        CALL HWCDEC
      ENDIF
  999 END
