use TMath::Abs() instead of ambiguous abs().

[u/mrichter/AliRoot.git] / HERWIG / src / hwbsud.f

CDECK  ID>, HWBSUD.

*CMZ :-        -14/07/92  13.28.23  by  Mike Seymour

*-- Author :    Bryan Webber

C-----------------------------------------------------------------------

      SUBROUTINE HWBSUD

C-----------------------------------------------------------------------

C     COMPUTES (OR READS) TABLES OF SUDAKOV FORM FACTORS

C-----------------------------------------------------------------------

      INCLUDE 'HERWIG61.INC'

      DOUBLE PRECISION HWUGAU,HWBVMC,HWBSUG,HWBSU1,HWBSU2,G1,G2,QRAT,

     & QLAM,POWER,AFAC,QMIN,QFAC,QNOW,ZMIN,ZMAX,Q1,QCOLD,VGOLD,VQOLD,

     & RMOLD(6),ACOLD,ZLO,ZHI

      INTEGER IQ,IS,L1,L2,L,LL,I,INOLD,NQOLD,NSOLD,NCOLD,NFOLD,SDOLD

      EXTERNAL HWUGAU,HWBVMC,HWBSUG,HWBSU1,HWBSU2

      SAVE NQOLD,NSOLD,NCOLD,NFOLD,SDOLD,QCOLD,VGOLD,VQOLD,RMOLD,ACOLD

      COMMON/HWSINT/QRAT,QLAM

      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) GOTO 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

CDECK  ID>, HWBSUG.

*CMZ :-        -13/07/92  20.15.54  by  Mike Seymour

*-- Author :    Bryan Webber, modified by Mike Seymour

C-----------------------------------------------------------------------

      FUNCTION HWBSUG(ZLOG)

C-----------------------------------------------------------------------

C     Z TIMES INTEGRAND IN EXPONENT OF GLUON SUDAKOV FORM FACTOR

C-----------------------------------------------------------------------

      DOUBLE PRECISION HWBSUG,HWBSUL,Z,ZLOG,W

      EXTERNAL HWBSUL

      Z=EXP(ZLOG)

      W=Z*(1.-Z)

      HWBSUG=HWBSUL(Z)*(W-2.+1./W)*Z

      END