From dd9a2777f709196b11c4daf93ca96085f8d8d531 Mon Sep 17 00:00:00 2001 From: Tobias Neumann Date: Thu, 29 Feb 2024 18:44:13 +0000 Subject: [PATCH] Initial merge with MCFM-Z to include Z @ N3LO; compiles but needs checks --- CMakeLists.txt | 4 +- README.md | 2 +- src/CuTe-MCFM/mod_Resummation.f90 | 38 ++- src/CuTe-MCFM/mod_Resummation_params.f90 | 3 +- src/Integrate/mcfm_vegas_adaptive.f | 7 - src/Mods/lhapdf_fortran.f90 | 66 ++++-- src/Mods/mod_Integration.f90 | 42 ++++ src/Mods/mod_SCET.f90 | 3 +- src/Need/mcfm_init.f | 45 ++-- src/Need/mcfmmain.f | 2 - src/Need/parseinput.f90 | 54 ++++- src/Need/setrunname.f | 13 +- src/Need/writereference.f | 8 + src/Parton/lhapdf_interface.cpp | 19 +- src/Phase/gen_njets.f | 4 +- src/Procdep/fragint.f | 6 +- src/Procdep/lowint.f | 8 +- src/Procdep/realint.f | 8 +- src/Procdep/resint.f90 | 280 +++++++++++++++++------ src/Procdep/scetint.f | 7 +- src/Procdep/virtint.f | 6 +- src/QT/qtint.f | 7 +- src/SCET/setupscet.f | 6 + src/SCET1j/mod_nnlo_w1jet.f90 | 3 + src/SCET1j/mod_nnlo_z1jet.f90 | 4 +- src/Singletop/boostx.f | 9 +- src/User/CMakeLists.txt | 84 +++---- src/User/mdata.f | 18 +- src/W2jet/qqb_w2jet_v.f | 2 + src/W2jet/qqb_wm2jet_g.f | 2 + src/W2jet/qqb_wp2jet_g.f | 2 + src/WpWp2j/CMakeLists.txt | 2 +- src/Z/AAAREADME.txt | 1 + src/Z2jet/qqb_z2jet_g.f | 2 + src/ptveto/ptint.f | 6 +- 35 files changed, 520 insertions(+), 253 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 3000af0..6318886 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -273,8 +273,8 @@ add_dependencies(objlib handyg_lib_static) if(use_internal_lhapdf) ExternalProject_Add(lhapdf - URL ${CMAKE_SOURCE_DIR}/lib/LHAPDF-6.5.1.tar.gz - SOURCE_DIR ${CMAKE_BINARY_DIR}/LHAPDF-6.5.1 + URL ${CMAKE_SOURCE_DIR}/lib/LHAPDF-6.5.4.tar.gz + SOURCE_DIR ${CMAKE_BINARY_DIR}/LHAPDF-6.5.4 PATCH_COMMAND mkdir -p ${CMAKE_BINARY_DIR}/local BUILD_IN_SOURCE ON CONFIGURE_COMMAND ./configure --with-pic --disable-lhaglue --disable-lhaglue-cxx --disable-python --prefix=${CMAKE_BINARY_DIR}/local CC=${CMAKE_C_COMPILER} CXX=${CMAKE_CXX_COMPILER} diff --git a/README.md b/README.md index 9f11b9a..4350afa 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ -# MCFM v10.3, CuTe-MCFM v1.2, January 2023 +# MCFM v10.4, CuTe-MCFM v1.3, March 2024 On the web: https://mcfm.fnal.gov/ diff --git a/src/CuTe-MCFM/mod_Resummation.f90 b/src/CuTe-MCFM/mod_Resummation.f90 index f121476..0897380 100644 --- a/src/CuTe-MCFM/mod_Resummation.f90 +++ b/src/CuTe-MCFM/mod_Resummation.f90 @@ -166,7 +166,7 @@ end function subroutine qtsubtraction(p, q2, x1, x2, f1, f2, hard, muH, & - order, muMult, resarray) + order, resarray, muMult) use LHAPDF use Beamfunctions3L use SCET @@ -193,13 +193,22 @@ real(dp), allocatable :: logMusqDivCutsq(:), taucuts(:) real(dp), allocatable :: outpiece(:) - allocate(resarray(size(tcutarray)+1)) - allocate(outpiece(size(tcutarray)+1)) - allocate(logMusqDivCutsq(size(tcutarray)+1)) - allocate(taucuts(size(tcutarray)+1)) + if (present(muMult)) then + allocate(resarray(1)) + allocate(outpiece(1)) + allocate(logMusqDivCutsq(1)) + allocate(taucuts(1)) + + taucuts(1) = taucut + else + allocate(resarray(size(tcutarray)+1)) + allocate(outpiece(size(tcutarray)+1)) + allocate(logMusqDivCutsq(size(tcutarray)+1)) + allocate(taucuts(size(tcutarray)+1)) - taucuts(1) = taucut - taucuts(2:size(tcutarray)+1) = tcutarray(:) + taucuts(1) = taucut + taucuts(2:size(tcutarray)+1) = tcutarray(:) + endif if (dynamictau) then do j=1,size(taucuts) @@ -672,12 +681,17 @@ initQuark = .true. endif - call find_root(qstar, 1._dp, 1.0e-3_dp, 100, root, ierr) - if (ierr .eqv. .false.) then - write (*,*) "WARNING: could not determine qstar, using default of 1.88" - root = 1.88_dp + + if (.not. present(resexp)) then + call find_root(qstar, 1._dp, 1.0e-3_dp, 100, root, ierr) + if (ierr .eqv. .false.) then + write (*,*) "WARNING: could not determine qstar, using default of 1.88" + root = 1.88_dp + endif + muStar = root + else + muStar = 1.88_dp endif - muStar = root ! Fourier integral is numerically unstable below mu~2GeV if (present(muMult)) then diff --git a/src/CuTe-MCFM/mod_Resummation_params.f90 b/src/CuTe-MCFM/mod_Resummation_params.f90 index b427118..409059e 100644 --- a/src/CuTe-MCFM/mod_Resummation_params.f90 +++ b/src/CuTe-MCFM/mod_Resummation_params.f90 @@ -15,7 +15,8 @@ module qtResummation_params real(dp), save, public :: qtmaxResexp = 80d0 real(dp), save, public :: qtminResexp = 1.0d0 - real(dp), save, public :: transitionSwitch = 0.4d0 + real(dp), save, public :: transitionSwitch = 0.4_dp + real(dp), save, allocatable, public :: transitionSwitches(:) logical, save, public :: enable_fixed_y = .false. real(dp), save, public :: fixed_y = 0.0d0 diff --git a/src/Integrate/mcfm_vegas_adaptive.f b/src/Integrate/mcfm_vegas_adaptive.f index ed35c7b..1ecb68d 100644 --- a/src/Integrate/mcfm_vegas_adaptive.f +++ b/src/Integrate/mcfm_vegas_adaptive.f @@ -115,9 +115,6 @@ c integer :: hist_readsnapshot, hist_writesnapshot logical :: bin common/bin/bin - logical :: dryrun - common/dryrun/dryrun - logical :: readin character(len=1024):: runname @@ -220,10 +217,6 @@ c integer :: hist_readsnapshot, hist_writesnapshot ! no binning for warumup stage bin = .false. - ! we don't need or want this here - dryrun = .false. - nprn = 0 - do part=1,maxParts; do ips=1,maxIps iterationStorage(part,ips)%used = .false. associate (info => iterationStorage(part,ips)%vinfo) diff --git a/src/Mods/lhapdf_fortran.f90 b/src/Mods/lhapdf_fortran.f90 index e867831..9d8e03e 100644 --- a/src/Mods/lhapdf_fortran.f90 +++ b/src/Mods/lhapdf_fortran.f90 @@ -26,11 +26,6 @@ module LHAPDF real(c_double) :: errplus real(c_double) :: errminus real(c_double) :: errsymm - real(c_double) :: scale - real(c_double) :: errplus_pdf - real(c_double) :: errminus_pdf - real(c_double) :: errsymm_pdf - real(c_double) :: err_par end type public :: pdf @@ -462,38 +457,51 @@ module LHAPDF integer :: accum integer, allocatable :: memberCounts(:) + integer, allocatable :: idx_j(:), idx_k(:) allocate(memberCounts(size(setnames))) do j=1,size(setnames) memberCounts(j) = lhapdf_number(trim(setnames(j))) enddo nmax = sum(memberCounts(:)) + allocate(idx_j(0:nmax-1)) + allocate(idx_k(0:nmax-1)) + if (.not. allocated(respdfs)) then allocate(respdfs(16,0:nmax-1)) accum = 0 do j=1,size(setnames) do k=0, memberCounts(j) - 1 - respdfs(1,accum) = lhapdf_loadMember(trim(setnames(j))//'_B10'//c_null_char, k) - respdfs(2,accum) = lhapdf_loadMember(trim(setnames(j))//'_B11'//c_null_char, k) - respdfs(3,accum) = lhapdf_loadMember(trim(setnames(j))//'_B20'//c_null_char, k) - respdfs(4,accum) = lhapdf_loadMember(trim(setnames(j))//'_B21'//c_null_char, k) - respdfs(5,accum) = lhapdf_loadMember(trim(setnames(j))//'_B22'//c_null_char, k) - respdfs(6,accum) = lhapdf_loadMember(trim(setnames(j))//'_B30'//c_null_char, k) - respdfs(7,accum) = lhapdf_loadMember(trim(setnames(j))//'_B31'//c_null_char, k) - respdfs(8,accum) = lhapdf_loadMember(trim(setnames(j))//'_B32'//c_null_char, k) - respdfs(9,accum) = lhapdf_loadMember(trim(setnames(j))//'_B33'//c_null_char, k) - respdfs(10,accum) = lhapdf_loadMember(trim(setnames(j))//'_G10'//c_null_char, k) - respdfs(11,accum) = lhapdf_loadMember(trim(setnames(j))//'_B44'//c_null_char, k) - respdfs(12,accum) = lhapdf_loadMember(trim(setnames(j))//'_B43'//c_null_char, k) - respdfs(13,accum) = lhapdf_loadMember(trim(setnames(j))//'_B55'//c_null_char, k) - - respdfs(14,accum) = lhapdf_loadMember(trim(setnames(j))//'_B42'//c_null_char, k) - respdfs(15,accum) = lhapdf_loadMember(trim(setnames(j))//'_B54'//c_null_char, k) - respdfs(16,accum) = lhapdf_loadMember(trim(setnames(j))//'_B66'//c_null_char, k) + idx_j(accum) = j + idx_k(accum) = k accum = accum + 1 enddo enddo +!$omp parallel do + do accum=0,nmax-1 + !if (rank == 0) then + ! write (*,*) "Allocating resummation member ", accum, idx_j(accum), idx_k(accum) + !endif + respdfs(1,accum) = lhapdf_loadMember(trim(setnames(idx_j(accum)))//'_B10'//c_null_char, idx_k(accum)) + respdfs(2,accum) = lhapdf_loadMember(trim(setnames(idx_j(accum)))//'_B11'//c_null_char, idx_k(accum)) + respdfs(3,accum) = lhapdf_loadMember(trim(setnames(idx_j(accum)))//'_B20'//c_null_char, idx_k(accum)) + respdfs(4,accum) = lhapdf_loadMember(trim(setnames(idx_j(accum)))//'_B21'//c_null_char, idx_k(accum)) + respdfs(5,accum) = lhapdf_loadMember(trim(setnames(idx_j(accum)))//'_B22'//c_null_char, idx_k(accum)) + respdfs(6,accum) = lhapdf_loadMember(trim(setnames(idx_j(accum)))//'_B30'//c_null_char, idx_k(accum)) + respdfs(7,accum) = lhapdf_loadMember(trim(setnames(idx_j(accum)))//'_B31'//c_null_char, idx_k(accum)) + respdfs(8,accum) = lhapdf_loadMember(trim(setnames(idx_j(accum)))//'_B32'//c_null_char, idx_k(accum)) + respdfs(9,accum) = lhapdf_loadMember(trim(setnames(idx_j(accum)))//'_B33'//c_null_char, idx_k(accum)) + respdfs(10,accum) = lhapdf_loadMember(trim(setnames(idx_j(accum)))//'_G10'//c_null_char, idx_k(accum)) + respdfs(11,accum) = lhapdf_loadMember(trim(setnames(idx_j(accum)))//'_B44'//c_null_char, idx_k(accum)) + respdfs(12,accum) = lhapdf_loadMember(trim(setnames(idx_j(accum)))//'_B43'//c_null_char, idx_k(accum)) + respdfs(13,accum) = lhapdf_loadMember(trim(setnames(idx_j(accum)))//'_B55'//c_null_char, idx_k(accum)) + + respdfs(14,accum) = lhapdf_loadMember(trim(setnames(idx_j(accum)))//'_B42'//c_null_char, idx_k(accum)) + respdfs(15,accum) = lhapdf_loadMember(trim(setnames(idx_j(accum)))//'_B54'//c_null_char, idx_k(accum)) + respdfs(16,accum) = lhapdf_loadMember(trim(setnames(idx_j(accum)))//'_B66'//c_null_char, idx_k(accum)) + enddo +!$omp end parallel do endif end subroutine initallResummation @@ -507,6 +515,7 @@ module LHAPDF integer :: j,k integer :: accum integer, allocatable :: memberCounts(:) + integer, allocatable :: idx_j(:), idx_k(:) allocate(memberCounts(size(setnames))) do j=1,size(setnames) @@ -514,16 +523,27 @@ module LHAPDF enddo nmax = sum(memberCounts(:)) + allocate(idx_j(0:nmax-1)) + allocate(idx_k(0:nmax-1)) if (.not. allocated(pdfs)) then allocate(pdfs(0:nmax-1)) accum = 0 do j=1,size(setnames) do k=0, memberCounts(j) - 1 - pdfs(accum) = lhapdf_loadMember(trim(setnames(j))//c_null_char, k) + idx_j(accum) = j + idx_k(accum) = k accum = accum + 1 enddo enddo +!$omp parallel do + do accum=0,nmax-1 + !if (rank == 0) then + ! write (*,*) "Allocating member ", accum, idx_j(accum), idx_k(accum) + !endif + pdfs(accum) = lhapdf_loadMember(trim(setnames(idx_j(accum)))//c_null_char, idx_k(accum)) + enddo +!$omp end parallel do endif end subroutine diff --git a/src/Mods/mod_Integration.f90 b/src/Mods/mod_Integration.f90 index e374a02..83b5a05 100644 --- a/src/Mods/mod_Integration.f90 +++ b/src/Mods/mod_Integration.f90 @@ -60,4 +60,46 @@ module Integration real(dp), save :: warmupChisqGoal + public :: mapPartStringToNumber + + contains + + function mapPartStringToNumber(str) result(ret) + implicit none + character(*), intent(in) :: str + integer :: ret + + select case (str) + case ("lo") + ret = lord + case ("nloVirt") + ret = nloVirt + case ("nloReal") + ret = nloReal + case ("nnloBelow") + ret = nnloBelow + case ("qtnloBelow") + ret = qtnloBelow + case ("qtnnloBelow") + ret = qtnnloBelow + case ("nnloVirtAbove") + ret = nnloVirtAbove + case ("nnloRealAbove") + ret = nnloRealAbove + case ("snloBelow") + ret = snloBelow + case ("snloAbove") + ret = snloAbove + case ("resexp") + ret = nloResexp + case ("resummed") + ret = nloResummed + case ("nloResAbove") + ret = nloResAbove + case default + error stop "Unknown part in mapPartStringToNumber: "//str + end select + + end function + end module diff --git a/src/Mods/mod_SCET.f90 b/src/Mods/mod_SCET.f90 index f207101..9c908ce 100644 --- a/src/Mods/mod_SCET.f90 +++ b/src/Mods/mod_SCET.f90 @@ -67,7 +67,8 @@ module SCET if (any(nprocbelow == [11,16,41]) .and. (pt34min > 0._dp)) then ! This branch for matching corrections in N4LL qt calculations, where pt34min > 0 ! The definition for W+jet and Z+jet divides by 1050 to match previous non-dynamic results. - getdynamictau = sqrt(pttwo(3,4,p)**2+puremass(p(3,:)+p(4,:))**2)/1050._dp + !getdynamictau = sqrt(pttwo(3,4,p)**2+puremass(p(3,:)+p(4,:))**2)/1050._dp + getdynamictau = pttwo(3,4,p)**2/400._dp elseif (any(nprocbelow == [11,16,41,204,210])) then getdynamictau = sqrt(pt(5,p)**2+puremass(p(3,:)+p(4,:))**2) elseif ( (kcase == kW_only) .or. (kcase == kW_1jet) & diff --git a/src/Need/mcfm_init.f b/src/Need/mcfm_init.f index 6597fac..3e35c27 100644 --- a/src/Need/mcfm_init.f +++ b/src/Need/mcfm_init.f @@ -58,7 +58,14 @@ c*********************************************************************** character(len=cfg_string_len) :: benchmark character(len=cfg_string_len) :: part - logical :: newPlotting + character(:), allocatable :: copt + integer :: io, ierr + character(len=255) :: imsg + + character(len=255) :: rundir + common/rundir/rundir + + include 'runstring.f' !$omp threadprivate(/pext/) @@ -178,33 +185,25 @@ c cutoff=cutoff*(sqrts/2000._dp)**2 c Set-up run name call setrunname(scale,facscale,cfg%config_directory) + open(newunit=io, file=trim(rundir)//"/"//trim(runname)//"_info.txt", status='replace', + & form='formatted', iostat=ierr, iomsg=imsg) + if (ierr == 0) then + allocate(character(len=2048) :: copt) + call get_command(copt) + write (io, '(A)') trim(copt) + + copt = compiler_version() + write (io, '(A)') "Compiler: " // copt + endif + close(unit=io) + c npart=9 is a dummy value, to ensure that all histograms are included npart=9 val=1.e-15_dp nprocbelow = nproc - call cfg_get(cfg,"histogram%newstyle",newPlotting) - - if (newPlotting .or. nproc == 1610 .or. nproc == 1650) then - call setup_plots() - call plots_allocate() - else - - ! This first call has two purposes: - ! 1. determine nplotmax for histogram allocation. - ! 2. must be run in parallel so that each thread initializes - ! in "tagbook" mode and takes runs fully in "tagplot" - ! for subsequent calls in the integration routines -!$omp parallel - call nplotter(p,val,val**2,1) -!$omp end parallel - -!$omp parallel - call initHistogramStorage(nplotmax) -!$omp end parallel - call initMasterStorage(nplotmax) - - endif + call setup_plots() + call plots_allocate() p(:,:) = 0._dp diff --git a/src/Need/mcfmmain.f b/src/Need/mcfmmain.f index b05594a..934703a 100644 --- a/src/Need/mcfmmain.f +++ b/src/Need/mcfmmain.f @@ -20,8 +20,6 @@ c * c*********************************************************************** real(dp):: integ,integ_err,r,er - logical:: dryrun - common/dryrun/dryrun integer*4 old_cw diff --git a/src/Need/parseinput.f90 b/src/Need/parseinput.f90 index 75d768d..78ef88d 100644 --- a/src/Need/parseinput.f90 +++ b/src/Need/parseinput.f90 @@ -30,6 +30,7 @@ contains call cfg_add(cfg, "general%nproc", 1, "process number") call cfg_add(cfg, "general%part", "nlo", "part: lo, nlo, nlocoeff, nnlocoeff") call cfg_add(cfg, "general%runstring", "run", "string identifying the run") + call cfg_add(cfg, "general%runname", "", "Filename prefix of snapshot to write out.") call cfg_add(cfg, "general%rundir", "./", "directory for output and snapshot files") call cfg_add(cfg, "general%sqrts", sqrts, "center of mass energy") call cfg_add(cfg, "general%ih1", +1, "ih1: +1 for proton, -1 for antiproton") @@ -57,8 +58,13 @@ contains call cfg_add(cfg, "resummation%res_range", [0._dp, 80._dp], "integration range of purely resummed contribution") call cfg_add(cfg, "resummation%resexp_range", [1._dp, 80._dp], "integration range of fixed-order expansion of resummed contribution") call cfg_add(cfg, "resummation%fo_cutoff", 1d0, "minimum qT for fixed-order calculation") - call cfg_add(cfg, "resummation%transitionswitch", 0.4d0, "transition switch to be used in plotting routines") + call cfg_add(cfg, "resummation%transitionswitch", 0.4_dp, "transition switches to be used in plotting routines") + call cfg_add(cfg, "resummation%transitionswitches", [real(dp) :: ], "transition switches to be used in plotting routines", dynamic_size=.true.) call cfg_add(cfg, "resummation%scalevar_rapidity", .false., "rapidity scale variation") + call cfg_add(cfg, "resummation%resexp_linPC", .false., "compute linear power corrections instead of resexp") + + call cfg_add(cfg, "resummation%LambdaNP_qq", 0._dp, "Lambda NP for qq") + call cfg_add(cfg, "resummation%LambdaNP_gg", 0._dp, "Lambda NP for qq") ! additional optional settings in [general] !call cfg_add(cfg, "general%nevtrequested", 0, "") ! JC: removed 3/8/22 (no longer operational) @@ -71,11 +77,19 @@ contains call cfg_add(cfg, "masses%mt", 173d0, "Top-quark mass") call cfg_add(cfg, "masses%mb", 4.66d0, "Bottom-quark mass") call cfg_add(cfg, "masses%mc", 1.275d0, "Charm-quark mass") + call cfg_add(cfg, "masses%wmass", 80.385_dp, "W-boson mass") + call cfg_add(cfg, "masses%zmass", 91.1876_dp, "Z-boson mass") + call cfg_add(cfg, "masses%wwidth", 2.091_dp, "W-boson width") + call cfg_add(cfg, "masses%zwidth", 2.4952_dp, "Z-boson width") call cfg_add(cfg, "masses%CKMdiag", .false., "Set CKM to diagonal and ignore CKMrotate") call cfg_add(cfg, "masses%CKMrotate", .false., "Implement VCKM by rotating pdfs") + call cfg_add(cfg, "couplings%ewscheme", 1, "1: Gmu scheme") + call cfg_add(cfg, "couplings%sin2th", 0.2234d0, "sin^2(theta) on-shell") + call cfg_add(cfg, "couplings%G_f", 1.16639e-5_dp, "Gf") + call cfg_add(cfg, "couplings%aemmz", 1._dp/128.89_dp, "alpha_EM, default at m_Z = 1/128.89") ! [scales] call cfg_add(cfg, "scales%renscale", 1d0, "Renormalization scale") call cfg_add(cfg, "scales%facscale", 1d0, "Factorization scale") @@ -260,6 +274,7 @@ contains use MCFMSettings, only: newStyleHistograms, resexp_linPC use qtResummation_params, only: qtcutoff, qtminRes, qtmaxRes, qtminResexp, & qtmaxResexp, transitionSwitch, & + transitionSwitches, & ! scalevar_rapidity, & scalevar_rapidity, & enable_fixed_y, fixed_y, enable_dsigma_dQ, generations, & @@ -271,6 +286,7 @@ contains use ggHwilson, only: expansionorder, Wilsonorder use SCET, only: useQT use ptveto + use iso_fortran_env implicit none include 'cplx.h' @@ -397,7 +413,6 @@ contains call cfg_get(cfg, "general%ih2", ih2) call cfg_get(cfg, "general%zerowidth", zerowidth) call cfg_get(cfg, "general%removebr", removebr) - call cfg_get_add(cfg, "general%ewscheme",ewscheme,ewscheme,"ew scheme") ! [nnlo] call cfg_get(cfg, "nnlo%dynamictau", dynamictau) @@ -420,16 +435,22 @@ contains ! [masses] call cfg_get(cfg, "masses%hmass", hmass) - call cfg_get_add(cfg, "masses%wmass", wmass_inp, wmass_inp, "W mass") - call cfg_get_add(cfg, "masses%zmass", zmass_inp, zmass_inp, "Z mass") - call cfg_get_add(cfg, "masses%wwidth", wwidth, wwidth, "W width") - call cfg_get_add(cfg, "masses%zwidth", zwidth, zwidth, "Z width") + call cfg_get(cfg, "masses%wmass", wmass_inp) + call cfg_get(cfg, "masses%zmass", zmass_inp) + call cfg_get(cfg, "masses%wwidth", wwidth) + call cfg_get(cfg, "masses%zwidth", zwidth) call cfg_get(cfg, "masses%mt", mt) call cfg_get(cfg, "masses%mb", mb) call cfg_get(cfg, "masses%mc", mc) call cfg_get(cfg, "masses%CKMdiag", CKMdiag) call cfg_get(cfg, "masses%CKMrotate", CKMrotate) + ! [couplings] + call cfg_get(cfg, "couplings%ewscheme", ewscheme) + call cfg_get(cfg, "couplings%sin2th", xw_inp) + call cfg_get(cfg, "couplings%G_f", gf_inp) + call cfg_get(cfg, "couplings%aemmz", aemmz_inp) + ! Yukawas mt_yuk=mt if (abs(mb) > 1.e-8_dp) then @@ -857,13 +878,30 @@ contains qtcutoff = res_fo_cutoff - call cfg_get(cfg, "resummation%transitionswitch", transitionSwitch) + if (cfg_var_configadded(cfg, "resummation%transitionswitches")) then + if (cfg_var_size(cfg, "resummation%transitionswitches") > 0) then + allocate(transitionswitches(cfg_var_size(cfg, "resummation%transitionswitches"))) + call cfg_get(cfg, "resummation%transitionswitches", transitionswitches) + + ! legacy plotting support + transitionswitch = transitionswitches(1) + else + allocate(transitionswitches(1)) + call cfg_get(cfg, "resummation%transitionswitch", transitionswitches(1)) + endif + else + allocate(transitionswitches(1)) + call cfg_get(cfg, "resummation%transitionswitch", transitionswitches(1)) + endif + call cfg_get_add(cfg, "benchmark%enable_fixed_y", enable_fixed_y, .false., "enable fixed y") call cfg_get_add(cfg, "benchmark%fixed_y", fixed_y, 0.0d0, "fixed y value") call cfg_get_add(cfg, "benchmark%enable_dsigma_dQ", enable_dsigma_dQ, .false., "calculate dsigma/dQ (mV)") call cfg_get_add(cfg, "benchmark%generations", generations, 5, "number nf to sum over") call cfg_get_add(cfg, "benchmark%fix_alphas_nf5", fix_alphas_nf5, .false., "fix alphas to 5-flavor 3-loop running") + call cfg_get(cfg, "resummation%resexp_linPC", resexp_linPC) + if (rank == 0) then if (resexp_linPC) then write (*,*) "******************************************" @@ -1060,7 +1098,7 @@ contains if (any(kpart==[ksnlo,knnlo,kn3lo]) & .or. ((kpart==kresummed) .and. (usept) .and. (kresorder==6))) then if (useqt .or. useqt_nnlo) then - if (any(nproc == [31,32])) then + if (any(nproc == [1,6,31,32])) then ! useqt is the new implementation up to N3LO from 2207.07056, but only set up for Z useqt = .true. useqt_nnlo = .false. diff --git a/src/Need/setrunname.f b/src/Need/setrunname.f index fd807c7..27327c3 100644 --- a/src/Need/setrunname.f +++ b/src/Need/setrunname.f @@ -4,6 +4,8 @@ ! subroutine setrunname(scalestart,fscalestart,werkdir) use singletop2_scale_m, only: use_DDIS + use m_config + use parseinput, only : cfg implicit none include 'types.f' include 'flags.f' @@ -18,7 +20,6 @@ include 'taucut.f' include 'ewcorr.f' real(dp):: scalestart,fscalestart - integer:: nlength character(len=255):: outlabel1 character(len=1024):: runname character(len=*), intent(in) :: werkdir @@ -28,7 +29,13 @@ character(len=20):: part,kpartstring character(len=6):: case,kcasestring common/runname/runname - common/nlength/nlength + character(len=cfg_string_len) :: runname_in + + call cfg_get(cfg, "general%runname", runname_in) + if (trim(runname_in) /= "") then + runname = trim(runname_in)//'_'//trim(runstring) + return + endif c if (abs(scalestart-fscalestart) < 1._dp) then c--- if the scales are the same, use this scale as the label @@ -116,8 +123,6 @@ c--- add working directory, if necessary runname=trim(werkdir)//'/'//trim(runname) endif - nlength=len(trim(runname)) - return end diff --git a/src/Need/writereference.f b/src/Need/writereference.f index f71adc8..d0cc99a 100644 --- a/src/Need/writereference.f +++ b/src/Need/writereference.f @@ -83,6 +83,14 @@ c common/writerefs/writerefs write (6,58) '* *' endif + if (any(nproc == [1,6,11,16])) then + write (6,58) '* Third order QCD predictions *' + write (6,58) '* for fiducial W-boson production *' + write (6,58) '* J.M. Campbell and T. Neumann *' + write (6,58) '* arXiv: 2308.15382 *' + write (6,58) '* *' + endif + if (any(nproc == [31,32]) .and. kpart == kresummed .or. kpart == kn3lo) then write (6,58) '* Fiducial Drell-Yan production at the LHC *' write (6,58) '* improved by qT-resummation at N4LL+N3LO *' diff --git a/src/Parton/lhapdf_interface.cpp b/src/Parton/lhapdf_interface.cpp index e6efe5b..7e0af3c 100644 --- a/src/Parton/lhapdf_interface.cpp +++ b/src/Parton/lhapdf_interface.cpp @@ -85,11 +85,26 @@ int lhapdf_orderqcd(LHAPDF::PDF* pdf) return pdf->orderQCD(); } -LHAPDF::PDFUncertainty lhapdf_computeUncertainty(LHAPDF::PDF *pdf, double *values, int nval) +struct mypdfuncertainty +{ + double central; + double errplus; + double errminus; + double errsymm; +}; + +struct mypdfuncertainty lhapdf_computeUncertainty(LHAPDF::PDF *pdf, double *values, int nval) { LHAPDF::PDFSet set = pdf->set(); + LHAPDF::PDFUncertainty uncert; + mypdfuncertainty myuncert; std::vector vals(&values[0], &values[nval]); - return set.uncertainty(vals); + uncert = set.uncertainty(vals); + myuncert.central = uncert.central; + myuncert.errplus = uncert.errplus; + myuncert.errminus = uncert.errminus; + myuncert.errsymm = uncert.errsymm; + return myuncert; } void lhapdf_getconfig(LHAPDF::PDF *pdf, char key[], char* value, int vlen) diff --git a/src/Phase/gen_njets.f b/src/Phase/gen_njets.f index 7dc6044..75c274c 100644 --- a/src/Phase/gen_njets.f +++ b/src/Phase/gen_njets.f @@ -65,8 +65,8 @@ c--- relax cuts for real and calculations with SCET above taucut if (origkpart == kresummed) then pt = r(ijet)**2*sqrts/2._dp xjac = pt*sqrts/2._dp*2._dp*r(ijet) - elseif (((usescet .eqv. .false.) .and. (kpart == kvirt)) .or. - & ((usescet) .and. (abovecut .eqv. .false.) .and. (origkpart /= kn3lo))) then + elseif (ptjetmin > 0._dp .and. (((usescet .eqv. .false.) .and. (kpart == kvirt)) .or. + & ((usescet) .and. (abovecut .eqv. .false.) .and. (origkpart /= kn3lo)))) then call genpt(r(ijet),ptjetmin,.true.,pt,xjac) else ! linear diff --git a/src/Procdep/fragint.f b/src/Procdep/fragint.f index fdcb1d4..4f95109 100644 --- a/src/Procdep/fragint.f +++ b/src/Procdep/fragint.f @@ -399,11 +399,7 @@ c--- loop over all PDF error sets, if necessary if (bin) then - if (newStyleHistograms) then - call nplotter_new(pjet,val) - else - call nplotter(pjet,val,val2,0) - endif + call nplotter_new(pjet,val) endif c --- Check weights : diff --git a/src/Procdep/lowint.f b/src/Procdep/lowint.f index ac52075..0203160 100644 --- a/src/Procdep/lowint.f +++ b/src/Procdep/lowint.f @@ -110,7 +110,7 @@ c--- ensure isolation code does not think this is fragmentation piece if (maxPDFsets > 0) pdfreweight(:) = 0._dp - if (usescet .and. kcase == kZ_2jet) then + if (usescet .and. (kcase == kZ_2jet .or. kcase == kW_2jet)) then if ( genps_z1jet_r(r,p,pswt) .eqv. .false. ) then goto 999 endif @@ -971,11 +971,7 @@ c--- for EW corrections, make additional weight available inside common block if (kewcorr /= knone) then wt_noew=lowint_noew*wgt endif - if (newStyleHistograms) then - call nplotter_new(pjet,val) - else - call nplotter(pjet,val,val2,0) - endif + call nplotter_new(pjet,val) endif if (includeTaucutgrid(0) .eqv. .false.) then diff --git a/src/Procdep/realint.f b/src/Procdep/realint.f index f72c304..e4dd363 100644 --- a/src/Procdep/realint.f +++ b/src/Procdep/realint.f @@ -299,7 +299,7 @@ c small safety cuts goto 999 endif else - if (.not. passed_smallnew(p,npart,1d-9)) then + if (.not. passed_smallnew(p,npart,1d-12)) then goto 999 endif endif @@ -1360,11 +1360,7 @@ c endif currentNd = nd call getptildejet(nd,pjet) call dotem(nvec,pjet,s) - if (newStyleHistograms) then - call nplotter_new(pjet,xmsq(nd)*wgt) - else - call nplotter(pjet,xmsq(nd)*wgt,(xmsq(nd)*wgt)**2,nd) - endif + call nplotter_new(pjet,xmsq(nd)*wgt) endif endif enddo ! loop over dipole contributions diff --git a/src/Procdep/resint.f90 b/src/Procdep/resint.f90 index f1bf434..f49c1ea 100644 --- a/src/Procdep/resint.f90 +++ b/src/Procdep/resint.f90 @@ -20,6 +20,156 @@ contains + subroutine computeqtsubint(p, xx, qtsubout, muMult_in) + use ieee_arithmetic + use qtResummation + use SCET + use qqb_z_hardfun + use Scalevar + use MCFMSetupPlots + use PDFerrors + implicit none + include 'src/Inc/mxdim.f' + include 'src/Inc/mxpart.f' + include 'src/Inc/npart.f' + include 'src/Inc/nproc.f' + include 'src/Inc/kpart.f' + include 'src/Inc/kprocess.f' + include 'src/Inc/nf.f' + include 'src/Inc/energy.f' + include 'src/Inc/dynamicscale.f' + include 'src/Inc/constants.f' + include 'src/Inc/initialscales.f' + include 'src/Inc/xmin.f' + include 'src/Inc/scale.f' + include 'src/Inc/qcdcouple.f' + include 'src/Inc/ewcharge.f' + include 'src/Inc/ewcouple.f' + + real(dp), intent(in) :: p(mxpart,4) + real(dp), optional :: muMult_in + real(dp), intent(out) :: qtsubout + real(dp), intent(in) :: xx(2) + + real(dp) :: muMult + real(dp), allocatable :: xmsq(:) + real(dp), allocatable :: results(:) + + real(dp) :: qsq + real(dp) :: msqall(-nf:nf,-nf:nf,0:3) + real(dp) :: msq(-nf:nf,-nf:nf) + real(dp) :: hard(0:3) + integer :: j,k, order + real(dp) :: nfv, hard3(3) + real(dp) :: twomass + + logical:: bin + common/bin/bin + + if (present(muMult_in)) then + muMult = muMult_in + else + muMult = 1._dp + endif + + if (dynamicscale) then + call scaleset(initscale*muMult, initfacscale*muMult, p) + else + call usescales(initscale*muMult, initfacscale*muMult) + endif + + if (origkpart == ksnlo) then + order = 1 + elseif (origkpart == knnlo) then + order = 2 + elseif (origkpart == kn3lo) then + order = 3 + else + write (*,*) __FILE__//": undefined kpart, line ", __LINE__ + error stop + endif + + if (kcase == kZ_only) then + qsq = twomass(1,2,p)**2 + + ! must not contain the axial singlet + ! contribution since it has no counterpart in the Z+jet + ! NNLO calculation + call qqb_z_loop(p, musq, msqall, withaxsinglet=.false.) + + msqall(:,:,1) = msqall(:,:,1) * (as/4._dp/pi) + msqall(:,:,2) = msqall(:,:,2) * (as/4._dp/pi)**2 + msqall(:,:,3) = msqall(:,:,3) * (as/4._dp/pi)**3 + elseif (kcase == kW_only) then + qsq = twomass(1,2,p)**2 + + call hardqq3(qsq,musq,NfV,hard3) + call qqb_w(p, msq) + + + where(msq /= 0._dp) msqall(:,:,0) = msq(:,:) + where(msq /= 0._dp) msqall(:,:,1) = hard3(1)*msq(:,:)*(as/2._dp/pi) + where(msq /= 0._dp) msqall(:,:,2) = hard3(2)*msq(:,:)*(as/2._dp/pi)**2 + where(msq /= 0._dp) msqall(:,:,3) = hard3(3)*msq(:,:)*(as/2._dp/pi)**3 + + else + write (*,*) __FILE__//": undefined kcase, line ", __LINE__ + error stop + endif + + if (muMult == 1._dp) then + + allocate(xmsq(size(tcutarray)+1), source=0._dp) + + do j=-5,5; do k =-5,5 + ! cycle if no contribution through NLO + if ( all(msqall(j,k,[0,1]) == 0._dp) ) cycle + + ! coefficient in as/4/pi + hard(0) = msqall(j,k,0) + hard(1) = msqall(j,k,1) / msqall(j,k,0) / (as/4._dp/pi) + hard(2) = msqall(j,k,2) / msqall(j,k,0) / (as/4._dp/pi)**2 + hard(3) = msqall(j,k,3) / msqall(j,k,0) / (as/4._dp/pi)**3 + + call qtsubtraction(p, qsq, xx(1), xx(2), j, k, hard, & + sqrt(musq), order, results) + + xmsq = xmsq + results + + enddo; enddo + + if (doMultitaucut) then + scetreweight(:) = 0._dp + if (xmsq(1) /= 0._dp) then + scetreweight(:) = xmsq(2:size(tcutarray)+1) / xmsq(1) + endif + endif + else ! we do scale variation + allocate(xmsq(1), source=0._dp) + + do j=-5,5; do k =-5,5 + ! cycle if no contribution through NLO + if ( all(msqall(j,k,[0,1]) == 0._dp) ) cycle + + ! coefficient in as/4/pi + hard(0) = msqall(j,k,0) + hard(1) = msqall(j,k,1) / msqall(j,k,0) / (as/4._dp/pi) + hard(2) = msqall(j,k,2) / msqall(j,k,0) / (as/4._dp/pi)**2 + hard(3) = msqall(j,k,3) / msqall(j,k,0) / (as/4._dp/pi)**3 + + call qtsubtraction(p, qsq, xx(1), xx(2), j, k, hard, & + sqrt(musq), order, results, muMult=muMult) + + xmsq = xmsq + results + + enddo; enddo + + endif + + qtsubout = xmsq(1) + + end subroutine computeqtsubint + function qtsubint(r,wgt) use ieee_arithmetic use qtResummation @@ -27,6 +177,7 @@ use qqb_z_hardfun use Scalevar use MCFMSetupPlots + use PDFerrors implicit none include 'src/Inc/mxdim.f' include 'src/Inc/mxpart.f' @@ -50,10 +201,11 @@ real(dp) :: p(mxpart,4), pjet(mxpart,4), pswt, xx(2) real(dp) :: val,flux,qsq real(dp) :: msqall(-nf:nf,-nf:nf,0:3) - real(dp), allocatable :: results(:) - real(dp), allocatable :: xmsq(:) + real(dp) :: msq(-nf:nf,-nf:nf) + real(dp) :: xmsq, xmsqMu real(dp) :: hard(0:3) integer :: j,k, order + real(dp) :: nfv, hard3(3) ! functions logical :: includedipole @@ -68,8 +220,13 @@ qtsubint = 0._dp p(:,:) = 0._dp pjet(:,:) = 0._dp + currentPDF = 0 + + if (doScalevar .and. bin) then + scalereweight(:) = 1._dp + endif - if (any(nproc == [31,32])) then ! Z -> e e + if (any(nproc == [1,6,31,32])) then ! Z -> e e or W -> e nu npart = 2 ! generate phase space at zero pT if (gen_res(r,0._dp,p,pswt) .eqv. .false. ) then @@ -107,70 +264,35 @@ call storeptilde(0,p) call getptildejet(0,pjet) - if (dynamicscale) then - call scaleset(initscale, initfacscale, p) - else - call usescales(initscale, initfacscale) - endif - - flux = fbGeV2/(two*xx(1)*xx(2)*sqrts**2) - - if (origkpart == ksnlo) then - order = 1 - elseif (origkpart == knnlo) then - order = 2 - elseif (origkpart == kn3lo) then - order = 3 - else - write (*,*) __FILE__//": undefined kpart, line ", __LINE__ - error stop - endif - - if (kcase == kZ_only) then - qsq = twomass(1,2,p)**2 - - ! must not contain the axial singlet - ! contribution since it has no counterpart in the Z+jet - ! NNLO calculation - call qqb_z_loop(p, musq, msqall, withaxsinglet=.false.) + ! central value + call computeqtsubint(p, xx, xmsq) - msqall(:,:,1) = msqall(:,:,1) * (as/4._dp/pi) - msqall(:,:,2) = msqall(:,:,2) * (as/4._dp/pi)**2 - msqall(:,:,3) = msqall(:,:,3) * (as/4._dp/pi)**3 - else - write (*,*) __FILE__//": undefined kcase, line ", __LINE__ - error stop + if (doScalevar .and. bin) then + if (xmsq /= 0._dp) then + ! multiplies by two + call computeqtsubint(p, xx, xmsqMu, muMult_in=scalevarmult(1)) + scalereweight(1) = xmsqMu / xmsq + ! divides by two + call computeqtsubint(p, xx, xmsqMu, muMult_in=scalevarmult(2)) + scalereweight(2) = xmsqMu / xmsq + endif endif - allocate(xmsq(size(tcutarray)+1), source=0._dp) - - do j=-5,5; do k =-5,5 - ! cycle if no contribution through NLO - if ( all(msqall(j,k,[0,1]) == 0._dp) ) cycle - - ! coefficient in as/4/pi - hard(0) = msqall(j,k,0) - hard(1) = msqall(j,k,1) / msqall(j,k,0) / (as/4._dp/pi) - hard(2) = msqall(j,k,2) / msqall(j,k,0) / (as/4._dp/pi)**2 - hard(3) = msqall(j,k,2) / msqall(j,k,0) / (as/4._dp/pi)**2 - - call qtsubtraction(p, qsq, xx(1), xx(2), j, k, hard, & - sqrt(musq), order, 1._dp, results) - - xmsq = xmsq + results - - enddo; enddo - + flux = fbGeV2/(two*xx(1)*xx(2)*sqrts**2) + qtsubint = flux*pswt*xmsq/BrnRat - if (doMultitaucut) then - scetreweight(:) = 0._dp - if (xmsq(1) /= 0._dp) then - scetreweight(:) = xmsq(2:size(tcutarray)+1) / xmsq(1) + if ((maxPDFsets > 0) .and. bin) then + pdfreweight(:) = 0._dp + do j=1,maxPDFsets + currentPDF = j + if (doPDFAlphas) then + call updateAlphas(scale) endif + call computeqtsubint(p, xx, xmsq) + pdfreweight(currentPDF) = (qtsubint - flux*pswt*xmsq/BrnRat)*wgt + enddo endif - qtsubint = flux*pswt*xmsq(1)/BrnRat - val=qtsubint*wgt if (ieee_is_nan(val) .OR. ( .NOT. ieee_is_finite(val))) then @@ -394,7 +516,6 @@ endif if (.not. passed_smallnew(p,npart,safety)) then - 999 continue resint = 0._dp return endif @@ -491,6 +612,12 @@ call nplotter_new(pjet,val) endif + return + +! alternate return + 999 continue + resint = 0._dp + return end function resint subroutine computeResint(r, p, p_born, qt, order, xx, resout, resexpout, muMultHard_in, muMultRes_in) @@ -581,7 +708,12 @@ if (kcase == kW_only) then qsq = twomass(1,2,p)**2 + if (order >= 7) then + call hardqq3(qsq,musq,NfV,hard3) + hard(1:2) = hard3(1:2) + else call hardqq(qsq, musq, hard) + endif if (boostBorn) then call qqb_w(p, msq) else @@ -592,6 +724,10 @@ where(msq /= 0._dp) msqall(:,:,1) = hard(1)*msq(:,:)*(as/2._dp/pi) where(msq /= 0._dp) msqall(:,:,2) = hard(2)*msq(:,:)*(as/2._dp/pi)**2 + if (order >= 7) then + call hardqq3(qsq,musq,NfV,hard3) + where(msq /= 0._dp) msqall(:,:,3) = hard3(3)*msq(:,:)*(as/2._dp/pi)**3 + endif elseif (.false. .and. kcase == kZ_only) then qsq = twomass(1,2,p)**2 call hardqq(qsq, musq, hard) @@ -1332,16 +1468,6 @@ if (pttwo(3,4,p) < 1.e-3_dp) goto 888 endif - - - !call smallnew(p,npart,*999) - if (.false.) then - 888 continue - resexpint = 0._dp - return - endif - - if (any(ieee_is_nan(p(1:npart+2,:)) .eqv. .true.)) then !write(6,*) 'Discarding NaN or infinite phase space point' resexpint = 0._dp @@ -1469,6 +1595,18 @@ call nplotter_new(pjet,val) endif + ! boost for calculation of linPC at small qT + if (resexp_linPC) then + val = val / qt**3 + endif + + return + +! alternate return + 888 continue + resexpint = 0._dp + return + end function resexpint end module diff --git a/src/Procdep/scetint.f b/src/Procdep/scetint.f index 3b610b3..a349e36 100644 --- a/src/Procdep/scetint.f +++ b/src/Procdep/scetint.f @@ -162,6 +162,7 @@ c--- bother calculating the matrix elements for it, instead bail out ! PDF VARIATION WITH CENTRAL SCALE if ((maxPDFsets > 0) .and. bin) then + pdfreweight(:) = 0._dp do j=1,maxPDFsets currentPDF = j if (doPDFAlphas) then @@ -178,11 +179,7 @@ c--- bother calculating the matrix elements for it, instead bail out if (bin) then includeTaucutgrid(0) = .true. - if (newStyleHistograms) then - call nplotter_new(pjet,val) - else - call nplotter(pjet,val,val2,0) - endif + call nplotter_new(pjet,val) endif if (enable_reweight_user) then diff --git a/src/Procdep/virtint.f b/src/Procdep/virtint.f index 52817f0..0f7d906 100644 --- a/src/Procdep/virtint.f +++ b/src/Procdep/virtint.f @@ -2644,11 +2644,7 @@ c--- for EW corrections, make additional weight available inside common block if (kewcorr /= knone) then wt_noew=virtint_noew*wgt endif - if (newStyleHistograms) then - call nplotter_new(pjet,val) - else - call nplotter(pjet,val,val2,0) - endif + call nplotter_new(pjet,val) endif diff --git a/src/QT/qtint.f b/src/QT/qtint.f index 22d8af7..2d97805 100644 --- a/src/QT/qtint.f +++ b/src/QT/qtint.f @@ -184,12 +184,7 @@ c pause if (bin) then includeTaucutgrid(0) = .true. - - if (newStyleHistograms) then - call nplotter_new(pjet,val) - else - call nplotter(pjet,val,val2,0) - endif + call nplotter_new(pjet,val) endif if (enable_reweight_user) then diff --git a/src/SCET/setupscet.f b/src/SCET/setupscet.f index 254ea60..8785f10 100644 --- a/src/SCET/setupscet.f +++ b/src/SCET/setupscet.f @@ -18,6 +18,12 @@ c---- additional jet; stops with error message if appropriate. & .or. (nproc == 31) .or. (nproc == 32) .or. (nproc == 33)) then nprocabove=nproc+10 ntau=0 + elseif (nproc == 11) then ! W+jet + nprocabove=22 + ntau = 1 + elseif (nproc == 16) then + nprocabove=27 + ntau = 1 elseif (nproc == 61) then nprocabove=461 ntau=0 diff --git a/src/SCET1j/mod_nnlo_w1jet.f90 b/src/SCET1j/mod_nnlo_w1jet.f90 index dda4343..9963497 100644 --- a/src/SCET1j/mod_nnlo_w1jet.f90 +++ b/src/SCET1j/mod_nnlo_w1jet.f90 @@ -428,6 +428,9 @@ module nnlo_w1jet else write(6,*) 'region not found:' write(6,*) 's12,s13,s23',s(q1,q2),s(q1,q3),s(q2,q3) + hard1 = 0._dp + hard2 = 0._dp + return endif ! write(6,*) 'iperm,region',iperm,region diff --git a/src/SCET1j/mod_nnlo_z1jet.f90 b/src/SCET1j/mod_nnlo_z1jet.f90 index 143e58d..f8f8633 100644 --- a/src/SCET1j/mod_nnlo_z1jet.f90 +++ b/src/SCET1j/mod_nnlo_z1jet.f90 @@ -209,7 +209,7 @@ module nnlo_z1jet #else write (*,*) "INVALID PSCHOICE" - pause + stop #endif #undef PSCHOICE @@ -509,7 +509,7 @@ module nnlo_z1jet #else write (*,*) "INVALID PSCHOICE" - pause + stop #endif genps_z1jet_rr = .true. diff --git a/src/Singletop/boostx.f b/src/Singletop/boostx.f index 39464cf..5ef31a5 100644 --- a/src/Singletop/boostx.f +++ b/src/Singletop/boostx.f @@ -10,13 +10,16 @@ c--- Boost input vector p_in to output vector p_out using the same c--- transformation as required to boost massive vector pt to ptt real(dp):: p_in(4),pt(4),ptt(4),p_out(4), - & p_tmp(4),beta(3),mass,gam,bdotp + & p_tmp(4),beta(3),mass,gam,bdotp,dotpr integer:: j - mass=pt(4)**2-pt(1)**2-pt(2)**2-pt(3)**2 + mass=dotpr(pt,pt) if (mass < 0._dp) then write(6,*) 'mass**2 < 0 in boostx.f, mass**2=',mass - stop + !stop + mass = 0._dp + else + mass=sqrt(mass) endif mass=sqrt(mass) diff --git a/src/User/CMakeLists.txt b/src/User/CMakeLists.txt index bafee97..05ac505 100644 --- a/src/User/CMakeLists.txt +++ b/src/User/CMakeLists.txt @@ -1,48 +1,48 @@ target_sources(objlib PRIVATE -autoplots.f -bookrelew.f +#autoplots.f +#bookrelew.f gencuts_user.f90 -genplots.f +#genplots.f mdata.f -nplotter_4ftwdk.f -nplotter_auto.f -nplotter_dirgam.f -nplotter_dm_mongam.f -nplotter_dm_monj.f -nplotter.f -nplotter_gamgam.f -nplotter_generic.f -nplotter_gmgmjt.f -nplotter_Hjets.f -nplotter_ktopanom.f -nplotter_qqZZqq.f -nplotter_tbbar.f -nplotter_trigam.f -nplotter_ttbar.f -nplotter_tt_tot.f -nplotter_ttw.f -nplotter_ttZ.f -nplotter_twojet.f -nplotter_user.f90 -nplotter_Vgamma.f -nplotter_VHbbarHXSWG.f -nplotter_VHgaga.f -nplotter_VHWW.f -nplotter_VV.f -nplotter_Wbbmas.f -nplotter_wgamgam.f -nplotter_Wjets.f -nplotter_W_only.f -nplotter_WWjet.f -nplotter_Zbbbar.f -nplotter_zgamgam.f -nplotter_zgamjet.f -nplotter_Zjets.f -nplotter_Z_only.f -nplotter_Ztjdk.f -nplotter_Ztj.f -nplotter_ZZlept.f -topreconstruct.f +#nplotter_4ftwdk.f +#nplotter_auto.f +#nplotter_dirgam.f +#nplotter_dm_mongam.f +#nplotter_dm_monj.f +#nplotter.f +#nplotter_gamgam.f +#nplotter_generic.f +#nplotter_gmgmjt.f +#nplotter_Hjets.f +#nplotter_ktopanom.f +#nplotter_qqZZqq.f +#nplotter_tbbar.f +#nplotter_trigam.f +#nplotter_ttbar.f +#nplotter_tt_tot.f +#nplotter_ttw.f +#nplotter_ttZ.f +#nplotter_twojet.f +#nplotter_user.f90 +#nplotter_Vgamma.f +#nplotter_VHbbarHXSWG.f +#nplotter_VHgaga.f +#nplotter_VHWW.f +#nplotter_VV.f +#nplotter_Wbbmas.f +#nplotter_wgamgam.f +#nplotter_Wjets.f +#nplotter_W_only.f +#nplotter_WWjet.f +#nplotter_Zbbbar.f +#nplotter_zgamgam.f +#nplotter_zgamjet.f +#nplotter_Zjets.f +#nplotter_Z_only.f +#nplotter_Ztjdk.f +#nplotter_Ztj.f +#nplotter_ZZlept.f +#topreconstruct.f nplotter_Z_new.f90 nplotter_W_new.f90 diff --git a/src/User/mdata.f b/src/User/mdata.f index 32248c1..59898cf 100644 --- a/src/User/mdata.f +++ b/src/User/mdata.f @@ -31,12 +31,15 @@ c input here. You have to know what you're doing. include 'ewinput.f' - data ewscheme / +1 / ! Chooses EW scheme - data Gf_inp / 1.16639e-5_dp / ! G_F - data aemmz_inp / 7.7585538055706e-03_dp / ! alpha_EM(m_Z)=1/128.89 - data xw_inp / 0.2223_dp / ! sin^2(theta_W) - data wmass_inp / 80.385_dp / ! W mass - data zmass_inp / 91.1876_dp / ! Z mass + !data ewscheme / +1 / ! Chooses EW scheme + !data Gf_inp / 1.16639e-5_dp / ! G_F + !data aemmz_inp / 7.7585538055706e-03_dp / ! alpha_EM(m_Z)=1/128.89 + !data xw_inp / 0.2312_dp / ! sin^2(theta_W) + !data xw_inp / 0.2234_dp / ! sin^2(theta_W) + !data GF_inp / 1.1663787e-05_dp / + !data aemmz_inp / 7.56239563671356e-3_dp / + !data wmass_inp / 80.385_dp / ! W mass + !data zmass_inp / 91.1876_dp / ! Z mass c*********************************************************************** @@ -72,7 +75,7 @@ c--- Widths: note that the top width is calculated in the program c--- The W width of 2.1054 is derived using the measured BR of c--- 10.86 +/- 0.09 % (PDG 2015) and the LO partial width calculation c--- for Mw=80.385 GeV - data wwidth,zwidth/2.085_dp,2.4952_dp/ + !data wwidth,zwidth/2.091_dp,2.4952_dp/ data tauwidth/2.269e-12_dp/ c--- Number of active flavours in the initial state: this parameter c--- may be changed in the program for some processes @@ -96,7 +99,6 @@ c data Vud , Vus , Vub , c & Vcd , Vcs , Vcb c & /1._dp,0._dp,0.000_dp, c & 0._dp,1._dp,0.000_dp/ -c Matches MATRIX use data Vud , Vus , Vub , & Vcd , Vcs , Vcb & /0.97417_dp,0.2248_dp,0.00409_dp, diff --git a/src/W2jet/qqb_w2jet_v.f b/src/W2jet/qqb_w2jet_v.f index b60ffbf..ec8cb42 100644 --- a/src/W2jet/qqb_w2jet_v.f +++ b/src/W2jet/qqb_w2jet_v.f @@ -75,6 +75,7 @@ c common/runstring/runstring if (first) then first=.false. +!$omp master if (rank == 0) then if ((Gflag) .or. (QandGflag)) then write(*,*) 'Using QQGG (VIRTUAL) matrix elements' @@ -100,6 +101,7 @@ c write(*,*) '[SLC is 1/N]' stop endif endif +!$omp end master endif if (useblha == 1) colourchoice=0 ! ensure calculation of full colour result diff --git a/src/W2jet/qqb_wm2jet_g.f b/src/W2jet/qqb_wm2jet_g.f index b34aa22..dd0a4a8 100644 --- a/src/W2jet/qqb_wm2jet_g.f +++ b/src/W2jet/qqb_wm2jet_g.f @@ -55,6 +55,7 @@ c ---> f(p5)+f(p6) if (first) then first=.false. +!$omp master if (rank == 0) then if ((Gflag) .or. (QandGflag)) then write(*,*) 'Using QQGG+G (REAL) matrix elements' @@ -80,6 +81,7 @@ c write(*,*) '[SLC is 1/N]' stop endif endif +!$omp end master endif msq(:,:)=0._dp diff --git a/src/W2jet/qqb_wp2jet_g.f b/src/W2jet/qqb_wp2jet_g.f index 6014361..cb54ef0 100644 --- a/src/W2jet/qqb_wp2jet_g.f +++ b/src/W2jet/qqb_wp2jet_g.f @@ -56,6 +56,7 @@ c ---> f(p5)+f(p6) if (first) then first=.false. +!$omp master if (rank == 0) then if ((Gflag) .or. (QandGflag)) then write(*,*) 'Using QQGG+G (REAL) matrix elements' @@ -81,6 +82,7 @@ c write(*,*) '[SLC is 1/N]' stop endif endif +!$omp end master endif msq(:,:)=0._dp diff --git a/src/WpWp2j/CMakeLists.txt b/src/WpWp2j/CMakeLists.txt index fe10c06..3ab282c 100644 --- a/src/WpWp2j/CMakeLists.txt +++ b/src/WpWp2j/CMakeLists.txt @@ -1,5 +1,5 @@ target_sources(objlib PRIVATE -nplotter_WPWP.f +#nplotter_WPWP.f qqb_wpwp_qqb.f qqb_wpwp_qqb_g.f qqqqampl.f90 diff --git a/src/Z/AAAREADME.txt b/src/Z/AAAREADME.txt index 01f4f16..5d4782d 100644 --- a/src/Z/AAAREADME.txt +++ b/src/Z/AAAREADME.txt @@ -1,2 +1,3 @@ The purpose of the files in this directory is the NLO calculation of Z production. +qqb_z_loop.f90: 1,2,3 loop corrections to Z production, see arXiv:2207.07056 diff --git a/src/Z2jet/qqb_z2jet_g.f b/src/Z2jet/qqb_z2jet_g.f index a51665a..3188280 100644 --- a/src/Z2jet/qqb_z2jet_g.f +++ b/src/Z2jet/qqb_z2jet_g.f @@ -52,6 +52,7 @@ c --> e^-(p3)+e^+(p4) logical :: filled(1:8) if (first) then first=.false. +!$omp master if (rank == 0) then if ((Gflag) .or. (QandGflag)) then write(*,*) 'Using QQGG+G (REAL) matrix elements' @@ -77,6 +78,7 @@ c write(*,*) '[SLC is 1/N]' stop endif endif +!$omp end master endif msq(:,:)=0._dp diff --git a/src/ptveto/ptint.f b/src/ptveto/ptint.f index b8ca9f9..ba58671 100644 --- a/src/ptveto/ptint.f +++ b/src/ptveto/ptint.f @@ -184,11 +184,7 @@ c pause if (bin) then includeTaucutgrid(0) = .true. - if (newStyleHistograms) then - call nplotter_new(pjet,val) - else - call nplotter(pjet,val,val2,0) - endif + call nplotter_new(pjet,val) endif if (enable_reweight_user) then