Off-shell singletop plotting using new infrastructure

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

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

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

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
 )
 

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