C-- Modified version of ggh.f originally taken from C-- http://durpdg.dur.ac.uk/hepdata/mrs.html C-- Use MSTW 2008 PDFs with varying heavy-quark masses and flavours. C-- Compile using e.g. gfortran ggh_hq_mstwpdf.f mstwpdf.f alphaS.f C-- Switch between NLO/NNLO PDFs or 68%/90% C.L. PDF uncertainties C-- by changing the couple of lines indicated in the code below. C-- Comments to Graeme Watt . C-- Paper reference: arXiv:1007.2624 [hep-ph]. PROGRAM GGH_HQ_MSTWPDF c This program computes the leading-order Higgs cross section c at a hadron-hadron colllider from the contribution g g --> H. C-- Note that it is not meaningful to combine a LO partonic cross C-- section with NLO/NNLO PDFs as done here. The only purpose of C-- this example program is to demonstrate how to interface to PDFs. IMPLICIT NONE INTEGER alphaSorder,alphaSnfmax,ieigen,neigen,imc,imb,nf,ifit, & lentrim PARAMETER(neigen=20) ! number of eigenvectors DOUBLE PRECISION GetOnePDF,xg,ALPHAS, & distance,tolerance,mCharm,mBottom,alphaSQ0,alphaSMZ, & RTS,MH,sigp,sigm,xsec,sig0,errp,errm,errsym, & imcSet,imbSet PARAMETER(MH=120.D0) ! Higgs mass in GeV CHARACTER prefix0*50,prefix*50,mcSet*4,mbSet*4 COMMON/mstwCommon/distance,tolerance, & mCharm,mBottom,alphaSQ0,alphaSMZ,alphaSorder,alphaSnfmax C-- Prefix for the PDF grid files with the default heavy-quark masses. C-- Uncomment ONE of the two lines below. C prefix0 = "Grids/mstw2008nlo" ! NLO prefix0 = "Grids/mstw2008nnlo" ! NNLO C-- Centre-of-mass energy in GeV. C RTS = 1.96D3 ! Tevatron at sqrt(s) = 1.96 TeV RTS = 7.D3 ! LHC at sqrt(s) = 7 TeV C 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---------------------------------------------------------------------- C-- First consider PDF sets with different m_c values. WRITE(6,*) WRITE(6,*) "First consider PDF sets with different m_c values." C-- Loop over fits with free/fixed alpha_S(M_Z), 5- or 3- flavours. DO ifit = 1, 4 WRITE(6,*) IF (ifit.EQ.1) THEN WRITE(6,*) "Fit with free alpha_S(M_Z), 5-flavours." ELSE IF (ifit.EQ.2) THEN WRITE(6,*) "Fit with fixed alpha_S(M_Z), 5-flavours." ELSE IF (ifit.EQ.3) THEN WRITE(6,*) "Fit with free alpha_S(M_Z), 3-flavours." ELSE IF (ifit.EQ.4) THEN WRITE(6,*) "Fit with fixed alpha_S(M_Z), 3-flavours." END IF DO imc = 0, 14 ! loop over m_c values IF (imc.EQ.0) THEN ! default m_c value IF (ifit.le.2) THEN ! 5-flavour evolution prefix = prefix0 ELSE ! 3-flavour evolution prefix = prefix0(1:lentrim(prefix0))//"_nf3" END IF ELSE imcSet = 1.00D0 + imc*0.05D0 IF (imc.GT.7) imcSet = imcSet + 0.05D0 WRITE(mcSet,'(F4.2)') imcSet ! convert to string IF (ifit.EQ.1) THEN C-- Fits with varying m_c and free alpha_S(M_Z), 5-flavour evolution. prefix = prefix0(1:lentrim(prefix0)) & //"_mc"//mcSet//"GeV" ELSE IF (ifit.EQ.2) THEN C-- Fits with varying m_c and fixed alpha_S(M_Z), 5-flavour evolution. prefix = prefix0(1:lentrim(prefix0)) & //"_mc"//mcSet//"GeV_fixasmz" ELSE IF (ifit.EQ.3) THEN C-- Fits with varying m_c and free alpha_S(M_Z), 3-flavour evolution. prefix = prefix0(1:lentrim(prefix0)) & //"_mc"//mcSet//"GeV_nf3" ELSE IF (ifit.EQ.4) THEN C-- Fits with varying m_c and fixed alpha_S(M_Z), 3-flavour evolution. prefix = prefix0(1:lentrim(prefix0)) & //"_mc"//mcSet//"GeV_fixasmz_nf3" ELSE WRITE(6,*) "Error: ifit = ",ifit STOP END IF END IF WRITE(6,*) C-- Initialise alpha_S automatically using value stored in PDF grids. xg = GetOnePDF(prefix,0,0.1D0,10.D0,0) ! dummy call to initialise IF (alphaSnfmax.EQ.5) THEN ! maximum of 5 flavours CALL INITALPHAS(alphaSorder,1.D0,1.D0,alphaSQ0, & mCharm,mBottom,1.D10) ! i.e. set m_t large ELSE IF (alphaSnfmax.EQ.4) THEN ! maximum of 4 flavours CALL INITALPHAS(alphaSorder,1.D0,1.D0,alphaSQ0, & mCharm,0.9D10,1.D10) ! i.e. set m_b and m_t large ELSE IF (alphaSnfmax.EQ.3) THEN ! maximum of 3 flavours CALL INITALPHAS(alphaSorder,1.D0,1.D0,alphaSQ0, & 0.8D10,0.9D10,1.D10) ! i.e. set m_c, m_b and m_t large ELSE WRITE(6,*) "Error: alphaSnfmax = ",alphaSnfmax STOP END IF C-- Print out values of heavy-quark masses. WRITE(6, & '(" mCharm = ",F4.2," GeV, mBottom = ",F4.2," GeV")') & mCharm,mBottom 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-- Check calculated value of alpha_S(Q_0) matches stored value. WRITE(6,'(" alphaS(Q0) = ",F7.5," = ",F7.5)') & ALPHAS(1.D0),alphaSQ0 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 ! end loop over m_c values END DO ! end loop over ifit C---------------------------------------------------------------------- C-- Now consider PDF sets with different m_b values. WRITE(6,*) WRITE(6,*) "Now consider PDF sets with different m_b values." C-- Loop over fits with free alpha_S(M_Z), 5-or 4-flavour evolution. DO ifit = 1, 2 WRITE(6,*) IF (ifit.EQ.1) THEN WRITE(6,*) "Fit with free alpha_S(M_Z), 5-flavours." ELSE IF (ifit.EQ.2) THEN WRITE(6,*) "Fit with free alpha_S(M_Z), 4-flavours." END IF DO imb = 0, 6 ! loop over m_b values IF (imb.EQ.0) THEN IF (ifit.eq.1) THEN ! 5-flavour evolution prefix = prefix0 ELSE ! 4-flavour evolution prefix = prefix0(1:lentrim(prefix0))//"_nf4" END IF ELSE imbSet = 4.00D0 + (imb-1)*0.25D0 IF (imb.GT.3) imbSet = imbSet + 0.25D0 WRITE(mbSet,'(F4.2)') imbSet ! convert to string IF (ifit.EQ.1) THEN C-- Fits with varying m_b and free alpha_S(M_Z), 5-flavour evolution. prefix = prefix0(1:lentrim(prefix0)) & //"_mb"//mbSet//"GeV" ELSE IF (ifit.EQ.2) THEN C-- Fits with varying m_b and free alpha_S(M_Z), 4-flavour evolution. prefix = prefix0(1:lentrim(prefix0)) & //"_mb"//mbSet//"GeV_nf4" C-- Fits with varying m_b and fixed alpha_S(M_Z) are not provided C-- as the correlation between m_b and alpha_S(M_Z) is negligible. ELSE WRITE(6,*) "Error: ifit = ",ifit STOP END IF END IF WRITE(6,*) C-- Initialise alpha_S automatically using value stored in PDF grids. xg = GetOnePDF(prefix,0,0.1D0,10.D0,0) ! dummy call to initialise IF (alphaSnfmax.EQ.5) THEN ! maximum of 5 flavours CALL INITALPHAS(alphaSorder,1.D0,1.D0,alphaSQ0, & mCharm,mBottom,1.D10) ! i.e. set m_t large ELSE IF (alphaSnfmax.EQ.4) THEN ! maximum of 4 flavours CALL INITALPHAS(alphaSorder,1.D0,1.D0,alphaSQ0, & mCharm,0.9D10,1.D10) ! i.e. set m_b and m_t large ELSE IF (alphaSnfmax.EQ.3) THEN ! maximum of 3 flavours CALL INITALPHAS(alphaSorder,1.D0,1.D0,alphaSQ0, & 0.8D10,0.9D10,1.D10) ! i.e. set m_c, m_b and m_t large ELSE WRITE(6,*) "Error: alphaSnfmax = ",alphaSnfmax END IF C-- Print out values of heavy-quark masses. WRITE(6, & '(" mCharm = ",F4.2," GeV, mBottom = ",F4.2," GeV")') & mCharm,mBottom 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-- Check calculated value of alpha_S(Q_0) matches stored value. WRITE(6,'(" alphaS(Q0) = ",F7.5," = ",F7.5)') & ALPHAS(1.D0),alphaSQ0 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 ! end loop over m_b values END DO ! end loop over ifit C---------------------------------------------------------------------- C-- Now consider the eigenvector PDF sets. WRITE(6,*) WRITE(6,*) "Now consider the eigenvector PDF sets." C-- Calculate the PDF uncertainty on the cross section using both C-- the formula for asymmetric errors [eqs.(51,52) of arXiv:0901.0002] C-- and the formula for symmetric errors [eq.(50) of same paper]. DO nf = 5, 3, -1 ! loop over number of flavours (5, 4, 3). WRITE(6,*) WRITE(6,'(" Evolve with maximum of ",I1," flavours.")') nf IF (nf.EQ.5) THEN ! default 5-flavour PDFs prefix = prefix0 ELSE IF (nf.EQ.4) THEN ! evolve with 4 flavours maximum prefix = prefix0(1:lentrim(prefix0))//"_nf4" ELSE IF (nf.EQ.3) THEN ! evolve with 3 flavours prefix = prefix0(1:lentrim(prefix0))//"_nf3" ELSE WRITE(6,*) "Error: nf = ",nf 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 IF (alphaSnfmax.EQ.5) THEN ! maximum of 5 flavours CALL INITALPHAS(alphaSorder,1.D0,1.D0,alphaSQ0, & mCharm,mBottom,1.D10) ! i.e. set m_t large ELSE IF (alphaSnfmax.EQ.4) THEN ! maximum of 4 flavours CALL INITALPHAS(alphaSorder,1.D0,1.D0,alphaSQ0, & mCharm,0.9D10,1.D10) ! i.e. set m_b and m_t large ELSE IF (alphaSnfmax.EQ.3) THEN ! maximum of 3 flavours CALL INITALPHAS(alphaSorder,1.D0,1.D0,alphaSQ0, & 0.8D10,0.9D10,1.D10) ! i.e. set m_c, m_b and m_t large ELSE WRITE(6,*) "Error: alphaSnfmax = ",alphaSnfmax END IF C-- Print out values of heavy-quark masses. WRITE(6,'(" mCharm = ",F4.2," GeV, mBottom = ",F4.2," GeV")') & mCharm,mBottom 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-- Check calculated value of alpha_S(Q_0) matches stored value. WRITE(6,'(" alphaS(Q0) = ",F7.5," = ",F7.5)') & ALPHAS(1.D0),alphaSQ0 C-- Get cross section for central set. CALL SIGMAH(prefix,RTS,MH,sig0,0) C-- Now calculate PDF uncertainty on this value. errp = 0.d0 errm = 0.d0 errsym = 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 = errp + & (max(sigp-sig0,sigm-sig0,0.d0))**2 errm = errm + & (max(sig0-sigp,sig0-sigm,0.d0))**2 errsym = errsym + (sigp-sigm)**2 END DO errp = sqrt(errp) ! positive PDF error errm = sqrt(errm) ! negative PDF error errsym = 0.5D0*sqrt(errsym) ! symmetric PDF error PRINT 98,sig0,errp,errm,errsym END DO 98 FORMAT(' xsec (pb) = ',F8.5,' +',F8.5,' -',F8.5, & ' ( +/- ',F8.5,' )') STOP END C---------------------------------------------------------------------- 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-- Uncomment ONE of the two lines below to choose 68% or 90% C.L. 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