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