C-- Modified version of ggh.f from C-- http://durpdg.dur.ac.uk/hepdata/mrs.html C-- Use MSTW 2008 PDFs with varying alpha_S sets. C-- Link to LHAPDF V5.7.0 and make sure grid files are in C-- the share/lhapdf/PDFsets/ directory. C-- Switch between NLO/NNLO PDFs or 68%/90% C.L. uncertainties C-- by changing the handful of lines indicated in the code below. C-- Comments to Graeme Watt . C-- Paper reference: arXiv:0905.3531 [hep-ph]. PROGRAM GGH_LHAPDF c This program computes the leading-order Higgs cross section c at a hadron-hadron colllider from the contribution g g --> H. IMPLICIT NONE INTEGER ieigen,neigen,ias PARAMETER(neigen=20) ! number of eigenvectors DOUBLE PRECISION MZ,alphasPDF, & RTS,MH,sigp,sigm,xsec,sigmax,sigmin, & sig0(0:4),errp(0:4),errm(0:4),errsym(0:4) PARAMETER(MH=120.D0) ! Higgs mass in GeV PARAMETER(MZ=91.1876D0) ! Z mass in GeV CHARACTER aslabel(0:4)*11 DATA aslabel/"best-fit aS","aS -1sigma ","aS -sigma/2", & "aS +sigma/2","aS +1sigma "/ C-- Centre-of-mass energy in GeV. C RTS = 1.96D3 ! Tevatron at sqrt(s) = 1.96 TeV C RTS = 10.D3 ! LHC at sqrt(s) = 10 TeV RTS = 14.D3 ! LHC at sqrt(s) = 14 TeV WRITE(6,*) "Leading-order estimates for sigma(gg --> H) in pb." WRITE(6,*) "Centre-of-mass energy = ",RTS WRITE(6,*) "Higgs mass = ",MH C-- First consider the PDF sets with different fixed alpha_S values. C call InitPDFsetByName("MSTW2008nlo_asmzrange.LHgrid") call InitPDFsetByName("MSTW2008nnlo_asmzrange.LHgrid") DO ias = 0, 21 C-- Write out value of alpha_S(M_Z). WRITE(6,'(" alphaS(MZ) = ",F7.5)') alphasPDF(MZ) C-- Get Higgs cross section for central PDF set and print out result. CALL SIGMAH(RTS,MH,xsec,ias) WRITE(6,'(" xsec = ",F8.5," pb")') xsec END DO C-- Now consider the eigenvector PDF sets. C-- Loop over the central set and the 4 sets with different alpha_S. C-- Store the results in an array. Use asymmetric PDF uncertainties. C-- See Eqs.(1,7,8) of arXiv:0905.3531 for relevant formulae. DO ias = 0, 4 C-- To use NLO PDFs and alpha_S, just replace "nnlo" with "nlo". C-- To use 68% C.L. uncertainties, just replace "90cl" with "68cl". IF (ias.EQ.0) THEN C call InitPDFsetByName( C & "MSTW2008nnlo68cl.LHgrid") ! best-fit alpha_S call InitPDFsetByName( & "MSTW2008nnlo90cl.LHgrid") ! best-fit alpha_S ELSE IF (ias.EQ.1) THEN C call InitPDFsetByName( C & "MSTW2008nnlo68cl_asmz-68cl.LHgrid") call InitPDFsetByName( & "MSTW2008nnlo90cl_asmz-90cl.LHgrid") ELSE IF (ias.EQ.2) THEN C call InitPDFsetByName( C & "MSTW2008nnlo68cl_asmz-68clhalf.LHgrid") call InitPDFsetByName( & "MSTW2008nnlo90cl_asmz-90clhalf.LHgrid") ELSE IF (ias.EQ.3) THEN C call InitPDFsetByName( C & "MSTW2008nnlo68cl_asmz+68clhalf.LHgrid") call InitPDFsetByName( & "MSTW2008nnlo90cl_asmz+90clhalf.LHgrid") ELSE IF (ias.EQ.4) THEN C call InitPDFsetByName( C & "MSTW2008nnlo68cl_asmz+68cl.LHgrid") call InitPDFsetByName( & "MSTW2008nnlo90cl_asmz+90cl.LHgrid") ELSE WRITE(6,*) "Error: ias = ",ias STOP END IF C-- Write out value of alpha_S(M_Z). WRITE(6,'(" alphaS(MZ) = ",F7.5)') alphasPDF(MZ) C-- Get cross section for central set. CALL SIGMAH(RTS,MH,sig0(ias),0) C-- Now calculate PDF uncertainty on this value. errp(ias) = 0.d0 errm(ias) = 0.d0 errsym(ias) = 0.d0 DO ieigen = 1, neigen ! loop over eigenvector sets CALL SIGMAH(RTS,MH,sigp,2*ieigen-1) ! "+" CALL SIGMAH(RTS,MH,sigm,2*ieigen) ! "-" errp(ias) = errp(ias) + & (max(sigp-sig0(ias),sigm-sig0(ias),0.d0))**2 errm(ias) = errm(ias) + & (max(sig0(ias)-sigp,sig0(ias)-sigm,0.d0))**2 errsym(ias) = errsym(ias) + (sigp-sigm)**2 END DO errp(ias) = sqrt(errp(ias)) ! positive PDF error errm(ias) = sqrt(errm(ias)) ! negative PDF error errsym(ias) = 0.5D0*sqrt(errsym(ias)) ! symmetric PDF error END DO C-- Calculate combined "PDF+alphaS" uncertainty. C-- See Eqs.(9,10) of arXiv:0905.3531 for relevant formulae. WRITE(6,*) "Leading-order estimates for sigma(gg --> H) in pb." WRITE(6,*) "Centre-of-mass energy = ",RTS WRITE(6,*) "Higgs mass = ",MH sigmax = sig0(0) sigmin = sig0(0) DO ias = 0, 4 PRINT 98,aslabel(ias),sig0(ias),errp(ias),errm(ias),errsym(ias) IF (sig0(ias)+errp(ias).GT.sigmax) sigmax = sig0(ias)+errp(ias) IF (sig0(ias)-errm(ias).LT.sigmin) sigmin = sig0(ias)-errm(ias) END DO PRINT 99,sig0(0),sigmax-sig0(0),sig0(0)-sigmin 98 FORMAT(1X,A11,', xsec = ',F8.5, & ' +',F8.5,' -',F8.5,' ( +/- ',F8.5,' )') 99 FORMAT(1X,'tot. spread',', xsec = ',F8.5, & ' +',F8.5,' -',F8.5) STOP END SUBROUTINE SIGMAH(W,HMASS,TOTAL,ISET) IMPLICIT REAL*8(A-H,O-Z) REAL*8 DSIG(0:100),XPDF(-6:6) DATA GF,TMASS/1.16639D-5,175D0/ DATA PI/3.1415927D0/ DATA NPTS/100/ Call InitPDF(ISET) ! LHAPDF X=(TMASS/HMASS)**2 C-------------------------------------------------- IF(X.GE.0.25D0) THEN FR=-2.D0*(DASIN(0.5D0/DSQRT(X)))**2 FI=0.D0 ELSEIF(X.LT.0.25D0) THEN ROOT=DSQRT(0.25D0-X) ETAP=0.5D0+ROOT ETAM=0.5D0-ROOT RLOG=DLOG(ETAP/ETAM) FR=0.5D0*(RLOG**2-PI**2) FI=PI*RLOG ENDIF AR=2.D0*X+X*(4.D0*X-1.D0)*FR AI= X*(4.D0*X-1.D0)*FI ABISQ=9.D0*(AR**2+AI**2) C-------------------------------------------------- CONS=GF/(288D0*DSQRT(2D0)*PI) RNORM=389379.66D3*CONS*alphasPDF(HMASS)**2*ABISQ RTAU=HMASS/W TAU=RTAU**2 YMAX=DLOG(W/HMASS) YMIN=0D0 YRANGE=YMAX-YMIN DY=YRANGE/NPTS DO 99 IY=0,NPTS Y=YMIN+IY*DY X1=RTAU*DEXP(Y) X2=TAU/X1 HGG=0D0 IF(DMAX1(X1,X2).GE.1D0) GOTO 99 CALL EvolvePDF(X1,HMASS,XPDF) G1 = XPDF(0) CALL EvolvePDF(X2,HMASS,XPDF) G2 = XPDF(0) HGG=G1*G2 99 DSIG(IY)=HGG C SIMPSONS RULE FOR THE TOTAL CROSS SECTION NN=NPTS/2 TOTAL=DSIG(0)-DSIG(NPTS) DO J=1,NN TOTAL=TOTAL+4D0*DSIG(2*J-1)+2D0*DSIG(2*J) ENDDO TOTAL=TOTAL*DY*RNORM/3D0*2D0 RETURN END