diff --git a/CHANGES b/CHANGES new file mode 100644 index 0000000..459d045 --- /dev/null +++ b/CHANGES @@ -0,0 +1,16 @@ +This is a list of changes to pphqq.f since June 2000: + +June 23 2015: version 1.10: + - Link to LHAPDF implemented + - input alpha_s(M_Z) instead of Lambda_n + +September 12 2007: - Link to MRST2001 PDFs corrected due to errorneous + interchange of u and d. (thanks to M. Ciccolini) + +April 12 2005: - Link to CTEQ6 and MRST2001 PDFs + (NGROUP = -1: CTEQ6 -2: MRST2001) + +June 2 2003: - Reading in the mu parameter which was missing! + +June 23 2000: - Insertion of missing scale variations for scales in units + of Q (thanks to S. Dawson) diff --git a/makefile b/makefile new file mode 100644 index 0000000..15a4e5c --- /dev/null +++ b/makefile @@ -0,0 +1,30 @@ +# Makro-Definitions: + +OBJS = pphqq.o Cteq61Pdf.o + +LIBS = /afs/psi.ch/user/s/spira/data_nobackup/LHAPDF/lib/libLHAPDF.so + +#FFLAGS= -fno-emulate-complex -fno-automatic -ffixed-line-length-none -ffast-math -march=pentiumpro -malign-double -Wall -fno-silent + +#LIBS = $(CERN)/libpdflib804.a $(CERN)/libpawlib.a $(CERN)/libgrafX11.a $(CERN)/libpacklib.a $(CERN)/libmathlib.a $(CERN)/libkernlib.a + +#LIBS = `cernlib -v pro pdflib804,pawlib,graflib/X11,packlib,mathlib,kernlib` + +#FC=f77 +#FC=g77 + +#FFLAGS= -pc 64 -g77libs + +#FC=pgf77 + +FC=gfortran + +# Commands: +.f.o: + $(FC) -c $(FFLAGS) $*.f + +pphqq: $(OBJS) + $(FC) $(FFLAGS) $(OBJS) $(LIBS) -o run + +clean: + rm -f $(OBJS) diff --git a/pphqq.f b/pphqq.f new file mode 100644 index 0000000..eb9bf79 --- /dev/null +++ b/pphqq.f @@ -0,0 +1,2267 @@ +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 diff --git a/pphqq.in b/pphqq.in new file mode 100644 index 0000000..09a54ed --- /dev/null +++ b/pphqq.in @@ -0,0 +1,44 @@ +MSSM = 0 +TGBET = 3.D0 +HIGGS = 1 +MA1 = 125.D0 +MA2 = 200.D0 +NA = 1 +PP/PPBAR = 0 +ENERGY = 14000.D0 + +MH=0/Q=1 = 1 +SCALE1 = 1.0D0 +SCALE2 = 1.0D0 +NSCALE = 1 +EPSILON = 1.D-8 + +NGROUP = 0 +NSET = 0 +MS0/DIS1 = 0 +PDFPATH = Grids +PDFNAME = MSTW2008nlo68cl.LHgrid + +IPOINT1 = 10000 +IPOINT = 100000 +ITERAT = 5 +NPRN = 10 +SUBST12Z = 1 + +MC = 1.40D0 +MB = 4.75D0 +MT = 173.2D0 +ALS(MZ) = 0.12018D0 +LOOP = 2 +MW = 80.41D0 +MZ = 91.187D0 +SW2 = 0.2315D0 +GF = 1.16639D-5 +MPROTON = 0.93827231D0 +Q0 = 2.D0 + +MSQ = 1000.D0 +MUR = 1000.D0 +AT = 0.D0 +AB = 0.D0 +MU = 0.D0 diff --git a/struc.tar.gz b/struc.tar.gz new file mode 100644 index 0000000..eaa2c5c Binary files /dev/null and b/struc.tar.gz differ