From 03f633b536c117a6a1b5caa6d08e67da59c5fcfe Mon Sep 17 00:00:00 2001 From: Tobias Neumann Date: Thu, 7 Mar 2024 20:51:18 +0000 Subject: [PATCH] Off-shell singletop plotting using new infrastructure --- src/Mods/bbfrac.f | 3 + src/Mods/mod_SetupPlots.f90 | 4 + src/Procdep/realint.f | 14 +- src/User/CMakeLists.txt | 1 + src/User/nplotter_ktopanom_new.f90 | 199 +++++++++++++++++++++++++++++ 5 files changed, 219 insertions(+), 2 deletions(-) create mode 100644 src/User/nplotter_ktopanom_new.f90 diff --git a/src/Mods/bbfrac.f b/src/Mods/bbfrac.f index 32e2f17..ea17487 100644 --- a/src/Mods/bbfrac.f +++ b/src/Mods/bbfrac.f @@ -14,4 +14,7 @@ real(dp), save, public :: bbfrac(0:maxd) !$omp threadprivate(bbfrac) + logical, save, public :: bbpart = .false. +!$omp threadprivate(bbpart) + end module diff --git a/src/Mods/mod_SetupPlots.f90 b/src/Mods/mod_SetupPlots.f90 index 9f24413..3ffb1e2 100644 --- a/src/Mods/mod_SetupPlots.f90 +++ b/src/Mods/mod_SetupPlots.f90 @@ -58,6 +58,7 @@ module MCFMSetupPlots use nplotter_ZH_bbar_tautau_gamgam, only: setup_ZH_bbar_tautau_gamgam => setup, & book_ZH_bbar_tautau_gamgam => book use nplotter_ZH_WW, only: setup_ZH_WW => setup, book_ZH_WW => book + use nplotter_ktopanom, only: setup_ktopanom => setup, book_ktopanom => book use MCFMPlotting use parseinput, only: cfg, cfg_get use types @@ -120,6 +121,9 @@ module MCFMSetupPlots elseif (kcase == kZHbbar .or. kcase == kZHgaga .or. kcase == kZH1jet) then call setup_ZH_bbar_tautau_gamgam() pbook => book_ZH_bbar_tautau_gamgam + elseif (kcase == ktopanom) then + call setup_ktopanom() + pbook => book_ktopanom endif ! could also add auto-histograms here diff --git a/src/Procdep/realint.f b/src/Procdep/realint.f index e4dd363..7ded2d2 100644 --- a/src/Procdep/realint.f +++ b/src/Procdep/realint.f @@ -9,7 +9,7 @@ use Scalevar use ieee_arithmetic use VVconfig_m - use bbfrac_m, only : bbfrac + use bbfrac_m, only : bbfrac, bbpart use MCFMStorage use SCET use m_gencuts, only : enable_reweight_user, reweight_user @@ -710,6 +710,7 @@ c call singcheck(qqb_QQb_g,qqb_QQb_gs,p) call qqb_QQb_gs(p,msqc) elseif (kcase==ktopanom) then bbfrac = 0._dp + bbpart = .false. msq_light = 0._dp msq_heavy = 0._dp @@ -1251,6 +1252,7 @@ c write (*,*) j,k,singletop2_dipole_pdfs(nd,1,j),singletop2_dipol c if (nd == 13 .or. nd == 14) the if ((j == 0 .and. abs(k) /= 5) .or. (abs(j) /= 5 .and. k ==0)) then bbfrac(nd) = bbfrac(nd) + xmsqjk + bbpart = .true. endif c endif endif @@ -1360,7 +1362,15 @@ c endif currentNd = nd call getptildejet(nd,pjet) call dotem(nvec,pjet,s) - call nplotter_new(pjet,xmsq(nd)*wgt) + if (bbpart) then + ! separately bin bb part and other part + bbpart = .true. + call nplotter_new(pjet, xmsq(nd)*wgt*bbfrac(nd)) + bbpart = .false. + call nplotter_new(pjet, xmsq(nd)*wgt*(1-bbfrac(nd))) + else + call nplotter_new(pjet,xmsq(nd)*wgt) + endif endif endif enddo ! loop over dipole contributions diff --git a/src/User/CMakeLists.txt b/src/User/CMakeLists.txt index 05ac505..3ac1366 100644 --- a/src/User/CMakeLists.txt +++ b/src/User/CMakeLists.txt @@ -58,6 +58,7 @@ nplotter_WW_new.f90 nplotter_WZ_new.f90 nplotter_ZZ_new.f90 nplotter_singletop_new.f90 +nplotter_ktopanom_new.f90 transition.f90 ) diff --git a/src/User/nplotter_ktopanom_new.f90 b/src/User/nplotter_ktopanom_new.f90 new file mode 100644 index 0000000..463e87f --- /dev/null +++ b/src/User/nplotter_ktopanom_new.f90 @@ -0,0 +1,199 @@ +! +! SPDX-License-Identifier: GPL-3.0-or-later +! Copyright (C) 2019-2022, respective authors of MCFM. +! + +module nplotter_ktopanom + use types + use MCFMPlotting + + implicit none + + public :: setup, book + private + + integer, save, allocatable :: histos(:) + + contains + + subroutine setup() + use types + use parseinput + implicit none + include 'mpicommon.f' + + allocate(histos(1)) + + histos(1) = plot_setup_uniform(0.0_dp, 200.0_dp, 5.0_dp, 'mbl') + + end subroutine + + subroutine book(p,wt,ids,vals,wts) + use types + use ResummationTransition + use ieee_arithmetic + use jettagging + use bbfrac_m, only: bbpart + use qsort_m + use mathfun + use SCET, only: currentNd + implicit none + include 'mxpart.f' + include 'kpart.f' + include 'jetlabel.f' + include 'plabel.f' + + real(dp), intent(in) :: p(mxpart,4) + real(dp), intent(in) :: wt + + integer, allocatable, intent(out) :: ids(:) + real(dp), allocatable, intent(out) :: vals(:) + real(dp), allocatable, intent(out) :: wts(:) + + integer :: jetindices(3) + + integer :: bjetindices(3) + integer :: nonbjetindices(3) + + integer :: bjetindices_bbfrac(3) + integer :: nonbjetindices_bbfrac(3) + + real(dp) :: plep(4), pnu(4) + + integer :: i,n + + integer :: imax, imaxb + + real(dp) :: puremass + real(dp) :: mbl + + jetindices(:) = [5,6,7] + + ! the 'bbpart' variable (logical, true or false) designates whether we + ! have two b-quarks in the final state or not + + bjetindices = 0 + nonbjetindices = 0 + + ! for main part with at most one b + if (bbpart .eqv. .false.) then + imax = 1 + imaxb = 1 + do i=1,jets + if (iand(bitflag(5),jetcontent(jetindices(i),currentNd)) > 0) then + bjetindices(imaxb) = jetindices(i) + imaxb = imaxb + 1 + else + nonbjetindices(imax) = jetindices(i) + imax = imax + 1 + endif + enddo + else + imax = 1 + imaxb = 1 + do i=1,jets + if ((iand(bitflag(5),jetcontent(jetindices(i),currentNd)) > 0) .or. & + (iand(bitflag(7),jetcontent(jetindices(i),currentNd)) > 0)) then + bjetindices(imaxb) = jetindices(i) + imaxb = imaxb + 1 + else + nonbjetindices(imax) = jetindices(i) + imax = imax + 1 + endif + enddo + + endif + + ! require at least one b-jet and one light jet + if ( bjetindices(1) == 0 .or. nonbjetindices(1) == 0 ) then + ids = histos + vals = [0._dp] + wts = [0._dp] + return + endif + + if (plabel(3) == 'el') then + plep(:) = p(3,:) + pnu(:) = p(4,:) + else + plep(:) = p(4,:) + pnu(:) = p(3,:) + endif + + mbl = puremass(p(bjetindices(1),:) + plep(:)) + + ids = histos + vals = [mbl] + wts = [wt] + + end subroutine book + + + subroutine wrecons(pl_in, pnu_in, pnu_out) + use types + implicit none + include 'masses.f' + + real(dp), intent(in) :: pl_in(4) + real(dp), intent(in) :: pnu_in(4) + real(dp), intent(out) :: pnu_out(4) + + integer, parameter :: ielectron = 1 + integer, parameter :: ineutrino = 2 + + real(dp) :: ps(0:3, 2) + real(dp) :: pn(0:3) + + real(dp) :: a,b,c,d, fix + logical :: good + + integer :: i,j + + ps(0:3,1) = [pl_in(4), pl_in(1), pl_in(2), pl_in(3)] + ps(0:3,2) = [pnu_in(4), pnu_in(1), pnu_in(2), pnu_in(3)] + + pn(:) = ps(:,2) + + + ! from MG 1 + + i = ielectron + j = ineutrino + a=wmass*wmass*.5d0+ps(1,i)*ps(1,j)+ps(2,i)*ps(2,j) + b = ps(0,i) + c = sqrt(ps(1,j)**2+ps(2,j)**2) + d = ps(3,i) + good =(a**2 - b**2*c**2 + c**2*d**2 > 0d0) + fix=1d0 + do while (.not. good) + fix=fix*0.9d0 + ps(1,j)=ps(1,j)*.9d0 + ps(2,j)=ps(2,j)*.9d0 + a = wmass*wmass*.5d0+ps(1,i)*ps(1,j)+ps(2,i)*ps(2,j) + b = ps(0,i) + c = sqrt(ps(1,j)**2+ps(2,j)**2) + d = ps(3,i) + good =(a**2 - b**2*c**2 + c**2*d**2 > 0d0) + enddo + ps(3,j) = (2d0*a*d/(b**2 - d**2) + & + 2d0*b*Sqrt(max(a**2 - b**2*c**2 + c**2*d**2,0d0))/ & + ((b - d)*(b + d)))/2d0 + + pn(3) = (2d0*a*d/(b**2 - d**2) - & + 2d0*b*Sqrt(max(a**2 - b**2*c**2 + c**2*d**2,0d0))/ & + ((b - d)*(b + d)))/2d0 + if (abs(ps(3,j)) > abs(pn(3))) then + a=ps(3,j) + ps(3,j)=pn(3) + pn(3)=a + endif + ps(0,j) = sqrt(ps(1,j)**2+ps(2,j)**2+ps(3,j)**2) + pn(0) = sqrt(pn(1)**2+pn(2)**2+pn(3)**2) + + pnu_out(4) = ps(0,j) + pnu_out(1:3) = ps(1:3,j) + + + end subroutine + +end module