tobias
/
MCFM-10.4
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

Implementation of table plotting

master
Tobias Neumann 3 years ago
parent
9337bfcf6f
commit
3d785371ef
4 changed files with 236 additions and 60 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 133
      src/Mods/mod_Plotting.f90
  2. 29
      src/Mods/mod_SetupPlots.f90
  3. 43
      src/Mods/mod_Superhisto.f90
  4. 91
      src/User/nplotter_Z_new.f90

133
src/Mods/mod_Plotting.f90
View File

@ -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

29
src/Mods/mod_SetupPlots.f90
View File

@ -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

43
src/Mods/mod_Superhisto.f90
View File

@ -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

91
src/User/nplotter_Z_new.f90
View File

@ -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

Write Preview
Loading…
Cancel
Save