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-- Compile using e.g. g77 ggh_mstwpdf.f mstwpdf.f alphaS.f 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]. C-- Updated 14/07/2010: Avoid using Fortran 90 intrinsic "len_trim". PROGRAM GGH_MSTWPDF 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 alphaSorder,alphaSnfmax,ieigen,neigen,ias,iasSet,lentrim PARAMETER(neigen=20) ! number of eigenvectors DOUBLE PRECISION GetOnePDF,xg,ALPHAS, & distance,tolerance,mCharm,mBottom,alphaSQ0,alphaSMZ, & 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 CHARACTER prefix0*50,prefix*50,aslabel(0:4)*11,asSet*3 COMMON/mstwCommon/distance,tolerance, & mCharm,mBottom,alphaSQ0,alphaSMZ,alphaSorder,alphaSnfmax DATA aslabel/"best-fit aS","aS -1sigma ","aS -sigma/2", & "aS +sigma/2","aS +1sigma "/ C-- Prefix for the PDF grid files with the best-fit alphaS. C prefix0 = "Grids/mstw2008nlo" prefix0 = "Grids/mstw2008nnlo" 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. DO ias = 0, 21 IF (ias.EQ.0) THEN prefix = prefix0 ! best-fit alpha_S ELSE C iasSet = 110 + (ias-1) ! = alpha_S(M_Z)*1000 (NLO) iasSet = 107 + (ias-1) ! = alpha_S(M_Z)*1000 (NNLO) WRITE(asSet,'(I3.3)') iasSet ! convert integer to string prefix = prefix0(1:lentrim(prefix0))//"_asmz"//asSet END IF C-- Initialise alpha_S automatically using value stored in PDF grids. xg = GetOnePDF(prefix,0,0.1D0,10.D0,0) ! dummy call to initialise CALL INITALPHAS(alphaSorder,1.D0,1.D0,alphaSQ0, & mCharm,mBottom,1.D10) C-- Check calculated value of alpha_S(M_Z) matches stored value. WRITE(6,'(" alphaS(MZ) = ",F7.5," = ",F7.5)') & ALPHAS(91.1876D0),alphaSMZ C-- Get Higgs cross section for central PDF set and print out result. CALL SIGMAH(prefix,RTS,MH,xsec,0) 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 IF (ias.EQ.0) THEN prefix = prefix0 ELSE IF (ias.EQ.1) THEN C prefix = prefix0(1:lentrim(prefix0))//"_asmz-68cl" prefix = prefix0(1:lentrim(prefix0))//"_asmz-90cl" ELSE IF (ias.EQ.2) THEN C prefix = prefix0(1:lentrim(prefix0))//"_asmz-68clhalf" prefix = prefix0(1:lentrim(prefix0))//"_asmz-90clhalf" ELSE IF (ias.EQ.3) THEN C prefix = prefix0(1:lentrim(prefix0))//"_asmz+68clhalf" prefix = prefix0(1:lentrim(prefix0))//"_asmz+90clhalf" ELSE IF (ias.EQ.4) THEN C prefix = prefix0(1:lentrim(prefix0))//"_asmz+68cl" prefix = prefix0(1:lentrim(prefix0))//"_asmz+90cl" ELSE WRITE(6,*) "Error: ias = ",ias STOP END IF C-- Initialise alpha_S automatically using value stored in PDF grids. xg = GetOnePDF(prefix,0,0.1D0,10.D0,0) ! dummy call to initialise CALL INITALPHAS(alphaSorder,1.D0,1.D0,alphaSQ0, & mCharm,mBottom,1.D10) C-- Check calculated value of alpha_S(M_Z) matches stored value. WRITE(6,'(" alphaS(MZ) = ",F7.5," = ",F7.5)') & ALPHAS(91.1876D0),alphaSMZ C-- Get cross section for central set. CALL SIGMAH(prefix,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(prefix,RTS,MH,sigp,2*ieigen-1) ! "+" CALL SIGMAH(prefix,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(prefix,W,HMASS,TOTAL,ISET) IMPLICIT REAL*8(A-H,O-Z) REAL*8 DSIG(0:100) DATA GF,TMASS/1.16639D-5,175D0/ DATA PI/3.1415927D0/ DATA NPTS/100/ CHARACTER prefix*50,prefix1*55 prefix1 = prefix C IF (ISET.GT.0) prefix1 = prefix(1:lentrim(prefix))//'.68cl' IF (ISET.GT.0) prefix1 = prefix(1:lentrim(prefix))//'.90cl' 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*ALPHAS(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 G1 = GetOnePDF(prefix1,ISET,X1,HMASS,0) G2 = GetOnePDF(prefix1,ISET,X2,HMASS,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