ExhumeInterface/src/HDecay.f

0001 C         Last modification on February 27th 1997 by M.S.
0002 C ==================================================================
0003 C ================= PROGRAM HDECAY: COMMENTS =======================
0004 C ==================================================================
0005 C
0006 C                       ***************
0007 C                       * VERSION 2.0 *
0008 C                       ***************
0009 C
0010 C
0011 C  This program calculates the total decay widths and the branching
0012 C  ratios of the C Standard Model Higgs boson (HSM) as well as those
0013 C  of the neutral (HL= the light CP-even, HH= the heavy CP-even, HA=
0014 C  the pseudoscalar) and the charged (HC) Higgs bosons of the Minimal
0015 C  Supersymmetric extension of the Standard Model (MSSM). It includes:
0016 C
0017 C - All the decay channels which are kinematically allowed and which
0018 C   have branching ratios larger than 10**(-4).
0019 C
0020 C - All QCD corrections to the fermionic and gluonic decay modes.
0021 C   Most of these corrections are mapped into running masses in a
0022 C   consistent way with some freedom for including high order terms.
0023 C
0024 C - Below--threshold three--body decays with off--shell top quarks
0025 C   or ONE off-shell gauge boson, as well as some decays with one
0026 C   off-shell Higgs boson in the MSSM.
0027 C
0028 C - Double off-shell decays: HSM,HL,HH --> W*W*,Z*Z* -->4 fermions,
0029 C   which could be important for Higgs masses close to MW or MZ.
0030 C
0031 C - In the MSSM, the radiative corrections with full squark mixing and
0032 C   uses the RG improved values of Higgs masses and couplings with the
0033 C   main NLO corrections implemented (based on M.Carena, M. Quiros and
0034 C   C.E.M. Wagner, hep-ph/9508343).
0035 C
0036 C - In the MSSM, all the decays into CHARGINOS, NEUTRALINOS, SLEPTONS
0037 C   and SQUARKS (with mixing in the stop and sbottom sectors).
0038 C
0039 C - Chargino, slepton and squark loops in the 2 photon decays and squark
0040 C   loops in the gluonic decays (including QCD corrections).
0041 C
0042 C  ===================================================================
0043 C  This peogram has been written by A.Djouadi, J.Kalinowski and M.Spira.
0044 C  For details on how to use the program see: hep-ph/96XXX. For any
0045 C  question, comment, suggestion or complaint, please contact us at:
0046
0047
0048 C ================ IT USES AS INPUT PARAMETERS:
0049 C
0050 C   IHIGGS: =0: CALCULATE BRANCHING RATIOS OF SM HIGGS BOSON
0051 C           =1: CALCULATE BRANCHING RATIOS OF MSSM h BOSON
0052 C           =2: CALCULATE BRANCHING RATIOS OF MSSM H BOSON
0053 C           =3: CALCULATE BRANCHING RATIOS OF MSSM A BOSON
0054 C           =4: CALCULATE BRANCHING RATIOS OF MSSM H+ BOSON
0055 C           =5: CALCULATE BRANCHING RATIOS OF ALL MSSM HIGGS BOSONS
0056 C
0057 C TGBET:    TAN(BETA) FOR MSSM
0058 C MABEG:    START VALUE OF M_A FOR MSSM AND M_H FOR SM
0059 C MAEND:    END VALUE OF M_A FOR MSSM AND M_H FOR SM
0060 C NMA:      NUMBER OF ITERATIONS FOR M_A
0061 C ALS(MZ):  VALUE FOR ALPHA_S(M_Z)
0062 C MSBAR(1): MSBAR MASS OF STRANGE QUARK AT SCALE Q=1 GEV
0063 C MC:       CHARM POLE MASS
0064 C MB:       BOTTOM POLE MASS
0065 C MT:       TOP POLE MASS
0066 C MTAU:     TAU MASS
0067 C MMUON:    MUON MASS
0068 C ALPH:     INVERSE QED COUPLING
0069 C GF:       FERMI CONSTANT
0070 C GAMW:     W WIDTH
0071 C GAMZ:     Z WIDTH
0072 C MZ:       Z MASS
0073 C MW:       W MASS
0074 C VUS:      CKM PARAMETER V_US
0075 C VCB:      CKM PARAMETER V_CB
0076 C VUB/VCB:  RATIO V_UB/V_CB
0077 C NNLO (M): =1: USE O(ALPHA_S) FORMULA FOR POLE MASS --> MSBAR MASS
0078 C           =2: USE O(ALPHA_S**2) FORMULA FOR POLE MASS --> MSBAR MASS
0079 C
0080 C ON-SHELL: =0: INCLUDE OFF_SHELL DECAYS H,A --> T*T*, A --> Z*H,
0081 C               H --> W*H+,Z*A, H+ --> W*A, W*H, T*B
0082 C           =1: EXCLUDE THE OFF-SHELL DECAYS ABOVE
0083 C
0084 C ON-SH-WZ: =0: INCLUDE DOUBLE OFF-SHELL PAIR DECAYS PHI --> W*W*,Z*Z*
0085 C           =1: INCLUDE ONLY SINGLE OFF-SHELL DECAYS PHI --> W*W,Z*Z
0086 C
0087 C IPOLE:    =0 COMPUTES RUNNING HIGGS MASSES (FASTER)
0088 C           =1 COMPUTES POLE HIGGS MASSES
0089 C
0090 C OFF-SUSY: =0: INCLUDE DECAYS (AND LOOPS) INTO SUPERSYMMETRIC PARTICLES
0091 C           =1: EXCLUDE DECAYS (AND LOOPS) INTO SUPERSYMMETRIC PARTICLES
0092 C
0093 C INIDEC:   =0: PRINT OUT SUMS OF CHARGINO/NEUTRALINO/SFERMION DECAYS
0094 C           =1: PRINT OUT INDIVIDUAL CHARGINO/NEUTRALINO/SFERMION DECAYS
0095 C
0096 C NF-GG:    NUMBER OF LIGHT FLAVORS INCLUDED IN THE GLUONIC DECAYS
0097 C            PHI --> GG* --> GQQ (3,4 OR 5)
0098
0099 C =====================================================================
0100 C =========== BEGINNING OF THE SUBROUTINE FOR THE DECAYS ==============
0101 C !!!!!!!!!!!!!! Any change below this line is at your own risk!!!!!!!!
0102 C =====================================================================
0103
0104 C fix unused TGBET warning
0105 C     SUBROUTINE HDEC(TGBET)
0106       SUBROUTINE HDEC()
0107       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
0108       DOUBLE PRECISION LAMB
0109       DIMENSION XX(4),YY(4)
0110       COMPLEX*16 CF,CG,CI1,CI2,CA,CB,CTT,CTB,CTC,CTW,CLT,CLB,CLW,
0111      .           CAT,CAB,CAC,CAW,CTL,CAL
0112       COMMON/HMASS/AMSM,AMA,AML,AMH,AMCH
0113       COMMON/MASSES/AMS,AMC,AMB,AMT
0114       COMMON/ALS/XLAMBDA,AMC0,AMB0,AMT0,N0
0115       COMMON/PARAM/GF,ALPH,AMTAU,AMMUON,AMZ,AMW
0116       COMMON/CKMPAR/VUS,VCB,VUB
0117       COMMON/WZWDTH/GAMC0,GAMT0,GAMT1,GAMW,GAMZ
0118       COMMON/ONSHELL/IONSH,IONWZ,IOFSUSY
0119       COMMON/OLDFASH/NFGG
0120       COMMON/FLAG/IHIGGS,NNLO,IPOLE
0121       COMMON/WIDTHSM/SMBRB,SMBRL,SMBRM,SMBRS,SMBRC,SMBRT,SMBRG,SMBRGA,
0122      .               SMBRZGA,SMBRW,SMBRZ,SMWDTH
0123       COMMON/COUP/GAT,GAB,GLT,GLB,GHT,GHB,GZAH,GZAL,
0124      .            GHHH,GLLL,GHLL,GLHH,GHAA,GLAA,GLVV,GHVV,
0125      .            GLPM,GHPM,B,A
0126       BETA(X)=DSQRT(1.D0-4.D0*X)
0127       LAMB(X,Y)=DSQRT((1.D0-X-Y)**2-4.D0*X*Y)
0128       HVV(X,Y)= GF/(4.D0*PI*DSQRT(2.D0))*X**3/2.D0*BETA(Y)
0129      .            *(1.D0-4.D0*Y+12.D0*Y**2)
0130       HFF(X,Y)= GF/(4*PI*DSQRT(2.D0))*X**3*Y*(BETA(Y))**3
0131       HV(V)=3.D0*(1.D0-8.D0*V+20.D0*V**2)/DSQRT((4.D0*V-1.D0))
0132      .      *DACOS((3.D0*V-1.D0)/2.D0/DSQRT(V**3))
0133      .      -(1.D0-V)*(47.D0/2.D0*V-13.D0/2.D0+1.D0/V)
0134      .      -3.D0/2.D0*(1.D0-6.D0*V+4.D0*V**2)*DLOG(V)
0135       QCD0(X) = (1+X**2)*(4*SP((1-X)/(1+X)) + 2*SP((X-1)/(X+1))
0136      .        - 3*DLOG((1+X)/(1-X))*DLOG(2/(1+X))
0137      .        - 2*DLOG((1+X)/(1-X))*DLOG(X))
0138      .        - 3*X*DLOG(4/(1-X**2)) - 4*X*DLOG(X)
0139       HQCDM(X)=QCD0(X)/X+(3+34*X**2-13*X**4)/16/X**3*DLOG((1+X)/(1-X))
0140      .        + 3.D0/8/X**2*(7*X**2-1)
0141       HQCD(X)=(4.D0/3*HQCDM(BETA(X))
0142      .        +2*(4.D0/3-DLOG(X))*(1-10*X)/(1-4*X))*ASH/PI
0143      .       + (29.14671D0 + RATCOUP*(1.570D0 - 2*DLOG(HIGTOP)/3
0144      .                                     + DLOG(X)**2/9))*(ASH/PI)**2
0145      .       + (164.14D0 - 25.77D0*5 + 0.259D0*5**2)*(ASH/PI)**3
0146       QCDH(X)=1.D0+HQCD(X)
0147       TQCDH(X)=1.D0+4.D0/3*HQCDM(BETA(X))*ASH/PI
0148 C
0149       XI(X,Y) = 2*X/(1-X-Y+LAMB(X,Y))
0150       BIJ(X,Y) = (1-X-Y)/LAMB(X,Y)*(
0151      .          4*SP(XI(X,Y)*XI(Y,X))
0152      .        - 2*SP(-XI(X,Y)) - 2*SP(-XI(Y,X))
0153      .        + 2*DLOG(XI(X,Y)*XI(Y,X))*DLOG(1-XI(X,Y)*XI(Y,X))
0154      .        - DLOG(XI(X,Y))*DLOG(1+XI(X,Y))
0155      .        - DLOG(XI(Y,X))*DLOG(1+XI(Y,X))
0156      .          )
0157      .        -4*(DLOG(1-XI(X,Y)*XI(Y,X))
0158      .      +XI(X,Y)*XI(Y,X)/(1-XI(X,Y)*XI(Y,X))*DLOG(XI(X,Y)*XI(Y,X)))
0159      .        + (LAMB(X,Y)+X-Y)/LAMB(X,Y)*(DLOG(1+XI(X,Y))
0160      .                 - XI(X,Y)/(1+XI(X,Y))*DLOG(XI(X,Y)))
0161      .        + (LAMB(X,Y)-X+Y)/LAMB(X,Y)*(DLOG(1+XI(Y,X))
0162      .                 - XI(Y,X)/(1+XI(Y,X))*DLOG(XI(Y,X)))
0163       ELW(AMH,AMF,QF,ACF)=ALPH/PI*3.D0/2*QF**2
0164      .                              *(3.D0/2-DLOG(AMH**2/AMF**2))
0165      .      +GF/8/DSQRT(2.D0)/PI**2*(ACF*AMT**2
0166      .        +AMW**2*(3*DLOG(CS)/SS-5)+AMZ**2*(0.5D0
0167      .          -3*(1-4*SS*DABS(QF))**2))
0168       CF(CA) = -CDLOG(-(1+CDSQRT(1-CA))/(1-CDSQRT(1-CA)))**2/4
0169       CG(CA) = CDSQRT(1-CA)/2*CDLOG(-(1+CDSQRT(1-CA))/(1-CDSQRT(1-CA)))
0170       CI1(CA,CB) = CA*CB/2/(CA-CB)
0171      .           + CA**2*CB**2/2/(CA-CB)**2*(CF(CA)-CF(CB))
0172      .           + CA**2*CB/(CA-CB)**2*(CG(CA)-CG(CB))
0173       CI2(CA,CB) = -CA*CB/2/(CA-CB)*(CF(CA)-CF(CB))
0174       HGGQCD(ASG,NF)=1.D0+ASG/PI*(95.D0/4.D0-NF*7.D0/6.D0)
0175       HFFSELF(AMH)=1.D0+GF*AMH**2/16.D0/PI**2/DSQRT(2.D0)*2.117203D0
0176      .            -(GF*AMH**2/16.D0/PI**2/DSQRT(2.D0))**2*32.6567D0
0177       HVVSELF(AMH)=1.D0+GF*AMH**2/16.D0/PI**2/DSQRT(2.D0)*2.800952D0
0178      .            +(GF*AMH**2/16.D0/PI**2/DSQRT(2.D0))**2*62.0308D0
0179
0180       PI=4D0*DATAN(1D0)
0181       SS=1.D0-(AMW/AMZ)**2
0182
0183       CS=1.D0-SS
0184
0185 C--DECOUPLING THE TOP QUARK FROM ALPHAS
0186       AMT0=3.D8
0187
0188 C        =========================================================
0189 C                              SM HIGGS DECAYS
0190 C        =========================================================
0191       AMXX=AMH
0192       AMH=AMSM
0193 C     =============  RUNNING MASSES
0194       RMS = RUNM(AMH,3)
0195       RMC = RUNM(AMH,4)
0196       RMB = RUNM(AMH,5)
0197       RMT = RUNM(AMH,6)
0198       RATCOUP = 1
0199       HIGTOP = AMH**2/AMT**2
0200
0201       ASH=ALPHAS(AMH,2)
0202       AMC0=1.D8
0203       AMB0=2.D8
0204       AS3=ALPHAS(AMH,2)
0205       AMC0=AMC
0206       AS4=ALPHAS(AMH,2)
0207       AMB0=AMB
0208 C     AMT0=AMT
0209 C     =============== PARTIAL WIDTHS
0210 C  H ---> G G
0211 C
0212        EPS=1.D-8
0213        NFEXT = 3
0214        ASG = AS3
0215        CTT = 4*AMT**2/AMH**2*DCMPLX(1D0,-EPS)
0216        CTB = 4*AMB**2/AMH**2*DCMPLX(1D0,-EPS)
0217        CAT = 2*CTT*(1+(1-CTT)*CF(CTT))
0218        CAB = 2*CTB*(1+(1-CTB)*CF(CTB))
0219        FQCD=HGGQCD(ASG,NFEXT)
0220        XFAC = CDABS(CAT+CAB)**2*FQCD
0221        HGG=HVV(AMH,0.D0)*(ASG/PI)**2*XFAC/8
0222 C  H ---> G G* ---> G CC   TO BE ADDED TO H ---> CC
0223        NFEXT = 4
0224        ASG = AS4
0225        FQCD=HGGQCD(ASG,NFEXT)
0226        XFAC = CDABS(CAT+CAB)**2*FQCD
0227        DCC=HVV(AMH,0.D0)*(ASG/PI)**2*XFAC/8 - HGG
0228
0229 C  H ---> G G* ---> G BB   TO BE ADDED TO H ---> BB
0230        NFEXT = 5
0231        ASG = ASH
0232        FQCD=HGGQCD(ASG,NFEXT)
0233        XFAC = CDABS(CAT+CAB)**2*FQCD
0234        DBB=HVV(AMH,0.D0)*(ASG/PI)**2*XFAC/8 - HGG - DCC
0235
0236       IF(NFGG.EQ.5)THEN
0237        HGG = HGG + DBB + DCC
0238        DBB = 0
0239        DCC = 0
0240       ELSEIF(NFGG.EQ.4)THEN
0241        HGG = HGG + DCC
0242        DCC = 0
0243       ENDIF
0244
0245 C  H ---> MU MU
0246       IF(AMH.LE.2*AMMUON) THEN
0247        HMM = 0
0248       ELSE
0249       HMM=HFF(AMH,(AMMUON/AMH)**2)
0250      .    *(1+ELW(AMH,AMMUON,-1.D0,7.D0))
0251      .    *HFFSELF(AMH)
0252       ENDIF
0253 C  H ---> TAU TAU
0254       IF(AMH.LE.2*AMTAU) THEN
0255        HLL = 0
0256       ELSE
0257       HLL=HFF(AMH,(AMTAU/AMH)**2)
0258      .    *(1+ELW(AMH,AMTAU,-1.D0,7.D0))
0259      .    *HFFSELF(AMH)
0260       ENDIF
0261 C  H --> SS
0262
0263       IF(AMH.LE.2*AMS) THEN
0264        HSS = 0
0265       ELSE
0266        HS2=3.D0*HFF(AMH,(RMS/AMH)**2)
0267      .    *QCDH(RMS**2/AMH**2)
0268      .    *(1+ELW(AMH,RMS,-1.D0/3.D0,7.D0))
0269      .    *HFFSELF(AMH)
0270        IF(HS2.LT.0.D0) HS2 = 0
0271        HS1=3.D0*HFF(AMH,(AMS/AMH)**2)
0272      .    *TQCDH(AMS**2/AMH**2)
0273      .    *HFFSELF(AMH)
0274        RAT = 2*AMS/AMH
0275        HSS = QQINT(RAT,HS1,HS2)
0276       ENDIF
0277 C  H --> CC
0278       IF(AMH.LE.2*AMC) THEN
0279        HCC = 0
0280       ELSE
0281        HC2=3.D0*HFF(AMH,(RMC/AMH)**2)
0282      .    *QCDH(RMC**2/AMH**2)
0283      .    *(1+ELW(AMH,RMC,2.D0/3.D0,7.D0))
0284      .    *HFFSELF(AMH)
0285      .   + DCC
0286        IF(HC2.LT.0.D0) HC2 = 0
0287        HC1=3.D0*HFF(AMH,(AMC/AMH)**2)
0288      .    *TQCDH(AMC**2/AMH**2)
0289      .    *HFFSELF(AMH)
0290        RAT = 2*AMC/AMH
0291        HCC = QQINT(RAT,HC1,HC2)
0292       ENDIF
0293 C  H --> BB :
0294       IF(AMH.LE.2*AMB) THEN
0295        HBB = 0
0296       ELSE
0297        HB2=3.D0*HFF(AMH,(RMB/AMH)**2)
0298      .    *QCDH(RMB**2/AMH**2)
0299      .    *(1+ELW(AMH,RMB,-1.D0/3.D0,1.D0))
0300      .    *HFFSELF(AMH)
0301      .   + DBB
0302        IF(HB2.LT.0.D0) HB2 = 0
0303        HB1=3.D0*HFF(AMH,(AMB/AMH)**2)
0304      .    *TQCDH(AMB**2/AMH**2)
0305      .    *HFFSELF(AMH)
0306        RAT = 2*AMB/AMH
0307        HBB = QQINT(RAT,HB1,HB2)
0308       ENDIF
0309 C  H ---> TT
0310       RATCOUP = 0
0311       IF(IONSH.EQ.0)THEN
0312        DLD=3D0
0313        DLU=5D0
0314        XM1 = 2D0*AMT-DLD
0315        XM2 = 2D0*AMT+DLU
0316        IF (AMH.LE.AMT+AMW+AMB) THEN
0317        HTT=0.D0
0318        ELSEIF (AMH.LE.XM1) THEN
0319         FACTT=6.D0*GF**2*AMH**3*AMT**2/2.D0/128.D0/PI**3
0320         CALL HTOTTS(AMH,AMT,AMB,AMW,HTTS)
0321         HTT=FACTT*HTTS
0322        ELSEIF (AMH.LE.XM2) THEN
0323         XX(1) = XM1-1D0
0324         XX(2) = XM1
0325         XX(3) = XM2
0326         XX(4) = XM2+1D0
0327         FACTT=6.D0*GF**2*XX(1)**3*AMT**2/2.D0/128.D0/PI**3
0328         CALL HTOTTS(XX(1),AMT,AMB,AMW,HTTS)
0329         YY(1)=FACTT*HTTS
0330         FACTT=6.D0*GF**2*XX(2)**3*AMT**2/2.D0/128.D0/PI**3
0331         CALL HTOTTS(XX(2),AMT,AMB,AMW,HTTS)
0332         YY(2)=FACTT*HTTS
0333         XMT = RUNM(XX(3),6)
0334         XY2=3.D0*HFF(XX(3),(XMT/XX(3))**2)
0335      .    *QCDH(XMT**2/XX(3)**2)
0336      .    *HFFSELF(XX(3))
0337         IF(XY2.LT.0.D0) XY2 = 0
0338         XY1=3.D0*HFF(XX(3),(AMT/XX(3))**2)
0339      .    *TQCDH(AMT**2/XX(3)**2)
0340      .    *HFFSELF(XX(3))
0341         RAT = 2*AMT/XX(3)
0342         YY(3) = QQINT(RAT,XY1,XY2)
0343         XMT = RUNM(XX(4),6)
0344         XY2=3.D0*HFF(XX(4),(XMT/XX(4))**2)
0345      .    *QCDH(XMT**2/XX(4)**2)
0346      .    *HFFSELF(XX(4))
0347         IF(XY2.LT.0.D0) XY2 = 0
0348         XY1=3.D0*HFF(XX(4),(AMT/XX(4))**2)
0349      .    *TQCDH(AMT**2/XX(4)**2)
0350      .    *HFFSELF(XX(4))
0351         RAT = 2*AMT/XX(4)
0352         YY(4) = QQINT(RAT,XY1,XY2)
0353         HTT = FINT_(AMH,XX,YY)
0354        ELSE
0355         HT2=3.D0*HFF(AMH,(RMT/AMH)**2)
0356      .    *QCDH(RMT**2/AMH**2)
0357      .    *HFFSELF(AMH)
0358         IF(HT2.LT.0.D0) HT2 = 0
0359         HT1=3.D0*HFF(AMH,(AMT/AMH)**2)
0360      .    *TQCDH(AMT**2/AMH**2)
0361      .    *HFFSELF(AMH)
0362         RAT = 2*AMT/AMH
0363         HTT = QQINT(RAT,HT1,HT2)
0364        ENDIF
0365       ELSE
0366        IF (AMH.LE.2.D0*AMT) THEN
0367         HTT=0.D0
0368        ELSE
0369         HT2=3.D0*HFF(AMH,(RMT/AMH)**2)
0370      .    *QCDH(RMT**2/AMH**2)
0371      .    *HFFSELF(AMH)
0372         IF(HT2.LT.0.D0) HT2 = 0
0373         HT1=3.D0*HFF(AMH,(AMT/AMH)**2)
0374      .    *TQCDH(AMT**2/AMH**2)
0375      .    *HFFSELF(AMH)
0376         RAT = 2*AMT/AMH
0377         HTT = QQINT(RAT,HT1,HT2)
0378        ENDIF
0379       ENDIF
0380 C  H ---> GAMMA GAMMA
0381        EPS=1.D-8
0382        CTT = 4*AMT**2/AMH**2*DCMPLX(1D0,-EPS)
0383        CTB = 4*AMB**2/AMH**2*DCMPLX(1D0,-EPS)
0384        CTC = 4*AMC**2/AMH**2*DCMPLX(1D0,-EPS)
0385        CTL = 4*AMTAU**2/AMH**2*DCMPLX(1D0,-EPS)
0386        CTW = 4*AMW**2/AMH**2*DCMPLX(1D0,-EPS)
0387        CAW = -(2+3*CTW+3*CTW*(2-CTW)*CF(CTW))
0388        CAT = 4/3D0 * 2*CTT*(1+(1-CTT)*CF(CTT))
0389        CAB = 1/3D0 * 2*CTB*(1+(1-CTB)*CF(CTB))
0390        CAC = 4/3D0 * 2*CTC*(1+(1-CTC)*CF(CTC))
0391        CAL =         2*CTL*(1+(1-CTL)*CF(CTL))
0392        XFAC = CDABS(CAT+CAB+CAC+CAL+CAW)**2
0393        HGA=HVV(AMH,0.D0)*(ALPH/PI)**2/16.D0*XFAC
0394
0395 C  H ---> Z GAMMA
0396       IF(AMH.LE.AMZ)THEN
0397        HZGA=0
0398       ELSE
0399        EPS=1.D-8
0400        TS = SS/CS
0401        FT = -3*2D0/3*(1-4*2D0/3*SS)/DSQRT(SS*CS)
0402        FB = 3*1D0/3*(-1+4*1D0/3*SS)/DSQRT(SS*CS)
0403        CTT = 4*AMT**2/AMH**2*DCMPLX(1D0,-EPS)
0404        CTB = 4*AMB**2/AMH**2*DCMPLX(1D0,-EPS)
0405        CTW = 4*AMW**2/AMH**2*DCMPLX(1D0,-EPS)
0406        CLT = 4*AMT**2/AMZ**2*DCMPLX(1D0,-EPS)
0407        CLB = 4*AMB**2/AMZ**2*DCMPLX(1D0,-EPS)
0408        CLW = 4*AMW**2/AMZ**2*DCMPLX(1D0,-EPS)
0409        CAT = FT*(CI1(CTT,CLT) - CI2(CTT,CLT))
0410        CAB = FB*(CI1(CTB,CLB) - CI2(CTB,CLB))
0411        CAW = -1/DSQRT(TS)*(4*(3-TS)*CI2(CTW,CLW)
0412      .     + ((1+2/CTW)*TS - (5+2/CTW))*CI1(CTW,CLW))
0413        XFAC = CDABS(CAT+CAB+CAW)**2
0414        ACOUP = DSQRT(2D0)*GF*AMZ**2*SS*CS/PI**2
0415        HZGA = GF/(4.D0*PI*DSQRT(2.D0))*AMH**3*(ALPH/PI)*ACOUP/16.D0
0416      .        *XFAC*(1-AMZ**2/AMH**2)**3
0417       ENDIF
0418
0419 C  H ---> W W
0420       IF(IONWZ.EQ.0)THEN
0421        DLD=2D0
0422        DLU=2D0
0423        XM1 = 2D0*AMW-DLD
0424        XM2 = 2D0*AMW+DLU
0425        IF (AMH.LE.XM1) THEN
0426         CALL HTOVV(AMH,AMW,GAMW,HTWW)
0427         HWW = 3D0/2D0*GF*AMW**4/DSQRT(2D0)/PI/AMH**3*HTWW
0428        ELSEIF (AMH.LE.XM2) THEN
0429         XX(1) = XM1-1D0
0430         XX(2) = XM1
0431         XX(3) = XM2
0432         XX(4) = XM2+1D0
0433         CALL HTOVV(XX(1),AMW,GAMW,HTWW)
0434         YY(1)=3D0/2D0*GF*AMW**4/DSQRT(2D0)/PI/XX(1)**3*HTWW
0435         CALL HTOVV(XX(2),AMW,GAMW,HTWW)
0436         YY(2)=3D0/2D0*GF*AMW**4/DSQRT(2D0)/PI/XX(2)**3*HTWW
0437         YY(3)=HVV(XX(3),AMW**2/XX(3)**2)
0438      .       *HVVSELF(XX(3))
0439         YY(4)=HVV(XX(4),AMW**2/XX(4)**2)
0440      .       *HVVSELF(XX(4))
0441         HWW = FINT_(AMH,XX,YY)
0442        ELSE
0443         HWW=HVV(AMH,AMW**2/AMH**2)
0444      .     *HVVSELF(AMH)
0445        ENDIF
0446       ELSE
0447        DLD=2D0
0448        DLU=2D0
0449        XM1 = 2D0*AMW-DLD
0450        XM2 = 2D0*AMW+DLU
0451       IF (AMH.LE.AMW) THEN
0452        HWW=0
0453       ELSE IF (AMH.LE.XM1) THEN
0454        CWW=3.D0*GF**2*AMW**4/16.D0/PI**3
0455        HWW=HV(AMW**2/AMH**2)*CWW*AMH
0456       ELSE IF (AMH.LT.XM2) THEN
0457        CWW=3.D0*GF**2*AMW**4/16.D0/PI**3
0458        XX(1) = XM1-1D0
0459        XX(2) = XM1
0460        XX(3) = XM2
0461        XX(4) = XM2+1D0
0462        YY(1)=HV(AMW**2/XX(1)**2)*CWW*XX(1)
0463        YY(2)=HV(AMW**2/XX(2)**2)*CWW*XX(2)
0464        YY(3)=HVV(XX(3),AMW**2/XX(3)**2)
0465      .      *HVVSELF(XX(3))
0466        YY(4)=HVV(XX(4),AMW**2/XX(4)**2)
0467      .      *HVVSELF(XX(4))
0468        HWW = FINT_(AMH,XX,YY)
0469       ELSE
0470        HWW=HVV(AMH,AMW**2/AMH**2)
0471      .     *HVVSELF(AMH)
0472       ENDIF
0473       ENDIF
0474
0475 C  H ---> Z Z
0476       IF(IONWZ.EQ.0)THEN
0477        DLD=2D0
0478        DLU=2D0
0479        XM1 = 2D0*AMZ-DLD
0480        XM2 = 2D0*AMZ+DLU
0481        IF (AMH.LE.XM1) THEN
0482         CALL HTOVV(AMH,AMZ,GAMZ,HTZZ)
0483         HZZ = 3D0/4D0*GF*AMZ**4/DSQRT(2D0)/PI/AMH**3*HTZZ
0484        ELSEIF (AMH.LE.XM2) THEN
0485         XX(1) = XM1-1D0
0486         XX(2) = XM1
0487         XX(3) = XM2
0488         XX(4) = XM2+1D0
0489         CALL HTOVV(XX(1),AMZ,GAMZ,HTZZ)
0490         YY(1)=3D0/4D0*GF*AMZ**4/DSQRT(2D0)/PI/XX(1)**3*HTZZ
0491         CALL HTOVV(XX(2),AMZ,GAMZ,HTZZ)
0492         YY(2)=3D0/4D0*GF*AMZ**4/DSQRT(2D0)/PI/XX(2)**3*HTZZ
0493         YY(3)=HVV(XX(3),AMZ**2/XX(3)**2)/2
0494      .       *HVVSELF(XX(3))
0495         YY(4)=HVV(XX(4),AMZ**2/XX(4)**2)/2
0496      .       *HVVSELF(XX(4))
0497         HZZ = FINT_(AMH,XX,YY)
0498        ELSE
0499         HZZ=HVV(AMH,AMZ**2/AMH**2)/2.D0
0500      .      *HVVSELF(AMH)
0501        ENDIF
0502       ELSE
0503       DLD=2D0
0504       DLU=2D0
0505       XM1 = 2D0*AMZ-DLD
0506       XM2 = 2D0*AMZ+DLU
0507       IF (AMH.LE.AMZ) THEN
0508       HZZ=0
0509       ELSE IF (AMH.LE.XM1) THEN
0510       CZZ=3.D0*GF**2*AMZ**4/192.D0/PI**3*(7-40/3.D0*SS+160/9.D0*SS**2)
0511       HZZ=HV(AMZ**2/AMH**2)*CZZ*AMH
0512       ELSE IF (AMH.LT.XM2) THEN
0513       CZZ=3.D0*GF**2*AMZ**4/192.D0/PI**3*(7-40/3.D0*SS+160/9.D0*SS**2)
0514       XX(1) = XM1-1D0
0515       XX(2) = XM1
0516       XX(3) = XM2
0517       XX(4) = XM2+1D0
0518       YY(1)=HV(AMZ**2/XX(1)**2)*CZZ*XX(1)
0519       YY(2)=HV(AMZ**2/XX(2)**2)*CZZ*XX(2)
0520       YY(3)=HVV(XX(3),AMZ**2/XX(3)**2)/2
0521      .      *HVVSELF(XX(3))
0522       YY(4)=HVV(XX(4),AMZ**2/XX(4)**2)/2
0523      .      *HVVSELF(XX(4))
0524       HZZ = FINT_(AMH,XX,YY)
0525       ELSE
0526       HZZ=HVV(AMH,AMZ**2/AMH**2)/2.D0
0527      .   *HVVSELF(AMH)
0528       ENDIF
0529       ENDIF
0530
0531 C
0532 C    ==========  TOTAL WIDTH AND BRANCHING RATIOS
0533 C
0534       WTOT=HLL+HMM+HSS+HCC+HBB+HTT+HGG+HGA+HZGA+HWW+HZZ
0535
0536       SMBRT=HTT/WTOT
0537       SMBRB=HBB/WTOT
0538       SMBRL=HLL/WTOT
0539       SMBRM=HMM/WTOT
0540       SMBRC=HCC/WTOT
0541       SMBRS=HSS/WTOT
0542       SMBRG=HGG/WTOT
0543       SMBRGA=HGA/WTOT
0544       SMBRZGA=HZGA/WTOT
0545       SMBRW=HWW/WTOT
0546       SMBRZ=HZZ/WTOT
0547       SMWDTH=WTOT
0548
0549       AMH=AMXX
0550
0551       RETURN
0552       END
0553 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
0554
0555 C     ************************************************************
0556 C     SUBROUTINE FOR HSM ---> V*V* ---> 4F
0557 C     ************************************************************
0558       SUBROUTINE HTOVV(AMH,AMV,GAMV,HTVV)
0559       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
0560       COMMON/VVOFF/AMH1,AMV1,GAMV1
0561       COMMON/PREC/IP
0562       EXTERNAL FTOVV1
0563       IP=20
0564       AMH1=AMH
0565       AMV1=AMV
0566       GAMV1=GAMV
0567       DLT=1D0/IP
0568       SUM=0D0
0569       DO 1 I=1,IP
0570        UU=DLT*I
0571        DD=UU-DLT
0572        CALL QGAUS1(FTOVV1,DD,UU,RES)
0573        SUM=SUM+RES
0574 1     CONTINUE
0575       HTVV=SUM
0576       RETURN
0577       END
0578
0579       DOUBLE PRECISION FUNCTION FTOVV1(XX)
0580       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
0581       COMMON/FIRST/X1
0582       COMMON/PREC/IP
0583       EXTERNAL FTOVV2
0584       X1=XX
0585       DLT=1D0/IP
0586       SUM=0D0
0587       DO 1 I=1,IP
0588        UU=DLT*I
0589        DD=UU-DLT
0590        CALL QGAUS2(FTOVV2,DD,UU,RES)
0591        SUM=SUM+RES
0592 1     CONTINUE
0593       FTOVV1=SUM
0594       RETURN
0595       END
0596
0597       DOUBLE PRECISION FUNCTION FTOVV2(XX)
0598       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
0599       DIMENSION YY(2)
0600       COMMON/FIRST/X1
0601       YY(1)=X1
0602       YY(2)=XX
0603       FTOVV2=FTOVV(YY)
0604       RETURN
0605       END
0606
0607       DOUBLE PRECISION FUNCTION FTOVV(XX)
0608       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
0609       DOUBLE PRECISION LAMB
0610       DIMENSION XX(2)
0611       COMMON/VVOFF/AMH,AMW,GAMW
0612       LAMB(X,Y)=DSQRT((1.D0-X-Y)**2-4.D0*X*Y)
0613       PI=4D0*DATAN(1D0)
0614       ICASE = 1
0615       IF(ICASE.EQ.0)THEN
0616        YY = AMH**2
0617        Y1 = DATAN((YY-AMW**2)/AMW/GAMW)
0618        Y2 = -DATAN((AMW**2)/AMW/GAMW)
0619        DJAC = Y1-Y2
0620        T1 = TAN(Y1*XX(1)+Y2*(1.D0-XX(1)))
0621        SP = AMW**2 + AMW*GAMW*T1
0622        YY = (AMH-DSQRT(SP))**2
0623        Y1 = DATAN((YY-AMW**2)/AMW/GAMW)
0624        Y2 = -DATAN((AMW**2)/AMW/GAMW)
0625        DJAC = DJAC*(Y1-Y2)
0626        T2 = TAN(Y1*XX(2)+Y2*(1.D0-XX(2)))
0627        SM = AMW**2 + AMW*GAMW*T2
0628        AM2=AMH**2
0629        GAM = AM2*LAMB(SP/AM2,SM/AM2)*(1+LAMB(SP/AM2,SM/AM2)**2*AMH**4
0630      .                               /SP/SM/12)
0631        PRO1 = SP/AMW**2
0632        PRO2 = SM/AMW**2
0633        FTOVV = PRO1*PRO2*GAM*DJAC/PI**2
0634       ELSE
0635        SP = AMH**2*XX(1)
0636        SM = (AMH-DSQRT(SP))**2*XX(2)
0637        DJAC = AMH**2*(AMH-DSQRT(SP))**2/PI**2
0638        AM2=AMH**2
0639        GAM = AM2*LAMB(SP/AM2,SM/AM2)*(1+LAMB(SP/AM2,SM/AM2)**2*AMH**4
0640      .                               /SP/SM/12)
0641        PRO1 = SP*GAMW/AMW/((SP-AMW**2)**2+AMW**2*GAMW**2)
0642        PRO2 = SM*GAMW/AMW/((SM-AMW**2)**2+AMW**2*GAMW**2)
0643        FTOVV = PRO1*PRO2*GAM*DJAC
0644       ENDIF
0645       RETURN
0646       END
0647
0648 C     ************************************************************
0649 C     SUBROUTINE FOR HSM ---> TT* ---> TBW
0650 C     ************************************************************
0651       SUBROUTINE HTOTTS(AMH,AMT,AMB,AMW,HTTS)
0652       IMPLICIT REAL*8(A-Z)
0653       INTEGER IP,K
0654       COMMON/PREC1/IP
0655       EXTERNAL FUNSTT1
0656       COMMON/IKSY0/X1,X2,M1,M2,M3,ECM,S
0657       COMMON/TOP0/AMH0,AMT0,AMB0,AMW0
0658       AMH0=AMH
0659       AMT0=AMT
0660       AMB0=AMB
0661       AMW0=AMW
0662       IP=5
0663       M1=AMB
0664       M2=AMT
0665       M3=AMW
0666 C     FIRST INTEGRATE OVER X2, i.e. (1+3) SYSTEM
0667 C        CHECK WHETHER ENOUGH PHASE SPACE
0668       MASTOT=M1+M2+M3
0669       IF(MASTOT.GE.AMH) GOTO 12
0670       ECM=AMH
0671       S=ECM**2
0672       U1=(ECM-M2)**2
0673       D1=(M1+M3)**2
0674       U=(S-D1+M2**2)/s
0675       D=(S-U1+M2**2)/s
0676       DEL=(U-D)/IP
0677       U=D+DEL
0678       XSEC=0.D0
0679       DO K=1,IP
0680       CALL QGAUS1(FUNSTT1,D,U,SS)
0681       D=U
0682       U=D+DEL
0683       XSEC=XSEC+SS
0684       ENDDO
0685       HTTS=XSEC
0686 12    CONTINUE
0687       RETURN
0688       END
0689
0690       DOUBLE PRECISION FUNCTION FUNSTT1(XL)
0691       IMPLICIT REAL*8(A-Z)
0692       INTEGER IP,I
0693       COMMON/IKSY0/X1,X2,M1,M2,M3,ECM,S
0694       COMMON/PREC1/IP
0695       EXTERNAL FUNSTT2
0696       X2=XL
0697       S13=S-S*X2+M2**2
0698       TEM=2.D0*DSQRT(S13)
0699       E2S=(S-S13-M2**2)/TEM
0700       E3S=(S13+M3**2-M1**2)/TEM
0701 C     SECOND INTEGRAL OVER X1, i.e. (2+3) SYSTEM
0702       U1=(E2S+E3S)**2-(DSQRT(E2S**2-M2**2)-DSQRT(E3S**2-M3**2))**2
0703       D1=(E2S+E3S)**2-(DSQRT(E2S**2-M2**2)+DSQRT(E3S**2-M3**2))**2
0704       U=(S-D1+M1**2)/s
0705       D=(S-U1+M1**2)/s
0706       DEL=(U-D)/IP
0707       FUNSTT1=0.d0
0708       U=D+DEL
0709       DO I=1,IP
0710       CALL QGAUS2(FUNSTT2,D,U,SS)
0711       FUNSTT1=FUNSTT1+SS
0712       D=U
0713       U=D+DEL
0714       ENDDO
0715       RETURN
0716       END
0717
0718       DOUBLE PRECISION FUNCTION FUNSTT2(XK)
0719       IMPLICIT REAL*8(A-Z)
0720       COMMON/IKSY0/X1,X2,M1,M2,M3,ECM,S
0721       X1=XK
0722       CALL ELEMSTT(SS)
0723       FUNSTT2=SS
0724       RETURN
0725       END
0726
0727       SUBROUTINE ELEMSTT(RES)
0728       IMPLICIT REAL*8(A-Z)
0729       COMMON/IKSY0/X1,X2,M1,M2,M3,ECM,S
0730       COMMON/TOP0/AMH,AMT,AMB,AMW
0731       COMMON/WZWDTH/GAMC0,GAMT0,GAMT1,GAMW0,GAMZ0
0732       GAMT=GAMT0**2*AMT**2/AMH**4
0733       GAMW=GAMW0**2*AMW**2/AMH**4
0734       W=AMW**2/AMH**2
0735       T=AMT**2/AMH**2
0736       Y1=1-X2
0737       Y2=1-X1
0738       X0=2.D0-X1-X2
0739       W1=(1.D0-X2)
0740       W3=(1.-X1-X2)
0741       W11=1.D0/((1.D0-X2)**2+GAMT)
0742       W33=1.D0/(W3**2+GAMW**2)
0743       W13=W1*W3*W11*W33
0744
0745       R11=4*T*W-16.*T*W*Y1-4.*T*Y2*Y1+8.*T*Y1+32.*T*W**2-20
0746      . .*T*Y1**2+8.*W*Y2*Y1+4.*W*Y1**2-4.*Y2*Y1**2-16.*T**2*W-
0747      .  32.*T**2*Y1+4.*T**2-16.*T**3-8.*W**2+4.*Y1**2-4.*Y1**3
0748       R33=-4.*T*W+4.*T*W*Y2-2.*T*W*Y2*Y1+4.*T*W*Y1+T*W*Y2**2-
0749      .  3.*T*W*Y1**2+2.*T*Y2*Y1-3.*T*Y2*Y1**2+4.*T*W**2-4.*T*W**3
0750      .  +T*Y2**2-3.*T*Y2**2*Y1-T*Y2**3+T*Y1**2-T*Y1**3+4.*T**2
0751      .  *W-4.*T**2*W*Y2-4.*T**2*W*Y1-2.*T**2*Y2*Y1-4.*T**2*W**2-
0752      .  T**2*Y2**2-T**2*Y1**2+4.*W**2*Y2*Y1-8.*W**3*Y2-8.*W**3*Y1
0753      .  +4.*W**3+8.*W**4
0754       R13=8.*W-24.*T*W+16.*T*W*Y1 -4.*T*Y2+16.*T*Y2*Y1-4.*T*
0755      .  Y1+16.*T*W**2+4.*T*Y2**2+12.*T*Y1**2-8.*W*Y2-12.*W*Y2*Y1
0756      .  -8.*W*Y1+4.*W*Y1**2-4.*Y2*Y1+8.*Y2*Y1**2+16.*T**2*W+8.
0757      .  *T**2*Y2+8.*T**2*Y1+16.*W**2*Y2+24.*W**2*Y1+4.*Y2**2*Y1-
0758      .  32.*W**3-4.*Y1**2+4.*Y1**3
0759       RES=R11*W11+4.D0*R33*W33/T-2.D0*R13*W13
0760       RETURN
0761       END
0762
0763       DOUBLE PRECISION FUNCTION FUNHTT2(XK)
0764       IMPLICIT REAL*8(A-Z)
0765       COMMON/IKSY2/X1,X2,M1,M2,M3,ECM,S
0766       X1=XK
0767       CALL ELEMHTT(SS)
0768       FUNHTT2=SS
0769       RETURN
0770       END
0771
0772       SUBROUTINE ELEMHTT(RES)
0773       IMPLICIT REAL*8(A-Z)
0774       COMMON/IKSY2/X1,X2,M1,M2,M3,ECM,S
0775       COMMON/TOP2/AMH,AMT,AMB,AMW,AMCH,TB,GHT,GAT,GHVV
0776       COMMON/WZWDTH/GAMC0,GAMT0,GAMT1,GAMW0,GAMZ0
0777       GAMT=GAMT1**2*AMT**2/AMH**4
0778       GAMC=GAMC0**2*AMCH**2/AMH**4
0779       GAMW=GAMW0**2*AMW**2/AMH**4
0780       CH=AMCH**2/AMH**2
0781       W=AMW**2/AMH**2
0782       T=AMT**2/AMH**2
0783       Y1=1-X2
0784       Y2=1-X1
0785       X0=2.D0-X1-X2
0786       W1=(1.D0-X2)
0787       W2=(1.D0-X0+W-CH)
0788       W3=-(1.-X1-X2)
0789       W22=1.D0/ ((1.D0-X0+W-CH)**2+GAMC)
0790       W11=1.D0/((1.D0-X2)**2+GAMT)
0791       W33=1.D0/(W3**2+GAMW**2)
0792       W12=W1*W2*W11*W22
0793       W13=W1*W3*W11*W33
0794       W23=W2*W3*W22*W33
0795
0796       R11=4*T*W-16.*T*W*Y1-4.*T*Y2*Y1+8.*T*Y1+32.*T*W**2-20
0797      . .*T*Y1**2+8.*W*Y2*Y1+4.*W*Y1**2-4.*Y2*Y1**2-16.*T**2*W-
0798      .  32.*T**2*Y1+4.*T**2-16.*T**3-8.*W**2+4.*Y1**2-4.*Y1**3
0799       R22=-16.*W+16.*T*W-8.*T*Y2*Y1-4.*T*Y2**2-4.*T*Y1**2+16
0800      .  .*W*Y2 + 8.*W*Y2*Y1 + 16.*W*Y1 + 4.*W*Y2**2 + 4.*W*Y1**2+8.*Y2*
0801      .  Y1-12.*Y2*Y1**2-12.*Y2**2*Y1-16.*W**2+4.*Y2**2-4.*Y2**3
0802      .  +4.*Y1**2-4.*Y1**3
0803       R33=-4.*T*W+4.*T*W*Y2-2.*T*W*Y2*Y1+4.*T*W*Y1+T*W*Y2**2-
0804      .  3.*T*W*Y1**2+2.*T*Y2*Y1-3.*T*Y2*Y1**2+4.*T*W**2-4.*T*W**3
0805      .  +T*Y2**2-3.*T*Y2**2*Y1-T*Y2**3+T*Y1**2-T*Y1**3+4.*T**2
0806      .  *W-4.*T**2*W*Y2-4.*T**2*W*Y1-2.*T**2*Y2*Y1-4.*T**2*W**2-
0807      .  T**2*Y2**2-T**2*Y1**2+4.*W**2*Y2*Y1-8.*W**3*Y2-8.*W**3*Y1
0808      .  +4.*W**3+8.*W**4
0809       R12=-16.*W+48.*T*W-16.*T*W*Y2+16.*T*W*Y1+8.*T*Y2-32.*T
0810      .  *Y2*Y1+8.*T*Y1-8.*T*Y2**2 - 24.*T*Y1**2+16.*W*Y2+8.*W*Y2*
0811      .  Y1+16.*W*Y1+8.*W*Y1**2+8.*Y2*Y1-16.*Y2*Y1**2-16.*T**2*Y2
0812      .  -16.*T**2*Y1-8.*Y2**2*Y1-16.*W**2+8.*Y1**2-8.*Y1**3
0813       R13=8.*W-24.*T*W+16.*T*W*Y1 -4.*T*Y2+16.*T*Y2*Y1-4.*T*
0814      .  Y1+16.*T*W**2+4.*T*Y2**2+12.*T*Y1**2-8.*W*Y2-12.*W*Y2*Y1
0815      .  -8.*W*Y1+4.*W*Y1**2-4.*Y2*Y1+8.*Y2*Y1**2+16.*T**2*W+8.
0816      .  *T**2*Y2+8.*T**2*Y1+16.*W**2*Y2+24.*W**2*Y1+4.*Y2**2*Y1-
0817      .  32.*W**3-4.*Y1**2+4.*Y1**3
0818       R23=16.*W-16.*T*W+8.*T*W*Y2+8.*T*W*Y1+8.*T*Y2*Y1+4.*T*
0819      .  Y2**2+4.*T*Y1**2-16.*W*Y2-16.*W*Y1-4.*W*Y2**2+4.*W*Y1**2
0820      .  -8.*Y2*Y1+12.*Y2*Y1**2+8.*W**2*Y2-8.*W**2*Y1+12.*Y2**2*
0821      .  Y1-4.*Y2**2+4.*Y2**3-4.*Y1**2+4.*Y1**3
0822       GLVV=DSQRT(1.D0-GHVV**2)
0823       RES=GHT**2*R11*W11+GLVV**2*GAT**2*R22*W22+
0824      .    4.D0*GHVV**2*R33*W33/T+2.D0*GHT*GLVV*GAT*R12*W12+
0825      .    2.D0*GHT*GHVV*R13*W13+2.D0*GHVV*GLVV*GAT*R23*W23
0826       RETURN
0827       END
0828
0829
0830 C   *****************  INTEGRATION ROUTINE ***********************
0831 C    Returns SS as integral of FUNC from A to B, by 10-point Gauss-
0832 C    Legendre integration
0833       SUBROUTINE QGAUS1(FUNC,A,B,SS)
0834       IMPLICIT REAL*8(A-Z)
0835       INTEGER J
0836       DIMENSION X(5),W(5)
0837       EXTERNAL FUNC
0838       DATA X/.1488743389D0,.4333953941D0,.6794095682D0
0839      .  ,.8650633666D0,.9739065285D0/
0840       DATA W/.2955242247D0,.2692667193D0,.2190863625D0
0841      .  ,.1494513491D0,.0666713443D0/
0842       XM=0.5D0*(B+A)
0843       XR=0.5D0*(B-A)
0844       SS=0.D0
0845       DO 11 J=1,5
0846         DX=XR*X(J)
0847         SS=SS+W(J)*(FUNC(XM+DX)+FUNC(XM-DX))
0848 11    CONTINUE
0849       SS=XR*SS
0850       RETURN
0851       END
0852
0853 C     Returns SS as integral of FUNC from A to B, by 10-point Gauss-
0854 C      Legendre integration
0855       SUBROUTINE QGAUS2(FUNC,A,B,SS)
0856       IMPLICIT REAL*8(A-Z)
0857       INTEGER J
0858       DIMENSION X(5),W(5)
0859       EXTERNAL FUNC
0860       DATA X/.1488743389D0,.4333953941D0,.6794095682D0
0861      .  ,.8650633666D0,.9739065285D0/
0862       DATA W/.2955242247D0,.2692667193D0,.2190863625D0
0863      .  ,.1494513491D0,.0666713443D0/
0864       XM=0.5D0*(B+A)
0865       XR=0.5D0*(B-A)
0866       SS=0.D0
0867       DO 11 J=1,5
0868         DX=XR*X(J)
0869         SS=SS+W(J)*(FUNC(XM+DX)+FUNC(XM-DX))
0870 11    CONTINUE
0871       SS=XR*SS
0872       RETURN
0873       END
0874
0875 C ******************************************************************
0876
0877       DOUBLE PRECISION FUNCTION RUNM(Q,NF)
0878       PARAMETER (NN=6)
0879       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
0880       PARAMETER (ZETA3 = 1.202056903159594D0)
0881       DIMENSION AM(NN),YMSB(NN)
0882       COMMON/ALS/XLAMBDA,AMCA,AMBA,AMTA,N0A
0883       COMMON/MASSES/AMS,AMC,AMB,AMT
0884       COMMON/STRANGE/AMSB
0885       COMMON/RUN/XMSB,XMHAT,XKFAC
0886       COMMON/FLAG/IHIGGS,NNLO,IPOLE
0887       SAVE ISTRANGE
0888       B0(NF)=(33.D0-2.D0*NF)/12D0
0889       B1(NF) = (102D0-38D0/3D0*NF)/16D0
0890       B2(NF) = (2857D0/2D0-5033D0/18D0*NF+325D0/54D0*NF**2)/64D0
0891       G0(NF) = 1D0
0892       G1(NF) = (202D0/3D0-20D0/9D0*NF)/16D0
0893       G2(NF) = (1249D0-(2216D0/27D0+160D0/3D0*ZETA3)*NF
0894      .       - 140D0/81D0*NF**2)/64D0
0895       C1(NF) = G1(NF)/B0(NF) - B1(NF)*G0(NF)/B0(NF)**2
0896       C2(NF) = ((G1(NF)/B0(NF) - B1(NF)*G0(NF)/B0(NF)**2)**2
0897      .       + G2(NF)/B0(NF) + B1(NF)**2*G0(NF)/B0(NF)**3
0898      .       - B1(NF)*G1(NF)/B0(NF)**2 - B2(NF)*G0(NF)/B0(NF)**2)/2D0
0899       TRAN(X,XK)=1D0+4D0/3D0*ALPHAS(X,2)/PI+XK*(ALPHAS(X,2)/PI)**2
0900       CQ(X,NF)=(2D0*B0(NF)*X)**(G0(NF)/B0(NF))
0901      .            *(1D0+C1(NF)*X+C2(NF)*X**2)
0902       DATA ISTRANGE/0/
0903       PI=4D0*DATAN(1D0)
0904       ACC = 1.D-8
0905       AM(1) = 0
0906       AM(2) = 0
0907 C--------------------------------------------
0908       IMSBAR = 0
0909       IF(IMSBAR.EQ.1)THEN
0910        IF(ISTRANGE.EQ.0)THEN
0911 C--STRANGE POLE MASS FROM MSBAR-MASS AT 1 GEV
0912         AMSD = XLAMBDA
0913         AMSU = 1.D8
0914 123     AMS  = (AMSU+AMSD)/2
0915         AM(3) = AMS
0916         XMSB = AMS/CQ(ALPHAS(AMS,2)/PI,3)
0917      .            *CQ(ALPHAS(1.D0,2)/PI,3)/TRAN(AMS,0D0)
0918         DD = (XMSB-AMSB)/AMSB
0919         IF(DABS(DD).GE.ACC)THEN
0920          IF(DD.LE.0.D0)THEN
0921           AMSD = AM(3)
0922          ELSE
0923           AMSU = AM(3)
0924          ENDIF
0925          GOTO 123
0926         ENDIF
0927         ISTRANGE=1
0928        ENDIF
0929        AM(3) = AMSB
0930       ELSE
0931        AMS=AMSB
0932        AM(3) = AMS
0933       ENDIF
0934 C--------------------------------------------
0935       AM(3) = AMSB
0936       AM(4) = AMC
0937       AM(5) = AMB
0938       AM(6) = AMT
0939       XK = 16.11D0
0940       DO 1 I=1,NF-1
0941        XK = XK - 1.04D0*(1.D0-AM(I)/AM(NF))
0942 1     CONTINUE
0943       IF(NF.GE.4)THEN
0944        XMSB = AM(NF)/TRAN(AM(NF),0D0)
0945        XMHAT = XMSB/CQ(ALPHAS(AM(NF),2)/PI,NF)
0946       ELSE
0947        XMSB = 0
0948        XMHAT = 0
0949       ENDIF
0950       YMSB(3) = AMSB
0951
0952       IF(NF.EQ.3)THEN
0953        YMSB(4) = YMSB(3)*CQ(ALPHAS(AM(4),2)/PI,3)/
0954      .                   CQ(ALPHAS(1.D0,2)/PI,3)
0955        YMSB(5) = YMSB(4)*CQ(ALPHAS(AM(5),2)/PI,4)/
0956      .                   CQ(ALPHAS(AM(4),2)/PI,4)
0957        YMSB(6) = YMSB(5)*CQ(ALPHAS(AM(6),2)/PI,5)/
0958      .                   CQ(ALPHAS(AM(5),2)/PI,5)
0959       ELSEIF(NF.EQ.4)THEN
0960        YMSB(4) = XMSB
0961        YMSB(5) = YMSB(4)*CQ(ALPHAS(AM(5),2)/PI,4)/
0962      .                   CQ(ALPHAS(AM(4),2)/PI,4)
0963        YMSB(6) = YMSB(5)*CQ(ALPHAS(AM(6),2)/PI,5)/
0964      .                   CQ(ALPHAS(AM(5),2)/PI,5)
0965       ELSEIF(NF.EQ.5)THEN
0966        YMSB(5) = XMSB
0967        YMSB(4) = YMSB(5)*CQ(ALPHAS(AM(4),2)/PI,4)/
0968      .                   CQ(ALPHAS(AM(5),2)/PI,4)
0969        YMSB(6) = YMSB(5)*CQ(ALPHAS(AM(6),2)/PI,5)/
0970      .                   CQ(ALPHAS(AM(5),2)/PI,5)
0971       ELSEIF(NF.EQ.6)THEN
0972        YMSB(6) = XMSB
0973        YMSB(5) = YMSB(6)*CQ(ALPHAS(AM(5),2)/PI,5)/
0974      .                   CQ(ALPHAS(AM(6),2)/PI,5)
0975        YMSB(4) = YMSB(5)*CQ(ALPHAS(AM(4),2)/PI,4)/
0976      .                   CQ(ALPHAS(AM(5),2)/PI,4)
0977       ENDIF
0978       IF(Q.LT.AMC)THEN
0979        N0=3
0980        Q0 = 1.D0
0981       ELSEIF(Q.LE.AMB)THEN
0982        N0=4
0983        Q0 = AMC
0984       ELSEIF(Q.LE.AMT)THEN
0985        N0=5
0986        Q0 = AMB
0987       ELSE
0988        N0=6
0989        Q0 = AMT
0990       ENDIF
0991       IF(NNLO.EQ.1.AND.NF.GT.3)THEN
0992        XKFAC = TRAN(AM(NF),0D0)/TRAN(AM(NF),XK)
0993       ELSE
0994        XKFAC = 1D0
0995       ENDIF
0996
0997       RUNM = YMSB(N0)*CQ(ALPHAS(Q,2)/PI,N0)/
0998      .               CQ(ALPHAS(Q0,2)/PI,N0)
0999      .       * XKFAC
1000       RETURN
1001       END
1002
1003       DOUBLE PRECISION FUNCTION ALPHAS(Q,N)
1004       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1005       DIMENSION XLB(6)
1006       COMMON/ALSLAM/XLB1(6),XLB2(6)
1007       COMMON/ALS/XLAMBDA,AMC,AMB,AMT,N0
1008       B0(NF)=33.D0-2.D0*NF
1009       B1(NF)=6.D0*(153.D0-19.D0*NF)/B0(NF)**2
1010       ALS1(NF,X)=12.D0*PI/(B0(NF)*DLOG(X**2/XLB(NF)**2))
1011       ALS2(NF,X)=12.D0*PI/(B0(NF)*DLOG(X**2/XLB(NF)**2))
1012      .          *(1.D0-B1(NF)*DLOG(DLOG(X**2/XLB(NF)**2))
1013      .           /DLOG(X**2/XLB(NF)**2))
1014       PI=4.D0*DATAN(1.D0)
1015       IF(N.EQ.1)THEN
1016        DO 1 I=1,6
1017         XLB(I)=XLB1(I)
1018 1      CONTINUE
1019       ELSE
1020        DO 2 I=1,6
1021         XLB(I)=XLB2(I)
1022 2      CONTINUE
1023       ENDIF
1024       IF(Q.LT.AMC)THEN
1025        NF=3
1026       ELSEIF(Q.LE.AMB)THEN
1027        NF=4
1028       ELSEIF(Q.LE.AMT)THEN
1029        NF=5
1030       ELSE
1031        NF=6
1032       ENDIF
1033       IF(N.EQ.1)THEN
1034         ALPHAS=ALS1(NF,Q)
1035       ELSE
1036         ALPHAS=ALS2(NF,Q)
1037
1038       ENDIF
1039       RETURN
1040       END
1041
1042       SUBROUTINE ALSINI(ACC)
1043       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1044       DIMENSION XLB(6)
1045       COMMON/ALSLAM/XLB1(6),XLB2(6)
1046       COMMON/ALS/XLAMBDA,AMC,AMB,AMT,N0
1047       PI=4.D0*DATAN(1.D0)
1048       XLB1(1)=0D0
1049       XLB1(2)=0D0
1050       XLB2(1)=0D0
1051       XLB2(2)=0D0
1052 C...  Initialize with zero
1053       DO I=1,6
1054          XLB(I) = 0.0
1055       ENDDO
1056       IF(N0.EQ.3)THEN
1057        XLB(3)=XLAMBDA
1058        XLB(4)=XLB(3)*(XLB(3)/AMC)**(2.D0/25.D0)
1059        XLB(5)=XLB(4)*(XLB(4)/AMB)**(2.D0/23.D0)
1060        XLB(6)=XLB(5)*(XLB(5)/AMT)**(2.D0/21.D0)
1061       ELSEIF(N0.EQ.4)THEN
1062        XLB(4)=XLAMBDA
1063        XLB(5)=XLB(4)*(XLB(4)/AMB)**(2.D0/23.D0)
1064        XLB(3)=XLB(4)*(XLB(4)/AMC)**(-2.D0/27.D0)
1065        XLB(6)=XLB(5)*(XLB(5)/AMT)**(2.D0/21.D0)
1066       ELSEIF(N0.EQ.5)THEN
1067        XLB(5)=XLAMBDA
1068        XLB(4)=XLB(5)*(XLB(5)/AMB)**(-2.D0/25.D0)
1069        XLB(3)=XLB(4)*(XLB(4)/AMC)**(-2.D0/27.D0)
1070        XLB(6)=XLB(5)*(XLB(5)/AMT)**(2.D0/21.D0)
1071       ELSEIF(N0.EQ.6)THEN
1072        XLB(6)=XLAMBDA
1073        XLB(5)=XLB(6)*(XLB(6)/AMT)**(-2.D0/23.D0)
1074        XLB(4)=XLB(5)*(XLB(5)/AMB)**(-2.D0/25.D0)
1075        XLB(3)=XLB(4)*(XLB(4)/AMC)**(-2.D0/27.D0)
1076       ENDIF
1077       DO 1 I=1,6
1078        XLB1(I)=XLB(I)
1079 1     CONTINUE
1080       IF(N0.EQ.3)THEN
1081        XLB(3)=XLAMBDA
1082        XLB(4)=XLB(3)*(XLB(3)/AMC)**(2.D0/25.D0)
1083      .             *(2.D0*DLOG(AMC/XLB(3)))**(-107.D0/1875.D0)
1084        XLB(4)=XITER(AMC,XLB(3),3,XLB(4),4,ACC)
1085        XLB(5)=XLB(4)*(XLB(4)/AMB)**(2.D0/23.D0)
1086      .             *(2.D0*DLOG(AMB/XLB(4)))**(-963.D0/13225.D0)
1087        XLB(5)=XITER(AMB,XLB(4),4,XLB(5),5,ACC)
1088        XLB(6)=XLB(5)*(XLB(5)/AMT)**(2.D0/21.D0)
1089      .            *(2.D0*DLOG(AMT/XLB(5)))**(-321.D0/3381.D0)
1090        XLB(6)=XITER(AMT,XLB(5),5,XLB(6),6,ACC)
1091       ELSEIF(N0.EQ.4)THEN
1092        XLB(4)=XLAMBDA
1093        XLB(5)=XLB(4)*(XLB(4)/AMB)**(2.D0/23.D0)
1094      .             *(2.D0*DLOG(AMB/XLB(4)))**(-963.D0/13225.D0)
1095        XLB(5)=XITER(AMB,XLB(4),4,XLB(5),5,ACC)
1096        XLB(3)=XLB(4)*(XLB(4)/AMC)**(-2.D0/27.D0)
1097      .             *(2.D0*DLOG(AMC/XLB(4)))**(107.D0/2025.D0)
1098        XLB(3)=XITER(AMC,XLB(4),4,XLB(3),3,ACC)
1099        XLB(6)=XLB(5)*(XLB(5)/AMT)**(2.D0/21.D0)
1100      .            *(2.D0*DLOG(AMT/XLB(5)))**(-321.D0/3381.D0)
1101        XLB(6)=XITER(AMT,XLB(5),5,XLB(6),6,ACC)
1102       ELSEIF(N0.EQ.5)THEN
1103        XLB(5)=XLAMBDA
1104        XLB(4)=XLB(5)*(XLB(5)/AMB)**(-2.D0/25.D0)
1105      .             *(2.D0*DLOG(AMB/XLB(5)))**(963.D0/14375.D0)
1106        XLB(4)=XITER(AMB,XLB(5),5,XLB(4),4,ACC)
1107        XLB(3)=XLB(4)*(XLB(4)/AMC)**(-2.D0/27.D0)
1108      .             *(2.D0*DLOG(AMC/XLB(4)))**(107.D0/2025.D0)
1109        XLB(3)=XITER(AMC,XLB(4),4,XLB(3),3,ACC)
1110        XLB(6)=XLB(5)*(XLB(5)/AMT)**(2.D0/21.D0)
1111      .            *(2.D0*DLOG(AMT/XLB(5)))**(-321.D0/3381.D0)
1112        XLB(6)=XITER(AMT,XLB(5),5,XLB(6),6,ACC)
1113       ELSEIF(N0.EQ.6)THEN
1114        XLB(6)=XLAMBDA
1115        XLB(5)=XLB(6)*(XLB(6)/AMT)**(-2.D0/23.D0)
1116      .            *(2.D0*DLOG(AMT/XLB(6)))**(321.D0/3703.D0)
1117        XLB(5)=XITER(AMT,XLB(6),6,XLB(5),5,ACC)
1118        XLB(4)=XLB(5)*(XLB(5)/AMB)**(-2.D0/25.D0)
1119      .             *(2.D0*DLOG(AMB/XLB(5)))**(963.D0/14375.D0)
1120        XLB(4)=XITER(AMB,XLB(5),5,XLB(4),4,ACC)
1121        XLB(3)=XLB(4)*(XLB(4)/AMC)**(-2.D0/27.D0)
1122      .             *(2.D0*DLOG(AMC/XLB(4)))**(107.D0/2025.D0)
1123        XLB(3)=XITER(AMC,XLB(4),4,XLB(3),3,ACC)
1124       ENDIF
1125       DO 2 I=1,6
1126        XLB2(I)=XLB(I)
1127 2     CONTINUE
1128       RETURN
1129       END
1130
1131       DOUBLE PRECISION FUNCTION XITER(Q,XLB1,NF1,XLB,NF2,ACC)
1132       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1133       B0(NF)=33.D0-2.D0*NF
1134       B1(NF)=6.D0*(153.D0-19.D0*NF)/B0(NF)**2
1135       ALS2(NF,X,XLB)=12.D0*PI/(B0(NF)*DLOG(X**2/XLB**2))
1136      .              *(1.D0-B1(NF)*DLOG(DLOG(X**2/XLB**2))
1137      .              /DLOG(X**2/XLB**2))
1138       AA(NF)=12D0*PI/B0(NF)
1139       BB(NF)=B1(NF)/AA(NF)
1140       XIT(A,B,X)=A/2.D0*(1D0+DSQRT(1D0-4D0*B*DLOG(X)))
1141       PI=4.D0*DATAN(1.D0)
1142       XLB2=XLB
1143       II=0
1144 1     II=II+1
1145       X=DLOG(Q**2/XLB2**2)
1146       ALP=ALS2(NF1,Q,XLB1)
1147       A=AA(NF2)/ALP
1148       B=BB(NF2)*ALP
1149       XX=XIT(A,B,X)
1150       XLB2=Q*DEXP(-XX/2.D0)
1151       Y1=ALS2(NF1,Q,XLB1)
1152       Y2=ALS2(NF2,Q,XLB2)
1153       DY=DABS(Y2-Y1)/Y1
1154       IF(DY.GE.ACC) GOTO 1
1155       XITER=XLB2
1156       RETURN
1157       END
1158
1159       DOUBLE PRECISION FUNCTION FINT_(Z,XX,YY)
1160 C--ONE-DIMENSIONAL CUBIC INTERPOLATION
1161 C--Z  = WANTED POINT
1162 C--XX = ARRAY OF 4 DISCRETE X-VALUES AROUND Z
1163 C--YY = ARRAY OF 4 DISCRETE FUNCTION-VALUES AROUND Z
1164       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1165       DIMENSION XX(4),YY(4)
1166       X = DLOG(Z)
1167       X0=DLOG(XX(1))
1168       X1=DLOG(XX(2))
1169       X2=DLOG(XX(3))
1170       X3=DLOG(XX(4))
1171       Y0=DLOG(YY(1))
1172       Y1=DLOG(YY(2))
1173       Y2=DLOG(YY(3))
1174       Y3=DLOG(YY(4))
1175       A0=(X-X1)*(X-X2)*(X-X3)/(X0-X1)/(X0-X2)/(X0-X3)
1176       A1=(X-X0)*(X-X2)*(X-X3)/(X1-X0)/(X1-X2)/(X1-X3)
1177       A2=(X-X0)*(X-X1)*(X-X3)/(X2-X0)/(X2-X1)/(X2-X3)
1178       A3=(X-X0)*(X-X1)*(X-X2)/(X3-X0)/(X3-X1)/(X3-X2)
1179       FINT_=DEXP(A0*Y0+A1*Y1+A2*Y2+A3*Y3)
1180       RETURN
1181       END
1182
1183       DOUBLE PRECISION FUNCTION SP(X)
1184 C--REAL DILOGARITHM (SPENCE-FUNCTION)
1185       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1186       COMPLEX*16 CX,LI2
1187       CX = DCMPLX(X,0.D0)
1188       SP = DREAL(LI2(CX))
1189       RETURN
1190       END
1191
1192       COMPLEX*16 FUNCTION LI2(X)
1193 C--COMPLEX DILOGARITHM (SPENCE-FUNCTION)
1194       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1195       COMPLEX*16 X,Y,CLI2
1196       COMMON/CONST/ZETA2,ZETA3
1197       ZERO=1.D-16
1198       XR=DREAL(X)
1199       XI=DIMAG(X)
1200       R2=XR*XR+XI*XI
1201       LI2=0
1202       IF(R2.LE.ZERO)THEN
1203         LI2=X
1204         RETURN
1205       ENDIF
1206       RR=XR/R2
1207       IF(R2.EQ.1.D0.AND.XI.EQ.0.D0)THEN
1208         IF(XR.EQ.1.D0)THEN
1209           LI2=DCMPLX(ZETA2)
1210         ELSE
1211           LI2=-DCMPLX(ZETA2/2.D0)
1212         ENDIF
1213         RETURN
1214       ELSEIF(R2.GT.1.D0.AND.RR.GT.0.5D0)THEN
1215         Y=(X-1.D0)/X
1216         LI2=CLI2(Y)+ZETA2-CDLOG(X)*CDLOG(1.D0-X)+0.5D0*CDLOG(X)**2
1217         RETURN
1218       ELSEIF(R2.GT.1.D0.AND.RR.LE.0.5D0)THEN
1219         Y=1.D0/X
1220         LI2=-CLI2(Y)-ZETA2-0.5D0*CDLOG(-X)**2
1221         RETURN
1222       ELSEIF(R2.LE.1.D0.AND.XR.GT.0.5D0)THEN
1223         Y=1.D0-X
1224         LI2=-CLI2(Y)+ZETA2-CDLOG(X)*CDLOG(1.D0-X)
1225        RETURN
1226       ELSEIF(R2.LE.1.D0.AND.XR.LE.0.5D0)THEN
1227         Y=X
1228         LI2=CLI2(Y)
1229         RETURN
1230       ENDIF
1231       END
1232
1233       COMPLEX*16 FUNCTION CLI2(X)
1234 C--TAYLOR-EXPANSION FOR COMPLEX DILOGARITHM (SPENCE-FUNCTION)
1235       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1236       COMPLEX*16 X,Z
1237       COMMON/BERNOULLI/B2(18),B12(18),B3(18)
1238       COMMON/POLY/NBER
1239       N=NBER-1
1240       Z=-CDLOG(1.D0-X)
1241       CLI2=B2(NBER)
1242       DO 111 I=N,1,-1
1243         CLI2=Z*CLI2+B2(I)
1244 111   CONTINUE
1245       CLI2=Z**2*CLI2+Z
1246       RETURN
1247       END
1248
1249       DOUBLE PRECISION FUNCTION FACULT(N)
1250 C--DOUBLE PRECISION VERSION OF FACULTY
1251       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1252       FACULT=1.D0
1253       IF(N.EQ.0)RETURN
1254       DO 999 I=1,N
1255         FACULT=FACULT*DFLOAT(I)
1256 999   CONTINUE
1257       RETURN
1258       END
1259
1260       SUBROUTINE BERNINI(N)
1261 C--INITIALIZATION OF COEFFICIENTS FOR POLYLOGARITHMS
1262       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1263       DIMENSION B(18),PB(19)
1264       COMMON/BERNOULLI/B2(18),B12(18),B3(18)
1265       COMMON/CONST/ZETA2,ZETA3
1266       COMMON/POLY/NBER
1267
1268       NBER=N
1269       PI=4.D0*DATAN(1.D0)
1270
1271       B(1)=-1.D0/2.D0
1272       B(2)=1.D0/6.D0
1273       B(3)=0.D0
1274       B(4)=-1.D0/30.D0
1275       B(5)=0.D0
1276       B(6)=1.D0/42.D0
1277       B(7)=0.D0
1278       B(8)=-1.D0/30.D0
1279       B(9)=0.D0
1280       B(10)=5.D0/66.D0
1281       B(11)=0.D0
1282       B(12)=-691.D0/2730.D0
1283       B(13)=0.D0
1284       B(14)=7.D0/6.D0
1285       B(15)=0.D0
1286       B(16)=-3617.D0/510.D0
1287       B(17)=0.D0
1288       B(18)=43867.D0/798.D0
1289       ZETA2=PI**2/6.D0
1290       ZETA3=1.202056903159594D0
1291
1292       DO 995 I=1,18
1293         B2(I)=B(I)/FACULT(I+1)
1294         B12(I)=DFLOAT(I+1)/FACULT(I+2)*B(I)/2.D0
1295         PB(I+1)=B(I)
1296         B3(I)=0.D0
1297 995   CONTINUE
1298       PB(1)=1.D0
1299       DO 996 I=1,18
1300       DO 996 J=0,I
1301         B3(I)=B3(I)+PB(J+1)*PB(I-J+1)/FACULT(I-J)/FACULT(J+1)
1302      .                                            /DFLOAT(I+1)
1303 996   CONTINUE
1304
1305       RETURN
1306       END
1307
1308       DOUBLE PRECISION FUNCTION QQINT(RAT,H1,H2)
1309       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1310       N = 2
1311       QQINT = RAT**N * H1 + (1-RAT**N) * H2
1312       RETURN
1313       END
1314
1315       DOUBLE PRECISION FUNCTION XITLA(NO,ALP,ACC)
1316 C--ITERATION ROUTINE TO DETERMINE IMPROVED LAMBDA'S
1317       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1318       COMMON/PARAM/GF,ALPH,AMTAU,AMMUON,AMZ,AMW
1319       B0(NF)=33.D0-2.D0*NF
1320       B1(NF)=6.D0*(153.D0-19.D0*NF)/B0(NF)**2
1321       ALS2(NF,X,XLB)=12.D0*PI/(B0(NF)*DLOG(X**2/XLB**2))
1322      .              *(1.D0-B1(NF)*DLOG(DLOG(X**2/XLB**2))
1323      .              /DLOG(X**2/XLB**2))
1324       AA(NF)=12D0*PI/B0(NF)
1325       BB(NF)=B1(NF)/AA(NF)
1326       XIT(A,B,X)=A/2.D0*(1D0+DSQRT(1D0-4D0*B*DLOG(X)))
1327       PI=4.D0*DATAN(1.D0)
1328       NF=5
1329       Q=AMZ
1330       XLB=Q*DEXP(-AA(NF)/ALP/2.D0)
1331       IF(NO.EQ.1)GOTO 111
1332       II=0
1333 1     II=II+1
1334       X=DLOG(Q**2/XLB**2)
1335       A=AA(NF)/ALP
1336       B=BB(NF)*ALP
1337       XX=XIT(A,B,X)
1338       XLB=Q*DEXP(-XX/2.D0)
1339       Y1=ALP
1340       Y2=ALS2(NF,Q,XLB)
1341       DY=DABS(Y2-Y1)/Y1
1342       IF(DY.GE.ACC) GOTO 1
1343 111   XITLA=XLB
1344       RETURN
1345       END