! ! run this with ! ! ./test_qed_evol_lhapdf -Qlo 3.0 -Qhi 100 -pdf CT14qed_proton -nqcdloop-qed 0 -x 1e-3 >! a ! ! !! TESTS !! * Test against plain qcd evolution, with QED alpha turned to zero !! -> works OK, but watch out for the fact that the tau introduces !! an extra threshold, so that numerically there are differences !! at the evolution precision; reducing du by a factor of 10 !! significantly reduces these differences. See basic test output with !! ./test_qed_evol_lhapdf -no-qed -x 0.5 | grep moment !! ./test_qed_evol_lhapdf -alpha-qed 0 -x 0.5 | grep moment !! !! * test momentum sum rule, with both nqcdloop_qed values !! -> this seems to work fine, e.g. test with !! ./test_qed_evol_lhapdf -ymax 25 -nqcdloop-qed 0 |& grep moment !! ./test_qed_evol_lhapdf -ymax 25 -nqcdloop-qed 1 |& grep moment !! !! * compare to photon evolution from another source? !! -> done in with-lhapdf/test_qed_evol_lhapdf.f90 !! generally, answers are to within O(1%) of the MRST2004qed and !! CT14qed_proton results for evolution from 3 to 100 GeV. E.g. !! ./test_qed_evol_lhapdf -Qlo 3.0 -Qhi 100 -pdf CT14qed_proton -nqcdloop-qed 0 -x 1e-3 >! a !! program test_qed_evol use hoppet_v1 use qed_evolution; use qed_objects; use qed_coupling_module use sub_defs_io implicit none character(len=200) :: pdfname type(grid_def) :: grid type(qed_split_mat) :: qed_split type(dglap_holder) :: dh real(dp) :: quark_masses(4:6) type(running_coupling) :: coupling type(qed_coupling) :: coupling_qed real(dp) :: ymax, dy, x, alphas_mz, lha_qmin, scale_qg real(dp), parameter :: mz = 91.1876_dp real(dp), pointer :: pdf_in(:,:), pdf_out(:,:), pdf_lhapdf_out(:,:) real(dp), pointer :: photon_lumi(:,:) real(dp) :: Qlo, Qhi, moments(ncompmin:ncompmaxLeptons) real(dp) :: alphaspdf ! from LHAPDF integer :: compact_output_n = -1 integer :: nloop_qcd = 3, nqcdloop_qed, imem, iflv, iunit logical :: no_qed, use_lhapdf if (log_val_opt("-out")) then iunit = idev_open_opt("-out") else iunit = 6 end if compact_output_n = int_val_opt("-compact-output", compact_output_n) write(iunit,'(a)') "# "//trim(command_line()) ymax = dble_val_opt("-ymax",20.0_dp) dy = dble_val_opt("-dy",0.1_dp) call InitGridDefDefault(grid, dy, ymax) write(iunit,'(a,2f10.6)') "# ymax, dy ", ymax, dy ! allow option to just write out the x values of the grid ! and exit if (log_val_opt("-xvalues")) then write(iunit,"(a)") "# next lines contain number of x values, and actual x values" write(iunit,*) grid%ny+1 write(iunit,*) xValues(grid) stop end if Qlo = sqrt(two) ! the initial scale Qhi = 40.0_dp Qlo = sqrt(dble_val_opt("-Q2lo", Qlo**2)) Qhi = sqrt(dble_val_opt("-Q2hi", Qhi**2)) Qlo = dble_val_opt("-Qlo",Qlo) Qhi = dble_val_opt("-Qhi",Qhi) write(iunit,'(a,es14.6)') "# Qlo = ", Qlo write(iunit,'(a,es14.6)') "# Qhi = ", Qhi scale_qg = dble_val_opt("-scale-qg",one) use_lhapdf = log_val_opt("-pdf") if (use_lhapdf) then pdfname = string_val_opt("-pdf") imem = int_val_opt("-imem",0) call InitPDFsetByName(trim(pdfname)) call InitPDFm(1,imem) write(iunit,*) "# using LHAPDF ", trim(pdfname), ' imem = ', imem quark_masses(4:6) = (/lhapdf_qmass(4), lhapdf_qmass(5), lhapdf_qmass(6)/) alphas_mz = alphaspdf(mz) call getorderas(nloop_qcd); nloop_qcd = nloop_qcd + 1 ! LHAPDF convention is that LO=0 ! the following lines are problematic for some reason ! if (Qlo < lhapdf_qmin()) & ! &write(0,*) "WARNING: Qlo = ", Qlo, " < lhapdf Qmin = ", lhapdf_qmin() else write(iunit,*) "# using toy pdf" quark_masses(4:6) = (/1.414213563_dp, 4.5_dp, 175.0_dp/) alphas_mz = 0.118_dp end if write(iunit,*) '# quark masses (c,b,t) = ', quark_masses nloop_qcd = int_val_opt("-nloop-qcd",nloop_qcd) write(iunit,*) "# nloop_qcd = ", nloop_qcd nqcdloop_qed = int_val_opt('-nqcdloop-qed',0) write(iunit,*) "# nqcdloop_qed = ", nqcdloop_qed no_qed = log_val_opt("-no-qed",.false.) x = dble_val_opt("-x",1e-3_dp) call SetDefaultEvolutionDu(dy/3.0_dp) ! generally a good choice call InitQEDSplitMat(grid, qed_split) call InitDglapHolder(grid,dh,factscheme=factscheme_MSbar,& & nloop=nloop_qcd,nflo=3,nfhi=6) call InitRunningCoupling(coupling,alfas=alphas_mz,Q=mz,nloop=nloop_qcd,& & quark_masses = quark_masses) if (log_val_opt("-alpha-qed")) then call InitQEDCoupling(coupling_qed, quark_masses(4), quark_masses(5), quark_masses(6), & & dble_val_opt("-alpha-qed")) else call InitQEDCoupling(coupling_qed, quark_masses(4), quark_masses(5), quark_masses(6)) end if write(iunit,*) "# 1/QED coupling at MZ = ", one/Value(coupling_qed,MZ) write(iunit,*) "# QCD coupling at MZ = ", Value(coupling,MZ) write(iunit,*) "# Qlo, Qhi ", Qlo, Qhi call AllocPDFWithLeptons(grid, pdf_in) call AllocPDFWithLeptons(grid, pdf_out) call AllocPDFWithLeptons(grid, pdf_lhapdf_out) call AllocGridQuant(grid, photon_lumi, 1, 1) if (use_lhapdf) then call fill_from_lhapdf(pdf_in, Qlo) call fill_from_lhapdf(pdf_lhapdf_out, Qhi) else pdf_in = zero ! fill in the QCD part pdf_in(:,:ncompmax) = unpolarized_dummy_pdf(xValues(grid)) ! then add in a photon of some kind pdf_in(:,iflv_photon) = alpha_qed_scale_0 * (one - xValues(grid))**4 end if if (scale_qg /= one) then pdf_in(:,:6) = pdf_in(:,:6) * scale_qg end if ! allow user to override the photon from a suitably formatted file if (log_val_opt("-read-photon")) pdf_in(:,iflv_photon) = read_photon(grid) if (.not.CheckAllArgsUsed(0)) stop pdf_out = pdf_in if (no_qed) then call EvolvePDF(dh, pdf_out(:,ncompmin:ncompmax), coupling, Qlo, Qhi, nloop=nloop_qcd) else call QEDQCDEvolvePDF(dh, qed_split, pdf_out, coupling, coupling_qed,& & Qlo, Qhi, nloop_qcd, nqcdloop_qed) end if call write_moments(pdf_in, ' in') call write_moments(pdf_out,'out') if (use_lhapdf) call write_moments(pdf_lhapdf_out,'lha') call write_pdf(pdf_in,'in') call write_pdf(pdf_out,'out') if (use_lhapdf) call write_pdf(pdf_lhapdf_out,'lhapdf out') photon_lumi(:,1) = PartonLuminosity(grid, pdf_in(:,8),pdf_in(:,8)) call write_pdf(photon_lumi,'photon-lumi-in') photon_lumi(:,1) = PartonLuminosity(grid, pdf_out(:,8),pdf_out(:,8)) call write_pdf(photon_lumi,'photon-lumi-out') if (use_lhapdf) then photon_lumi(:,1) = PartonLuminosity(grid, pdf_lhapdf_out(:,8),pdf_lhapdf_out(:,8)) call write_pdf(photon_lumi,'photon-lumi-lhapdf-out') end if write(0,*) 'x = ', x, 'Q = ', Qhi write(0,'(a5,4a17)') 'iflv','lhapdf','hoppet' do iflv = -6, 8 write(0,'(i5,4f17.7)') iflv, pdf_lhapdf_out(:,iflv).atx.(x.with.grid),& & pdf_out (:,iflv).atx.(x.with.grid) end do contains function lhapdf_qmin() result(qmin) real(dp) :: qmin, q2min call getq2min(q2min) write(0,*) q2min qmin = sqrt(q2min) end function lhapdf_qmin function lhapdf_qmass(iflv) result(res) integer, intent(in) :: iflv real(dp) :: res call getqmass(iflv,res) end function lhapdf_qmass !---------------------------------------------------------------------- subroutine fill_from_lhapdf(this_pdf, Q) real(dp), intent(out) :: this_pdf(0:,-6:) real(dp), intent(in) :: Q !--------- real(dp) :: xv(0:grid%ny) integer :: i xv = xValues(grid) this_pdf = zero do i = 0, grid%ny call EvolvePDFPhoton(xv(i), Q, this_pdf(i,-6:6), this_pdf(i,8)) !write(iunit,*) xv(i), this_pdf(i,8) end do end subroutine fill_from_lhapdf !---------------------------------------------------------------------- ! this assumes a 6-column format where column 1 is x and column 6 is ! xf_{gamma/p}(x) function read_photon(grid) result(photon) type(grid_def), intent(in) :: grid !--------------------------------- integer :: iy, iread character(len=1000) :: line real(dp) :: photon(0:grid%ny), xVals(0:grid%ny), numbers(6) real(dp), parameter :: tolerance = 1e-7_dp iread = idev_open_opt("-read-photon",status='OLD') xVals = xValues(grid) write(0,*) 'reading photon from file',trim(string_val_opt("-read-photon")) do iy = 0, grid%ny ! skip comment lines do while (.true.) read(iread,'(a)') line if (line(1:1) /= '#') exit end do ! then get & check the numbers read(line,*) numbers if (abs(log(numbers(1)/xVals(iy))) > tolerance) then write(0,*) 'grid x value = ', xVals(iy) write(0,*) 'input x value = ', numbers(1) call wae_error('read_photon', 'x value is inconsistent with grid') end if photon(iy) = numbers(6) end do end function read_photon !---------------------------------------------------------------------- subroutine write_moments(pdf,label) real(dp), intent(in) :: pdf(0:,ncompmin:) character(len=*), intent(in) :: label !----------------- real(dp) :: moments(ncompmin:ubound(pdf,dim=2)) write(iunit,"(a)", advance="no"), "# total momentum "//trim(label)//" (& components)" moments = TruncatedMoment(grid, pdf, one) write(iunit,"(40f10.7)") sum(moments), moments write(iunit,"(a)", advance="no"), "# total number "//trim(label)//" (& components)" moments(1:6) = TruncatedMoment(grid, pdf(:,1:6)-pdf(:,-1:-6:-1), zero) write(iunit,"(40f11.7)") sum(moments(1:6)), moments(1:6) end subroutine write_moments !---------------------------------------------------------------------- subroutine write_pdf(pdf,label) real(dp), intent(in) :: pdf(0:,ncompmin:) character(len=*), intent(in) :: label !---------------------- real(dp) :: y, yVals(0:grid%ny), last_y real(dp) :: zeta, zeta_max, a_stretch=7.0_dp, pdf_val integer :: i,j write(iunit,"(a)"), "# pdf "//trim(label) !do i = 0, 100 ! y = 0.1_dp * i ! write(iunit,*) y, pdf.aty.(y.with.grid) !end do if (compact_output_n > 0) then ! value of compact_output_n sets number of points being written zeta_max = zeta_of_y(grid%ymax, a_stretch) do i = 0, compact_output_n zeta = (i*zeta_max)/compact_output_n y = y_of_zeta(zeta, a_stretch) write(iunit,'(f10.6)',advance="no") y do j = lbound(pdf,dim=2), ubound(pdf,dim=2) pdf_val = pdf(:,j).aty.(y.with.grid) if (abs(pdf_val) < 1e-100_dp) then write(iunit,'(a)',advance='no') ' 0 ' else write(iunit,'(es14.6)',advance='no') pdf_val end if end do write(iunit,*) end do else last_y = -one yVals = yValues(grid) do i = 0, grid%ny y = yVals(i) if (y < 1.000000001_dp * last_y) cycle write(iunit,*) y, pdf(i,:) last_y = y end do end if write(iunit,*) write(iunit,*) end subroutine write_pdf !====================================================================== !! The dummy PDF suggested by Vogt as the initial condition for the !! unpolarized evolution (as used in hep-ph/0511119). function unpolarized_dummy_pdf(xvals) result(pdf) real(dp), intent(in) :: xvals(:) real(dp) :: pdf(size(xvals),ncompmin:ncompmax) real(dp) :: uv(size(xvals)), dv(size(xvals)) real(dp) :: ubar(size(xvals)), dbar(size(xvals)) !--------------------- real(dp), parameter :: N_g = 1.7_dp, N_ls = 0.387975_dp real(dp), parameter :: N_uv=5.107200_dp, N_dv = 3.064320_dp real(dp), parameter :: N_db = half*N_ls pdf = zero ! clean method for labelling as PDF as being in the human representation ! (not actually needed after setting pdf=0 call LabelPdfAsHuman(pdf) !-- remember that these are all xvals*q(xvals) uv = N_uv * xvals**0.8_dp * (1-xvals)**3 dv = N_dv * xvals**0.8_dp * (1-xvals)**4 dbar = N_db * xvals**(-0.1_dp) * (1-xvals)**6 ubar = dbar * (1-xvals) ! labels iflv_g, etc., come from the hoppet_v1 module, inherited ! from the main program pdf(:, iflv_g) = N_g * xvals**(-0.1_dp) * (1-xvals)**5 pdf(:,-iflv_s) = 0.2_dp*(dbar + ubar) pdf(:, iflv_s) = pdf(:,-iflv_s) pdf(:, iflv_u) = uv + ubar pdf(:,-iflv_u) = ubar pdf(:, iflv_d) = dv + dbar pdf(:,-iflv_d) = dbar end function unpolarized_dummy_pdf !----------------------------------------------------------------- !! return zeta = ln 1/x + a*(1-x) (x = exp(-y)) function zeta_of_y(y, a) result(zeta) real(dp), intent(in) :: y, a real(dp) :: zeta zeta = y + a*(one - exp(-y)) end function zeta_of_y !----------------------------------------------------------------- !! return inversion of zeta = ln 1/x + a*(1-x) (x = exp(-y)) function y_of_zeta(zeta, a) result(y) real(dp), intent(in) :: zeta, a real(dp) :: y, x, diff_from_zero, deriv integer :: iter real(dp), parameter :: eps = 1e-12_dp integer, parameter :: maxiter = 100 ! starting condition (and soln if a = 0) y = zeta if (a /= zero) then do iter = 0, maxiter x = exp(-y); diff_from_zero = zeta - y - a*(one-x); ! we have found good solution if (abs(diff_from_zero) < eps) exit deriv = -one - a*x; y = y - diff_from_zero / deriv; end do end if if (iter > maxiter) write(0,*) "y_of_zeta reached maxiter" end function y_of_zeta end program test_qed_evol