ltpth/VV2H
0
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
main
VV2H/pphqq.f
Basil Bruhn 6c1ead93f4 initial 
Signed-off-by: Basil Bruhn <basil.bruhn@psi.ch>
2025-11-17 10:14:53 +01:00

2268 lines
67 KiB
Fortran
Raw Permalink Blame History

C**********************************************************************
C
C ****************
C * VERSION 1.10 *
C ****************
C
C--This program calculates the production cross section of Higgs
C bosons via WW/ZZ -> h,H at hadron colliders at NLO QCD
C according to the formulae presented in
C
C T. Han, G. Valencia and S. Willenbrock, Phys. Rev. Lett. 69
C (1992) 3274.
C
C Interference effects between W- and Z-boson fusion are neglected,
C since they are less than 1%. The program allows to calculate the
C total production cross section in the DIS as well as the MSbar scheme
C for the parton densities for the scalar Higgs bosons of the SM and
C MSSM. The MSSM Higgs sector is implemented in the approximate
C 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 pphqq.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 pp (0) or ppbar (1) collider
C PP/PPBAR = 0
C
C--Energy of hadron collider in GeV
C ENERGY = 14000.D0
C
C--Choice of scale for parton densities: 0 -> mu_0 = M_H
C 1 -> mu_0 = Q_1,Q_2
C MH=0/Q=1 = 1
C
C--Scale loop: mu = xi * mu_0, xi = SCALE1 -> SCALE2 with NSCALE
C equidistant points in total
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--Parton densities: 0 = MSbar scheme, 1 = DIS scheme
C MS0/DIS1 = 0
C
C--LHAPDF: PDF path and name
C PDFPATH = Grids
C PDFNAME = MSTW2008nlo68cl.LHgrid
C
C--VEGAS: IPOINT1 = points for Born term
C IPOINT = points for NLO terms
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 IPOINT1 = 10000
C IPOINT = 100000
C ITERAT = 5
C NPRN = 10
C
C--Substitution for integration over z in NLO (should be left as it is)
C SUBST12Z = 1
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, proton mass
C MW = 80.41D0
C MZ = 91.187D0
C SW2 = 0.2315D0
C GF = 1.16639D-5
C MPROTON = 0.93827231D0
C
C--Lower cut in Q_1, Q_2: -Q_i**2 > Q0**2
C Q0 = 2.D0
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 PPHQQ
C--PROGRAM FOR WW/ZZ -> H AT HADRON COLLIDERS
PARAMETER(NIN=95,NOUT=96)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CHARACTER*100 PDFNAME,PATHNAME
COMMON/PARAM/S,AM,QQ,EPS,GF,AMW,AMZ,SW2,XMP,GEVPB
COMMON/ALS/XLAMBDA,AMC,AMB,AMT,N0
COMMON/RESULT/RES,ERR,DCHI2,DUM
COMMON/PDFLIB/NGROUP,NSET
COMMON/QCUT/Q0
COMMON/SUBSTIT/NSUBST
COMMON/MSDIS/NDIS
COMMON/COLL/ICOLL
COMMON/UNIT/NUNIT
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='pphqq.in')
OPEN(NOUT,FILE='pphqq.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)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,100)W
READ(NIN,*)
READ(NIN,101)NUNIT
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,101)NDIS
READ(NIN,102)PATHNAME
READ(NIN,102)PDFNAME
READ(NIN,*)
READ(NIN,101)IPOINT1
READ(NIN,101)IPOINT
READ(NIN,101)ITERAT
READ(NIN,101)NPRN
READ(NIN,101)NSUBST
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,100)XMP
READ(NIN,100)Q0
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--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
WRITE(NOUT,*)' PP ---> HQQ + X'
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
QQ = SCALE*AM
C--INTEGRATION OF BORN TERM
IDIM = 6
CALL VEGASN(SIG,ABSERR,IDIM,IPOINT1,ITERAT,NPRN,IGRAPH)
SIGB = RES*COUPLING
DSIGB = ERR*COUPLING
IF(LOOP.NE.1)THEN
C--INTEGRATION OF NLO QQBAR INITIAL STATE
IDIM = 7
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,QQ,EPS,GF,AMW,AMZ,SW2,XMP,GEVPB
DIMENSION XX(6),Y(6)
PI = 4*DATAN(1.D0)
DO I=1,6
Y(I) = EPS + (1-2*EPS)*XX(I)
ENDDO
SIG = GEVPB*XMAT(Y)*(1-2*EPS)**6
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,QQ,EPS,GF,AMW,AMZ,SW2,XMP,GEVPB
DIMENSION XX(7),Y(7)
EXTERNAL QF1,QF2,QF3
DO I=1,7
Y(I) = EPS + (1-2*EPS)*XX(I)
ENDDO
Y(7) = XX(7)
DSIGQQ = GEVPB*DMAT(Y,QF1,QF2,QF3)*(1-2*EPS)**6
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,QQ,EPS,GF,AMW,AMZ,SW2,XMP,GEVPB
DIMENSION XX(7),Y(7)
EXTERNAL GF1,GF2,GF3
DO I=1,7
Y(I) = EPS + (1-2*EPS)*XX(I)
ENDDO
Y(7) = XX(7)
DSIGGQ = GEVPB*DMAT(Y,GF1,GF2,GF3)*(1-2*EPS)**6
RETURN
END
DOUBLE PRECISION FUNCTION ZMAT(F11,F12,F21,F22,F31,F32)
C--GENERIC FORM OF MATRIX ELEMENTS
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
COMMON/VMAT/Q1Q1,Q1Q2,Q2Q2,P1Q1,P1Q2,P2Q1,P2Q2,P1P2,XMP
ZMAT = F11*F12*(2+Q1Q2**2/Q1Q1/Q2Q2)
. + F11*F22/P2Q2*(P2Q2**2/Q2Q2 - XMP**2
. + (P2Q1-P2Q2*Q1Q2/Q2Q2)**2/Q1Q1)
. + F21*F12/P1Q1*(P1Q1**2/Q1Q1 - XMP**2
. + (P1Q2-P1Q1*Q1Q2/Q1Q1)**2/Q2Q2)
. + F21*F22/P1Q1/P2Q2*(P1P2 - P1Q1*P2Q1/Q1Q1
. - P2Q2*P1Q2/Q2Q2 + P1Q1*P2Q2*Q1Q2/Q1Q1/Q2Q2)**2
. + F31*F32/2/P1Q1/P2Q2*(P1P2*Q1Q2 - P1Q2*P2Q1)
RETURN
END
DOUBLE PRECISION FUNCTION XMAT(Y)
C--BORN MATRIX ELEMENT AND KINEMATICS
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION PDF(-6:6), Y(6)
COMMON/COLL/ICOLL
COMMON/PARAM/S,AM,QQ,EPS,GF,AMW,AMZ,SW2,XMP,GEVPB
COMMON/VMAT/Q1Q1,Q1Q2,Q2Q2,P1Q1,P1Q2,P2Q1,P2Q2,P1P2,XXMP
COMMON/UNIT/NUNIT
COMMON/QCUT/Q0
COMMON/SUBSTIT/NSUBST
PI = 4*DATAN(1.D0)
XXMP = XMP
RES = 0
DO NPROC=-1,1
IF(NPROC.EQ.0)THEN
AMV = AMZ
I1=0
I2=0
ELSEIF(NPROC.EQ.-1)THEN
AMV = AMW
I1=-1
I2=1
ELSEIF(NPROC.EQ.1)THEN
AMV = AMW
I1=1
I2=-1
ENDIF
TAUH = AM**2/S
TAU = TAUH + (1-TAUH)*Y(1)
Y1 = -DLOG(TAU)*Y(2)
X1 = DEXP(-Y1)
X2 = TAU/X1
SHAT = X1*X2*S
XMU = AM**2/SHAT
XX1 = (1-XMU)*Y(3)
XX2 = (1-XX1-XMU)*(1-Y(4)) + (1-XMU/(1-XX1))*Y(4)
XX3 = 2-XX1-XX2
CT12 = 2/XX1/XX2*(1-XX1-XX2+XX1*XX2/2-XMU)
ST12 = DSQRT(1-CT12**2)
IF(NSUBST.NE.0)THEN
SQ10 = 1/(AMV**2+XX1*SHAT)
SQ11 = 1/(AMV**2+Q0**2)
SQ1Q1 = SQ10 + (SQ11-SQ10)*Y(5)
Q1Q1 = -1/SQ1Q1 + AMV**2
PRO1 = Q1Q1 - AMV**2
CTH = 1+Q1Q1*2/SHAT/XX1
STH = DSQRT(1-CTH**2)
Q2M = -XX2*SHAT/2*(1+CT12*CTH-ST12*STH)
Q2P = -XX2*SHAT/2*(1+CT12*CTH+ST12*STH)
SQ20 = 1/(AMV**2-Q2P)
SQ21 = 1/(AMV**2-Q2M)
SQ2Q2 = SQ20 + (SQ21-SQ20)*Y(6)
Q2Q2 = -1/SQ2Q2 + AMV**2
PRO2 = Q2Q2 - AMV**2
CCHI = -(2*Q2Q2/XX2/SHAT+1+CT12*CTH)/ST12/STH
SCHI = DSQRT(1-CCHI**2)
FAC = -(1-TAUH)*DLOG(TAU)*(1-XMU)
. * (XX1+XMU-XMU/(1-XX1))
. * 2/SHAT/XX1*PRO1**2*(SQ11-SQ10)
. * 4/SHAT/XX2*PRO2**2/ST12/STH/SCHI*(SQ21-SQ20)
ELSE
CTH = -1 + 2*Y(5)
CHI = PI*Y(6)
FAC = -(1-TAUH)*DLOG(TAU)*(1-XMU)
. * (XX1+XMU-XMU/(1-XX1))*2*2*PI
CCHI = DCOS(CHI)
STH = DSQRT(1-CTH**2)
ENDIF
I1I2 = SHAT/2
F1I1 = SHAT/4*XX1*(1-CTH)
F1I2 = SHAT/4*XX1*(1+CTH)
F2I1 = SHAT/4*XX2*(1-ST12*CCHI*STH-CT12*CTH)
F2I2 = SHAT/4*XX2*(1+ST12*CCHI*STH+CT12*CTH)
F3I1 = SHAT/2 - F1I1 - F2I1
F3I2 = SHAT/2 - F1I2 - F2I2
F1F2 = SHAT/2*(1-XX3+XMU)
F1F3 = SHAT/2*(1-XX2-XMU)
F2F3 = SHAT/2*(1-XX1-XMU)
Q1Q1 = -2*F1I1
Q2Q2 = -2*F2I2
Q1Q2 = F1F2 + I1I2 - F1I2 - F2I1
P1Q1 = F1I1/X1
P2Q2 = F2I2/X2
P1Q2 = (F2I1-I1I2)/X1
P2Q1 = (F1I2-I1I2)/X2
P1P2 = S/2
IF(NUNIT.EQ.1)THEN
SCALE = QQ/AM
Q1 = SCALE*DSQRT(-Q1Q1)
Q2 = SCALE*DSQRT(-Q2Q2)
ELSE
Q1 = QQ
Q2 = QQ
ENDIF
IF((-Q1Q1).GT.Q0**2.AND.(-Q2Q2).GT.Q0**2)THEN
PRO1 = Q1Q1 - AMV**2
PRO2 = Q2Q2 - AMV**2
F11 = F1(I1,X1,Q1,1)
F12 = F1(I2,X2,Q2,ICOLL)
F21 = F2(I1,X1,Q1,1)
F22 = F2(I2,X2,Q2,ICOLL)
F31 = F3(I1,X1,Q1,1)
F32 = F3(I2,X2,Q2,ICOLL)
ZZ = ZMAT(F11,F12,F21,F22,F31,F32)
RES = RES + DSQRT(2.D0)*GF**3*AMV**8/PRO1**2/PRO2**2*ZZ
. * F1I1*F2I2
. * FAC / 512/PI**4
ENDIF
ENDDO
XMAT = RES
RETURN
END
DOUBLE PRECISION FUNCTION DMAT(Y,DF1,DF2,DF3)
C--NLO MATRIX ELEMENT AND KINEMATICS
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION PDF(-6:6), Y(7)
COMMON/COLL/ICOLL
COMMON/PARAM/S,AM,QQ,EPS,GF,AMW,AMZ,SW2,XMP,GEVPB
COMMON/VMAT/Q1Q1,Q1Q2,Q2Q2,P1Q1,P1Q2,P2Q1,P2Q2,P1P2,XXMP
COMMON/UNIT/NUNIT
COMMON/QCUT/Q0
COMMON/SUBSTIT/NSUBST
EXTERNAL DF1,DF2,DF3
PI = 4*DATAN(1.D0)
XXMP = XMP
RES = 0
DO NPROC=-1,1
IF(NPROC.EQ.0)THEN
AMV = AMZ
I1=0
I2=0
ELSEIF(NPROC.EQ.-1)THEN
AMV = AMW
I1=-1
I2=1
ELSEIF(NPROC.EQ.1)THEN
AMV = AMW
I1=1
I2=-1
ENDIF
TAUH = AM**2/S
TAU = TAUH + (1-TAUH)*Y(1)
Y1 = -DLOG(TAU)*Y(2)
X1 = DEXP(-Y1)
X2 = TAU/X1
SHAT = X1*X2*S
XMU = AM**2/SHAT
XX1 = (1-XMU)*Y(3)
XX2 = (1-XX1-XMU)*(1-Y(4)) + (1-XMU/(1-XX1))*Y(4)
XX3 = 2-XX1-XX2
CT12 = 2/XX1/XX2*(1-XX1-XX2+XX1*XX2/2-XMU)
ST12 = DSQRT(1-CT12**2)
IF(NSUBST.NE.0)THEN
SQ10 = 1/(AMV**2+XX1*SHAT)
SQ11 = 1/(AMV**2+Q0**2)
SQ1Q1 = SQ10 + (SQ11-SQ10)*Y(5)
Q1Q1 = -1/SQ1Q1 + AMV**2
PRO1 = Q1Q1 - AMV**2
CTH = 1+Q1Q1*2/SHAT/XX1
STH = DSQRT(1-CTH**2)
Q2M = -XX2*SHAT/2*(1+CT12*CTH-ST12*STH)
Q2P = -XX2*SHAT/2*(1+CT12*CTH+ST12*STH)
SQ20 = 1/(AMV**2-Q2P)
SQ21 = 1/(AMV**2-Q2M)
SQ2Q2 = SQ20 + (SQ21-SQ20)*Y(6)
Q2Q2 = -1/SQ2Q2 + AMV**2
PRO2 = Q2Q2 - AMV**2
CCHI = -(2*Q2Q2/XX2/SHAT+1+CT12*CTH)/ST12/STH
SCHI = DSQRT(1-CCHI**2)
FAC = -(1-TAUH)*DLOG(TAU)*(1-XMU)
. * (XX1+XMU-XMU/(1-XX1))
. * 2/SHAT/XX1*PRO1**2*(SQ11-SQ10)
. * 4/SHAT/XX2*PRO2**2/ST12/STH/SCHI*(SQ21-SQ20)
ELSE
CTH = -1 + 2*Y(5)
CHI = PI*Y(6)
FAC = -(1-TAUH)*DLOG(TAU)*(1-XMU)
. * (XX1+XMU-XMU/(1-XX1))*2*2*PI
CCHI = DCOS(CHI)
STH = DSQRT(1-CTH**2)
ENDIF
I1I2 = SHAT/2
F1I1 = SHAT/4*XX1*(1-CTH)
F1I2 = SHAT/4*XX1*(1+CTH)
F2I1 = SHAT/4*XX2*(1-ST12*CCHI*STH-CT12*CTH)
F2I2 = SHAT/4*XX2*(1+ST12*CCHI*STH+CT12*CTH)
F3I1 = SHAT/2 - F1I1 - F2I1
F3I2 = SHAT/2 - F1I2 - F2I2
F1F2 = SHAT/2*(1-XX3+XMU)
F1F3 = SHAT/2*(1-XX2-XMU)
F2F3 = SHAT/2*(1-XX1-XMU)
Q1Q1 = -2*F1I1
Q2Q2 = -2*F2I2
Q1Q2 = F1F2 + I1I2 - F1I2 - F2I1
P1Q1 = F1I1/X1
P2Q2 = F2I2/X2
P1Q2 = (F2I1-I1I2)/X1
P2Q1 = (F1I2-I1I2)/X2
P1P2 = S/2
IF(NUNIT.EQ.1)THEN
SCALE = QQ/AM
Q1 = SCALE*DSQRT(-Q1Q1)
Q2 = SCALE*DSQRT(-Q2Q2)
ELSE
Q1 = QQ
Q2 = QQ
ENDIF
Z = Y(7)
IF((-Q1Q1).GT.Q0**2.AND.(-Q2Q2).GT.Q0**2)THEN
PRO1 = Q1Q1 - AMV**2
PRO2 = Q2Q2 - AMV**2
F11 = F1(I1,X1,Q1,1)
F12 = F1(I2,X2,Q2,ICOLL)
F21 = F2(I1,X1,Q1,1)
F22 = F2(I2,X2,Q2,ICOLL)
F31 = F3(I1,X1,Q1,1)
F32 = F3(I2,X2,Q2,ICOLL)
XF11 = DF1(I1,Z,X1,Q1,-Q1Q1,1)
XF12 = DF1(I2,Z,X2,Q2,-Q2Q2,ICOLL)
XF21 = DF2(I1,Z,X1,Q1,-Q1Q1,1)
XF22 = DF2(I2,Z,X2,Q2,-Q2Q2,ICOLL)
XF31 = DF3(I1,Z,X1,Q1,-Q1Q1,1)
XF32 = DF3(I2,Z,X2,Q2,-Q2Q2,ICOLL)
ZZ = ZMAT(XF11,F12,XF21,F22,XF31,F32)
. + ZMAT(F11,XF12,F21,XF22,F31,XF32)
RES = RES + DSQRT(2.D0)*GF**3*AMV**8/PRO1**2/PRO2**2*ZZ
. * F1I1*F2I2
. * FAC / 512/PI**4
ENDIF
ENDDO
DMAT = RES
RETURN
END
DOUBLE PRECISION FUNCTION F1(IV,X,Q,ICOLL)
C--STRUCTURE FUNCTION F_1
C--IV = 0: Z = +-1: W+-
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION PDF(-6:6)
COMMON/PARAM/S,AM,QQ,EPS,GF,AMW,AMZ,SW2,XMP,GEVPB
COMMON/QCUT/Q0
F1 = 0
IF(Q.GT.Q0)THEN
CALL STRUC(X,Q,PDF)
IF(IV.EQ.0)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 + PDF(I) + PDF(-I)
J = I+1
IF(J.LT.6) DLU = DLU + PDF(J) + PDF(-J)
ENDDO
F1 = (CU*DLU+CD*DLD)/X
ELSEIF(IV.EQ.1)THEN
DL = PDF(2*ICOLL) + PDF(4*ICOLL) + PDF(-1*ICOLL) + PDF(-3*ICOLL)
F1 = 4*DL/X
ELSEIF(IV.EQ.-1)THEN
DL = PDF(-2*ICOLL) + PDF(-4*ICOLL) + PDF(1*ICOLL) + PDF(3*ICOLL)
F1 = 4*DL/X
ENDIF
ENDIF
RETURN
END
DOUBLE PRECISION FUNCTION F2(IV,X,Q,ICOLL)
C--STRUCTURE FUNCTION F_2
C--IV = 0: Z = +-1: W+-
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION PDF(-6:6)
COMMON/PARAM/S,AM,QQ,EPS,GF,AMW,AMZ,SW2,XMP,GEVPB
COMMON/QCUT/Q0
F2 = 0
IF(Q.GT.Q0)THEN
CALL STRUC(X,Q,PDF)
IF(IV.EQ.0)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 + PDF(I) + PDF(-I)
J = I+1
IF(J.LT.6) DLU = DLU + PDF(J) + PDF(-J)
ENDDO
F2 = 2*(CU*DLU+CD*DLD)
ELSEIF(IV.EQ.1)THEN
DL = PDF(2*ICOLL) + PDF(4*ICOLL) + PDF(-1*ICOLL) + PDF(-3*ICOLL)
F2 = 8*DL
ELSEIF(IV.EQ.-1)THEN
DL = PDF(-2*ICOLL) + PDF(-4*ICOLL) + PDF(1*ICOLL) + PDF(3*ICOLL)
F2 = 8*DL
ENDIF
ENDIF
RETURN
END
DOUBLE PRECISION FUNCTION F3(IV,X,Q,ICOLL)
C--STRUCTURE FUNCTION F_3
C--IV = 0: Z = +-1: W+-
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION PDF(-6:6)
COMMON/PARAM/S,AM,QQ,EPS,GF,AMW,AMZ,SW2,XMP,GEVPB
COMMON/QCUT/Q0
F3 = 0
IF(Q.GT.Q0)THEN
CALL STRUC(X,Q,PDF)
IF(IV.EQ.0)THEN
CU = (1.D0-8.D0/3.D0*SW2)
CD = (1.D0-4.D0/3.D0*SW2)
DLU = 0
DLD = 0
DO I=1,5,2
DLD = DLD + PDF(I*ICOLL) - PDF(-I*ICOLL)
J = I+1
IF(J.LT.6) DLU = DLU + PDF(J*ICOLL) - PDF(-J*ICOLL)
ENDDO
F3 = 4*(CU*DLU+CD*DLD)/X
ELSEIF(IV.EQ.1)THEN
DL = PDF(2*ICOLL) + PDF(4*ICOLL) - PDF(-1*ICOLL) - PDF(-3*ICOLL)
F3 = 8*DL/X
ELSEIF(IV.EQ.-1)THEN
DL = -PDF(-2*ICOLL) - PDF(-4*ICOLL) + PDF(1*ICOLL) + PDF(3*ICOLL)
F3 = 8*DL/X
ENDIF
ENDIF
RETURN
END
DOUBLE PRECISION FUNCTION QF1(I,XX,X,Q,Q2,ICOLL)
C--STRUCTURE FUNCTION F_1: NLO QUARK CONTRIBUTION
C--I = 0: Z = +-1: W+-
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
COMMON/PARAM/S,AM,QQ,EPS,GF,AMW,AMZ,SW2,XMP,GEVPB
COMMON/QCUT/Q0
COMMON/MSDIS/NDIS
COMMON/SUBSTIT/NSUBST
QF1 = 0
IF(Q.GT.Q0)THEN
PI = 4*DATAN(1.D0)
IF(NSUBST.EQ.2)THEN
XL = EPS/(1-EPS)
AL = -DLOG((1-X)/X) + (-DLOG(XL)+DLOG((1-X)/X))*XX
Z = 1/(1+DEXP(-AL))
FAC = (-DLOG(XL)+DLOG((1-X)/X))*(1-Z)
ELSE
AL = -DLOG(X)*XX
Z = DEXP(-AL)
FAC = -DLOG(X)
ENDIF
Y = X/Z
SLG = DLOG(Q**2/Q2)
NALP = 2
ALP = ALPHAS(Q,NALP)/PI
IF(NDIS.NE.1)THEN
W1 = D0(I,X,Z,Q,ICOLL)*FAC*(-3/2.D0-2*SLG)
W2 = D1(I,X,Z,Q,ICOLL)*FAC
W3 = -(1+Z**2)/(1-Z)*DLOG(Z)+3+2*Z - 2*Z + (1+Z)*SLG
W4 = -9/2.D0-PI**2/3 - 3/2.D0*SLG
W3 = W3*F1(I,Y,Q,ICOLL)*FAC
W4 = W4*F1(I,X,Q,ICOLL)
QF1 = 2/3.D0*ALP*(W1+W2+W3+W4)
ELSE
WP1 = D0(I,X,Z,Q,ICOLL)*FAC*(-2*SLG)
WP3 = - 2*Z + (1+Z)*SLG
WP4 = - 3/2.D0*SLG
WP3 = WP3*F1(I,Y,Q,ICOLL)*FAC
WP4 = WP4*F1(I,X,Q,ICOLL)
QF1 = 2/3.D0*ALP*(WP1+WP3+WP4)
ENDIF
ENDIF
RETURN
END
DOUBLE PRECISION FUNCTION QF2(I,XX,X,Q,Q2,ICOLL)
C--STRUCTURE FUNCTION F_2: NLO QUARK CONTRIBUTION
C--I = 0: Z = +-1: W+-
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
COMMON/PARAM/S,AM,QQ,EPS,GF,AMW,AMZ,SW2,XMP,GEVPB
COMMON/QCUT/Q0
COMMON/MSDIS/NDIS
COMMON/SUBSTIT/NSUBST
QF2 = 0
IF(Q.GT.Q0)THEN
PI = 4*DATAN(1.D0)
IF(NSUBST.EQ.2)THEN
XL = EPS/(1-EPS)
AL = -DLOG((1-X)/X) + (-DLOG(XL)+DLOG((1-X)/X))*XX
Z = 1/(1+DEXP(-AL))
FAC = (-DLOG(XL)+DLOG((1-X)/X))*(1-Z)
ELSE
AL = -DLOG(X)*XX
Z = DEXP(-AL)
FAC = -DLOG(X)
ENDIF
Y = X/Z
SLG = DLOG(Q**2/Q2)
NALP = 2
ALP = ALPHAS(Q,NALP)/PI
IF(NDIS.NE.1)THEN
W1 = D0(I,X,Z,Q,ICOLL)*FAC*(-3/2.D0-2*SLG)
W2 = D1(I,X,Z,Q,ICOLL)*FAC
W3 = -(1+Z**2)/(1-Z)*DLOG(Z)+3+2*Z + (1+Z)*SLG
W4 = -9/2.D0-PI**2/3 - 3/2.D0*SLG
W3 = W3*F1(I,Y,Q,ICOLL)*FAC
W4 = W4*F1(I,X,Q,ICOLL)
QF2 = 2*X*2/3.D0*ALP*(W1+W2+W3+W4)
ELSE
WP1 = D0(I,X,Z,Q,ICOLL)*FAC*(-2*SLG)
WP3 = (1+Z)*SLG
WP4 = - 3/2.D0*SLG
WP3 = WP3*F1(I,Y,Q,ICOLL)*FAC
WP4 = WP4*F1(I,X,Q,ICOLL)
QF2 = 2*X*2/3.D0*ALP*(WP1+WP3+WP4)
ENDIF
ENDIF
RETURN
END
DOUBLE PRECISION FUNCTION QF3(I,XX,X,Q,Q2,ICOLL)
C--STRUCTURE FUNCTION F_3: NLO QUARK CONTRIBUTION
C--I = 0: Z = +-1: W+-
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
COMMON/PARAM/S,AM,QQ,EPS,GF,AMW,AMZ,SW2,XMP,GEVPB
COMMON/QCUT/Q0
COMMON/MSDIS/NDIS
COMMON/SUBSTIT/NSUBST
QF3 = 0
IF(Q.GT.Q0)THEN
PI = 4*DATAN(1.D0)
IF(NSUBST.EQ.2)THEN
XL = EPS/(1-EPS)
AL = -DLOG((1-X)/X) + (-DLOG(XL)+DLOG((1-X)/X))*XX
Z = 1/(1+DEXP(-AL))
FAC = (-DLOG(XL)+DLOG((1-X)/X))*(1-Z)
ELSE
AL = -DLOG(X)*XX
Z = DEXP(-AL)
FAC = -DLOG(X)
ENDIF
Y = X/Z
SLG = DLOG(Q**2/Q2)
NALP = 2
ALP = ALPHAS(Q,NALP)/PI
IF(NDIS.NE.1)THEN
W1 = D30(I,X,Z,Q,ICOLL)*FAC*(-3/2.D0-2*SLG)
W2 = D31(I,X,Z,Q,ICOLL)*FAC
W3 = -(1+Z**2)/(1-Z)*DLOG(Z)+3+2*Z - 1-Z + (1+Z)*SLG
W4 = -9/2.D0-PI**2/3 - 3/2.D0*SLG
W3 = W3*F3(I,Y,Q,ICOLL)*FAC
W4 = W4*F3(I,X,Q,ICOLL)
QF3 = 2/3.D0*ALP*(W1+W2+W3+W4)
ELSE
WP1 = D30(I,X,Z,Q,ICOLL)*FAC*(-2*SLG)
WP3 = (1+Z)*(-1+SLG)
WP4 = - 3/2.D0*SLG
WP3 = WP3*F3(I,Y,Q,ICOLL)*FAC
WP4 = WP4*F3(I,X,Q,ICOLL)
QF3 = 2/3.D0*ALP*(WP1+WP3+WP4)
ENDIF
ENDIF
RETURN
END
DOUBLE PRECISION FUNCTION GF1(I,XX,X,Q,Q2,ICOLL)
C--STRUCTURE FUNCTION F_1: NLO GLUON CONTRIBUTION
C--I = 0: Z = +-1: W+-
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION PDF(-6:6)
COMMON/PARAM/S,AM,QQ,EPS,GF,AMW,AMZ,SW2,XMP,GEVPB
COMMON/QCUT/Q0
COMMON/MSDIS/NDIS
COMMON/SUBSTIT/NSUBST
GF1 = 0
IF(Q.GT.Q0)THEN
IF(I.EQ.0)THEN
CU = 1.D0 + (1.D0-8.D0/3.D0*SW2)**2
CD = 1.D0 + (1.D0-4.D0/3.D0*SW2)**2
FACTOR = 4*CU + 6*CD
ELSE
FACTOR = 16
ENDIF
PI = 4*DATAN(1.D0)
IF(NSUBST.EQ.2)THEN
XL = EPS/(1-EPS)
AL = -DLOG((1-X)/X) + (-DLOG(XL)+DLOG((1-X)/X))*XX
Z = 1/(1+DEXP(-AL))
FAC = (-DLOG(XL)+DLOG((1-X)/X))*(1-Z)
ELSE
AL = -DLOG(X)*XX
Z = DEXP(-AL)
FAC = -DLOG(X)
ENDIF
Y = X/Z
SLG = DLOG(Q**2/Q2)
CALL STRUC(Y,Q,PDF)
NALP = 2
ALP = ALPHAS(Q,NALP)/PI
IF(NDIS.NE.1)THEN
WW = (Z**2+(1-Z)**2)*DLOG((1-Z)/Z)+8*Z*(1-Z)-1 - 4*Z*(1-Z)
. - (Z**2+(1-Z)**2)*SLG
GF1 = FACTOR/4.D0*ALP*WW*FAC*PDF(0)/Y
ELSE
WP = - 4*Z*(1-Z)
. - (Z**2+(1-Z)**2)*SLG
GF1 = FACTOR/4.D0*ALP*WP*FAC*PDF(0)/Y
ENDIF
ENDIF
RETURN
END
DOUBLE PRECISION FUNCTION GF2(I,XX,X,Q,Q2,ICOLL)
C--STRUCTURE FUNCTION F_2: NLO GLUON CONTRIBUTION
C--I = 0: Z = +-1: W+-
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION PDF(-6:6)
COMMON/PARAM/S,AM,QQ,EPS,GF,AMW,AMZ,SW2,XMP,GEVPB
COMMON/QCUT/Q0
COMMON/MSDIS/NDIS
COMMON/SUBSTIT/NSUBST
GF2 = 0
IF(Q.GT.Q0)THEN
IF(I.EQ.0)THEN
CU = 1.D0 + (1.D0-8.D0/3.D0*SW2)**2
CD = 1.D0 + (1.D0-4.D0/3.D0*SW2)**2
FACTOR = 4*CU + 6*CD
ELSE
FACTOR = 16
ENDIF
PI = 4*DATAN(1.D0)
IF(NSUBST.EQ.2)THEN
XL = EPS/(1-EPS)
AL = -DLOG((1-X)/X) + (-DLOG(XL)+DLOG((1-X)/X))*XX
Z = 1/(1+DEXP(-AL))
FAC = (-DLOG(XL)+DLOG((1-X)/X))*(1-Z)
ELSE
AL = -DLOG(X)*XX
Z = DEXP(-AL)
FAC = -DLOG(X)
ENDIF
Y = X/Z
SLG = DLOG(Q**2/Q2)
CALL STRUC(Y,Q,PDF)
NALP = 2
ALP = ALPHAS(Q,NALP)/PI
IF(NDIS.NE.1)THEN
WW = (Z**2+(1-Z)**2)*DLOG((1-Z)/Z)+8*Z*(1-Z)-1
. - (Z**2+(1-Z)**2)*SLG
GF2 = 2*X*FACTOR/4.D0*ALP*WW*FAC*PDF(0)/Y
ELSE
WP = - (Z**2+(1-Z)**2)*SLG
GF2 = 2*X*FACTOR/4.D0*ALP*WP*FAC*PDF(0)/Y
ENDIF
ENDIF
RETURN
END
DOUBLE PRECISION FUNCTION GF3(I,XX,X,Q,Q2,ICOLL)
C--STRUCTURE FUNCTION F_3: NLO GLUON CONTRIBUTION
C--I = 0: Z = +-1: W+-
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
GF3 = 0.D0
RETURN
END
DOUBLE PRECISION FUNCTION D0(I,TAU,X,Q,ICOLL)
C--STRUCTURE FUNCTION F_1: PLUS DISTRIBUTION 1/(1-Z)_+
C--I = 0: Z = +-1: W+-
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
FFZ = G0(I,TAU,X,Q,ICOLL)
FF1 = X*G0(I,TAU,1.D0,Q,ICOLL)
D0 = 1/(1-X)*(FFZ - FF1) + DLOG(1-TAU)/(1-TAU)*FF1
RETURN
END
DOUBLE PRECISION FUNCTION D1(I,TAU,X,Q,ICOLL)
C--STRUCTURE FUNCTION F_1: PLUS DISTRIBUTION (LOG(1-Z)/(1-Z))_+
C--I = 0: Z = +-1: W+-
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
FFZ = G1(I,TAU,X,Q,ICOLL)
FF1 = X*G1(I,TAU,1.D0,Q,ICOLL)
D1 = DLOG(1-X)/(1-X)*(FFZ - FF1)
. + DLOG(1-TAU)**2/2/(1-TAU)*FF1
RETURN
END
DOUBLE PRECISION FUNCTION D30(I,TAU,X,Q,ICOLL)
C--STRUCTURE FUNCTION F_3: PLUS DISTRIBUTION 1/(1-Z)_+
C--I = 0: Z = +-1: W+-
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
FFZ = G30(I,TAU,X,Q,ICOLL)
FF1 = X*G30(I,TAU,1.D0,Q,ICOLL)
D30 = 1/(1-X)*(FFZ - FF1) + DLOG(1-TAU)/(1-TAU)*FF1
RETURN
END
DOUBLE PRECISION FUNCTION D31(I,TAU,X,Q,ICOLL)
C--STRUCTURE FUNCTION F_3: PLUS DISTRIBUTION (LOG(1-Z)/(1-Z))_+
C--I = 0: Z = +-1: W+-
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
FFZ = G31(I,TAU,X,Q,ICOLL)
FF1 = X*G31(I,TAU,1.D0,Q,ICOLL)
D31 = DLOG(1-X)/(1-X)*(FFZ - FF1)
. + DLOG(1-TAU)**2/2/(1-TAU)*FF1
RETURN
END
DOUBLE PRECISION FUNCTION G1(I,TAU,Z,Q,ICOLL)
C--FUNCTION FOR D1(I,TAU,X,Q,ICOLL)
C--I = 0: Z = +-1: W+-
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
Y = TAU/Z
G1 = (1+Z**2)*F1(I,Y,Q,ICOLL)
RETURN
END
DOUBLE PRECISION FUNCTION G0(I,TAU,Z,Q,ICOLL)
C--FUNCTION FOR D0(I,TAU,X,Q,ICOLL)
C--I = 0: Z = +-1: W+-
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
Y = TAU/Z
G0 = F1(I,Y,Q,ICOLL)
RETURN
END
DOUBLE PRECISION FUNCTION G31(I,TAU,Z,Q,ICOLL)
C--FUNCTION FOR D31(I,TAU,X,Q,ICOLL)
C--I = 0: Z = +-1: W+-
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
Y = TAU/Z
G31 = (1+Z**2)*F3(I,Y,Q,ICOLL)
RETURN
END
DOUBLE PRECISION FUNCTION G30(I,TAU,Z,Q,ICOLL)
C--FUNCTION FOR D30(I,TAU,X,Q,ICOLL)
C--I = 0: Z = +-1: W+-
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
Y = TAU/Z
G30 = F3(I,Y,Q,ICOLL)
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,QQ,EPS,GF,AMW,AMZ,SW2,XMP,GEVPB
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,QQ,EPS,GF,AMW,AMZ,SW2,XMP,GEVPB
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,QQ,EPS,GF,AMW,AMZ,SW2,XMP,GEVPB
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,QQ,EPS,GF,AMW,AMZ,SW2,XMP,GEVPB
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
Reference in New Issue View Git Blame Copy Permalink