C**********************************************************************
C
C                       ****************
C                       * VERSION 1.10 *
C                       ****************
C
C--This program calculates the production cross section of Higgs
C  bosons via W*/Z* -> h,H + W/Z at hadron colliders at NLO QCD
C  according to the formulae presented in
C
C  T. Han and S. Willenbrock, Phys. Lett. B273 (1991) 167.
C
C  The program allows to calculate the total production cross sections
C  for the scalar Higgs bosons of the SM and MSSM. The MSSM Higgs sector
C  is implemented in the approximate two-loop RGE approach of
C
C  M. Carena, J. Espinosa, M. Quiros and C.E.M. Wagner, Phys. Lett. B355
C  (1995) 209.
C
C  The parton densities are defined in the subroutine STRUC at the end
C  of the program and can be changed there. The default is the use of
C  the PDFLIB.
C
C--Description of the input file pphv.in with default values:
C  ==========================================================
C
C--Model: MSSM = 0 -> SM
C         MSSM = 1 -> MSSM
C  MSSM     = 0
C
C--tan(beta) for MSSM
C  TGBET    = 3.D0
C
C--Higgs type: (SM: HIGGS = 2)
C              HIGGS = 1 -> h
C              HIGGS = 2 -> H
C  HIGGS    = 1
C
C--Higgs mass loop: MA = MA1 -> MA2 with NA equidistant points in total
C                   (SM: MH = MA)
C  MA1      = 100.D0
C  MA2      = 1000.D0
C  NA       = 10
C
C--Choice of vector boson: 1 = Z, 2 = W
C  Z=1/W=2  = 2
C
C--Energy of hadron collider in GeV
C  ENERGY   = 14000.D0
C
C--Choice pp (0) or ppbar (1) collider
C  PP/PPBAR = 0
C
C--Scale loop: mu = xi * mu_0, xi = SCALE1 -> SCALE2 with NSCALE
C              equidistant points in total
C              mu_0 = Q = M_{HV}
C  SCALE1   = 1.0D0
C  SCALE2   = 1.0D0
C  NSCALE   = 1
C
C--Integration cut (should be left as it is)
C  EPSILON  = 1.D-8
C
C--Parton densities from LHAPDF.
C  NGROUP   = 0
C  NSET     = 0
C  Special cases: NGROUP = 1 --> CTEQ6
C
C--LHAPDF: PDF path and name
C  PDFPATH  =  Grids
C  PDFNAME  =  MSTW2008nlo68cl.LHgrid
C
C--VEGAS: IPOINT  = number of points
C         ITERAT  = number of iterations
C         NPRN: 0  = no output of individual iterations
C               1  = full output of individual iterations
C               10 = table output of individual iterations
C  IPOINT   = 10000
C  ITERAT   = 5
C  NPRN     = 10
C
C--Heavy quark masses for running masses and alpha_s (top quark
C                                                     decoupled!!!)
C  MC       = 1.40D0
C  MB       = 4.75D0
C  MT       = 173.2D0
C
C--ALSMZ = alpha_s(M_Z)
C  ALSMZ    = 0.12018D0
C
C--Order of calculation: LOOP = 1 -> Born cross section
C                        LOOP = 2 -> NLO cross section
C  LOOP     = 2
C
C--W mass, Z mass, Weinberg angle, Fermi constant
C  MW       = 80.41D0
C  MZ       = 91.187D0
C  SW2      = 0.2315D0
C  GF       = 1.16639D-5
C
C--MSSM parameters: in GeV
C                   MSQ = common squark mass
C                   MUR = right-handed up-type squark mass
C                   AT  = trilinear coupling A_t
C                   AB  = trilinear coupling A_b
C                   MU  = Higgsino mass paramater mu
C  MSQ      = 1000.D0
C  MUR      = 1000.D0
C  AT       = 0.D0
C  AB       = 0.D0
C  MU       = 0.D0
C
C
C  written by Michael Spira, e-mail: michael.spira@psi.ch.
C
C  June 23, 2015
C=======================================================
C
      PROGRAM PPHV
C--PROGRAM FOR W*/Z* -> H + W/Z AT HADRON COLLIDERS
      PARAMETER(NIN=95,NOUT=96)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      CHARACTER*100 PDFNAME,PATHNAME
      COMMON/PARAM/S,AM,SCF,EPS,GF,AMW,AMZ,SW2,GEVPB,NPROC
      COMMON/ALS/XLAMBDA,AMC,AMB,AMT,N0
      COMMON/RESULT/RES,ERR,DCHI2,DUM
      COMMON/PDFLIB/NGROUP,NSET
      COMMON/COLLIDER/ICOLL
      COMMON/MASSES/AMS0,AMC0,AMB0,AMT0
      COMMON/STRANGE/AMSB
      COMMON/HMASS/AMSM,AMA,AML,AMH,AMCH
      COMMON/BREAK/AMSQ,AMUR,AU,AD,AMU,AM2,ASL
      COMMON/COUP/GAT,GAB,GLT,GLB,GHT,GHB,GZAH,GZAL,
     .            GHHH,GLLL,GHLL,GLHH,GHAA,GLAA,GLVV,GHVV,
     .            GLPM,GHPM,B,A

      EXTERNAL SIG,DSIGQQ,DSIGGQ

      PI = 4*DATAN(1.D0)
C--CONVERSION FACTOR 1/GEV**2 --> PB
      GEVPB = 389379660.D0

      OPEN(NIN,FILE='pphv.in')
      OPEN(NOUT,FILE='pphv.out')

C--READING INPUT FILE
      READ(NIN,101)IMODEL
      READ(NIN,100)TGBET
      READ(NIN,101)IHIGGS
      IF(IMODEL.EQ.0) IHIGGS = 2
      READ(NIN,100)AM1
      READ(NIN,100)AM2
      READ(NIN,101)NAM
      READ(NIN,101)NPROC
      READ(NIN,100)W
      READ(NIN,101)NCOLL
C--NCOLL = 0 -> ICOLL =  1 -> PP COLLIDER
C--NCOLL = 1 -> ICOLL = -1 -> PPBAR COLLIDER
      ICOLL=1
      IF(NCOLL.EQ.1)ICOLL=-1
      READ(NIN,*)
      READ(NIN,100)SCALE1
      READ(NIN,100)SCALE2
      READ(NIN,101)NSCALE
      READ(NIN,100)EPS
      READ(NIN,*)
      READ(NIN,101)NGROUP
      READ(NIN,101)NSET
      READ(NIN,102)PATHNAME
      READ(NIN,102)PDFNAME
      READ(NIN,*)
      READ(NIN,101)IPOINT
      READ(NIN,101)ITERAT
      READ(NIN,101)NPRN
      READ(NIN,*)
      READ(NIN,100)AMC0
      READ(NIN,100)AMB0
      READ(NIN,100)AMT0
      READ(NIN,100)ALSMZ
      READ(NIN,101)LOOP
      READ(NIN,100)AMW
      READ(NIN,100)AMZ
      READ(NIN,100)SW2
      READ(NIN,100)GF
      READ(NIN,*)
      READ(NIN,100)AMSQ
      READ(NIN,100)AMUR
      READ(NIN,100)AU
      READ(NIN,100)AD
      READ(NIN,100)AMU
C--Parameters not needed
      ASL = 0
      AMSB = 0.190D0
C--Parameters for COMMON block MASSES and alpha_s (top quark
C                                                  decoupled!!!)
      AMS0 = AMSB
      AMC = AMC0
      AMB = AMB0
      AMT = AMT0*1.D8

C--Initializing vector boson mass
      IF(NPROC.EQ.1)THEN
       AMV = AMZ
      ELSE
       AMV = AMW
      ENDIF

C--VEGAS PARAMETERS
      ABSERR = 0.D0
      IGRAPH = 0

C--INTIALIZING ALPHA_S
      ACC = 1.D-10
      N0 = 5
      XLAMBDA = XITLA(LOOP,ALSMZ,ACC)
      CALL ALSINI(ACC)
C--INTIALIZING RANDOM NUMBER GENERATOR
      CALL RSTART(12,34,56,78)

C--INITIALIZE PDFSET
      CALL PDFSET(PATHNAME,PDFNAME)

      S = W**2

C--WRITING HEADER OF OUTPUT FILE
      IF(NPROC.EQ.1)THEN
       WRITE(NOUT,*)' PP ---> HZ + X'
      ELSE
       WRITE(NOUT,*)' PP ---> HW + X'
      ENDIF
      WRITE(NOUT,*)' =============='
      WRITE(NOUT,*)
      WRITE(NOUT,*)' NGROUP   = ',NGROUP,'   NSET = ',NSET
      WRITE(NOUT,*)' W        = ',W,' GEV'
      WRITE(NOUT,*)' LOOP     = ',LOOP
      IF(N0.EQ.4)THEN
       WRITE(NOUT,*)' LAMBDA_4 = ',XLAMBDA
      ELSE
       WRITE(NOUT,*)' LAMBDA_5 = ',XLAMBDA
      ENDIF
      WRITE(NOUT,*)
     .         ' ----------------------------------------------------'
      WRITE(NOUT,*)

      DO 9999 I = 1,NAM
       IF(NAM.EQ.1)THEN
        AM = AM1
       ELSE
        AM = AM1 + (AM2-AM1)/(NAM-1.D0)*(I-1)
       ENDIF
       AMA = AM
       COUPLING = 1.D0

C--Call MSSM couplings, if MSSM = 1
       IF(IMODEL.EQ.1) THEN
        CALL SUSYCP(TGBET)
        IF(IHIGGS.EQ.1) THEN
         AM = AML
         COUPLING = GLVV**2
        ELSE
         AM = AMH
         COUPLING = GHVV**2
        ENDIF
       ENDIF
      DO 9999 J = 1,NSCALE
       IF(NSCALE.EQ.1)THEN
        SCALE = SCALE1
       ELSE
        SCALE = SCALE1 + (SCALE2-SCALE1)/(NSCALE-1.D0)*(J-1)
       ENDIF

       SCF = SCALE
C--INTEGRATION OF BORN TERM
       IDIM = 2
       CALL VEGASN(SIG,ABSERR,IDIM,IPOINT,ITERAT,NPRN,IGRAPH)
       SIGB = RES*COUPLING
       DSIGB = ERR*COUPLING
       IF(LOOP.NE.1)THEN
C--INTEGRATION OF NLO QQBAR INITIAL STATE
        IDIM = 3
        CALL VEGASN(DSIGQQ,ABSERR,IDIM,IPOINT,ITERAT,NPRN,IGRAPH)
        SIGQ = RES*COUPLING
        DSIGQ = ERR*COUPLING
C--INTEGRATION OF NLO GQ INITIAL STATE
        CALL VEGASN(DSIGGQ,ABSERR,IDIM,IPOINT,ITERAT,NPRN,IGRAPH)
        SIGG = RES*COUPLING
        DSIGG = ERR*COUPLING
 
C--CALCULATION OF RELATIVE CORRECTIONS
        XKQ    = SIGQ/SIGB
        DKQ    = DSIGQ/SIGB
        XKG    = SIGG/SIGB
        DKG    = DSIGG/SIGB
        XKTOT  = 1 + XKQ + XKG
        DKTOT  = DSQRT(DKQ**2 + DKG**2)
        SIGH   = SIGB + SIGQ + SIGG
        DSIGH  = DSQRT(DSIGB**2 + DSIGQ**2 + DSIGG**2)
       ENDIF
       IF(IMODEL.EQ.0)THEN
        WRITE(NOUT,*)' MH_SM    = ',AM,'      GEV'
       ELSE
        WRITE(NOUT,*)' MA_MSSM  = ',AMA,'      GEV'
        WRITE(NOUT,*)' TG(BETA) = ',TGBET
        IF(IHIGGS.EQ.1)THEN
         WRITE(NOUT,*)' Mh_MSSM  = ',AM,'      GEV'
         WRITE(NOUT,*)' G_V^h**2 = ',COUPLING
        ELSE
         WRITE(NOUT,*)' MH_MSSM  = ',AM,'      GEV'
         WRITE(NOUT,*)' G_V^H**2 = ',COUPLING
        ENDIF
       ENDIF
       WRITE(NOUT,*)' SCALE    = ',SCALE
       IF(LOOP.NE.1)THEN
        WRITE(NOUT,*)' KQQ      = ',XKQ,'      +- ',DKQ
        WRITE(NOUT,*)' KGQ      = ',XKG,'      +- ',DKG
        WRITE(NOUT,*)' KTOT     = ',XKTOT,'      +- ',DKTOT
       ENDIF
        WRITE(NOUT,*)' SIG_LO   = (',SIGB,'      +- ',DSIGB,') PB'
       IF(LOOP.NE.1)THEN
        WRITE(NOUT,*)' SIG_NLO  = (',SIGH,'      +- ',DSIGH,') PB'
       ENDIF
       WRITE(NOUT,*)

9999  CONTINUE

100   FORMAT(10X,G30.20)
101   FORMAT(10X,I30)
102   FORMAT(11X,A100)
      STOP
      END

      DOUBLE PRECISION FUNCTION SIG(XX)
C--FUNCTION FOR BORN TERM
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON/PARAM/S,AM,SCF,EPS,GF,AMW,AMZ,SW2,GEVPB,NPROC
      DIMENSION XX(2)
      PI = 4*DATAN(1.D0)
      IF(NPROC.EQ.1)THEN
       AMV = AMZ
      ELSE
       AMV = AMW
      ENDIF
      TAU = (AM+AMV)**2/S
      X0 = EPS + (1-2*EPS)*XX(1)
      Y0 = EPS + (1-2*EPS)*XX(2)
c     X = TAU + (1-TAU)*X0
      X = DEXP(DLOG(TAU)*X0)
      Y = -DLOG(X)*Y0
c     FAC = -(1-TAU)*DLOG(X)*(1-2*EPS)**2
      FAC = DLOG(TAU)*X*DLOG(X)*(1-2*EPS)**2
      QSQ = X*S
      SCA = SCF*DSQRT(QSQ)
      PD = DLUMQQ(X,Y,SCA)
      COUP = SIG0(QSQ)
      SIG = GEVPB*COUP*PD*FAC
      RETURN
      END

      DOUBLE PRECISION FUNCTION DSIGQQ(XX)
C--FUNCTION FOR NLO QQBAR INITIAL STATE
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON/PARAM/S,AM,SCF,EPS,GF,AMW,AMZ,SW2,GEVPB,NPROC
      DIMENSION XX(3)
      EXTERNAL F0,F1
      PI = 4*DATAN(1.D0)
      ZETA2 = PI**2/6
      IF(NPROC.EQ.1)THEN
       AMV = AMZ
      ELSE
       AMV = AMW
      ENDIF
      TAU = (AM+AMV)**2/S
      X0 = EPS + (1-2*EPS)*XX(1)
      Y0 = EPS + (1-2*EPS)*XX(2)
      Z0 = EPS + (1-2*EPS)*XX(3)
c     X = TAU + (1-TAU)*X0
      X = DEXP(DLOG(TAU)*X0)
      Y = -DLOG(X)*Y0
      Z = TAU/X + (1-TAU/X)*Z0
c     FAC = -(1-TAU)*DLOG(X)*(1-TAU/X)*(1-2*EPS)**3
      FAC = DLOG(TAU)*X*DLOG(X)*(1-TAU/X)*(1-2*EPS)**3
      QSQ = X*S
      SCA = SCF*DSQRT(QSQ)
      ALP = ALPHAS(SCA,2)/PI
      PD0 = DLUMQQ(X,Y,SCA)
      COUP = SIG0(QSQ)/(1-TAU/X)*ALP*PD0
      W0 = (-2*DLOG(SCA**2/X/S) + 8.D0/3.D0*(ZETA2-2))*COUP
      RQQ = (-1-Z)
      QSQ = X*Z*S
      SCA = SCF*DSQRT(QSQ)
      ALP = ALPHAS(SCA,2)/PI
      PD = DLUMQQ(X,Y,SCA)
      COUP = SIG0(QSQ)*ALP*PD
      W1 = 4.D0/3.D0*(-RQQ*DLOG(SCA**2/X/S)-2*(1+Z)*DLOG(1-Z))*COUP
      RES0 = D0(TAU,X,Z,PD,PD0)
      RES1 = D1(TAU,X,Z,PD,PD0)
      W2 = RES0 + RES1
      WW = W0 + W1 + W2
      DSIGQQ = GEVPB*WW*FAC
      RETURN
      END

      DOUBLE PRECISION FUNCTION DSIGGQ(XX)
C--FUNCTION FOR NLO GQ INITIAL STATE
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON/PARAM/S,AM,SCF,EPS,GF,AMW,AMZ,SW2,GEVPB,NPROC
      DIMENSION XX(3)
      PI = 4*DATAN(1.D0)
      IF(NPROC.EQ.1)THEN
       AMV = AMZ
      ELSE
       AMV = AMW
      ENDIF
      TAU = (AM+AMV)**2/S
      X0 = EPS + (1-2*EPS)*XX(1)
      Y0 = EPS + (1-2*EPS)*XX(2)
      Z0 = EPS + (1-2*EPS)*XX(3)
c     X = TAU + (1-TAU)*X0
      X = DEXP(DLOG(TAU)*X0)
      Y = -DLOG(X)*Y0
      Z = TAU/X + (1-TAU/X)*Z0
c     FAC = -(1-TAU)*DLOG(X)*(1-TAU/X)*(1-2*EPS)**3
      FAC = DLOG(TAU)*X*DLOG(X)*(1-TAU/X)*(1-2*EPS)**3
      QSQ = X*Z*S
      COUP = SIG0(QSQ)
      PGQ = (Z**2 + (1-Z)**2)/2
      SCA = SCF*DSQRT(QSQ)
      ALP = ALPHAS(SCA,2)/PI
      PD = DLUMGQ(X,Y,SCA)
      WW = (-PGQ/2*DLOG(SCA**2/X/(1-Z)**2/S)+(1+6*Z-7*Z**2)/8)*COUP
      DSIGGQ = GEVPB*WW*ALP*PD*FAC
      RETURN
      END

      DOUBLE PRECISION FUNCTION SIG0(QSQ)
C--KERNEL FOR BORN TERM
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON/PARAM/S,AM,SCF,EPS,GF,AMW,AMZ,SW2,GEVPB,NPROC
      PI = 4*DATAN(1.D0)
      IF(NPROC.EQ.1)THEN
       AMV = AMZ
      ELSE
       AMV = AMW
      ENDIF
      SIG0 = GF**2*AMV**4/288.D0/PI/QSQ*DSQRT(DLAM(AMV**2,AM**2,QSQ))
     .     * (DLAM(AMV**2,AM**2,QSQ)+12*AMV**2/QSQ)/(1-AMV**2/QSQ)**2
      RETURN
      END

      DOUBLE PRECISION FUNCTION DLAM(X,Y,Z)
C--KAELLEN FUNCTION
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DLAM = (1-X/Z-Y/Z)**2 - 4*X*Y/Z**2
      RETURN
      END

      DOUBLE PRECISION FUNCTION D0(TAU0,TAU,X,PD,PD0)
C--PLUS DISTRIBUTION 1/(1-Z)_+
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      FFZ = F0(TAU,X)*PD
      FF1 = F0(TAU,1.D0)*PD0
      D0 = 1/(1-X)*(FFZ - FF1) + DLOG(1-TAU0/TAU)/(1-TAU0/TAU)*FF1
      RETURN
      END

      DOUBLE PRECISION FUNCTION D1(TAU0,TAU,X,PD,PD0)
C--PLUS DISTRIBUTION (LOG(1-Z)/(1-Z))_+
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      FFZ = F1(TAU,X)*PD
      FF1 = F1(TAU,1.D0)*PD0
      D1 = DLOG(1-X)/(1-X)*(FFZ - FF1)
     .   + DLOG(1-TAU0/TAU)**2/2/(1-TAU0/TAU)*FF1
      RETURN
      END

      DOUBLE PRECISION FUNCTION F0(TAU,Z)
C--FUNCTION FOR D0(TAU0,TAU,X,PD,PD0)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON/PARAM/S,AM,SCF,EPS,GF,AMW,AMZ,SW2,GEVPB,NPROC
      QSQ = TAU*Z*S
      PI = 4*DATAN(1.D0)
      SCA = SCF*DSQRT(QSQ)
      ALP = ALPHAS(SCA,2)/PI
      SLG = DLOG(SCA**2/TAU/S)
      F0 = -2*4.D0/3.D0*SIG0(QSQ)*SLG*ALP
      RETURN
      END

      DOUBLE PRECISION FUNCTION F1(TAU,Z)
C--FUNCTION FOR D1(TAU0,TAU,X,PD,PD0)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON/PARAM/S,AM,SCF,EPS,GF,AMW,AMZ,SW2,GEVPB,NPROC
      QSQ = TAU*Z*S
      PI = 4*DATAN(1.D0)
      SCA = SCF*DSQRT(QSQ)
      ALP = ALPHAS(SCA,2)/PI
      F1 = 4*4.D0/3.D0*SIG0(QSQ)*ALP
      RETURN
      END

      DOUBLE PRECISION FUNCTION DLUMQQ(TAU,Y,QQ)
C--Q QBAR-LUMINOSITY WITH Z/W-CHARGES
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION PDF1(-6:6),PDF2(-6:6)
      COMMON/PARAM/S,AM,SCF,EPS,GF,AMW,AMZ,SW2,GEVPB,NPROC
      COMMON/COLLIDER/ICOLL
      X = DEXP(-Y)
      CALL STRUC(X,QQ,PDF1)
      CALL STRUC(TAU/X,QQ,PDF2)
      DLUMQQ = 0
      IF(NPROC.EQ.1)THEN
       CU = 1.D0 + (1.D0-8.D0/3.D0*SW2)**2
       CD = 1.D0 + (1.D0-4.D0/3.D0*SW2)**2
       DLU = 0
       DLD = 0
       DO I=1,5,2
        DLD = DLD + PDF1(I)*PDF2(-ICOLL*I) + PDF2(ICOLL*I)*PDF1(-I)
        J = I+1
        IF(J.LT.6) DLU = DLU + PDF1(J)*PDF2(-ICOLL*J)
     .                       + PDF2(ICOLL*J)*PDF1(-J)
       ENDDO
       DLUMQQ = (CU*DLU+CD*DLD)/TAU
      ELSE
       DL = 0
       DO I=1,3,2
         J = I+1
         DL = DL + PDF1(I)*PDF2(-ICOLL*J) + PDF1(J)*PDF2(-ICOLL*I)
     .           + PDF2(ICOLL*I)*PDF1(-J) + PDF2(ICOLL*J)*PDF1(-I)
       ENDDO
       DLUMQQ = 4*DL/TAU
      ENDIF
      RETURN
      END

      DOUBLE PRECISION FUNCTION DLUMGQ(TAU,Y,QQ)
C--G Q-LUMINOSITY WITH Z/W-CHARGES
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION PDF1(-6:6),PDF2(-6:6)
      COMMON/PARAM/S,AM,SCF,EPS,GF,AMW,AMZ,SW2,GEVPB,NPROC
      COMMON/COLLIDER/ICOLL
      X = DEXP(-Y)
      IF(NPROC.EQ.1)THEN
       CU = 1.D0 + (1.D0-8.D0/3.D0*SW2)**2
       CD = 1.D0 + (1.D0-4.D0/3.D0*SW2)**2
       IMAX = 5
      ELSE
       CU = 4.D0
       CD = 4.D0
       IMAX = 3
      ENDIF
      CALL STRUC(X,QQ,PDF1)
      CALL STRUC(TAU/X,QQ,PDF2)
      DLU = 0
      DLD = 0
      DO I=1,IMAX,2
       DLD = DLD + PDF1(0)*(PDF2(I) + PDF2(-I))
     .           + PDF2(0)*(PDF1(I) + PDF1(-I))
       J = I+1
       IF(J.LT.6)
     . DLU = DLU + PDF1(0)*(PDF2(J) + PDF2(-J))
     .           + PDF2(0)*(PDF1(J) + PDF1(-J))
      ENDDO
      DLUMGQ = (CU*DLU+CD*DLD)/TAU
      RETURN
      END

      DOUBLE PRECISION FUNCTION ALPHAS(Q,N)
C--ALPHA_S: Q = SCALE, N = 1 -> LO, N = 2 -> NLO
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION XLB(6)
      COMMON/ALSLAM/XLB1(6),XLB2(6)
      COMMON/ALS/XLAMBDA,AMC,AMB,AMT,N0
      B0(NF)=33.D0-2.D0*NF
      B1(NF)=6.D0*(153.D0-19.D0*NF)/B0(NF)**2
      ALS1(NF,X)=12.D0*PI/(B0(NF)*DLOG(X**2/XLB(NF)**2))
      ALS2(NF,X)=12.D0*PI/(B0(NF)*DLOG(X**2/XLB(NF)**2))
     .          *(1.D0-B1(NF)*DLOG(DLOG(X**2/XLB(NF)**2))
     .           /DLOG(X**2/XLB(NF)**2))
      PI=4.D0*DATAN(1.D0)
      IF(N.EQ.1)THEN
       DO 1 I=1,6
        XLB(I)=XLB1(I)
1      CONTINUE
      ELSE
       DO 2 I=1,6
        XLB(I)=XLB2(I)
2      CONTINUE
      ENDIF

      IF(Q.LT.AMC)THEN
       NF=3
      ELSEIF(Q.LE.AMB)THEN
       NF=4
      ELSEIF(Q.LE.AMT)THEN
       NF=5
      ELSE
       NF=6
      ENDIF
      IF(N.EQ.1)THEN
        ALPHAS=ALS1(NF,Q)
      ELSE
        ALPHAS=ALS2(NF,Q)
      ENDIF
      RETURN
      END
 
      SUBROUTINE ALSINI(ACC)
C--ALPHA_S: INITIALIZATION OF LAMBDA_NF, ACC = ACCURAY AT THRESHOLDS
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION XLB(6)
      COMMON/ALSLAM/XLB1(6),XLB2(6)
      COMMON/ALS/XLAMBDA,AMC,AMB,AMT,N0
      B0(NF)=33.D0-2.D0*NF
      B1(NF)=6.D0*(153.D0-19.D0*NF)/B0(NF)**2
      ALS2(NF,X)=12.D0*PI/(B0(NF)*DLOG(X**2/XLB(NF)**2))
     .          *(1.D0-B1(NF)*DLOG(DLOG(X**2/XLB(NF)**2))
     .           /DLOG(X**2/XLB(NF)**2))
      PI=4.D0*DATAN(1.D0)
      XLB1(1)=0D0
      XLB1(2)=0D0
      XLB2(1)=0D0
      XLB2(2)=0D0
      IF(N0.EQ.4)THEN
       XLB(4)=XLAMBDA
       XLB(5)=XLB(4)*(XLB(4)/AMB)**(2.D0/23.D0)
      ELSEIF(N0.EQ.5)THEN
       XLB(5)=XLAMBDA
       XLB(4)=XLB(5)*(XLB(5)/AMB)**(-2.D0/25.D0)
      ENDIF
      XLB(3)=XLB(4)*(XLB(4)/AMC)**(-2.D0/27.D0)
      XLB(6)=XLB(5)*(XLB(5)/AMT)**(2.D0/21.D0)
      DO 1 I=1,6
       XLB1(I)=XLB(I)
1     CONTINUE
      IF(N0.EQ.4)THEN
       XLB(4)=XLAMBDA
       XLB(5)=XLB(4)*(XLB(4)/AMB)**(2.D0/23.D0)
     .             *(2.D0*DLOG(AMB/XLB(4)))**(-963.D0/13225.D0)
       XLB(5)=XITER(AMB,XLB(4),4,XLB(5),5,ACC)
      ELSEIF(N0.EQ.5)THEN
       XLB(5)=XLAMBDA
       XLB(4)=XLB(5)*(XLB(5)/AMB)**(-2.D0/25.D0)
     .             *(2.D0*DLOG(AMB/XLB(5)))**(963.D0/14375.D0)
       XLB(4)=XITER(AMB,XLB(5),5,XLB(4),4,ACC)
      ENDIF
      XLB(3)=XLB(4)*(XLB(4)/AMC)**(-2.D0/27.D0)
     .            *(2.D0*DLOG(AMC/XLB(4)))**(107.D0/2025.D0)
      XLB(3)=XITER(AMC,XLB(4),4,XLB(3),3,ACC)
      XLB(6)=XLB(5)*(XLB(5)/AMT)**(2.D0/21.D0)
     .           *(2.D0*DLOG(AMT/XLB(5)))**(-321.D0/3381.D0)
      XLB(6)=XITER(AMT,XLB(5),5,XLB(6),6,ACC)
      DO 2 I=1,6
       XLB2(I)=XLB(I)
2     CONTINUE
      RETURN
      END
 
      DOUBLE PRECISION FUNCTION XITER(Q,XLB1,NF1,XLB,NF2,ACC)
C--ALPHA_S: ITERATION FOR ALPHA_S(M_Z) -> LAMBDA_5
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      B0(NF)=33.D0-2.D0*NF
      B1(NF)=6.D0*(153.D0-19.D0*NF)/B0(NF)**2
      ALS2(NF,X,XLB)=12.D0*PI/(B0(NF)*DLOG(X**2/XLB**2))
     .              *(1.D0-B1(NF)*DLOG(DLOG(X**2/XLB**2))
     .              /DLOG(X**2/XLB**2))
      AA(NF)=12D0*PI/B0(NF)
      BB(NF)=B1(NF)/AA(NF)
      XIT(A,B,X)=A/2.D0*(1D0+DSQRT(1D0-4D0*B*DLOG(X)))
      PI=4.D0*DATAN(1.D0)
      XLB2=XLB
      IF(ACC.GE.1.D0) GOTO 111
      II=0
1     II=II+1
      X=DLOG(Q**2/XLB2**2)
      ALP=ALS2(NF1,Q,XLB1)
      A=AA(NF2)/ALP
      B=BB(NF2)*ALP
      XX=XIT(A,B,X)
      XLB2=Q*DEXP(-XX/2.D0)
      Y1=ALS2(NF1,Q,XLB1)
      Y2=ALS2(NF2,Q,XLB2)
      DY=DABS(Y2-Y1)/Y1
      IF(DY.GE.ACC) GOTO 1
111   XITER=XLB2
      RETURN
      END

      DOUBLE PRECISION FUNCTION XITLA(NO,ALP,ACC)
C--ITERATION ROUTINE TO DETERMINE IMPROVED LAMBDAS
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON/PARAM/S,AM,SCF,EPS,GF,AMW,AMZ,SW2,GEVPB,NPROC
      B0(NF)=33.D0-2.D0*NF
      B1(NF)=6.D0*(153.D0-19.D0*NF)/B0(NF)**2
      B2(NF)=27/2.D0*(2857-5033/9.D0*NF+325/27.D0*NF**2)/B0(NF)**3
      B3(NF)= 81*(149753/6.d0+3564*zeta3-(1078361/162.d0+6508*zeta3/27)
     .        *nf+(50065/162.d0+6472*zeta3/81)*nf**2+1093/729.d0*nf**3)
     .      / B0(NF)**4
      ALS2(NF,X,XLB)=12.D0*PI/(B0(NF)*DLOG(X**2/XLB**2))
     .              *(1.D0-B1(NF)*DLOG(DLOG(X**2/XLB**2))
     .              /DLOG(X**2/XLB**2))
      ALS3(NF,X,XLB)=12.D0*PI/(B0(NF)*DLOG(X**2/XLB**2))
     .          *(1.D0-B1(NF)*DLOG(DLOG(X**2/XLB**2))
     .           /DLOG(X**2/XLB**2)
     .           +(B1(NF)**2*(DLOG(DLOG(X**2/XLB**2))**2
     .                      -DLOG(DLOG(X**2/XLB**2))-1)+B2(NF))
     .           /DLOG(X**2/XLB**2)**2)
      ALS4(NF,X,XLB)=12.D0*PI/(B0(NF)*DLOG(X**2/XLB**2))
     .          *(1.D0-B1(NF)*DLOG(DLOG(X**2/XLB**2))
     .           /DLOG(X**2/XLB**2)
     .           +(B1(NF)**2*(DLOG(DLOG(X**2/XLB**2))**2
     .                      -DLOG(DLOG(X**2/XLB**2))-1)+B2(NF))
     .           /DLOG(X**2/XLB**2)**2
     .           -(B1(NF)**3*(DLOG(DLOG(X**2/XLB**2))**3
     .                      -5*DLOG(DLOG(X**2/XLB**2))**2/2
     .                      -2*DLOG(DLOG(X**2/XLB**2))+1/2.d0)
     .            +3*B1(NF)*B2(NF)*DLOG(DLOG(X**2/XLB**2))
     .            -B3(NF)/2)/DLOG(X**2/XLB**2)**3)
      AA(NF)=12D0*PI/B0(NF)
      BB(NF)=B1(NF)/AA(NF)
      CC(NF)=B2(NF)/AA(NF)
      DD(NF)=B3(NF)/AA(NF)
      XIT(A,B,X)=A/2.D0*(1D0+DSQRT(1D0-4D0*B*DLOG(X)))
      XIT3(A,B,C,X)=A/2.D0*(1D0+DSQRT(1D0-4D0*B*DLOG(X)
     .          *(1-(A*B*(DLOG(X)**2-DLOG(X)-1)+C/B)/X/DLOG(X))))
      XIT4(A,B,C,D,X)=A/2.D0*(1D0+DSQRT(1D0-4D0*B*DLOG(X)
     .          *(1-(A*B*(DLOG(X)**2-DLOG(X)-1)+C/B)/X/DLOG(X)
     .          +(A**2*B**2*(DLOG(X)**3-5*DLOG(X)**2/2-2*DLOG(X)+1/2.D0)
     .           +3*A*C*DLOG(X)-D/B/2)/X**2/DLOG(X))))
      PI=4.D0*DATAN(1.D0)
      N3LO = 0
      ZETA2 = PI**2/6
      ZETA3 = 1.2020569031595942853997381D0
      NF=5
      Q=AMZ
      XLB=Q*DEXP(-AA(NF)/ALP/2.D0)
      IF(NO.EQ.1)GOTO 111
      II=0
1     II=II+1
      X=DLOG(Q**2/XLB**2)
      A=AA(NF)/ALP
      B=BB(NF)*ALP
      C=CC(NF)*ALP
      D=DD(NF)*ALP
      IF(NO.EQ.2)THEN
       XX=XIT(A,B,X)
      ELSEIF(NO.EQ.3)THEN
       XX=XIT3(A,B,C,X)
      ELSE
       XX=XIT4(A,B,C,D,X)
      ENDIF
      IF(N3LO.NE.0) XX=XIT4(A,B,C,D,X)
      XLB=Q*DEXP(-XX/2.D0)
      Y1=ALP
      IF(NO.EQ.2)THEN
       Y2=ALS2(NF,Q,XLB)
      ELSEIF(NO.EQ.3)THEN
       Y2=ALS3(NF,Q,XLB)
      ELSE
       Y2=ALS4(NF,Q,XLB)
      ENDIF
      IF(N3LO.NE.0) Y2=ALS4(NF,Q,XLB)
      DY=DABS(Y2-Y1)/Y1
      IF(DY.GE.ACC) GOTO 1
111   XITLA=XLB
      RETURN
      END

      SUBROUTINE VEGASN(FXN,ACC,NDIM,NPOINT,NITT,NPRN,IGRAPH)
C--VEGAS: WARMUP AND MAIN INTEGRATION
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      EXTERNAL FXN
      NPOINT0=NPOINT/10
      NITT0=2*NITT
      CALL VEGAS(FXN,ACC,NDIM,NPOINT0,NITT0,NPRN,IGRAPH)
      CALL VEGAS1(FXN,ACC,NDIM,NPOINT,NITT,NPRN,IGRAPH)
      RETURN
      END

C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      SUBROUTINE VEGAS(FXN,BCC,NDIM,NCALL,ITMX,NPRN,IGRAPH)
C--VEGAS SUBROUTINE
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON/RESU/RES
      COMMON/VEGOUT/NV
      COMMON/BVEG2/NDO,IT,SI,SI2,SWGT,SCHI,XI(50,10),SCALLS
     1,D(50,10),DI(50,10),NXI(50,10)
      DIMENSION XIN(50),R(50),DX(10),IA(10),KG(10),DT(10)
      DIMENSION XL(10),XU(10),QRAN(10),X(10)
      COMMON/RESULT/S1,S2,S3,S4
       EXTERNAL FXN
      DATA XL,XU/10*0.D0,10*1.D0/
      DATA NDMX/50/,ALPH/1.5D0/,ONE/1.D0/,MDS/1/
      IPR=1
      IF(NPRN.GT.0)IPR=0
      NDO=1
      DO 1 J=1,NDIM
1     XI(1,J)=ONE
      ENTRY VEGAS1(FXN,BCC,NDIM,NCALL,ITMX,NPRN,IGRAPH)
      NOW=IGRAPH
CS    IF(IGRAPH.GT.0)CALL INPLOT(NOW,F1,W)
      IT=0
      SI=0.D0
      SI2=SI
      SWGT=SI
      SCHI=SI
      SCALLS=SI
      ENTRY VEGAS2(FXN,BCC,NDIM,NCALL,ITMX,NPRN,IGRAPH)
      ND=NDMX
      NG=1
      IF(MDS.EQ.0) GO TO 2
      NG=(NCALL*0.5)**(1./NDIM)
      MDS=1
      IF((2*NG-NDMX).LT.0) GO TO 2
      MDS=-1
      NPG=NG/NDMX+1
      ND=NG/NPG
      NG=NPG*ND
2     K=NG**NDIM
      NPG=NCALL/K
      IF(NPG.LT.2)NPG=2
      CALLS=NPG*K
      DXG=ONE/NG
      DV2G=DXG**(2*NDIM)/NPG/NPG/(NPG-ONE)
      XND=ND
      NDM=ND-1
      DXG=DXG*XND
      XJAC=ONE
      DO 3 J=1,NDIM
      DX(J)=XU(J)-XL(J)
3     XJAC=XJAC*DX(J)
      IF(ND.EQ.NDO) GO TO 8
      RC=NDO/XND
      DO 7 J=1,NDIM
      K=0
      XN=0.D0
      DR=XN
      I=K
4     K=K+1
      DR=DR+ONE
      XO=XN
      XN=XI(K,J)
5     IF(RC.GT.DR) GO TO 4
      I=I+1
      DR=DR-RC
      XIN(I)=XN-(XN-XO)*DR
      IF(I.LT.NDM) GO TO 5
      DO 6  I=1,NDM
6     XI(I,J)=XIN(I)
7     XI(ND,J)=ONE
      NDO=ND
      ACC=BCC
      IF(NPRN.NE.0.AND.NPRN.NE.10)WRITE(NV,200)NDIM,CALLS,IT,ITMX
     1,ACC,MDS,ND
8     CONTINUE
      IF(NPRN.EQ.10)WRITE(NV,290)NDIM,CALLS,ITMX,ACC,MDS,ND
      ENTRY VEGAS3(FXN,BCC,NDIM,NCALL,ITMX,NPRN,IGRAPH)
9     IT=IT+1
      TI=0.D0
      TSI=TI
CS    IF(IGRAPH.GT.0)CALL REPLOT(NOW,F1,W)
      DO 10 J=1,NDIM
      KG(J)=1
      DO 10 I=1,ND
      NXI(I,J)=0
      D(I,J)=TI
10    DI(I,J)=TI
11    FB=0.D0
      F2B=FB
      K=0
12    K=K+1
      DO 121 J=1,NDIM
121   QRAN(J)=RANDM(0)
      WGT=XJAC
      DO 15 J=1,NDIM
      XN=(KG(J)-QRAN(J))*DXG+ONE
      IA(J)=XN
      IAJ=IA(J)
      IAJ1=IAJ-1
      IF(IAJ.GT.1) GO TO 13
      XO=XI(IAJ,J)
      RC=(XN-IAJ)*XO
      GO TO 14
13    XO=XI(IAJ,J)-XI(IAJ1,J)
      RC=XI(IAJ1,J)+(XN-IAJ)*XO
14    X(J)=XL(J)+RC*DX(J)
15    WGT=WGT*XO*XND
      F=FXN(X)*WGT
      F1=F/CALLS
      W=WGT/CALLS
CS    IF(IGRAPH.GT.0)CALL XPLOT(NOW,F1,W)
      F2=F**2
      FB=FB+F
      F2B=F2B+F2
      DO 16 J=1,NDIM
      IAJ=IA(J)
      NXI(IAJ,J)=NXI(IAJ,J)+1
      DI(IAJ,J)=DI(IAJ,J)+F/CALLS
16    IF(MDS.GE.0)  D(IAJ,J)=D(IAJ,J)+F2
      IF(K.LT.NPG) GO TO 12
      F2B=F2B*NPG
      F2B=SQRT(F2B)
      F2B=(F2B-FB)*(F2B+FB)
      TI=TI+FB
      TSI=TSI+F2B
      IF(MDS.GE.0) GO TO 18
      DO 17 J=1,NDIM
      IAJ=IA(J)
17    D(IAJ,J)=D(IAJ,J)+F2B
18    K=NDIM
19    KG(K)=MOD(KG(K),NG)+1
      IF(KG(K).NE.1) GO TO 11
      K=K-1
      IF(K.GT.0) GO TO 19
      TI=TI/CALLS
      TSI=TSI*DV2G
      TI2=TI*TI
      WGT=TI2/TSI
      SI=SI+TI*WGT
      SI2=SI2+TI2
      SWGT=SWGT+WGT
      SCHI=SCHI+TI2*WGT
      SCALLS=SCALLS+CALLS
      AVGI=SI/SWGT
      SD=SWGT*IT/SI2
      CHI2A=0.D0
      IF(IT.GT.1)CHI2A=SD*(SCHI/SWGT-AVGI*AVGI)/(IT-1)
      SD=ONE/SD
      SD=SQRT(SD)
      IF(NPRN.EQ.0) GO TO 21
      TSI=SQRT(TSI)
      IF(NPRN.NE.10)WRITE(NV,201)IPR,IT,TI,TSI,AVGI,SD,CHI2A
      RES=AVGI
      IF(NPRN.EQ.10)WRITE(NV,203)IT,TI,TSI,AVGI,SD,CHI2A
      IF(NPRN.GE.0) GO TO 21
      DO 20 J=1,NDIM
      WRITE(NV,202)J
20    WRITE(NV,204)(XI(I,J),DI(I,J),D(I,J),I=1,ND)
21    IF(ABS(SD/AVGI).LE.ABS(ACC).OR.IT.GE.ITMX)NOW=2
      S1=AVGI
      S2=SD
      S3=TI
      S4=TSI
CS    IF(IGRAPH.GT.0)CALL PLOTIT(NOW,F1,W)
C      DO 23 J=1,NDIM
C      XO=D(1,J)
C      XN=D(2,J)
C      D(1,J)=(XO+XN)*0.5D0
C      DT(J)=D(1,J)
C      DO 22 I=2,NDM
C      D(I,J)=XO+XN
C      XO=XN
C      XN=D(I+1,J)
C      D(I,J)=(D(I,J)+XN)/3.D0
C22    DT(J)=DT(J)+D(I,J)
C      D(ND,J)=(XN+XO)*0.5D0
C23    DT(J)=DT(J)+D(ND,J)
C-----THIS PART OF THE VEGAS-ALGORITHM IS UNSTABLE
C-----IT SHOULD BE REPLACED BY
      DO 23 J=1,NDIM
      DT(J)=0.D0
      DO 23 I=1,ND
      IF(NXI(I,J).GT.0)D(I,J)=D(I,J)/NXI(I,J)
23    DT(J)=DT(J)+D(I,J)
      DO 28 J=1,NDIM
      RC=0.D0
      DO 24 I=1,ND
      R(I)=0.D0
      IF(D(I,J).LE.0.D0)GO TO 24
      XO=DT(J)/D(I,J)
      R(I)=((XO-ONE)/XO/LOG(XO))**ALPH
24    RC=RC+R(I)
      RC=RC/XND
      K=0
      XN=0.D0
      DR=XN
      I=K
25    K=K+1
      DR=DR+R(K)
      XO=XN
      XN=XI(K,J)
26    IF(RC.GT.DR) GO TO 25
      I=I+1
      DR=DR-RC
      XIN(I)=XN-(XN-XO)*DR/R(K)
      IF(I.LT.NDM) GO TO 26
      DO 27 I=1,NDM
27    XI(I,J)=XIN(I)
28    XI(ND,J)=ONE
      IF(IT.LT.ITMX.AND.ABS(ACC).LT.ABS(SD/AVGI))GO TO 9
200   FORMAT(35H0INPUT PARAMETERS FOR VEGAS   NDIM=,I3
     1,8H  NCALL=,F8.0/28X,5H  IT=,I5,8H  ITMX =,I5/28X
     2,6H  ACC=,G9.3/28X,6H  MDS=,I3,6H   ND=,I4//)
290   FORMAT(13H0VEGAS  NDIM=,I3,8H  NCALL=,F8.0,8H  ITMX =,I5
     1,6H  ACC=,G9.3,6H  MDS=,I3,6H   ND=,I4)
201   FORMAT(/I1,20HINTEGRATION BY VEGAS/13H0ITERATION NO,I3,
     114H.   INTEGRAL =,G14.8/20X,10HSTD DEV  =,G10.4/
     234H ACCUMULATED RESULTS.   INTEGRAL =,G14.8/
     324X,10HSTD DEV  =,G10.4 / 24X,18HCHI**2 PER ITN   =,G10.4)
202   FORMAT(14H0DATA FOR AXIS,I2 / 7X,1HX,7X,10H  DELT I  ,
     12X,11H CONVCE    ,11X,1HX,7X,10H  DELT I  ,2X,11H CONVCE
     2,11X,1HX,7X,10H  DELT I  ,2X,11H CONVCE     /)
204   FORMAT(1X,3G12.4,5X,3G12.4,5X,3G12.4)
203   FORMAT(1H ,I3,G20.8,G12.4,G20.8,G12.4,G12.4)
      S1=AVGI
      S2=SD
      S3=CHI2A
      RETURN
      END

C----------------------------------------------------------------------
C  A UNIVERSAL RANDOM NUMBER GENERATOR
 
        DOUBLE PRECISION FUNCTION RANDM(IDMY)
C--RANDOM NUMBER GENERATOR
        IMPLICIT REAL*8(A-H,O-Z)
        REAL*4 UNIV
        RANDM=DBLE(UNIV(1))
        RETURN 
        END

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

        FUNCTION UNIV(IDUM)
C--FUNCTION FOR RANDOM NUMBER GENERATOR
        REAL U(97)
        COMMON /SET1/ U,C,CD,CM,I,J
        UNIV=U(I)-U(J)
        IF(UNIV.LT.0.) UNIV=UNIV+1.
        U(I)=UNIV
        I=I-1
        IF(I.EQ.0) I=97
        J=J-1
        IF(J.EQ.0) J=97
        C=C-CD
        IF(C.LT.0.) C=C+CM
        UNIV=UNIV-C
        IF(UNIV.LT.0.) UNIV=UNIV+1
        RETURN
        END
 
C----------------------------------------------------------------------
C INITIALIZING THE RANDOM NUMBER GENERATOR
C TO INITIALIZE CALL RSTART(12,34,56,78)


        SUBROUTINE RSTART(I,J,K,L)
C--INITIALIZATION ROUTINE FOR RANDOM NUMBER GENERATOR
        REAL U(97)
        COMMON /SET1/ U,C,CD,CM,ISTART,JSTART
        IF ((I.LT.0).OR.(I.GT.178 )) STOP 'FIRST SEED .LT.0 OR .GT.178'
        IF ((J.LT.0).OR.(J.GT.178 )) STOP 'SECOND SEED .LT.0 OR .GT.178'
        IF ((K.LT.0).OR.(K.GT.178 )) STOP 'THIRD SEED .LT.0 OR .GT.178'
        IF ((L.LT.0).OR.(L.GT.168 )) STOP 'FOURTH SEED .LT.0 OR .GT.168'
        IF ( (I.EQ.1).AND.(J.EQ.1).AND.(K.EQ.1) ) STOP
     &     'FIRST, SECOND AND THIRD SEEDS ARE ALL EQUAL TO 1'
        ISTART=97
        JSTART=33
        IDUM=I
        JDUM=J
        KDUM=K
        LDUM=L
        DO 2 II=1,97
        S=0.
        T=.5
        DO 3 JJ=1,24
          M=MOD(MOD(IDUM*JDUM,179)*K,179)
          IDUM=JDUM
          JDUM=KDUM
          KDUM=M
          LDUM=MOD(53*LDUM+1,169)
          IF(MOD(LDUM*M,64).GE.32) S=S+T
3         T=.5*T
2         U(II)=S
        C=362436./16777216.
        CD=7654321./16777216.
        CM=16777213./16777216.
        RETURN
        END

C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

C *****************************************************************
C ************* SUBROUTINE FOR THE SUSY COUPLINGS *****************
C *****************************************************************
      SUBROUTINE SUSYCP(TGBET)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION LA1,LA2,LA3,LA4,LA5,LA6,LA7,LA3T
      COMMON/MASSES/AMS,AMC,AMB,AMT
      COMMON/HMASS/AMSM,AMA,AML,AMH,AMCH
      COMMON/HSELF/LA1,LA2,LA3,LA4,LA5,LA6,LA7
      COMMON/BREAK/AMSQ,AMUR,AU,AD,AMU,AM2,ASL
      COMMON/PARAM/S,AM,SCF,EPS,GF,AMW,AMZ,SW2,GEVPB,NPROC
      COMMON/COUP/GAT,GAB,GLT,GLB,GHT,GHB,GZAH,GZAL,
     .            GHHH,GLLL,GHLL,GLHH,GHAA,GLAA,GLVV,GHVV,
     .            GLPM,GHPM,B,A
      PI=4*DATAN(1D0)
      V=1.D0/DSQRT(DSQRT(2.D0)*GF)
      BET=DATAN(TGBET)
      SB = DSIN(BET)
      CB = DCOS(BET)
      CALL SUBH(AMA,TGBET,AMSQ,AMUR,AMT,AU,AD,AMU,AML,AMH,AMCH,SA,CA,
     .          TBA)
      LA3T=LA3+LA4+LA5
      AMA2=AMA**2
      AML2=AML**2
      AMH2=AMH**2
      AMP2=AMCH**2
      SBMA = SB*CA-CB*SA
      CBMA = CB*CA+SB*SA
      SBPA = SB*CA+CB*SA
      CBPA = CB*CA-SB*SA
      S2A = 2*SA*CA
      C2A = CA**2-SA**2
      S2B = 2*SB*CB
      C2B = CB**2-SB**2
      GLZZ = 1/V/2*AML2*SBMA
      GHZZ = 1/V/2*AMH2*CBMA
      GLWW = 2*GLZZ
      GHWW = 2*GHZZ
      GLAZ = 1/V*(AML2-AMA2)*CBMA
      GHAZ = -1/V*(AMH2-AMA2)*SBMA
      GLPW = -1/V*(AMP2-AML2)*CBMA
      GLMW = GLPW
      GHPW = 1/V*(AMP2-AMH2)*SBMA
      GHMW = GHPW
      GAPW = 1/V*(AMP2-AMA2)
      GAMW = -GAPW
      GHHH = V/2*(LA1*CA**3*CB + LA2*SA**3*SB + LA3T*SA*CA*SBPA
     .     + LA6*CA**2*(3*SA*CB+CA*SB) + LA7*SA**2*(3*CA*SB+SA*CB))
      GLLL = -V/2*(LA1*SA**3*CB - LA2*CA**3*SB + LA3T*SA*CA*CBPA
     .     - LA6*SA**2*(3*CA*CB-SA*SB) + LA7*CA**2*(3*SA*SB-CA*CB))
      GLHH = -3*V/2*(LA1*CA**2*CB*SA - LA2*SA**2*SB*CA
     .     + LA3T*(SA**3*CB-CA**3*SB+2*SBMA/3)
     .     - LA6*CA*(CB*C2A-SA*SBPA) - LA7*SA*(C2A*SB+CA*SBPA))
      GHLL = 3*V/2*(LA1*SA**2*CB*CA + LA2*CA**2*SB*SA
     .     + LA3T*(SA**3*SB+CA**3*CB-2*CBMA/3)
     .     - LA6*SA*(CB*C2A+CA*CBPA) + LA7*CA*(C2A*SB+SA*CBPA))
      GLAA = -V/2*(LA1*SB**2*CB*SA - LA2*CB**2*SB*CA
     .     - LA3T*(SB**3*CA-CB**3*SA) + 2*LA5*SBMA
     .     - LA6*SB*(CB*SBPA+SA*C2B) - LA7*CB*(C2B*CA-SB*SBPA))
      GHAA = V/2*(LA1*SB**2*CB*CA + LA2*CB**2*SB*SA
     .     + LA3T*(SB**3*SA+CB**3*CA) - 2*LA5*CBMA
     .     - LA6*SB*(CB*CBPA+CA*C2B) + LA7*CB*(SB*CBPA+SA*C2B))
      GLPM = 2*GLAA + V*(LA5 - LA4)*SBMA
      GHPM = 2*GHAA + V*(LA5 - LA4)*CBMA
      GLZZ = 2*GLZZ
      GHZZ = 2*GHZZ
      GLLL = 6*GLLL
      GHHH = 6*GHHH
      GLHH = 2*GLHH
      GHLL = 2*GHLL
      GLAA = 2*GLAA
      GHAA = 2*GHAA
      XNORM = AMZ**2/V
      GLLL = GLLL/XNORM
      GHLL = GHLL/XNORM
      GLHH = GLHH/XNORM
      GHHH = GHHH/XNORM
      GHAA = GHAA/XNORM
      GLAA = GLAA/XNORM
      GLPM = GLPM/XNORM
      GHPM = GHPM/XNORM
      GAT=1.D0/TGBET
      GAB=TGBET
      GLT=CA/SB
      GLB=-SA/CB
      GHT=SA/SB
      GHB=CA/CB
      GZAL=-CBMA
      GZAH=SBMA
      GLVV=SBMA
      GHVV=CBMA
      B=BET
      A=DATAN(SA/CA)
      IF(CA.LT.0D0)THEN
       IF(SA.LT.0D0)THEN
        A = A-PI
       ELSE
        A = A+PI
       ENDIF
      ENDIF
      RETURN
      END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     THIS PROGRAM COMPUTES THE RENORMALIZATION GROUP IMPROVED
C     VALUES OF HIGGS MASSES AND COUPLINGS IN THE MSSM.
C
C     INPUT: MA,TANB = TAN(BETA),MQ,MUR,MTOP,AU,AD,MU.
C
C     ALL MASSES IN GEV UNITS. MA IS THE CP-ODD HIGGS MASS,
C     MTOP IS THE PHYSICAL TOP MASS, MQ AND MUR ARE THE SOFT
C     SUPERSYMMETRY BREAKING MASS PARAMETERS OF LEFT HANDED
C     AND RIGHT HANDED STOPS RESPECTIVELY, AU AND AD ARE THE
C     STOP AND SBOTTOM TRILINEAR SOFT BREAKING TERMS,
C     RESPECTIVELY,  AND MU IS THE SUPERSYMMETRIC
C     HIGGS MASS PARAMETER. WE USE THE  CONVENTIONS FROM
C     THE PHYSICS REPORT OF HABER AND KANE: LEFT RIGHT
C     STOP MIXING TERM PROPORTIONAL TO (AU - MU/TANB).
C
C     WE USE AS INPUT TANB DEFINED AT THE SCALE MTOP.
C
C     OUTPUT: MH,HM,MHCH, SA = SIN(ALPHA), CA= COS(ALPHA), TANBA
C
C     WHERE MH AND HM ARE THE LIGHTEST AND HEAVIEST CP-EVEN
C     HIGGS MASSES, MHCH IS THE CHARGED HIGGS MASS AND
C     ALPHA IS THE HIGGS MIXING ANGLE.
C
C     TANBA IS THE ANGLE TANB AT THE CP-ODD HIGGS MASS SCALE.
C
C     RANGE OF VALIDITY:
C
C    (STOP1**2 - STOP2**2)/(STOP2**2 + STOP1**2) < 0.5
C    (SBOT1**2 - SBOT2**2)/(SBOT2**2 + SBOT2**2) < 0.5
C
C     WHERE STOP1, STOP2, SBOT1 AND SBOT2 ARE THE STOP AND
C     ARE THE SBOTTOM  MASS EIGENVALUES, RESPECTIVELY. THIS
C     RANGE AUTOMATICALLY EXCLUDES THE EXISTENCE OF TACHYONS.
C
C
C     FOR THE CHARGED HIGGS MASS COMPUTATION, THE METHOD IS
C     VALID IF
C
C     2 * |MB * AD* TANB|  < M_SUSY**2,  2 * |MTOP * AU| < M_SUSY**2
C
C     2 * |MB * MU * TANB| < M_SUSY**2,  2 * |MTOP * MU| < M_SUSY**2
C
C     WHERE M_SUSY**2 IS THE AVERAGE OF THE SQUARED STOP MASS
C     EIGENVALUES, M_SUSY**2 = (STOP1**2 + STOP2**2)/2. THE SBOTTOM
C     MASSES HAVE BEEN ASSUMED TO BE OF ORDER OF THE STOP ONES.
C
C     M_SUSY**2 = (MQ**2 + MUR**2)*0.5 + MTOP**2
C
C     PROGRAM BASED ON THE WORK BY M. CARENA, J.R. ESPINOSA,
C     M. QUIROS AND C.E.M. WAGNER, CERN-PREPRINT CERN-TH/95-45.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      SUBROUTINE SUBH (MA,TANB,MQ,MUR,MTOP,AU,AD,MU,MH,HM,
     * MHCH,SA,CA,TANBA)
      IMPLICIT REAL*8(A-H,L,M,O-Z)
      COMMON/PARAM/S,AM,SCF,EPS,GF,AMW,AMZ,SW2,GEVPB,NPROC
      COMMON/HSELF/LAMBDA1,LAMBDA2,LAMBDA3,LAMBDA4,LAMBDA5,
     .             LAMBDA6,LAMBDA7
C     MZ = 91.18
C     ALPHA1 = 0.0101
C     ALPHA2 = 0.0337
C     ALPHA3Z = 0.12
C     V = 174.1
C     PI = 3.14159
      TANBA = TANB
      TANBT = TANB

C     MBOTTOM(MTOP) = 3. GEV
C     MB = 3.
C     ALPHA3 = ALPHA3Z/(1. +(11. - 10./3.)/4./PI*ALPHA3Z*
C    *LOG(MTOP**2/MZ**2))

C     RMTOP= RUNNING TOP QUARK MASS
C     RMTOP = MTOP/(1.+4.*ALPHA3/3./PI)
C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
      MB = RUNM(MTOP,5)
      PI = 4*DATAN(1D0)
      MZ = AMZ
      V  = 1/DSQRT(2*DSQRT(2D0)*GF)
      CW = AMW**2/AMZ**2
      SW = 1-CW
      ALPHA2  = (2*AMW/V/DSQRT(2D0))**2/4/PI
      ALPHA1  = ALPHA2*SW/CW
      ALPHA3Z = ALPHAS(AMZ,2)
      ALPHA3  = ALPHAS(MTOP,2)
      RMTOP   = RUNM(MTOP,6)
C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C      RMTOP=MTOP
C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
      MS = ((MQ**2 + MUR**2)/2. + MTOP**2)**.5
      T = LOG(MS**2/MTOP**2)
      SINB = TANB/((1. + TANB**2)**.5)
      COSB = SINB/TANB
C      IF(MA.LE.MTOP) TANBA = TANBT
      IF(MA.GT.MTOP)
     *TANBA = TANBT*(1.-3./32./PI**2*
     *(RMTOP**2/V**2/SINB**2-MB**2/V**2/COSB**2)*
     *LOG(MA**2/MTOP**2))

      SINBT = TANBT/((1. + TANBT**2)**.5)
      COSBT = 1./((1. + TANBT**2)**.5)
      COS2BT = (TANBT**2 - 1.)/(TANBT**2 + 1.)
      G1 = (ALPHA1*4.*PI)**.5
      G2 = (ALPHA2*4.*PI)**.5
      G3 = (ALPHA3*4.*PI)**.5
      HU = RMTOP/V/SINBT
      HD =  MB/V/COSBT

C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C      G3=0
C      HU=0
C      HD=0
C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

      XAU = (2.*AU**2/MS**2)*(1. - AU**2/12./MS**2)
      XAD = (2.*AD**2/MS**2)*(1. - AD**2/12./MS**2)
      AUD = (-6.*MU**2/MS**2 - ( MU**2- AD*AU)**2/MS**4.
     *+ 3.*(AU + AD)**2/MS**2)/6.
      LAMBDA1 = ((G1**2 + G2**2)/4.)*(1.-3.*HD**2*T/8./PI**2)
     *+(3.*HD**4/8./PI**2) * (T + XAD/2. + (3.*HD**2/2. + HU**2/2.
     *- 8.*G3**2) * (XAD*T + T**2)/16./PI**2)
     *-(3.*HU**4* MU**4/96./PI**2/MS**4) * (1+ (9.*HU**2 -5.* HD**2
     *-  16.*G3**2) *T/16./PI**2)
      LAMBDA2 = ((G1**2 + G2**2)/4.)*(1.-3.*HU**2*T/8./PI**2)
     *+(3.*HU**4/8./PI**2) * (T + XAU/2. + (3.*HU**2/2. + HD**2/2.
     *- 8.*G3**2) * (XAU*T + T**2)/16./PI**2)
     *-(3.*HD**4* MU**4/96./PI**2/MS**4) * (1+ (9.*HD**2 -5.* HU**2
     *-  16.*G3**2) *T/16./PI**2)
      LAMBDA3 = ((G2**2 - G1**2)/4.)*(1.-3.*
     *(HU**2 + HD**2)*T/16./PI**2)
     *+(6.*HU**2*HD**2/16./PI**2) * (T + AUD/2. + (HU**2 + HD**2
     *- 8.*G3**2) * (AUD*T + T**2)/16./PI**2)
     *+(3.*HU**4/96./PI**2) * (3.*MU**2/MS**2 - MU**2*AU**2/
     *MS**4)* (1.+ (6.*HU**2 -2.* HD**2/2.
     *-  16.*G3**2) *T/16./PI**2)
     *+(3.*HD**4/96./PI**2) * (3.*MU**2/MS**2 - MU**2*AD**2/
     *MS**4)*(1.+ (6.*HD**2 -2.* HU**2
     *-  16.*G3**2) *T/16./PI**2)
      LAMBDA4 = (- G2**2/2.)*(1.-3.*(HU**2 + HD**2)*T/16./PI**2)
     *-(6.*HU**2*HD**2/16./PI**2) * (T + AUD/2. + (HU**2 + HD**2
     *- 8.*G3**2) * (AUD*T + T**2)/16./PI**2)
     *+(3.*HU**4/96./PI**2) * (3.*MU**2/MS**2 - MU**2*AU**2/
     *MS**4)*
     *(1+ (6.*HU**2 -2.* HD**2
     *-  16.*G3**2) *T/16./PI**2)
     *+(3.*HD**4/96./PI**2) * (3.*MU**2/MS**2 - MU**2*AD**2/
     *MS**4)*
     *(1+ (6.*HD**2 -2.* HU**2/2.
     *-  16.*G3**2) *T/16./PI**2)
      LAMBDA5 = -(3.*HU**4* MU**2*AU**2/96./PI**2/MS**4) *
     * (1- (2.*HD**2 -6.* HU**2 + 16.*G3**2) *T/16./PI**2)
     *-(3.*HD**4* MU**2*AD**2/96./PI**2/MS**4) *
     * (1- (2.*HU**2 -6.* HD**2 + 16.*G3**2) *T/16./PI**2)
      LAMBDA6 = (3.*HU**4* MU**3*AU/96./PI**2/MS**4) *
     * (1- (7.*HD**2/2. -15.* HU**2/2. + 16.*G3**2) *T/16./PI**2)
     *+(3.*HD**4* MU *(AD**3/MS**3 - 6.*AD/MS )/96./PI**2/MS) *
     * (1- (HU**2/2. -9.* HD**2/2. + 16.*G3**2) *T/16./PI**2)
      LAMBDA7 = (3.*HD**4* MU**3*AD/96./PI**2/MS**4) *
     * (1- (7.*HU**2/2. -15.* HD**2/2. + 16.*G3**2) *T/16./PI**2)
     *+(3.*HU**4* MU *(AU**3/MS**3 - 6.*AU/MS )/96./PI**2/MS) *
     * (1- (HD**2/2. -9.* HU**2/2. + 16.*G3**2) *T/16./PI**2)
      TRM2 = MA**2 + 2.*V**2* (LAMBDA1* COSBT**2 +
     *2.* LAMBDA6*SINBT*COSBT
     *+ LAMBDA5*SINBT**2 + LAMBDA2* SINBT**2 + 2.* LAMBDA7*SINBT*COSBT
     *+ LAMBDA5*COSBT**2)
      DETM2 = 4.*V**4*(-(SINBT*COSBT*(LAMBDA3 + LAMBDA4) +
     *LAMBDA6*COSBT**2
     *+ LAMBDA7* SINBT**2)**2 + (LAMBDA1* COSBT**2 +
     *2.* LAMBDA6* COSBT*SINBT
     *+ LAMBDA5*SINBT**2)*(LAMBDA2* SINBT**2 +2.* LAMBDA7* COSBT*SINBT
     *+ LAMBDA5*COSBT**2)) + MA**2*2.*V**2 *
     *((LAMBDA1* COSBT**2 +2.*
     *LAMBDA6* COSBT*SINBT + LAMBDA5*SINBT**2)*COSBT**2 +
     *(LAMBDA2* SINBT**2 +2.* LAMBDA7* COSBT*SINBT + LAMBDA5*COSBT**2)
     **SINBT**2
     * +2.*SINBT*COSBT* (SINBT*COSBT*(LAMBDA3
     * + LAMBDA4) + LAMBDA6*COSBT**2
     *+ LAMBDA7* SINBT**2))

      MH2 = (TRM2 - (TRM2**2 - 4.* DETM2)**.5)/2.
      HM2 = (TRM2 + (TRM2**2 - 4.* DETM2)**.5)/2.
      HM = HM2**.5
      MH = MH2**.5
      MHCH2 = MA**2 + (LAMBDA5 - LAMBDA4)* V**2
      MHCH = MHCH2**.5
      MHCH = MHCH2**.5

      SINALPHA = SQRT(((TRM2**2 - 4.* DETM2)**.5) -
     * ((2.*V**2*(LAMBDA1* COSBT**2 + 2.*
     *LAMBDA6* COSBT*SINBT
     *+ LAMBDA5*SINBT**2) + MA**2*SINBT**2)
     *- (2.*V**2*(LAMBDA2* SINBT**2 +2.* LAMBDA7* COSBT*SINBT
     *+ LAMBDA5*COSBT**2) + MA**2*COSBT**2)))/
     *SQRT(((TRM2**2 - 4.* DETM2)**.5))/2.**.5

      COSALPHA = (2.*(2.*V**2*(SINBT*COSBT*(LAMBDA3 + LAMBDA4) +
     *LAMBDA6*COSBT**2 + LAMBDA7* SINBT**2) -
     *MA**2*SINBT*COSBT))/2.**.5/
     *SQRT(((TRM2**2 - 4.* DETM2)**.5)*
     *(((TRM2**2 - 4.* DETM2)**.5) -
     * ((2.*V**2*(LAMBDA1* COSBT**2 + 2.*
     *LAMBDA6* COSBT*SINBT
     *+ LAMBDA5*SINBT**2) + MA**2*SINBT**2)
     *- (2.*V**2*(LAMBDA2* SINBT**2 +2.* LAMBDA7* COSBT*SINBT
     *+ LAMBDA5*COSBT**2) + MA**2*COSBT**2))))

      SA = -SINALPHA
      CA = -COSALPHA

 2242 RETURN
      END

      SUBROUTINE AMHAMA (ICASE,MH,MA,TANB,MQ,MUR,MTOP,AU,AD,MU)
C--CALCULATION OF PSEUDOSCALAR HIGGS MASS FROM HIGGS MASS MH
C--ICASE=0: MH=PSEUDOSCALAR MASS
C--ICASE=1: MH=LIGHT SCALAR MASS
C--ICASE=2: MH=HEAVY SCALAR MASS
C--ICASE=3: MH=CHARGED HIGGS MASS
      IMPLICIT REAL*8(A-H,L,M,O-Z)
      COMMON/PARAM/S,AM,SCF,EPS,GF,AMW,AMZ,SW2,GEVPB,NPROC
C     MZ = 91.18
C     ALPHA1 = 0.0101
C     ALPHA2 = 0.0337
C     ALPHA3Z = 0.12
C     V = 174.1
C     PI = 3.14159
      TANBI = TANB
      TANBA = TANB
      TANBT = TANB

C     MBOTTOM(MTOP) = 3. GEV
C     MB = 3.

C     ALPHA3 = ALPHA3Z/(1. +(11. - 10./3.)/4./PI*ALPHA3Z*
C    *LOG(MTOP**2./MZ**2.))

C     RMTOP= RUNNING TOP QUARK MASS
C     RMTOP = MTOP/(1.+4.*ALPHA3/3./PI)

C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
      MB = RUNM(MTOP,5)
      PI = 4*DATAN(1D0)
      MZ = AMZ
      V  = 1/DSQRT(2*DSQRT(2D0)*GF)
      CW = AMW**2/AMZ**2
      SW = 1-CW
      ALPHA2  = (2*AMW/V/DSQRT(2D0))**2/4/PI
      ALPHA1  = ALPHA2*SW/CW
      ALPHA3Z = ALPHAS(AMZ,2)
      ALPHA3  = ALPHAS(MTOP,2)
      RMTOP   = RUNM(MTOP,6)
C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C      RMTOP=MTOP
C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

      MS = ((MQ**2. + MUR**2.)/2. + MTOP**2.)**.5
      T = LOG(MS**2./MTOP**2.)
      SINB = TANBI/((1. + TANBI**2.)**.5)
      COSB = SINB/TANBI
C     IF(MA.LE.MTOP) TANBT = TANBI
      IF(MA.GT.MTOP)
     *TANBA = TANBT*(1.-3./32./PI**2.*
     *(RMTOP**2./V**2./SINB**2.-MB**2./V**2./COSB**2.)*
     *LOG(MA**2./MTOP**2.))

      SINBT = TANBT/((1. + TANBT**2.)**.5)
      COSBT = 1./((1. + TANBT**2.)**.5)
      COS2BT = (TANBT**2. - 1.)/(TANBT**2. + 1.)
      G1 = (ALPHA1*4.*PI)**.5
      G2 = (ALPHA2*4.*PI)**.5
      G3 = (ALPHA3*4.*PI)**.5
      HU = RMTOP/V/SINBT
      HD =  MB/V/COSBT

C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C      G3=0
C      HU=0
C      HD=0
C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

      XAU = (2.*AU**2./MS**2.)*(1. - AU**2./12./MS**2.)
      XAD = (2.*AD**2./MS**2.)*(1. - AD**2./12./MS**2.)
      AUD = (-6.*MU**2/MS**2. - ( MU**2- AD*AU)**2./MS**4.
     *+ 3.*(AU + AD)**2./MS**2.)/6.
      LAMBDA1 = ((G1**2. + G2**2.)/4.)*(1.-3.*HD**2.*T/8./PI**2.)
     *+(3.*HD**4./8./PI**2.) * (T + XAD/2. + (3.*HD**2./2. + HU**2./2.
     *- 8.*G3**2.) * (XAD*T + T**2.)/16./PI**2.)
     *-(3.*HU**4.* MU**4/96./PI**2./MS**4.) * (1+ (9.*HU**2. -5.* HD**2.
     *-  16.*G3**2.) *T/16./PI**2.)
      LAMBDA2 = ((G1**2. + G2**2.)/4.)*(1.-3.*HU**2.*T/8./PI**2.)
     *+(3.*HU**4./8./PI**2.) * (T + XAU/2. + (3.*HU**2./2. + HD**2./2.
     *- 8.*G3**2.) * (XAU*T + T**2.)/16./PI**2.)
     *-(3.*HD**4.* MU**4/96./PI**2./MS**4.) * (1+ (9.*HD**2. -5.* HU**2.
     *-  16.*G3**2.) *T/16./PI**2.)
      LAMBDA3 = ((G2**2. - G1**2.)/4.)*(1.-3.*
     *(HU**2. + HD**2.)*T/16./PI**2.)
     *+(6.*HU**2.*HD**2./16./PI**2.) * (T + AUD/2. + (HU**2. + HD**2.
     *- 8.*G3**2.) * (AUD*T + T**2.)/16./PI**2.)
     *+(3.*HU**4./96./PI**2.) * (3.*MU**2/MS**2. - MU**2*AU**2./
     *MS**4.)* (1.+ (6.*HU**2. -2.* HD**2./2.
     *-  16.*G3**2.) *T/16./PI**2.)
     *+(3.*HD**4./96./PI**2.) * (3.*MU**2/MS**2. - MU**2*AD**2./
     *MS**4.)*(1.+ (6.*HD**2. -2.* HU**2./2.
     *-  16.*G3**2.) *T/16./PI**2.)
      LAMBDA4 = (- G2**2./2.)*(1.-3.*(HU**2. + HD**2.)*T/16./PI**2.)
     *-(6.*HU**2.*HD**2./16./PI**2.) * (T + AUD/2. + (HU**2. + HD**2.
     *- 8.*G3**2.) * (AUD*T + T**2.)/16./PI**2.)
     *+(3.*HU**4./96./PI**2.) * (3.*MU**2/MS**2. - MU**2*AU**2./
     *MS**4.)*
     *(1+ (6.*HU**2. -2.* HD**2./2.
     *-  16.*G3**2.) *T/16./PI**2.)
     *+(3.*HD**4./96./PI**2.) * (3.*MU**2/MS**2. - MU**2*AD**2./
     *MS**4.)*
     *(1+ (6.*HD**2. -2.* HU**2./2.
     *-  16.*G3**2.) *T/16./PI**2.)
      LAMBDA5 = -(3.*HU**4.* MU**2*AU**2./96./PI**2./MS**4.) *
     * (1- (2.*HD**2. -6.* HU**2. + 16.*G3**2.) *T/16./PI**2.)
     *-(3.*HD**4.* MU**2*AD**2./96./PI**2./MS**4.) *
     * (1- (2.*HU**2. -6.* HD**2. + 16.*G3**2.) *T/16./PI**2.)
      LAMBDA6 = (3.*HU**4.* MU**3*AU/96./PI**2./MS**4.) *
     * (1- (7.*HD**2./2. -15.* HU**2./2. + 16.*G3**2.) *T/16./PI**2.)
     *+(3.*HD**4.* MU *(AD**3./MS**3. - 6.*AD/MS )/96./PI**2./MS) *
     * (1- (HU**2./2. -9.* HD**2./2. + 16.*G3**2.) *T/16./PI**2.)
      LAMBDA7 = (3.*HD**4.* MU**3*AD/96./PI**2./MS**4.) *
     * (1- (7.*HU**2./2. -15.* HD**2./2. + 16.*G3**2.) *T/16./PI**2.)
     *+(3.*HU**4.* MU *(AU**3./MS**3. - 6.*AU/MS )/96./PI**2./MS) *
     * (1- (HD**2./2. -9.* HU**2./2. + 16.*G3**2.) *T/16./PI**2.)

      DEL1 = 2.*V**2.* (LAMBDA1* COSBT**2 +
     *2.* LAMBDA6*SINBT*COSBT
     *+ LAMBDA5*SINBT**2. + LAMBDA2* SINBT**2 + 2.* LAMBDA7*SINBT*COSBT
     *+ LAMBDA5*COSBT**2.)
      DEL2 = 4.*V**4.*(-(SINBT*COSBT*(LAMBDA3 + LAMBDA4) +
     *LAMBDA6*COSBT**2.
     *+ LAMBDA7* SINBT**2.)**2. + (LAMBDA1* COSBT**2. +
     *2.* LAMBDA6* COSBT*SINBT
     *+ LAMBDA5*SINBT**2.)*(LAMBDA2* SINBT**2. +2.* LAMBDA7* COSBT*SINBT
     *+ LAMBDA5*COSBT**2.))
      DELA = 2.*V**2. *
     *((LAMBDA1* COSBT**2. +2.*
     *LAMBDA6* COSBT*SINBT + LAMBDA5*SINBT**2.)*COSBT**2. +
     *(LAMBDA2* SINBT**2. +2.* LAMBDA7* COSBT*SINBT + LAMBDA5*COSBT**2.)
     **SINBT**2.
     * +2.*SINBT*COSBT* (SINBT*COSBT*(LAMBDA3
     * + LAMBDA4) + LAMBDA6*COSBT**2.
     *+ LAMBDA7* SINBT**2.))

      IF(ICASE.EQ.1.OR.ICASE.EQ.2)THEN
       XX = (MH**4 - MH**2*DEL1 + DEL2)/(MH**2 - DELA)
       IF(XX.GT.0D0)THEN
        MA = DSQRT(XX)
       ELSE
        MA = -1.D8
       ENDIF
      ELSEIF(ICASE.EQ.3)THEN
       XX = MH**2 - (LAMBDA5 - LAMBDA4)* V**2.
       IF(XX.GT.0D0)THEN
        MA = DSQRT(MH**2 - (LAMBDA5 - LAMBDA4)* V**2.)
       ELSE
        MA = -1.D8
       ENDIF
      ELSE
       MA = MH
      ENDIF

      RETURN
      END

      DOUBLE PRECISION FUNCTION RUNM(Q,NF)
C--RUNNING QUARK MASSES: Q = SCALE, NF = 3,..,6 QUARK TYPE
      PARAMETER (NN=6)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (ZETA3 = 1.202056903159594D0)
      DIMENSION AM(NN),YMSB(NN)
      COMMON/ALS/XLAMBDA,AMCA,AMBA,AMTA,N0A
      COMMON/MASSES/AMS,AMC,AMB,AMT
      COMMON/STRANGE/AMSB
      COMMON/RUN/XMSB,XMHAT,XKFAC
      COMMON/FLAG/IHIGGS,NNLO
      SAVE ISTRANGE
      B0(NF)=(33.D0-2.D0*NF)/12D0
      B1(NF) = (102D0-38D0/3D0*NF)/16D0
      B2(NF) = (2857D0/2D0-5033D0/18D0*NF+325D0/54D0*NF**2)/64D0
      G0(NF) = 1D0
      G1(NF) = (202D0/3D0-20D0/9D0*NF)/16D0
      G2(NF) = (1249D0-(2216D0/27D0+160D0/3D0*ZETA3)*NF
     .       - 140D0/81D0*NF**2)/64D0
      C1(NF) = G1(NF)/B0(NF) - B1(NF)*G0(NF)/B0(NF)**2
      C2(NF) = ((G1(NF)/B0(NF) - B1(NF)*G0(NF)/B0(NF)**2)**2
     .       + G2(NF)/B0(NF) + B1(NF)**2*G0(NF)/B0(NF)**3
     .       - B1(NF)*G1(NF)/B0(NF)**2 - B2(NF)*G0(NF)/B0(NF)**2)/2D0
      TRAN(X,XK)=1D0+4D0/3D0*ALPHAS(X,2)/PI+XK*(ALPHAS(X,2)/PI)**2
      CQ(X,NF)=(2D0*B0(NF)*X)**(G0(NF)/B0(NF))
     .            *(1D0+C1(NF)*X+C2(NF)*X**2)
      DATA ISTRANGE/0/
      PI=4D0*DATAN(1D0)
      ACC = 1.D-8
      AM(1) = 0
      AM(2) = 0
C--------------------------------------------
      IMSBAR = 0
      NNLO = 0
      IF(IMSBAR.EQ.1)THEN
       IF(ISTRANGE.EQ.0)THEN
C--STRANGE POLE MASS FROM MSBAR-MASS AT 1 GEV
        AMSD = XLAMBDA
        AMSU = 1.D8
123     AMS  = (AMSU+AMSD)/2
        AM(3) = AMS
        XMSB = AMS/CQ(ALPHAS(AMS,2)/PI,3)
     .            *CQ(ALPHAS(1.D0,2)/PI,3)/TRAN(AMS,0D0)
        DD = (XMSB-AMSB)/AMSB
        IF(DABS(DD).GE.ACC)THEN
         IF(DD.LE.0.D0)THEN
          AMSD = AM(3)
         ELSE
          AMSU = AM(3)
         ENDIF
         GOTO 123
        ENDIF
        ISTRANGE=1
       ENDIF
       AM(3) = AMSB
      ELSE
       AMS=AMSB
       AM(3) = AMS
      ENDIF
C--------------------------------------------
      AM(3) = AMSB
      AM(4) = AMC
      AM(5) = AMB
      AM(6) = AMT
      XK = 16.11D0
      DO 1 I=1,NF-1
       XK = XK - 1.04D0*(1.D0-AM(I)/AM(NF))
1     CONTINUE
      IF(NF.GE.4)THEN
       XMSB = AM(NF)/TRAN(AM(NF),0D0)
       XMHAT = XMSB/CQ(ALPHAS(AM(NF),2)/PI,NF)
      ELSE
       XMSB = 0
       XMHAT = 0
      ENDIF
      YMSB(3) = AMSB
      IF(NF.EQ.3)THEN
       YMSB(4) = YMSB(3)*CQ(ALPHAS(AM(4),2)/PI,3)/
     .                   CQ(ALPHAS(1.D0,2)/PI,3)
       YMSB(5) = YMSB(4)*CQ(ALPHAS(AM(5),2)/PI,4)/
     .                   CQ(ALPHAS(AM(4),2)/PI,4)
       YMSB(6) = YMSB(5)*CQ(ALPHAS(AM(6),2)/PI,5)/
     .                   CQ(ALPHAS(AM(5),2)/PI,5)
      ELSEIF(NF.EQ.4)THEN
       YMSB(4) = XMSB
       YMSB(5) = YMSB(4)*CQ(ALPHAS(AM(5),2)/PI,4)/
     .                   CQ(ALPHAS(AM(4),2)/PI,4)
       YMSB(6) = YMSB(5)*CQ(ALPHAS(AM(6),2)/PI,5)/
     .                   CQ(ALPHAS(AM(5),2)/PI,5)
      ELSEIF(NF.EQ.5)THEN
       YMSB(5) = XMSB
       YMSB(4) = YMSB(5)*CQ(ALPHAS(AM(4),2)/PI,4)/
     .                   CQ(ALPHAS(AM(5),2)/PI,4)
       YMSB(6) = YMSB(5)*CQ(ALPHAS(AM(6),2)/PI,5)/
     .                   CQ(ALPHAS(AM(5),2)/PI,5)
      ELSEIF(NF.EQ.6)THEN
       YMSB(6) = XMSB
       YMSB(5) = YMSB(6)*CQ(ALPHAS(AM(5),2)/PI,5)/
     .                   CQ(ALPHAS(AM(6),2)/PI,5)
       YMSB(4) = YMSB(5)*CQ(ALPHAS(AM(4),2)/PI,4)/
     .                   CQ(ALPHAS(AM(5),2)/PI,4)
      ENDIF
      IF(Q.LT.AMC)THEN
       N0=3
       Q0 = 1.D0
      ELSEIF(Q.LE.AMB)THEN
       N0=4
       Q0 = AMC
      ELSEIF(Q.LE.AMT)THEN
       N0=5
       Q0 = AMB
      ELSE
       N0=6
       Q0 = AMT
      ENDIF
      IF(NNLO.EQ.1.AND.NF.GT.3)THEN
       XKFAC = TRAN(AM(NF),0D0)/TRAN(AM(NF),XK)
      ELSE
       XKFAC = 1D0
      ENDIF
      RUNM = YMSB(N0)*CQ(ALPHAS(Q,2)/PI,N0)/
     .               CQ(ALPHAS(Q0,2)/PI,N0)
     .       * XKFAC
      RETURN
      END

C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
C      SUBROUTINE STRUC(X,Q,PDF)
CC--PARTON DENSITIES
C      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      DIMENSION PDF(-6:6), VALUE(20)
C      CHARACTER*20 PARM(20)
C      COMMON/PDFLIB/NGROUP,NSET
C      QQ=Q
C      PARM(1)='NPTYPE'
C      PARM(2)='NGROUP'
C      PARM(3)='NSET'
C      VALUE(1)=1.D0
C      VALUE(2)=DFLOAT(NGROUP)
C      VALUE(3)=DFLOAT(NSET)
C      CALL PDFSET(PARM,VALUE)
C      CALL PFTOPDG(X,QQ,PDF)
C      RETURN
C      END
 
c     subroutine struc(x,q,pdf)
c     implicit double precision (a-h,o-z)
c     dimension pdf(-6:6), value(20)
c     character*20 parm(20)
c     common/pdflib/ngroup,nset
c     if(ngroup.gt.0)then
c      parm(1)='nptype'
c      parm(2)='ngroup'
c      parm(3)='nset'
c      value(1)=1.d0
c      value(2)=dfloat(ngroup)
c      value(3)=dfloat(nset)
c      call pdfset(parm,value)
c      call pftopdg(x,q,pdf)
c     elseif(ngroup.eq.-1)then
c      call SetCtq6(nset)
c      pdf(6)  = 0
c      pdf(-6) = 0
c      do i=-5,5
c       j = i
c       if(i.eq.1)j=2
c       if(i.eq.2)j=1
c       if(i.eq.-1)j=-2
c       if(i.eq.-2)j=-1
c       pdf(j) = x*Ctq6Pdf(i,x,q)
c      enddo
c     else
c      mode = nset
c      call mrst2001(x,q,mode,upv,dnv,usea,dsea,str,chm,bot,glu)
c      pdf(-6) = 0
c      pdf(-5) = bot
c      pdf(-4) = chm
c      pdf(-3) = str
c      pdf(-2) = usea
c      pdf(-1) = dsea
c      pdf(0) = glu
c      pdf(1) = dnv + dsea
c      pdf(2) = upv + usea
c      pdf(3) = str
c      pdf(4) = chm
c      pdf(5) = bot
c      pdf(6) = 0
c     endif
c     pdf( 6) = 0
c     pdf(-6) = 0
c     return
c     end

      subroutine struc(x,q,pdf)
      implicit double precision (a-h,o-z)
      dimension pdf(-6:6), value(20)
      common/pdflib/ngroup,nset
      ipdflib = ngroup
      if(ipdflib.eq.1)then
       call SetCtq6(nset)
       pdf(6)  = 0
       pdf(-6) = 0
       do i=-5,5
        j = i
        if(i.eq.1)j=2
        if(i.eq.2)j=1
        if(i.eq.-1)j=-2
        if(i.eq.-2)j=-1
        pdf(j) = x*Ctq6Pdf(i,x,q)
       enddo
      else
       call evolvePDF(x,q,pdf)
      endif
      pdf( 6) = 0
      pdf(-6) = 0
      return
      end

      subroutine pdfset(pathname,pdfname)
      implicit double precision (a-h,o-z)
      character*100 pdfname, pathname
      common/pdflib/ngroup,nset
      if(ngroup.eq.0)then
       call SetPDFpath(pathname)
       call InitPDFsetByName(pdfname)
       call InitPDF(nset)
      endif
      return
      end
