Files
V2HV/zzz/pphv.f
Basil Bruhn 2d2fbd40e1 initial
Signed-off-by: Basil Bruhn <basil.bruhn@psi.ch>
2025-11-17 10:17:15 +01:00

1792 lines
53 KiB
Fortran

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)
c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
SCA = SCF*(AM+AMV)/2
c write(6,*)SCA
c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
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