*
* $Id: hwbsud.F,v 1.1.1.1 1996/03/08 17:02:10 mclareni Exp $
*
* $Log: hwbsud.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>, HWBSUD.
*CMZ :-        -14/07/92  13.28.23  by  Mike Seymour
*-- Author :    Bryan Webber
C-----------------------------------------------------------------------
      SUBROUTINE HWBSUD
C     COMPUTES (OR READS) TABLES OF SUDAKOV FORM FACTORS
C-----------------------------------------------------------------------
C     14/07/92: MODIFIED TO USE JACOBEAN TRANSFORM INSTEAD OF
C               SUBTRACTION TO REMOVE DIVERGENCE IN HWBSU1 AND HWBSUG.
C               SUM OVER FLAVOUR THRESHOLDS TO IMPROVE CONVERGENCE.
C-----------------------------------------------------------------------
#include "herwig58/herwig58.inc"
      INTEGER IQ,IS,L1,L2,L,LL,I,INOLD,NQOLD,NSOLD,NCOLD,NFOLD,SDOLD
      DOUBLE PRECISION HWUGAU,HWBVMC,G1,G2,QRAT,QLAM,POWER,AFAC,QMIN,
     & QFAC,QNOW,ZMIN,ZMAX,Q1,QCOLD,VGOLD,VQOLD,RMOLD(6),
     & ACOLD,HWBSUG,HWBSU1,HWBSU2,ZLO,ZHI
      SAVE NQOLD,NSOLD,NCOLD,NFOLD,SDOLD,QCOLD,VGOLD,VQOLD,RMOLD,ACOLD
      COMMON/HWSINT/QRAT,QLAM
      EXTERNAL HWBSUG,HWBSU1,HWBSU2
      IF (LRSUD.EQ.0) THEN
        POWER=1./FLOAT(NQEV-1)
        AFAC=6.*CAFAC/BETAF
        QMIN=QG+QG
        QFAC=(1.1*QLIM/QMIN)**POWER
        SUD(1,1)=1.
        QEV(1,1)=QMIN
C--IS=1 FOR GLUON->GLUON+GLUON FORM FACTOR
        DO 10 IQ=2,NQEV
        QNOW=QFAC*QEV(IQ-1,1)
        QLAM=QNOW/QCDL3
        ZMIN=QG/QNOW
        QRAT=1./ZMIN
        G1=0
        DO 5 I=3,6
          ZLO=ZMIN
          ZHI=HALF
          IF (I.NE.6) ZLO=MAX(ZLO,QG/RMASS(I+1))
          IF (I.NE.3) ZHI=MIN(ZHI,QG/RMASS(I))
          IF (ZHI.GT.ZLO) G1=G1+HWUGAU(HWBSUG,LOG(ZLO),LOG(ZHI),ACCUR)
    5   CONTINUE
        SUD(IQ,1)=EXP(AFAC*G1)
   10   QEV(IQ,1)=QNOW
        AFAC=3.*CFFAC/BETAF
C--QUARK FORM FACTORS.
C--IS=2,3,4,5,6,7 FOR U/D,S,C,B,T,V
        DO 15 IS=2,NSUD
        Q1=HWBVMC(IS)
        IF (IS.EQ.7) Q1=HWBVMC(209)
        QMIN=Q1+QG
        IF (QMIN.GT.QLIM) GO TO 15
        QFAC=(1.1*QLIM/QMIN)**POWER
        SUD(1,IS)=1.
        QEV(1,IS)=QMIN
        DO 14 IQ=2,NQEV
        QNOW=QFAC*QEV(IQ-1,IS)
        QLAM=QNOW/QCDL3
        ZMIN=QG/QNOW
        QRAT=1./ZMIN
        ZMAX=QG/QMIN
        G1=0
        DO 12 I=3,6
          ZLO=ZMIN
          ZHI=ZMAX
          IF (I.NE.6) ZLO=MAX(ZLO,QG/RMASS(I+1))
          IF (I.NE.3) ZHI=MIN(ZHI,QG/RMASS(I))
          IF (ZHI.GT.ZLO) G1=G1+HWUGAU(HWBSU1,LOG(ZLO),LOG(ZHI),ACCUR)
   12   CONTINUE
        ZMIN=Q1/QNOW
        QRAT=1./ZMIN
        ZMAX=Q1/QMIN
        G2=0
        DO 13 I=3,6
          ZLO=ZMIN
          ZHI=ZMAX
          IF (I.NE.6) ZLO=MAX(ZLO,Q1/RMASS(I+1))
          IF (I.NE.3) ZHI=MIN(ZHI,Q1/RMASS(I))
          IF (ZHI.GT.ZLO) G2=G2+HWUGAU(HWBSU2,ZLO,ZHI,ACCUR)
   13   CONTINUE
        SUD(IQ,IS)=EXP(AFAC*(G1+G2))
   14   QEV(IQ,IS)=QNOW
   15   CONTINUE
        QCOLD=QCDLAM
        VGOLD=VGCUT
        VQOLD=VQCUT
        ACOLD=ACCUR
        INOLD=INTER
        NQOLD=NQEV
        NSOLD=NSUD
        NCOLD=NCOLO
        NFOLD=NFLAV
        SDOLD=SUDORD
        DO 16 IS=1,NSUD
   16   RMOLD(IS)=RMASS(IS)
      ELSE
        IF (LRSUD.GT.0) THEN
          IF (IPRINT.NE.0) WRITE (6,17) LRSUD
   17     FORMAT(10X,'READING SUDAKOV TABLE ON UNIT',I4)
          OPEN(UNIT=LRSUD,FORM='UNFORMATTED',STATUS='UNKNOWN')
          READ(UNIT=LRSUD) QCOLD,VGOLD,VQOLD,RMOLD,
     &       ACOLD,QEV,SUD,INOLD,NQOLD,NSOLD,NCOLD,NFOLD,SDOLD
          CLOSE(UNIT=LRSUD)
        ENDIF
C---CHECK THAT RELEVANT PARAMETERS ARE UNCHANGED
        IF (QCDLAM.NE.QCOLD) CALL HWWARN('HWBSUD',501,*999)
        IF (VGCUT .NE.VGOLD) CALL HWWARN('HWBSUD',502,*999)
        IF (VQCUT .NE.VQOLD) CALL HWWARN('HWBSUD',503,*999)
        IF (ACCUR .NE.ACOLD) CALL HWWARN('HWBSUD',504,*999)
        IF (INTER .NE.INOLD) CALL HWWARN('HWBSUD',505,*999)
        IF (NQEV  .NE.NQOLD) CALL HWWARN('HWBSUD',506,*999)
        IF (NSUD  .NE.NSOLD) CALL HWWARN('HWBSUD',507,*999)
        IF (NCOLO .NE.NCOLD) CALL HWWARN('HWBSUD',508,*999)
        IF (NFLAV .NE.NFOLD) CALL HWWARN('HWBSUD',509,*999)
        IF (SUDORD.NE.SDOLD) CALL HWWARN('HWBSUD',510,*999)
C---CHECK MASSES AND THAT TABLES ARE BIG ENOUGH FOR THIS RUN
        DO 18 IS=1,NSUD
          IF (RMASS(IS).NE.RMOLD(IS))
     &      CALL HWWARN('HWBSUD',510+IS,*999)
          IF (QEV(NQEV,IS).LT.QLIM.AND.HWBVMC(IS)+QG.LT.QLIM)
     &      CALL HWWARN('HWBSUD',500,*999)
   18   CONTINUE
      ENDIF
      IF (LWSUD.GT.0) THEN
        IF (IPRINT.NE.0) WRITE (6,19) LWSUD
   19   FORMAT(10X,'WRITING SUDAKOV TABLE ON UNIT',I4)
        OPEN (UNIT=LWSUD,FORM='UNFORMATTED',STATUS='UNKNOWN')
        WRITE(UNIT=LWSUD)  QCDLAM,VGCUT,VQCUT,(RMASS(I),I=1,6),
     &     ACCUR,QEV,SUD,INTER,NQEV,NSUD,NCOLO,NFLAV,SUDORD
        CLOSE(UNIT=LWSUD)
      ENDIF
      IF (IPRINT.GT.2) THEN
C--PRINT EXTRACTS FROM TABLES OF FORM FACTORS
        DO 40 IS=1,NSUD
        WRITE(6,20) IS,NQEV
   20   FORMAT(1H1//10X,'EXTRACT FROM TABLE OF SUDAKOV FORM FACTOR NO.',
     &  I2,' (',I5,' ACTUAL ENTRIES)'//10X,'SUD IS PROBABILITY THAT',
     &  ' PARTON WITH GIVEN UPPER LIMIT ON Q WILL REACH THRESHOLD',
     &  ' WITHOUT BRANCHING'///2X,8('      Q     SUD ')/)
        L2=NQEV/8
        L1=L2/32
        IF (L1.LT.1) L1=1
        DO 40 L=L1,L2,L1
        LL=L+7*L2
        WRITE(6,30) (QEV(I,IS),SUD(I,IS),I=L,LL,L2)
   30   FORMAT(2X,8(F9.2,F7.4))
   40   CONTINUE
        WRITE(6,50)
   50   FORMAT(1H1)
      ENDIF
  999 END