diff --git a/src/Mods/mod_Plotting.f90 b/src/Mods/mod_Plotting.f90 index 5ae3424..4639bc2 100644 --- a/src/Mods/mod_Plotting.f90 +++ b/src/Mods/mod_Plotting.f90 @@ -9,6 +9,7 @@ module MCFMPlotting public :: plot_setup_uniform public :: plot_setup_custom + public :: plot_setup_table public :: plots_allocate public :: plot_book @@ -19,6 +20,24 @@ module MCFMPlotting contains + function plot_setup_table(nbins,title) + use types + implicit none + + integer :: plot_setup_table + integer, intent(in) :: nbins + character(len=*), intent(in) :: title + + include 'mpicommon.f' + + if (rank == 0) then + write (*,'(A,A,A)') "Table histogram initialized for '", title, "'" + endif + + plot_setup_table = plot_setup_internal() + call histobuf(plot_setup_table)%init(title,0._dp,real(nbins,dp),1._dp,.true.) + end function + ! called once to initialize individual plots function plot_setup_uniform(xmin,xmax,dx,title) use types @@ -173,13 +192,14 @@ module MCFMPlotting end subroutine ! replaces bookplot - subroutine plot_book(ids, vals, wt0, wts) + subroutine plot_book(ids, vals, wt0, wts, twts) use types use MCFMStorage use Scalevar use PDFerrors use SCET use Integration + use Superhisto, only: TableEntries implicit none include 'kpart.f' include 'nproc.f' @@ -189,10 +209,27 @@ module MCFMPlotting real(dp), intent(in) :: wt0 real(dp), intent(in) :: wts(:) + type(TableEntries), intent(in), optional :: twts(:) + + real(dp), allocatable :: oneArray(:) + real(dp) :: nomTau integer :: i,n,nplots + integer :: numTables + integer, allocatable :: tablesMapping(:) nplots = size(ids) + allocate(tablesMapping(nplots)) + + numTables = 0 + if (present(twts)) then + do i=1,nplots + if (threadStorage(currentPart,currentIps)%histCentral%histos(i)%customTable) then + numTables = numTables + 1 + tablesMapping(i) = numTables + endif + enddo + endif ! for 161 and 165 we always call tmpbook, since these ! modified routines always call threadStorageTmpCommit to properly book everything @@ -203,11 +240,19 @@ module MCFMPlotting if (includeTaucutgrid(currentNd)) then if (kpart == kreal .or. nproc == 1610 .or. nproc == 1650) then do n=1,nplots - call threadStorage(currentPart,currentIps)%histCentral%histos(n)%tmpbook(vals(n),wts(n)) + if (threadStorage(currentPart,currentIps)%histCentral%histos(n)%customTable) then + call threadStorage(currentPart,currentIps)%histCentral%histos(n)%tmpbookTable(twts(tablesMapping(n))%entries) + else + call threadStorage(currentPart,currentIps)%histCentral%histos(n)%tmpbook(vals(n),wts(n)) + endif enddo else do n=1,nplots - call threadStorage(currentPart,currentIps)%histCentral%histos(n)%book(vals(n),wts(n)) + if (threadStorage(currentPart,currentIps)%histCentral%histos(n)%customTable) then + call threadStorage(currentPart,currentIps)%histCentral%histos(n)%bookTable(twts(tablesMapping(n))%entries) + else + call threadStorage(currentPart,currentIps)%histCentral%histos(n)%book(vals(n),wts(n)) + endif enddo endif @@ -215,15 +260,29 @@ module MCFMPlotting if (kpart == kreal .or. nproc == 1610 .or. nproc == 1650) then do i=1,maxscalevar+extrascalevar do n=1,nplots - call threadStorage(currentPart,currentIps)%histScalevar(i)% & - histos(n)%tmpbook(vals(n),wts(n) * (scalereweight(i) - 1._dp)) + if (threadStorage(currentPart,currentIps)%histScalevar(i)%histos(n)%customTable) then + !call threadStorage(currentPart,currentIps)%histScalevar(i)% & + ! histos(n)%tmpbookTable(twts(n) * (scalereweight(i) - 1._dp)) + error stop "to implement: tables with doScalevar" + else + call threadStorage(currentPart,currentIps)%histScalevar(i)% & + histos(n)%tmpbook(vals(n),wts(n) * (scalereweight(i) - 1._dp)) + endif enddo + !if (present(twts)) then + ! call threadStorage(currentPart,currentIps)%histScalevar(i)% & + ! histos(nplots+1)%tmpbookTable(twts * (scalereweight(i) - 1._dp)) + !endif enddo else do i=1,maxscalevar+extrascalevar do n=1,nplots - call threadStorage(currentPart,currentIps)%histScalevar(i)% & - histos(n)%book(vals(n),wts(n) * (scalereweight(i) - 1._dp)) + if (threadStorage(currentPart,currentIps)%histScalevar(i)%histos(n)%customTable) then + error stop "to implement: tables with doScalevar" + else + call threadStorage(currentPart,currentIps)%histScalevar(i)% & + histos(n)%book(vals(n),wts(n) * (scalereweight(i) - 1._dp)) + endif enddo enddo endif @@ -243,6 +302,10 @@ module MCFMPlotting call threadStorage(currentPart,currentIps)%histTaucut(i)% & histos(n)%tmpbook(vals(n), wts(n) * (nomTau - scetreweight(i))) enddo + !if (present(twts)) then + ! call threadStorage(currentPart,currentIps)%histTaucut(i)% & + ! histos(nplots+1)%tmpbookTable(twts * (nomTau - scetreweight(i))) + !endif enddo else do i=1,size(tcutarray) @@ -250,6 +313,10 @@ module MCFMPlotting call threadStorage(currentPart,currentIps)%histTaucut(i)% & histos(n)%book(vals(n), wts(n) * (nomTau - scetreweight(i))) enddo + !if (present(twts)) then + ! call threadStorage(currentPart,currentIps)%histTaucut(i)% & + ! histos(nplots+1)%bookTable(twts * (nomTau - scetreweight(i))) + !endif enddo endif @@ -263,13 +330,22 @@ module MCFMPlotting do i=1,maxPDFsets do n=1,nplots if (wt0 /= 0._dp) then - call threadStorage(currentPart,currentIps)%histPDFerrors(i)% & - histos(n)%tmpbook(vals(n), wts(n)/wt0 * pdfreweight(i)) + call threadStorage(currentPart,currentIps)%histPDFerrors(i)% & + histos(n)%tmpbook(vals(n), wts(n)/wt0 * pdfreweight(i)) else call threadStorage(currentPart,currentIps)%histPDFerrors(i)% & histos(n)%tmpbook(vals(n), pdfreweight(i)) endif enddo + !if (present(twts)) then + ! if (wt0 /= 0._dp) then + ! call threadStorage(currentPart,currentIps)%histPDFerrors(i)% & + ! histos(nplots+1)%tmpbookTable(twts/wt0 * pdfreweight(i)) + ! else + ! call threadStorage(currentPart,currentIps)%histPDFerrors(i)% & + ! histos(nplots+1)%tmpbookTable(oneArray*pdfreweight(i)) + ! endif + !endif enddo endif else @@ -277,13 +353,22 @@ module MCFMPlotting if ((maxPDFsets > 0) .and. kpart == kreal .and. currentPDF > 0 .and. includeTaucutgrid(currentNd)) then do n=1,nplots if (wt0 /= 0._dp) then - call threadStorage(currentPart,currentIps)%histPDFerrors(currentPDF)% & - histos(n)%tmpbook(vals(n), wts(n)/wt0 * pdfreweight(currentPDF)) + call threadStorage(currentPart,currentIps)%histPDFerrors(currentPDF)% & + histos(n)%tmpbook(vals(n), wts(n)/wt0 * pdfreweight(currentPDF)) else call threadStorage(currentPart,currentIps)%histPDFerrors(currentPDF)% & histos(n)%tmpbook(vals(n), pdfreweight(currentPDF)) endif enddo + !if (present(twts)) then + ! if (wt0 /= 0._dp) then + ! call threadStorage(currentPart,currentIps)%histPDFerrors(currentPDF)% & + ! histos(nplots+1)%tmpbookTable(twts/wt0 * pdfreweight(currentPDF)) + ! else + ! call threadStorage(currentPart,currentIps)%histPDFerrors(currentPDF)% & + ! histos(nplots+1)%tmpbookTable(oneArray * pdfreweight(currentPDF)) + ! endif + !endif endif endif @@ -293,25 +378,43 @@ module MCFMPlotting do i=1,maxPDFsets do n=1,nplots if (wt0 /= 0._dp) then - call threadStorage(currentPart,currentIps)%histPDFerrors(i)% & - histos(n)%tmpbook(vals(n), wts(n)/wt0 * pdfreweight(i)) + call threadStorage(currentPart,currentIps)%histPDFerrors(i)% & + histos(n)%tmpbook(vals(n), wts(n)/wt0 * pdfreweight(i)) else call threadStorage(currentPart,currentIps)%histPDFerrors(i)% & histos(n)%tmpbook(vals(n), pdfreweight(i)) endif enddo + !if (present(twts)) then + ! if (wt0 /= 0._dp) then + ! call threadStorage(currentPart,currentIps)%histPDFerrors(i)% & + ! histos(nplots+1)%tmpbookTable(twts/wt0 * pdfreweight(i)) + ! else + ! call threadStorage(currentPart,currentIps)%histPDFerrors(i)% & + ! histos(nplots+1)%tmpbookTable(oneArray * pdfreweight(i)) + ! endif + !endif enddo else do i=1,maxPDFsets do n=1,nplots if (wt0 /= 0._dp) then - call threadStorage(currentPart,currentIps)%histPDFerrors(i)% & - histos(n)%book(vals(n), wts(n)/wt0 * pdfreweight(i)) + call threadStorage(currentPart,currentIps)%histPDFerrors(i)% & + histos(n)%book(vals(n), wts(n)/wt0 * pdfreweight(i)) else call threadStorage(currentPart,currentIps)%histPDFerrors(i)% & histos(n)%book(vals(n), pdfreweight(i)) endif enddo + !if (present(twts)) then + ! if (wt0 /= 0._dp) then + ! call threadStorage(currentPart,currentIps)%histPDFerrors(i)% & + ! histos(nplots+1)%bookTable(twts/wt0 * pdfreweight(i)) + ! else + ! call threadStorage(currentPart,currentIps)%histPDFerrors(i)% & + ! histos(nplots+1)%bookTable(oneArray * pdfreweight(i)) + ! endif + !endif enddo endif endif diff --git a/src/Mods/mod_SetupPlots.f90 b/src/Mods/mod_SetupPlots.f90 index 8e9e4e6..9f24413 100644 --- a/src/Mods/mod_SetupPlots.f90 +++ b/src/Mods/mod_SetupPlots.f90 @@ -4,11 +4,11 @@ ! module MCFMSetupPlots + use types implicit none public :: setup_plots public :: nplotter_new - private interface @@ -21,16 +21,27 @@ module MCFMSetupPlots real(dp), allocatable, intent(out) :: vals(:) real(dp), allocatable, intent(out) :: wts(:) end subroutine + + subroutine tbook(p,wt,wts) + use types + use Superhisto, only: TableEntries + implicit none + include 'mxpart.f' + real(dp), intent(in) :: p(mxpart,4), wt + type(TableEntries), allocatable, intent(out) :: wts(:) + end subroutine end interface procedure (ibook), pointer :: pbook => null() + + procedure (tbook), pointer :: ptbook => null() logical, save :: disableHistograms = .false. contains subroutine setup_plots() - use nplotter_Z, only: setup_Z => setup, book_Z => book + use nplotter_Z, only: setup_Z => setup, book_Z => book, bookTable_Z => bookTable use nplotter_W, only: setup_W => setup, book_W => book use nplotter_singletop, only: setup_singletop => setup, book_singletop => book use nplotter_Higgs, only: setup_Higgs => setup, book_Higgs => book @@ -70,6 +81,7 @@ module MCFMSetupPlots elseif (kcase==kZ_only .or. kcase==kgg2lep .or. kcase == kZ_2jet .or. kcase == kZ_1jet) then call setup_Z() pbook => book_Z + ptbook => bookTable_Z elseif (kcase == kbq_tpq .or. kcase == kbq_tpq_jet) then call setup_singletop() pbook => book_singletop @@ -117,6 +129,7 @@ module MCFMSetupPlots subroutine nplotter_new(p,wt) use types use MCFMPlotting + use Superhisto, only: TableEntries implicit none include 'mxpart.f' real(dp), intent(in) :: p(mxpart,4), wt @@ -125,6 +138,8 @@ module MCFMSetupPlots real(dp), allocatable :: vals(:) real(dp), allocatable :: wts(:) + type(TableEntries), allocatable :: twts(:) + ! could call plotting for auto histograms here ! call process specific plotting @@ -137,7 +152,15 @@ module MCFMSetupPlots if (.not. allocated(wts)) then allocate(wts(1+size(ids)), source=wt) endif - call plot_book([1,ids],[0.5_dp,vals],wt,[wt,wts]) + + if (associated(ptbook)) then + ! we use a table here + call ptbook(p,wt,twts) + call plot_book([1,ids],[0.5_dp,vals],wt,[wt,wts],twts) + else + call plot_book([1,ids],[0.5_dp,vals],wt,[wt,wts]) + endif + endif end subroutine diff --git a/src/Mods/mod_Superhisto.f90 b/src/Mods/mod_Superhisto.f90 index 287759f..9d7fc76 100644 --- a/src/Mods/mod_Superhisto.f90 +++ b/src/Mods/mod_Superhisto.f90 @@ -8,6 +8,12 @@ module Superhisto use types implicit none + type TableEntries + real(dp), allocatable :: entries(:) + end type TableEntries + + public :: TableEntries + private type sh_histogram @@ -27,13 +33,17 @@ module Superhisto logical :: customBinning real(dp), allocatable :: xarr(:) + logical :: customTable + contains procedure :: init_custom => shinit_custom procedure :: init => shinit procedure :: initialized => shinitialized procedure :: book => shbook + procedure :: bookTable => shbookTable procedure :: tmpbook => shtmpbook + procedure :: tmpbookTable => shtmpbookTable procedure :: tmpcommit => shtmpcommit procedure :: reset => shreset procedure :: tmpreset => shtmpreset @@ -46,6 +56,7 @@ module Superhisto public :: sh_histogram public :: shbook + public :: shbookTable public :: shreset public :: shinit public :: shinitialized @@ -205,6 +216,7 @@ module Superhisto newsh%xmax = xarr(n) newsh%delx = 0._dp newsh%customBinning = .true. + newsh%customTable = .false. newsh%xarr = xarr newsh%nbins = n-1 @@ -222,24 +234,34 @@ module Superhisto end subroutine shinit_custom - subroutine shinit(newsh, title, xmin, xmax, delx) + subroutine shinit(newsh, title, xmin, xmax, delx, cTable_in) implicit none class(sh_histogram), intent(inout) :: newsh character(len=*), intent(in) :: title real(dp), intent(in) :: xmin, xmax, delx integer :: nbins + logical, optional, intent(in) :: cTable_in + logical :: cTable + + if (present(cTable_in)) then + cTable = cTable_in + else + cTable = .false. + endif + newsh%title = title newsh%xmin = xmin newsh%xmax = xmax newsh%delx = delx newsh%customBinning = .false. + newsh%customTable = cTable nbins = ceiling((xmax-xmin)/delx) newsh%nbins = nbins !write (*,*) "initializing histo ", title, xmin, xmax, delx, nbins, & - !"for thread", omp_get_thread_num() + ! "for thread", omp_get_thread_num() if (allocated(newsh%xx) .eqv. .false.) then allocate(newsh%xx(0:nbins+1), source=0._dp) @@ -253,6 +275,23 @@ module Superhisto end subroutine + subroutine shbookTable(hist, wgts) + use ieee_arithmetic + class(sh_histogram), intent(inout) :: hist + real(dp), intent(in) :: wgts(:) + + hist%xx(1:hist%nbins) = hist%xx(1:hist%nbins) + wgts + hist%xxsq(1:hist%nbins) = hist%xxsq(1:hist%nbins) + wgts**2 + end subroutine + + subroutine shtmpbookTable(hist, wgts) + use ieee_arithmetic + class(sh_histogram), intent(inout) :: hist + real(dp), intent(in) :: wgts(:) + + hist%tmp(1:hist%nbins) = hist%tmp(1:hist%nbins) + wgts + end subroutine + subroutine shbook(hist, x, wgt) use ieee_arithmetic use omp_lib diff --git a/src/User/nplotter_Z_new.f90 b/src/User/nplotter_Z_new.f90 index 8cd27db..281c23a 100644 --- a/src/User/nplotter_Z_new.f90 +++ b/src/User/nplotter_Z_new.f90 @@ -10,7 +10,7 @@ module nplotter_Z use qtResummation_params, only: transitionSwitch implicit none - public :: setup, book + public :: setup, book, bookTable private integer, save, allocatable :: histos(:) @@ -23,12 +23,7 @@ module nplotter_Z implicit none include 'mpicommon.f' - include 'first.f' - if (first .and. rank == 0) then - write(6,*) 'Using plotting routine nplotter_Z_new.f90' - first=.false. - endif - allocate(histos(5)) + allocate(histos(8)) if (rank == 0) then write (*,*) "RESUMMATION: Using transition with switch ", transitionSwitch @@ -43,32 +38,31 @@ module nplotter_Z 2.5119d0,3.1623d0,3.9811d0,5.0119d0,6.3096d0,7.9433d0, & 10.0000d0,12.5893d0,15.8489d0,19.9526d0,25.1189d0, & 31.6228d0,39.8107d0,50.1187d0,63.0957d0,79.4328d0,100.0000d0], & - 'pt34_fine_notrans') - - histos(2) = plot_setup_custom([0d0,2d0,4d0,6d0,8d0, & + 'pt34_fine') + histos(2) = plot_setup_uniform(0.0_dp,60._dp,1.0_dp,'pt34') + histos(3) = plot_setup_custom([0d0,2d0,4d0,6d0,8d0, & 10d0,12d0,14d0,16d0,18d0,20d0,22.5d0,25d0,27.5d0,30d0, & 33d0,36d0,39d0,42d0,45d0,48d0,51d0,54d0,57d0,61d0,65d0, & 70d0,75d0,80d0,85d0,95d0,105d0,125d0,150d0,175d0,200d0, & - 250d0,300d0,350d0,400d0,470d0,550d0,650d0,900d0],'pt34_big_trans04') - histos(3) = plot_setup_custom([0d0,0.004d0,0.008d0,0.012d0, & + 250d0,300d0,350d0,400d0,470d0,550d0,650d0,900d0],'pt34_atlas') + histos(4) = plot_setup_custom([0d0,0.004d0,0.008d0,0.012d0, & 0.016d0,0.02d0,0.024d0,0.029d0,0.034d0,0.039d0,0.045d0, & 0.051d0,0.057d0,0.064d0,0.072d0,0.081d0,0.091d0,0.102d0, & 0.114d0,0.128d0,0.145d0,0.165d0,0.189d0,0.219d0,0.258d0, & 0.312d0,0.391d0,0.524d0,0.695d0,0.918d0,1.153d0,1.496d0, & - 1.947d0,2.522d0,3.277d0,5d0,10d0],'phistar_trans04') + 1.947d0,2.522d0,3.277d0,5d0,10d0],'phistar_atlas') - histos(4) = plot_setup_custom([0d0,2d0,4d0,6d0,8d0, & - 10d0,12d0,14d0,16d0,18d0,20d0,22.5d0,25d0,27.5d0,30d0, & - 33d0,36d0,39d0,42d0,45d0,48d0,51d0,54d0,57d0,61d0,65d0, & - 70d0,75d0,80d0,85d0,95d0,105d0,125d0,150d0,175d0,200d0, & - 250d0,300d0,350d0,400d0,470d0,550d0,650d0,900d0],'pt34_big_trans06') - histos(5) = plot_setup_custom([0d0,0.004d0,0.008d0,0.012d0, & - 0.016d0,0.02d0,0.024d0,0.029d0,0.034d0,0.039d0,0.045d0, & - 0.051d0,0.057d0,0.064d0,0.072d0,0.081d0,0.091d0,0.102d0, & - 0.114d0,0.128d0,0.145d0,0.165d0,0.189d0,0.219d0,0.258d0, & - 0.312d0,0.391d0,0.524d0,0.695d0,0.918d0,1.153d0,1.496d0, & - 1.947d0,2.522d0,3.277d0,5d0,10d0],'phistar_trans06') + histos(5) = plot_setup_custom([0d0,1d0,2d0,3d0,4d0,5d0,6d0,7d0, & + 8d0,9d0,10d0,11d0,12d0,13d0,14d0,16d0,18d0,20d0,22d0,25d0, & + 28d0,32d0,37d0,43d0,52d0,65d0,85d0,120d0,160d0,190d0,220d0, & + 250d0,300d0,400d0,500d0,800d0,1650d0],'pt34_cms13') + histos(6) = plot_setup_custom([0d0,2.5d0,5d0,7.5d0,10d0,12.5d0, & + 15d0,17.5d0,20d0,30d0,40d0,50d0,70d0,90d0,110d0,150d0,190d0, & + 250d0, 600d0], 'pt34_cms7') + + histos(7) = plot_setup_table(5,'table_for_observable_one') + histos(8) = plot_setup_table(10,'another_table') end subroutine @@ -89,7 +83,7 @@ module nplotter_Z real(dp), allocatable, intent(out) :: wts(:) real(dp) :: pttwo, twomass, delphi, etarap - real(dp) :: pt34, trans04, trans06 + real(dp) :: pt34, trans real(dp) :: phistar, phiacop, costhetastar, delphi34 pt34 = pttwo(3,4,p) @@ -99,21 +93,15 @@ module nplotter_Z costhetastar = tanh((etarap(3,p)-etarap(4,p))/2._dp) phistar = tan(phiacop/2._dp)*sin(acos(costhetastar)) - ! the variable transitionSwitch is taken from the input file and can be used here - ! instead of the hardcoded 0.4 and 0.6 - if (origKpart == kresummed) then if (abovecut .eqv. .false.) then - trans04 = transition((pt34/twomass(3,4,p))**2d0,0.001d0, 0.4d0 ,0.001d0) - trans06 = transition((pt34/twomass(3,4,p))**2d0,0.001d0, 0.6d0 ,0.001d0) + trans = transition((pt34/twomass(3,4,p))**2d0,0.001d0,transitionSwitch,0.001d0) else ! fo piece without transition - trans04 = 1._dp - trans06 = 1._dp + trans = 1._dp endif else - trans04 = 1._dp - trans06 = 1._dp + trans = 1._dp endif if (ieee_is_nan(pt34)) then @@ -124,14 +112,37 @@ module nplotter_Z phistar = -1._dp endif - ! fill histograms: first without transition function - ! then with 0.4 transition function, then with 0.6 transition function - ! for estimating matching uncertainty - ids = histos - vals = [pt34, pt34,phistar, pt34,phistar] - wts = [wt, wt*trans04,wt*trans06, wt*trans04,wt*trans06] + vals = [pt34,pt34,pt34,phistar,pt34,pt34] + wts = [wt*trans,wt*trans,wt*trans,wt*trans,wt*trans,wt*trans] + + end subroutine + + subroutine bookTable(p,wt,wts) + use types + use Superhisto, only: TableEntries + implicit none + include 'mxpart.f' + real(dp), intent(in) :: p(mxpart,4), wt + type(TableEntries), allocatable, intent(out) :: wts(:) + + integer :: i + + ! calculate coefficients here based on p + + ! two tables, allocate more for additional tables + allocate(wts(2)) + + ! allocate five entries for the first table + allocate(wts(1)%entries(5)) + + ! allocate ten entries for the first table + allocate(wts(2)%entries(10)) + ! return array of five coefficients for first table + wts(1)%entries = [1._dp, 0.5_dp, 2._dp, wt, 1._dp] + ! ten for second + wts(2)%entries = [(wt, i=1,10)] end subroutine end module