Index: trunk/src/matching/Makefile.am =================================================================== --- trunk/src/matching/Makefile.am (revision 8803) +++ trunk/src/matching/Makefile.am (revision 8804) @@ -1,254 +1,273 @@ ## Makefile.am -- Makefile for WHIZARD ## ## Process this file with automake to produce Makefile.in # # Copyright (C) 1999-2022 by # Wolfgang Kilian # Thorsten Ohl # Juergen Reuter # with contributions from # cf. main AUTHORS file # # WHIZARD is free software; you can redistribute it and/or modify it # under the terms of the GNU General Public License as published by # the Free Software Foundation; either version 2, or (at your option) # any later version. # # WHIZARD is distributed in the hope that it will be useful, but # WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # # You should have received a copy of the GNU General Public License # along with this program; if not, write to the Free Software # Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. # ######################################################################## ## The files in this directory implement different matchings ## We create a library which is still to be combined with auxiliary libs. noinst_LTLIBRARIES = libmatching.la check_LTLIBRARIES = libmatching_ut.la COMMON_F90 = \ matching_base.f90 \ + mlm_matching.f90 \ ckkw_matching.f90 \ - mlm_matching.f90 + powheg_matching.f90 MPI_F90 = \ - powheg_matching.f90_mpi + powheg_matching_sub.f90_mpi SERIAL_F90 = \ - powheg_matching.f90_serial + powheg_matching_sub.f90_serial +MATCHING_SUBMODULES = \ + matching_base_sub.f90 \ + mlm_matching_sub.f90 \ + ckkw_matching_sub.f90 + +MATCHING_MODULES = \ + $(COMMON_F90) EXTRA_DIST = \ $(COMMON_F90) \ + $(MATCHING_SUBMODULES) \ $(SERIAL_F90) \ $(MPI_F90) nodist_libmatching_la_SOURCES = \ - $(COMMON_F90) \ - powheg_matching.f90 + $(MATCHING_MODULES) \ + $(MATCHING_SUBMODULES) \ + powheg_matching_sub.f90 -DISTCLEANFILES = powheg_matching.f90 +DISTCLEANFILES = powheg_matching_sub.f90 if FC_USE_MPI -powheg_matching.f90: powheg_matching.f90_mpi +powheg_matching_sub.f90: powheg_matching_sub.f90_mpi -cp -f $< $@ else -powheg_matching.f90: powheg_matching.f90_serial +powheg_matching_sub.f90: powheg_matching_sub.f90_serial -cp -f $< $@ endif libmatching_ut_la_SOURCES = \ powheg_matching_uti.f90 powheg_matching_ut.f90 ## Omitting this would exclude it from the distribution dist_noinst_DATA = matching.nw # Modules and installation # Dump module names into file Modules execmoddir = $(fmoddir)/whizard nodist_execmod_HEADERS = \ - ${nodist_libmatching_la_SOURCES:.f90=.$(FCMOD)} + ${MATCHING_MODULES:.f90=.$(FCMOD)} # Dump module names into file Modules libmatching_Modules = \ - $(nodist_libmatching_la_SOURCES:.f90=) \ + $(MATCHING_MODULES:.f90=) \ $(libmatching_ut_la_SOURCES:.f90=) Modules: Makefile @for module in $(libmatching_Modules); do \ echo $$module >> $@.new; \ done @if diff $@ $@.new -q >/dev/null; then \ rm $@.new; \ else \ mv $@.new $@; echo "Modules updated"; \ fi BUILT_SOURCES = Modules ## Fortran module dependencies # Get module lists from other directories module_lists = \ ../basics/Modules \ ../utilities/Modules \ ../testing/Modules \ ../system/Modules \ ../combinatorics/Modules \ ../parsing/Modules \ ../rng/Modules \ ../physics/Modules \ ../qft/Modules \ ../expr_base/Modules \ ../types/Modules \ ../particles/Modules \ ../beams/Modules \ ../matrix_elements/Modules \ ../me_methods/Modules \ ../events/Modules \ ../muli/Modules \ ../phase_space/Modules \ ../mci/Modules \ ../blha/Modules \ ../gosam/Modules \ ../openloops/Modules \ ../recola/Modules \ ../fks/Modules \ ../model_features/Modules \ ../variables/Modules \ ../shower/Modules \ ../threshold/Modules \ ../process_integration/Modules $(module_lists): $(MAKE) -C `dirname $@` Modules Module_dependencies.sed: $(nodist_libmatching_la_SOURCES) $(libmatching_ut_la_SOURCES) Module_dependencies.sed: $(module_lists) @rm -f $@ echo 's/, *only:.*//' >> $@ echo 's/, *&//' >> $@ echo 's/, *.*=>.*//' >> $@ echo 's/$$/.lo/' >> $@ for list in $(module_lists); do \ dir="`dirname $$list`"; \ for mod in `cat $$list`; do \ echo 's!: '$$mod'.lo$$!': $$dir/$$mod'.lo!' >> $@; \ done \ done DISTCLEANFILES += Module_dependencies.sed # The following line just says # include Makefile.depend # but in a portable fashion (depending on automake's AM_MAKE_INCLUDE @am__include@ @am__quote@Makefile.depend@am__quote@ Makefile.depend: Module_dependencies.sed Makefile.depend: $(nodist_libmatching_la_SOURCES) $(libmatching_ut_la_SOURCES) @rm -f $@ for src in $^; do \ module="`basename $$src | sed 's/\.f[90][0358]//'`"; \ grep '^ *use ' $$src \ | grep -v '!NODEP!' \ | sed -e 's/^ *use */'$$module'.lo: /' \ -f Module_dependencies.sed; \ done > $@ DISTCLEANFILES += Makefile.depend # Fortran90 module files are generated at the same time as object files .lo.$(FCMOD): @: # touch $@ AM_FCFLAGS = -I../basics -I../utilities -I../testing -I../system -I../combinatorics -I../parsing -I../rng -I../physics -I../qft -I../expr_base -I../types -I../particles -I../beams -I../matrix_elements -I../me_methods -I../events -I../muli -I../phase_space -I../mci -I../blha -I../gosam -I../openloops -I../recola -I../fks -I../model_features -I../variables -I../tauola -I../shower -I../threshold -I../process_integration -I../../circe1/src -I../../circe2/src -I../pdf_builtin -I../lhapdf -I../qed_pdf -I../fastjet -I../../vamp/src ######################################################################## +matching_base_sub.lo: matching_base.lo +mlm_matching_sub.lo: mlm_matching.lo +ckkw_matching_sub.lo: ckkw_matching.lo +powheg_matching_sub.lo: powheg_matching.lo + +######################################################################## ## Default Fortran compiler options ## Profiling if FC_USE_PROFILING AM_FCFLAGS += $(FCFLAGS_PROFILING) endif ## OpenMP if FC_USE_OPENMP AM_FCFLAGS += $(FCFLAGS_OPENMP) endif # MPI if FC_USE_MPI AM_FCFLAGS += $(FCFLAGS_MPI) endif if RECOLA_AVAILABLE AM_FCFLAGS += $(RECOLA_INCLUDES) endif ######################################################################## ## Non-standard targets and dependencies ## (Re)create F90 sources from NOWEB source. if NOWEB_AVAILABLE FILTER = -filter "sed 's/defn MPI:/defn/'" PRELUDE = $(top_srcdir)/src/noweb-frame/whizard-prelude.nw POSTLUDE = $(top_srcdir)/src/noweb-frame/whizard-postlude.nw matching.stamp: $(PRELUDE) $(srcdir)/matching.nw $(POSTLUDE) @rm -f matching.tmp @touch matching.tmp for src in \ $(COMMON_F90) $(libmatching_ut_la_SOURCES); do \ $(NOTANGLE) -R[[$$src]] $^ | $(CPIF) $$src; \ done + for src in $(MATCHING_SUBMODULES); do \ + $(NOTANGLE) -R[[$$src]] $^ | $(CPIF) $$src; \ + done for src in $(SERIAL_F90:.f90_serial=.f90); do \ $(NOTANGLE) -R[[$$src]] $^ | $(CPIF) $$src'_serial'; \ done for src in $(MPI_F90:.f90_mpi=.f90); do \ $(NOTANGLE) -R[[$$src]] $(FILTER) $^ | $(CPIF) $$src'_mpi'; \ done @mv -f matching.tmp matching.stamp -$(COMMON_F90) $(SERIAL_F90) $(MPI_F90) $(libmatching_ut_la_SOURCES): matching.stamp +$(COMMON_F90) $(MATCHING_SUBMODULES) $(SERIAL_F90) $(MPI_F90) $(libmatching_ut_la_SOURCES): matching.stamp ## Recover from the removal of $@ @if test -f $@; then :; else \ rm -f matching.stamp; \ $(MAKE) $(AM_MAKEFLAGS) matching.stamp; \ fi endif ######################################################################## ## Non-standard cleanup tasks ## Remove sources that can be recreated using NOWEB if NOWEB_AVAILABLE maintainer-clean-noweb: -rm -f *.f90 *.f90_serial *.f90_mpi *.c endif .PHONY: maintainer-clean-noweb ## Remove those sources also if builddir and srcdir are different if NOWEB_AVAILABLE clean-noweb: test "$(srcdir)" != "." && rm -f *.f90 *.f90_serial *.f90_mpi *.c || true endif .PHONY: clean-noweb ## Remove F90 module files clean-local: clean-noweb -rm -f matching.stamp matching.tmp -rm -f *.$(FCMOD) if FC_SUBMODULES - -rm -f *.smod + -rm -f *.smod *.sub endif ## Remove backup files maintainer-clean-backup: -rm -f *~ .PHONY: maintainer-clean-backup ## Register additional clean targets maintainer-clean-local: maintainer-clean-noweb maintainer-clean-backup Index: trunk/src/matching/matching.nw =================================================================== --- trunk/src/matching/matching.nw (revision 8803) +++ trunk/src/matching/matching.nw (revision 8804) @@ -1,4727 +1,5547 @@ % -*- ess-noweb-default-code-mode: f90-mode; noweb-default-code-mode: f90-mode; -*- % WHIZARD code as NOWEB source: Matching and Merging %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{Matching} \includemodulegraph{matching} <<[[matching_base.f90]]>>= <> module matching_base <> -<> - use diagnostics use sm_qcd use model_data use particles use variables use shower_base use instances, only: process_instance_t use rng_base <> <> <> <> <> + interface +<> + end interface + +end module matching_base +@ %def matching_base +@ +<<[[matching_base_sub.f90]]>>= +<> + +submodule (matching_base) matching_base_s + +<> + use diagnostics + + implicit none + contains <> -end module matching_base -@ %def matching_base +end submodule matching_base_s + +@ %def matching_base_s @ \section{Abstract Matching Type} A matching will need access to the [[shower]] as well as matrix elements that we currently get over [[process_instace]]. The [[model]] is intended for the backup [[model_hadrons]]. <>= public :: matching_t <>= type, abstract :: matching_t logical :: is_hadron_collision = .false. type(qcd_t) :: qcd class(shower_base_t), pointer :: shower => null () type(process_instance_t), pointer :: process_instance => null () class(model_data_t), pointer :: model => null () class(rng_t), allocatable :: rng type(string_t) :: process_name contains <> end type matching_t @ %def matching_t @ <>= procedure (matching_init), deferred :: init <>= abstract interface subroutine matching_init (matching, var_list, process_name) import class(matching_t), intent(out) :: matching type(var_list_t), intent(in) :: var_list type(string_t), intent(in) :: process_name end subroutine matching_init end interface @ %def matching_init @ If we use a polymorphic settings type, this boilerplate wouldn't be necessary but then we introduce [[select type]] statements all over the place. <>= type(var_list_t), intent(in) :: var_list type(string_t), intent(in) :: process_name if (debug_on) call msg_debug (D_MATCHING, "matching_init") call matching%settings%init (var_list) matching%process_name = process_name @ <>= procedure (matching_write), deferred :: write <>= abstract interface subroutine matching_write (matching, unit) import class(matching_t), intent(in) :: matching integer, intent(in), optional :: unit end subroutine matching_write end interface @ %def matching_write @ <>= procedure :: import_rng => matching_import_rng +<>= + pure module subroutine matching_import_rng (matching, rng) + class(matching_t), intent(inout) :: matching + class(rng_t), allocatable, intent(inout) :: rng + end subroutine matching_import_rng <>= - pure subroutine matching_import_rng (matching, rng) + pure module subroutine matching_import_rng (matching, rng) class(matching_t), intent(inout) :: matching class(rng_t), allocatable, intent(inout) :: rng call move_alloc (from = rng, to = matching%rng) end subroutine matching_import_rng @ %def matching_import_rng @ <>= procedure :: connect => matching_connect procedure :: base_connect => matching_connect +<>= + module subroutine matching_connect & + (matching, process_instance, model, shower) + class(matching_t), intent(inout) :: matching + type(process_instance_t), intent(in), target :: process_instance + class(model_data_t), intent(in), target, optional :: model + class(shower_base_t), intent(in), target, optional :: shower + end subroutine matching_connect <>= - subroutine matching_connect (matching, process_instance, model, shower) + module subroutine matching_connect (matching, process_instance, model, shower) class(matching_t), intent(inout) :: matching type(process_instance_t), intent(in), target :: process_instance class(model_data_t), intent(in), target, optional :: model class(shower_base_t), intent(in), target, optional :: shower if (debug_on) call msg_debug (D_MATCHING, "matching_connect") matching%process_instance => process_instance if (present (model)) matching%model => model if (present (shower)) matching%shower => shower end subroutine matching_connect @ %def matching_base_connect @ <>= procedure (matching_before_shower), deferred :: before_shower <>= abstract interface subroutine matching_before_shower (matching, particle_set, vetoed) import class(matching_t), intent(inout) :: matching type(particle_set_t), intent(inout) :: particle_set logical, intent(out) :: vetoed end subroutine matching_before_shower end interface @ %def matching_before_shower @ <>= procedure (matching_after_shower), deferred :: after_shower <>= abstract interface subroutine matching_after_shower (matching, particle_set, vetoed) import class(matching_t), intent(inout) :: matching type(particle_set_t), intent(inout) :: particle_set logical, intent(out) :: vetoed end subroutine matching_after_shower end interface @ %def matching_after_shower @ Per default, do nothing here. <>= procedure :: prepare_for_events => matching_prepare_for_events +<>= + module subroutine matching_prepare_for_events (matching) + class(matching_t), intent(inout), target :: matching + end subroutine matching_prepare_for_events <>= - subroutine matching_prepare_for_events (matching) + module subroutine matching_prepare_for_events (matching) class(matching_t), intent(inout), target :: matching end subroutine matching_prepare_for_events @ %def matching_prepare_for_events @ <>= procedure :: first_event => matching_first_event +<>= + module subroutine matching_first_event (matching) + class(matching_t), intent(inout), target :: matching + end subroutine matching_first_event <>= - subroutine matching_first_event (matching) + module subroutine matching_first_event (matching) class(matching_t), intent(inout), target :: matching end subroutine matching_first_event @ %def matching_first_event @ <>= procedure (matching_get_method), deferred :: get_method <>= abstract interface function matching_get_method (matching) result (method) import type(string_t) :: method class(matching_t), intent(in) :: matching end function matching_get_method end interface @ %def matching_after_shower @ <>= procedure :: final => matching_final +<>= + module subroutine matching_final (matching) + class(matching_t), intent(in) :: matching + end subroutine matching_final <>= - subroutine matching_final (matching) + module subroutine matching_final (matching) class(matching_t), intent(in) :: matching end subroutine matching_final @ %def matching_final @ \subsection{Matching implementations} <>= public :: MATCH_MLM, MATCH_CKKW, MATCH_POWHEG <>= integer, parameter :: MATCH_MLM = 1 integer, parameter :: MATCH_CKKW = 2 integer, parameter :: MATCH_POWHEG = 3 integer, parameter :: MATCH_UNDEFINED = 17 @ %def MATCH_MLM MATCH_CKKW MATCH_POWHEG MATCH_UNDEFINED @ A dictionary <>= public :: matching_method <>= interface matching_method module procedure matching_method_of_string module procedure matching_method_to_string end interface +<>= + elemental module function matching_method_of_string (string) result (i) + integer :: i + type(string_t), intent(in) :: string + end function matching_method_of_string + elemental module function matching_method_to_string (i) result (string) + type(string_t) :: string + integer, intent(in) :: i + end function matching_method_to_string <>= - elemental function matching_method_of_string (string) result (i) + elemental module function matching_method_of_string (string) result (i) integer :: i type(string_t), intent(in) :: string select case (char (string)) case ("MLM") i = MATCH_MLM case ("CKKW") i = MATCH_CKKW case ("POWHEG") i = MATCH_POWHEG case default i = MATCH_UNDEFINED end select end function matching_method_of_string - elemental function matching_method_to_string (i) result (string) + elemental module function matching_method_to_string (i) result (string) type(string_t) :: string integer, intent(in) :: i select case (i) case (MATCH_MLM) string = "MLM" case (MATCH_CKKW) string = "CKKW" case (MATCH_POWHEG) string = "POWHEG" case default string = "UNDEFINED" end select end function matching_method_to_string @ %def matching_method @ \section{MLM Matching} <<[[mlm_matching.f90]]>>= <> module mlm_matching <> <> -<> - use io_units use constants - use format_utils, only: write_separator - use diagnostics - use file_utils use lorentz - use subevents, only: PRT_OUTGOING use particles use variables - use shower_base - use ktclus use matching_base <> <> <> + interface +<> + end interface + +end module mlm_matching +@ %def mlm_matching +@ +<<[[mlm_matching_sub.f90]]>>= +<> + +submodule (mlm_matching) mlm_matching_s + +<> + use io_units + use format_utils, only: write_separator + use diagnostics + use file_utils + use subevents, only: PRT_OUTGOING + use shower_base + use ktclus + + implicit none + contains <> -end module mlm_matching -@ %def mlm_matching +end submodule mlm_matching_s + +@ %def mlm_matching_s @ <>= public :: mlm_matching_settings_t <>= type :: mlm_matching_settings_t real(default) :: mlm_Qcut_ME = one real(default) :: mlm_Qcut_PS = one real(default) :: mlm_ptmin, mlm_etamax, mlm_Rmin, mlm_Emin real(default) :: mlm_ETclusfactor = 0.2_default real(default) :: mlm_ETclusminE = five real(default) :: mlm_etaclusfactor = one real(default) :: mlm_Rclusfactor = one real(default) :: mlm_Eclusfactor = one integer :: kt_imode_hadronic = 4313 integer :: kt_imode_leptonic = 1111 integer :: mlm_nmaxMEjets = 0 contains <> end type mlm_matching_settings_t @ %def mlm_matching_settings_t @ <>= procedure :: init => mlm_matching_settings_init +<>= + module subroutine mlm_matching_settings_init (settings, var_list) + class(mlm_matching_settings_t), intent(out) :: settings + type(var_list_t), intent(in) :: var_list + end subroutine mlm_matching_settings_init <>= - subroutine mlm_matching_settings_init (settings, var_list) + module subroutine mlm_matching_settings_init (settings, var_list) class(mlm_matching_settings_t), intent(out) :: settings type(var_list_t), intent(in) :: var_list settings%mlm_Qcut_ME = & var_list%get_rval (var_str ("mlm_Qcut_ME")) settings%mlm_Qcut_PS = & var_list%get_rval (var_str ("mlm_Qcut_PS")) settings%mlm_ptmin = & var_list%get_rval (var_str ("mlm_ptmin")) settings%mlm_etamax = & var_list%get_rval (var_str ("mlm_etamax")) settings%mlm_Rmin = & var_list%get_rval (var_str ("mlm_Rmin")) settings%mlm_Emin = & var_list%get_rval (var_str ("mlm_Emin")) settings%mlm_nmaxMEjets = & var_list%get_ival (var_str ("mlm_nmaxMEjets")) settings%mlm_ETclusfactor = & var_list%get_rval (var_str ("mlm_ETclusfactor")) settings%mlm_ETclusminE = & var_list%get_rval (var_str ("mlm_ETclusminE")) settings%mlm_etaclusfactor = & var_list%get_rval (var_str ("mlm_etaclusfactor")) settings%mlm_Rclusfactor = & var_list%get_rval (var_str ("mlm_Rclusfactor")) settings%mlm_Eclusfactor = & var_list%get_rval (var_str ("mlm_Eclusfactor")) end subroutine mlm_matching_settings_init @ %def mlm_matching_settings_init @ <>= procedure :: write => mlm_matching_settings_write +<>= + module subroutine mlm_matching_settings_write (settings, unit) + class(mlm_matching_settings_t), intent(in) :: settings + integer, intent(in), optional :: unit + end subroutine mlm_matching_settings_write <>= - subroutine mlm_matching_settings_write (settings, unit) + module subroutine mlm_matching_settings_write (settings, unit) class(mlm_matching_settings_t), intent(in) :: settings integer, intent(in), optional :: unit integer :: u u = given_output_unit (unit); if (u < 0) return write (u, "(3x,A,ES19.12)") & "mlm_Qcut_ME = ", settings%mlm_Qcut_ME write (u, "(3x,A,ES19.12)") & "mlm_Qcut_PS = ", settings%mlm_Qcut_PS write (u, "(3x,A,ES19.12)") & "mlm_ptmin = ", settings%mlm_ptmin write (u, "(3x,A,ES19.12)") & "mlm_etamax = ", settings%mlm_etamax write (u, "(3x,A,ES19.12)") & "mlm_Rmin = ", settings%mlm_Rmin write (u, "(3x,A,ES19.12)") & "mlm_Emin = ", settings%mlm_Emin write (u, "(3x,A,1x,I0)") & "mlm_nmaxMEjets = ", settings%mlm_nmaxMEjets write (u, "(3x,A,ES19.12)") & "mlm_ETclusfactor (D=0.2) = ", settings%mlm_ETclusfactor write (u, "(3x,A,ES19.12)") & "mlm_ETclusminE (D=5.0) = ", settings%mlm_ETclusminE write (u, "(3x,A,ES19.12)") & "mlm_etaclusfactor (D=1.0) = ", settings%mlm_etaClusfactor write (u, "(3x,A,ES19.12)") & "mlm_Rclusfactor (D=1.0) = ", settings%mlm_RClusfactor write (u, "(3x,A,ES19.12)") & "mlm_Eclusfactor (D=1.0) = ", settings%mlm_EClusfactor end subroutine mlm_matching_settings_write @ %def mlm_matching_settings_write @ This is a container for the (colored) parton momenta as well as the jet momenta. <>= public :: mlm_matching_t <>= type, extends (matching_t) :: mlm_matching_t type(vector4_t), dimension(:), allocatable, public :: P_ME type(vector4_t), dimension(:), allocatable, public :: P_PS type(vector4_t), dimension(:), allocatable, private :: JETS_ME type(vector4_t), dimension(:), allocatable, private :: JETS_PS type(mlm_matching_settings_t) :: settings contains <> end type mlm_matching_t @ %def mlm_matching_t @ <>= procedure :: init => mlm_matching_init +<>= + module subroutine mlm_matching_init (matching, var_list, process_name) + class(mlm_matching_t), intent(out) :: matching + type(var_list_t), intent(in) :: var_list + type(string_t), intent(in) :: process_name + end subroutine mlm_matching_init <>= - subroutine mlm_matching_init (matching, var_list, process_name) + module subroutine mlm_matching_init (matching, var_list, process_name) class(mlm_matching_t), intent(out) :: matching <> end subroutine mlm_matching_init @ %def mlm_matching_init @ <>= procedure :: write => mlm_matching_write +<>= + module subroutine mlm_matching_write (matching, unit) + class(mlm_matching_t), intent(in) :: matching + integer, intent(in), optional :: unit + end subroutine mlm_matching_write <>= - subroutine mlm_matching_write (matching, unit) + module subroutine mlm_matching_write (matching, unit) class(mlm_matching_t), intent(in) :: matching integer, intent(in), optional :: unit integer :: i, u u = given_output_unit (unit); if (u < 0) return write (u, "(1x,A)") "MLM matching:" call matching%settings%write (u) write (u, "(3x,A)") "Momenta of ME partons:" if (allocated (matching%P_ME)) then do i = 1, size (matching%P_ME) write (u, "(4x)", advance = "no") call vector4_write (matching%P_ME(i), unit = u) end do else write (u, "(5x,A)") "[empty]" end if call write_separator (u) write (u, "(3x,A)") "Momenta of ME jets:" if (allocated (matching%JETS_ME)) then do i = 1, size (matching%JETS_ME) write (u, "(4x)", advance = "no") call vector4_write (matching%JETS_ME(i), unit = u) end do else write (u, "(5x,A)") "[empty]" end if call write_separator (u) write(u, "(3x,A)") "Momenta of shower partons:" if (allocated (matching%P_PS)) then do i = 1, size (matching%P_PS) write (u, "(4x)", advance = "no") call vector4_write (matching%P_PS(i), unit = u) end do else write (u, "(5x,A)") "[empty]" end if call write_separator (u) write (u, "(3x,A)") "Momenta of shower jets:" if (allocated (matching%JETS_PS)) then do i = 1, size (matching%JETS_PS) write (u, "(4x)", advance = "no") call vector4_write (matching%JETS_PS(i), unit = u) end do else write (u, "(5x,A)") "[empty]" end if call write_separator (u) end subroutine mlm_matching_write @ %def mlm_matching_write @ <>= procedure :: get_method => mlm_matching_get_method +<>= + module function mlm_matching_get_method (matching) result (method) + type(string_t) :: method + class(mlm_matching_t), intent(in) :: matching + end function mlm_matching_get_method <>= - function mlm_matching_get_method (matching) result (method) + module function mlm_matching_get_method (matching) result (method) type(string_t) :: method class(mlm_matching_t), intent(in) :: matching method = matching_method (MATCH_MLM) end function mlm_matching_get_method @ %def mlm_matching_get_method @ <>= procedure :: before_shower => mlm_matching_before_shower -<>= - subroutine mlm_matching_before_shower & +<>= + module subroutine mlm_matching_before_shower & (matching, particle_set, vetoed) + class(mlm_matching_t), intent(inout) :: matching + type(particle_set_t), intent(inout) :: particle_set + logical, intent(out) :: vetoed + end subroutine mlm_matching_before_shower +<>= + module subroutine mlm_matching_before_shower & + (matching, particle_set, vetoed) class(mlm_matching_t), intent(inout) :: matching type(particle_set_t), intent(inout) :: particle_set logical, intent(out) :: vetoed vetoed = .false. end subroutine mlm_matching_before_shower @ %def mlm_matching_before_shower @ <>= procedure :: after_shower => mlm_matching_after_shower +<>= + module subroutine mlm_matching_after_shower (matching, particle_set, vetoed) + class(mlm_matching_t), intent(inout) :: matching + type(particle_set_t), intent(inout) :: particle_set + logical, intent(out) :: vetoed + end subroutine mlm_matching_after_shower <>= - subroutine mlm_matching_after_shower (matching, particle_set, vetoed) + module subroutine mlm_matching_after_shower (matching, particle_set, vetoed) class(mlm_matching_t), intent(inout) :: matching type(particle_set_t), intent(inout) :: particle_set logical, intent(out) :: vetoed if (debug_on) call msg_debug (D_MATCHING, "mlm_matching_after_shower") call matching%shower%get_final_colored_ME_momenta (matching%P_ME) call matching%fill_P_PS (particle_set) !!! MLM stage 3 -> reconstruct and possibly reject call matching%apply (vetoed) if (debug_active (D_MATCHING)) call matching%write () if (allocated (matching%P_ME)) deallocate (matching%P_ME) if (allocated (matching%P_PS)) deallocate (matching%P_PS) if (allocated (matching%JETS_ME)) deallocate (matching%JETS_ME) if (allocated (matching%JETS_PS)) deallocate (matching%JETS_PS) end subroutine mlm_matching_after_shower @ %def mlm_matching_after_shower @ Transfer partons after parton shower to [[matching%P_PS]] <>= procedure :: fill_P_PS => mlm_matching_fill_P_PS +<>= + module subroutine mlm_matching_fill_P_PS (matching, particle_set) + class(mlm_matching_t), intent(inout) :: matching + type(particle_set_t), intent(in) :: particle_set + end subroutine mlm_matching_fill_P_PS <>= - subroutine mlm_matching_fill_P_PS (matching, particle_set) + module subroutine mlm_matching_fill_P_PS (matching, particle_set) class(mlm_matching_t), intent(inout) :: matching type(particle_set_t), intent(in) :: particle_set integer :: i, j, n_jets_PS integer, dimension(2) :: col type(particle_t) :: tempprt real(double) :: eta type(vector4_t) :: p_tmp !!! loop over particles and extract final colored ones with eta matching%settings%mlm_etaClusfactor * & matching%settings%mlm_etamax) then if (debug_active (D_MATCHING)) then call msg_debug (D_MATCHING, "Rejecting this particle") call tempprt%write () end if cycle end if end if n_jets_PS = n_jets_PS + 1 end do allocate (matching%P_PS(1:n_jets_PS)) if (debug_on) call msg_debug (D_MATCHING, "n_jets_ps", n_jets_ps) j = 1 do i = 1, particle_set%get_n_tot () tempprt = particle_set%get_particle (i) if (tempprt%get_status () /= PRT_OUTGOING) cycle col = tempprt%get_color () if (all(col == 0)) cycle ! TODO: (bcn 2015-04-28) where is the corresponding part for lepton colliders? if (matching%is_hadron_collision) then p_tmp = tempprt%get_momentum () if (energy (p_tmp) - longitudinal_part (p_tmp) < 1.E-10_default .or. & energy (p_tmp) + longitudinal_part (p_tmp) < 1.E-10_default) then eta = pseudorapidity (p_tmp) else eta = rapidity (p_tmp) end if if (eta > matching%settings%mlm_etaClusfactor * & matching%settings%mlm_etamax) cycle end if matching%P_PS(j) = tempprt%get_momentum () j = j + 1 end do end subroutine mlm_matching_fill_P_PS @ %def mlm_matching_fill_P_PS @ <>= procedure :: apply => mlm_matching_apply +<>= + module subroutine mlm_matching_apply (matching, vetoed) + class(mlm_matching_t), intent(inout) :: matching + logical, intent(out) :: vetoed + end subroutine mlm_matching_apply <>= - subroutine mlm_matching_apply (matching, vetoed) + module subroutine mlm_matching_apply (matching, vetoed) class(mlm_matching_t), intent(inout) :: matching logical, intent(out) :: vetoed integer :: i, j integer :: n_jets_ME, n_jets_PS, n_jets_PS_atycut real(double) :: ycut real(double), dimension(:, :), allocatable :: PP real(double), dimension(:), allocatable :: Y real(double), dimension(:,:), allocatable :: P_JETS real(double), dimension(:,:), allocatable :: P_ME integer, dimension(:), allocatable :: JET integer :: NJET, NSUB integer :: imode !!! TODO: (bcn 2014-03-26) Why is ECUT hard coded to 1? !!! It is the denominator of the KT measure. Candidate for removal real(double) :: ECUT = 1._double integer :: ip1,ip2 ! KTCLUS COMMON BLOCK INTEGER NMAX,NUM,HIST PARAMETER (NMAX=512) DOUBLE PRECISION P,KT,KTP,KTS,ETOT,RSQ,KTLAST COMMON /KTCOMM/ETOT,RSQ,P(9,NMAX),KTP(NMAX,NMAX),KTS(NMAX), & KT(NMAX),KTLAST(NMAX),HIST(NMAX),NUM vetoed = .true. if (signal_is_pending ()) return <> <> <> <> <> vetoed = .false. 999 continue end subroutine mlm_matching_apply @ %def mlm_matching_apply @ <>= if (allocated (matching%P_ME)) then ! print *, "number of partons after ME: ", size(matching%P_ME) n_jets_ME = size (matching%P_ME) else n_jets_ME = 0 end if if (allocated (matching%p_PS)) then ! print *, "number of partons after PS: ", size(matching%p_PS) n_jets_PS = size (matching%p_PS) else n_jets_PS = 0 end if @ <>= if (n_jets_ME > 0) then ycut = (matching%settings%mlm_ptmin)**2 allocate (PP(1:4, 1:N_jets_ME)) do i = 1, n_jets_ME PP(1:3,i) = matching%p_ME(i)%p(1:3) PP(4,i) = matching%p_ME(i)%p(0) end do <> allocate (P_ME(1:4,1:n_jets_ME)) allocate (JET(1:n_jets_ME)) allocate (Y(1:n_jets_ME)) if (signal_is_pending ()) return call KTCLUR (imode, PP, n_jets_ME, & dble (matching%settings%mlm_Rclusfactor * matching%settings%mlm_Rmin), ECUT, y, *999) call KTRECO (1, PP, n_jets_ME, ECUT, ycut, ycut, P_ME, JET, & NJET, NSUB, *999) n_jets_ME = NJET if (NJET > 0) then allocate (matching%JETS_ME (1:NJET)) do i = 1, NJET matching%JETS_ME(i) = vector4_moving (REAL(P_ME(4,i), default), & vector3_moving([REAL(P_ME(1,i), default), & REAL(P_ME(2,i), default), REAL(P_ME(3,i), default)])) end do end if deallocate (P_ME) deallocate (JET) deallocate (Y) deallocate (PP) end if @ <>= if (n_jets_PS > 0) then ycut = (matching%settings%mlm_ptmin + max (matching%settings%mlm_ETclusminE, & matching%settings%mlm_ETclusfactor * matching%settings%mlm_ptmin))**2 allocate (PP(1:4, 1:n_jets_PS)) do i = 1, n_jets_PS PP(1:3,i) = matching%p_PS(i)%p(1:3) PP(4,i) = matching%p_PS(i)%p(0) end do <> allocate (P_JETS(1:4,1:n_jets_PS)) allocate (JET(1:n_jets_PS)) allocate (Y(1:n_jets_PS)) if (signal_is_pending ()) return call KTCLUR (imode, PP, n_jets_PS, & dble (matching%settings%mlm_Rclusfactor * matching%settings%mlm_Rmin), & ECUT, y, *999) call KTRECO (1, PP, n_jets_PS, ECUT, ycut, ycut, P_JETS, JET, & NJET, NSUB, *999) n_jets_PS_atycut = NJET if (n_jets_ME == matching%settings%mlm_nmaxMEjets .and. NJET > 0) then ! print *, " resetting ycut to ", Y(matching%settings%mlm_nmaxMEjets) ycut = y(matching%settings%mlm_nmaxMEjets) call KTRECO (1, PP, n_jets_PS, ECUT, ycut, ycut, P_JETS, JET, & NJET, NSUB, *999) end if ! !Sample of code for a FastJet interface ! palg = 1d0 ! 1.0d0 = kt, 0.0d0 = Cam/Aachen, -1.0d0 = anti-kt ! R = 0.7_double ! radius parameter ! f = 0.75_double ! overlap threshold ! !call fastjetppgenkt(PP,n,R,palg,P_JETS,NJET) ! KT-Algorithm ! !call fastjetsiscone(PP,n,R,f,P_JETS,NJET) ! SiSCone-Algorithm if (NJET > 0) then allocate (matching%JETS_PS(1:NJET)) do i = 1, NJET matching%JETS_PS(i) = vector4_moving (REAL(P_JETS(4,i), default), & vector3_moving([REAL(P_JETS(1,i), default), & REAL(P_JETS(2,i), default), REAL(P_JETS(3,i), default)])) end do end if deallocate (P_JETS) deallocate (JET) deallocate (Y) else n_jets_PS_atycut = 0 end if @ <>= if (matching%is_hadron_collision) then imode = matching%settings%kt_imode_hadronic else imode = matching%settings%kt_imode_leptonic end if @ <>= if (n_jets_PS_atycut < n_jets_ME) then ! print *, "DISCARDING: Not enough PS jets: ", n_jets_PS_atycut return end if if (n_jets_PS_atycut > n_jets_ME .and. n_jets_ME /= matching%settings%mlm_nmaxMEjets) then ! print *, "DISCARDING: Too many PS jets: ", n_jets_PS_atycut return end if @ <>= if (allocated(matching%JETS_PS)) then ! print *, "number of jets after PS: ", size(matching%JETS_PS) n_jets_PS = size (matching%JETS_PS) else n_jets_PS = 0 end if if (n_jets_ME > 0 .and. n_jets_PS > 0) then n_jets_PS = size (matching%JETS_PS) if (allocated (PP)) deallocate(PP) allocate (PP(1:4, 1:n_jets_PS + 1)) do i = 1, n_jets_PS if (signal_is_pending ()) return PP(1:3,i) = matching%JETS_PS(i)%p(1:3) PP(4,i) = matching%JETS_PS(i)%p(0) end do if (allocated (Y)) deallocate(Y) allocate (Y(1:n_jets_PS + 1)) y = zero do i = 1, n_jets_ME PP(1:3,n_jets_PS + 2 - i) = matching%JETS_ME(i)%p(1:3) PP(4,n_jets_PS + 2 - i) = matching%JETS_ME(i)%p(0) !!! This makes more sense than hardcoding ! call KTCLUS (4313, PP, (n_jets_PS + 2 - i), 1.0_double, Y, *999) call KTCLUR (imode, PP, (n_jets_PS + 2 - i), & dble (matching%settings%mlm_Rclusfactor * matching%settings%mlm_Rmin), & ECUT, y, *999) if (0.99 * y(n_jets_PS + 1 - (i - 1)).gt.ycut) then ! print *, "DISCARDING: Jet ", i, " not clusterd" return end if !!! search for and remove PS jet clustered with ME Jet ip1 = HIST(n_jets_PS + 2 - i) / NMAX ip2 = mod(hist(n_jets_PS + 2 - i), NMAX) if ((ip2 /= n_jets_PS + 2 - i) .or. (ip1 <= 0)) then ! print *, "DISCARDING: Jet ", i, " not clustered ", ip1, ip2, & ! hist(n_jets_PS + 2 - i) return else ! print *, "PARTON clustered", ip1, ip2, hist(n_jets_PS + 2 - i) PP(:,IP1) = zero do j = IP1, n_jets_PS - i PP(:, j) = PP(:,j + 1) end do end if end do end if @ \section{CKKW matching} This module contains the CKKW matching. The type [[ckkw_pseudo_shower_weights_t]] gives the (relative) weights for different clusterings of the final particles, as given in Eq.~(2.7) of hep-ph/0503281v1. Each particle has a binary labelling (power of 2) (first particle = 1, second particle = 2, third particle = 4, ...). Each recombination therefore corresponds to an integer, that is not a power of 2. Fur multiple subsequent recombinations, no different weights for different sequences of clustering are stored. It is assumed that the weight of a multiply recombined state is a combination of the states with one fewer recombination and that these states' contributions are proportional to their weights. For a $2->n$ event, the weights array thus has the size $2^{(2 + n) - 1}$. The [[weights_by_type]] array gives the weights depending on the type of the particle, the first index is the same as for weights, the second index gives the type of the new mother particle: \begin{itemize} \item[0:] uncolored ($\gamma$, $Z$, $W$, Higgs) \item[1:] colored (quark) \item[2:] gluon \item[3:] squark \item[4:] gluino \end{itemize} [[alphaS]] gives the value for $alpha_s$ used in the generation of the matrix element. This is needed for the reweighting using the values for a running $alpha_s$ at the scales of the clusterings. <<[[ckkw_matching.f90]]>>= <> module ckkw_matching <> <> -<> - use io_units use constants - use format_utils, only: write_separator - use diagnostics - use physics_defs use lorentz use particles use rng_base use shower_base use shower_partons - use shower_core use variables use matching_base <> <> <> + interface +<> + end interface + +end module ckkw_matching +@ %def ckkw_matching +@ +<<[[ckkw_matching_sub.f90]]>>= +<> + +submodule (ckkw_matching) ckkw_matching_s + +<> + use io_units + use format_utils, only: write_separator + use diagnostics + use physics_defs + use shower_core + + implicit none + contains <> -end module ckkw_matching -@ %def ckkw_matching +end submodule ckkw_matching_s + +@ %def ckkw_matching_s @ The fundamental CKKW matching parameter are defined here: <>= public :: ckkw_matching_settings_t <>= type :: ckkw_matching_settings_t real(default) :: alphaS = 0.118_default real(default) :: Qmin = one integer :: n_max_jets = 0 contains <> end type ckkw_matching_settings_t @ %def ckkw_matching_settings_t @ This is empty for the moment. <>= procedure :: init => ckkw_matching_settings_init +<>= + module subroutine ckkw_matching_settings_init (settings, var_list) + class(ckkw_matching_settings_t), intent(out) :: settings + type(var_list_t), intent(in) :: var_list + end subroutine ckkw_matching_settings_init <>= - subroutine ckkw_matching_settings_init (settings, var_list) + module subroutine ckkw_matching_settings_init (settings, var_list) class(ckkw_matching_settings_t), intent(out) :: settings type(var_list_t), intent(in) :: var_list settings%alphaS = 1.0_default settings%Qmin = 1.0_default settings%n_max_jets = 3 end subroutine ckkw_matching_settings_init @ %def ckkw_matching_settings_init @ <>= procedure :: write => ckkw_matching_settings_write +<>= + module subroutine ckkw_matching_settings_write (settings, unit) + class(ckkw_matching_settings_t), intent(in) :: settings + integer, intent(in), optional :: unit + end subroutine ckkw_matching_settings_write <>= - subroutine ckkw_matching_settings_write (settings, unit) + module subroutine ckkw_matching_settings_write (settings, unit) class(ckkw_matching_settings_t), intent(in) :: settings integer, intent(in), optional :: unit integer :: u u = given_output_unit (unit); if (u < 0) return write (u, "(1x,A)") "CKKW matching settings:" call write_separator (u) write (u, "(3x,A,1x,ES19.12)") & "alphaS = ", settings%alphaS write (u, "(3x,A,1x,ES19.12)") & "Qmin = ", settings%Qmin write (u, "(3x,A,1x,I0)") & "n_max_jets = ", settings%n_max_jets end subroutine ckkw_matching_settings_write @ %def ckkw_matching_settings_write @ <>= public :: ckkw_pseudo_shower_weights_t <>= type :: ckkw_pseudo_shower_weights_t real(default) :: alphaS real(default), dimension(:), allocatable :: weights real(default), dimension(:,:), allocatable :: weights_by_type contains <> end type ckkw_pseudo_shower_weights_t @ %def ckkw_pseudo_shower_weights_t @ <>= procedure :: init => ckkw_pseudo_shower_weights_init +<>= + module subroutine ckkw_pseudo_shower_weights_init (weights) + class(ckkw_pseudo_shower_weights_t), intent(out) :: weights + end subroutine ckkw_pseudo_shower_weights_init <>= - subroutine ckkw_pseudo_shower_weights_init (weights) + module subroutine ckkw_pseudo_shower_weights_init (weights) class(ckkw_pseudo_shower_weights_t), intent(out) :: weights weights%alphaS = zero end subroutine ckkw_pseudo_shower_weights_init @ %def ckkw_pseudo_shower_weights_init @ <>= procedure :: write => ckkw_pseudo_shower_weights_write +<>= + module subroutine ckkw_pseudo_shower_weights_write (weights, unit) + class(ckkw_pseudo_shower_weights_t), intent(in) :: weights + integer, intent(in), optional :: unit + end subroutine ckkw_pseudo_shower_weights_write <>= - subroutine ckkw_pseudo_shower_weights_write (weights, unit) + module subroutine ckkw_pseudo_shower_weights_write (weights, unit) class(ckkw_pseudo_shower_weights_t), intent(in) :: weights integer, intent(in), optional :: unit integer :: s, i, u u = given_output_unit (unit); if (u < 0) return s = size (weights%weights) write (u, "(1x,A)") "CKKW (pseudo) shower weights: " do i = 1, s write (u, "(3x,I0,2(ES19.12))") i, weights%weights(i), & weights%weights_by_type(i,:) end do write (u, "(3x,A,1x,I0)") "alphaS =", weights%alphaS end subroutine ckkw_pseudo_shower_weights_write @ %def ckkw_pseudo_shower_weights_write @ Generate fake ckkw weights. This can be dropped, once information from the matrix element generation is available. <>= procedure :: fake => ckkw_pseudo_shower_weights_fake +<>= + pure module subroutine ckkw_pseudo_shower_weights_fake & + (weights, particle_set) + class(ckkw_pseudo_shower_weights_t), intent(inout) :: weights + type(particle_set_t), intent(in) :: particle_set + end subroutine ckkw_pseudo_shower_weights_fake <>= - pure subroutine ckkw_pseudo_shower_weights_fake (weights, particle_set) + pure module subroutine ckkw_pseudo_shower_weights_fake (weights, particle_set) class(ckkw_pseudo_shower_weights_t), intent(inout) :: weights type(particle_set_t), intent(in) :: particle_set integer :: i, j, n type(vector4_t) :: momentum n = 2**particle_set%n_tot if (allocated (weights%weights)) then deallocate (weights%weights) end if allocate (weights%weights (1:n)) do i = 1, n momentum = vector4_null do j = 1, particle_set%n_tot if (btest (i,j-1)) then momentum = momentum + particle_set%prt(j)%p end if end do if (momentum**1 > 0.0) then weights%weights(i) = 1.0 / (momentum**2) end if end do ! equally distribute the weights by type if (allocated (weights%weights_by_type)) then deallocate (weights%weights_by_type) end if allocate (weights%weights_by_type (1:n, 0:4)) do i = 1, n do j = 0, 4 weights%weights_by_type(i,j) = 0.2 * weights%weights(i) end do end do end subroutine ckkw_pseudo_shower_weights_fake @ %def ckkw_pseudo_shower_weights_fake @ <>= public :: ckkw_matching_t <>= type, extends (matching_t) :: ckkw_matching_t type(ckkw_matching_settings_t) :: settings type(ckkw_pseudo_shower_weights_t) :: weights contains <> end type ckkw_matching_t @ %def ckkw_matching_t @ <>= procedure :: init => ckkw_matching_init +<>= + module subroutine ckkw_matching_init (matching, var_list, process_name) + class(ckkw_matching_t), intent(out) :: matching + type(var_list_t), intent(in) :: var_list + type(string_t), intent(in) :: process_name + end subroutine ckkw_matching_init <>= - subroutine ckkw_matching_init (matching, var_list, process_name) + module subroutine ckkw_matching_init (matching, var_list, process_name) class(ckkw_matching_t), intent(out) :: matching <> end subroutine ckkw_matching_init @ %def ckkw_matching_init @ <>= procedure :: write => ckkw_matching_write +<>= + module subroutine ckkw_matching_write (matching, unit) + class(ckkw_matching_t), intent(in) :: matching + integer, intent(in), optional :: unit + end subroutine ckkw_matching_write <>= - subroutine ckkw_matching_write (matching, unit) + module subroutine ckkw_matching_write (matching, unit) class(ckkw_matching_t), intent(in) :: matching integer, intent(in), optional :: unit call matching%settings%write (unit) call matching%weights%write (unit) end subroutine ckkw_matching_write @ %def ckkw_matching_write @ <>= procedure :: get_method => ckkw_matching_get_method +<>= + module function ckkw_matching_get_method (matching) result (method) + type(string_t) :: method + class(ckkw_matching_t), intent(in) :: matching + end function ckkw_matching_get_method <>= - function ckkw_matching_get_method (matching) result (method) - type(string_t) :: method - class(ckkw_matching_t), intent(in) :: matching - method = matching_method (MATCH_CKKW) + module function ckkw_matching_get_method (matching) result (method) + type(string_t) :: method + class(ckkw_matching_t), intent(in) :: matching + method = matching_method (MATCH_CKKW) end function ckkw_matching_get_method @ %def ckkw_matching_get_method @ <>= procedure :: before_shower => ckkw_matching_before_shower -<>= - subroutine ckkw_matching_before_shower & +<>= + module subroutine ckkw_matching_before_shower & (matching, particle_set, vetoed) + class(ckkw_matching_t), intent(inout) :: matching + type(particle_set_t), intent(inout) :: particle_set + logical, intent(out) :: vetoed + end subroutine ckkw_matching_before_shower +<>= + module subroutine ckkw_matching_before_shower & + (matching, particle_set, vetoed) class(ckkw_matching_t), intent(inout) :: matching type(particle_set_t), intent(inout) :: particle_set logical, intent(out) :: vetoed call matching%weights%init () call matching%weights%fake (particle_set) select type (shower => matching%shower) type is (shower_t) call ckkw_matching_apply (shower%partons, & matching%settings, & matching%weights, matching%rng, vetoed) class default call msg_bug ("CKKW matching only works with WHIZARD shower.") end select end subroutine ckkw_matching_before_shower @ %def ckkw_matching_before_shower @ <>= public :: ckkw_matching_apply +<>= + module subroutine ckkw_matching_apply & + (partons, settings, weights, rng, vetoed) + type(parton_pointer_t), dimension(:), intent(inout), allocatable :: & + partons + type(ckkw_matching_settings_t), intent(in) :: settings + type(ckkw_pseudo_shower_weights_t), intent(in) :: weights + class(rng_t), intent(inout), allocatable :: rng + logical, intent(out) :: vetoed + end subroutine ckkw_matching_apply <>= - subroutine ckkw_matching_apply (partons, settings, weights, rng, vetoed) + module subroutine ckkw_matching_apply & + (partons, settings, weights, rng, vetoed) type(parton_pointer_t), dimension(:), intent(inout), allocatable :: & partons type(ckkw_matching_settings_t), intent(in) :: settings type(ckkw_pseudo_shower_weights_t), intent(in) :: weights class(rng_t), intent(inout), allocatable :: rng logical, intent(out) :: vetoed real(default), dimension(:), allocatable :: scales real(double) :: weight, sf real(default) :: rand integer :: i, n_partons if (signal_is_pending ()) return weight = one n_partons = size (partons) do i = 1, n_partons call partons(i)%p%write () end do !!! the pseudo parton shower is already simulated by shower_add_interaction !!! get the respective clustering scales allocate (scales (1:n_partons)) do i = 1, n_partons if (.not. associated (partons(i)%p)) cycle if (partons(i)%p%type == INTERNAL) then scales(i) = two * min (partons(i)%p%child1%momentum%p(0), & partons(i)%p%child2%momentum%p(0))**2 * & (1.0 - (space_part (partons(i)%p%child1%momentum) * & space_part (partons(i)%p%child2%momentum)) / & (space_part (partons(i)%p%child1%momentum)**1 * & space_part (partons(i)%p%child2%momentum)**1)) scales(i) = sqrt (scales(i)) partons(i)%p%ckkwscale = scales(i) print *, scales(i) end if end do print *, " scales finished" !!! if (highest multiplicity) -> reweight with PDF(mu_F) / PDF(mu_cut) do i = 1, n_partons call partons(i)%p%write () end do !!! Reweight and possibly veto the whole event !!! calculate the relative alpha_S weight !! calculate the Sudakov weights for internal lines !! calculate the Sudakov weights for external lines do i = 1, n_partons if (signal_is_pending ()) return if (.not. associated (partons(i)%p)) cycle if (partons(i)%p%type == INTERNAL) then !!! get type !!! check that all particles involved are colored if ((partons(i)%p%is_colored () .or. & partons(i)%p%ckkwtype > 0) .and. & (partons(i)%p%child1%is_colored () .or. & partons(i)%p%child1%ckkwtype > 0) .and. & (partons(i)%p%child1%is_colored () .or. & partons(i)%p%child1%ckkwtype > 0)) then print *, "reweight with alphaS(" , partons(i)%p%ckkwscale, & ") for particle ", partons(i)%p%nr if (partons(i)%p%belongstoFSR) then print *, "FSR" weight = weight * D_alpha_s_fsr (partons(i)%p%ckkwscale**2, & partons(i)%p%settings) / settings%alphas else print *, "ISR" weight = weight * & D_alpha_s_isr (partons(i)%p%ckkwscale**2, & partons(i)%p%settings) / settings%alphas end if else print *, "no reweight with alphaS for ", partons(i)%p%nr end if if (partons(i)%p%child1%type == INTERNAL) then print *, "internal line from ", & partons(i)%p%child1%ckkwscale, & " to ", partons(i)%p%ckkwscale, & " for type ", partons(i)%p%child1%ckkwtype if (partons(i)%p%child1%ckkwtype == 0) then sf = 1.0 else if (partons(i)%p%child1%ckkwtype == 1) then sf = SudakovQ (partons(i)%p%child1%ckkwscale, & partons(i)%p%ckkwscale, & partons(i)%p%settings, .true., rng) print *, "SFQ = ", sf else if (partons(i)%p%child1%ckkwtype == 2) then sf = SudakovG (partons(i)%p%child1%ckkwscale, & partons(i)%p%ckkwscale, & partons(i)%p%settings, .true., rng) print *, "SFG = ", sf else print *, "SUSY not yet implemented" end if weight = weight * min (one, sf) else print *, "external line from ", settings%Qmin, & partons(i)%p%ckkwscale if (partons(i)%p%child1%is_quark ()) then sf = SudakovQ (settings%Qmin, & partons(i)%p%ckkwscale, & partons(i)%p%settings, .true., rng) print *, "SFQ = ", sf else if (partons(i)%p%child1%is_gluon ()) then sf = SudakovG (settings%Qmin, & partons(i)%p%ckkwscale, & partons(i)%p%settings, .true., rng) print *, "SFG = ", sf else print *, "not yet implemented (", & partons(i)%p%child2%type, ")" sf = one end if weight = weight * min (one, sf) end if if (partons(i)%p%child2%type == INTERNAL) then print *, "internal line from ", partons(i)%p%child2%ckkwscale, & " to ", partons(i)%p%ckkwscale, & " for type ", partons(i)%p%child2%ckkwtype if (partons(i)%p%child2%ckkwtype == 0) then sf = 1.0 else if (partons(i)%p%child2%ckkwtype == 1) then sf = SudakovQ (partons(i)%p%child2%ckkwscale, & partons(i)%p%ckkwscale, & partons(i)%p%settings, .true., rng) print *, "SFQ = ", sf else if (partons(i)%p%child2%ckkwtype == 2) then sf = SudakovG (partons(i)%p%child2%ckkwscale, & partons(i)%p%ckkwscale, & partons(i)%p%settings, .true., rng) print *, "SFG = ", sf else print *, "SUSY not yet implemented" end if weight = weight * min (one, sf) else print *, "external line from ", settings%Qmin, & partons(i)%p%ckkwscale if (partons(i)%p%child2%is_quark ()) then sf = SudakovQ (settings%Qmin, & partons(i)%p%ckkwscale, & partons(i)%p%settings, .true., rng) print *, "SFQ = ", sf else if (partons(i)%p%child2%is_gluon ()) then sf = SudakovG (settings%Qmin, & partons(i)%p%ckkwscale, & partons(i)%p%settings, .true., rng) print *, "SFG = ", sf else print *, "not yet implemented (", & partons(i)%p%child2%type, ")" sf = one end if weight = weight * min (one, sf) end if end if end do call rng%generate (rand) print *, "final weight: ", weight !!!!!!! WRONG vetoed = .false. ! vetoed = (rand > weight) if (vetoed) then return end if !!! finally perform the parton shower !!! veto emissions that are too hard deallocate (scales) end subroutine ckkw_matching_apply @ %def ckkw_matching_apply @ @ <>= procedure :: after_shower => ckkw_matching_after_shower +<>= + module subroutine ckkw_matching_after_shower & + (matching, particle_set, vetoed) + class(ckkw_matching_t), intent(inout) :: matching + type(particle_set_t), intent(inout) :: particle_set + logical, intent(out) :: vetoed + end subroutine ckkw_matching_after_shower <>= - subroutine ckkw_matching_after_shower (matching, particle_set, vetoed) + module subroutine ckkw_matching_after_shower (matching, particle_set, vetoed) class(ckkw_matching_t), intent(inout) :: matching - type(particle_set_t), intent(inout) :: particle_set - logical, intent(out) :: vetoed - vetoed = .false. + type(particle_set_t), intent(inout) :: particle_set + logical, intent(out) :: vetoed + vetoed = .false. end subroutine ckkw_matching_after_shower @ %def ckkw_matching_after_shower <>= function GammaQ (smallq, largeq, settings, fsr) result (gamma) real(default), intent(in) :: smallq, largeq type(shower_settings_t), intent(in) :: settings logical, intent(in) :: fsr real(default) :: gamma gamma = (8._default / three) / (pi * smallq) gamma = gamma * (log(largeq / smallq) - 0.75) if (fsr) then gamma = gamma * D_alpha_s_fsr (smallq**2, settings) else gamma = gamma * D_alpha_s_isr (smallq**2, settings) end if end function GammaQ @ %def GammaQ @ <>= function GammaG (smallq, largeq, settings, fsr) result (gamma) real(default), intent(in) :: smallq, largeq type(shower_settings_t), intent(in) :: settings logical, intent(in) :: fsr real(default) :: gamma gamma = 6._default / (pi * smallq) gamma = gamma *( log(largeq / smallq) - 11.0 / 12.0) if (fsr) then gamma = gamma * D_alpha_s_fsr (smallq**2, settings) else gamma = gamma * D_alpha_s_isr (smallq**2, settings) end if end function GammaG @ %def GammaG @ <>= function GammaF (smallq, settings, fsr) result (gamma) real(default), intent(in) :: smallq type(shower_settings_t), intent(in) :: settings logical, intent(in) :: fsr real(default) :: gamma gamma = number_of_flavors (smallq, settings%max_n_flavors, & settings%min_virtuality) / (three * pi * smallq) if (fsr) then gamma = gamma * D_alpha_s_fsr (smallq**2, settings) else gamma = gamma * D_alpha_s_isr (smallq**2, settings) end if end function GammaF @ %def GammaF @ <>= function SudakovQ (Q1, Q, settings, fsr, rng) result (sf) real(default), intent(in) :: Q1, Q type(shower_settings_t), intent(in) :: settings class(rng_t), intent(inout), allocatable :: rng logical, intent(in) :: fsr real(default) :: sf real(default) :: integral integer, parameter :: NTRIES = 100 integer :: i real(default) :: rand integral = zero do i = 1, NTRIES call rng%generate (rand) integral = integral + GammaQ (Q1 + rand * (Q - Q1), Q, settings, fsr) end do integral = integral / NTRIES sf = exp (-integral) end function SudakovQ @ %def SudakovQ @ <>= function SudakovG (Q1, Q, settings, fsr, rng) result (sf) real(default), intent(in) :: Q1, Q type(shower_settings_t), intent(in) :: settings logical, intent(in) :: fsr real(default) :: sf real(default) :: integral class(rng_t), intent(inout), allocatable :: rng integer, parameter :: NTRIES = 100 integer :: i real(default) :: rand integral = zero do i = 1, NTRIES call rng%generate (rand) integral = integral + & GammaG (Q1 + rand * (Q - Q1), Q, settings, fsr) + & GammaF (Q1 + rand * (Q - Q1), settings, fsr) end do integral = integral / NTRIES sf = exp (-integral) end function SudakovG @ %def SudakovG @ \section{POWHEG} This module generates radiation according to the POWHEG Sudakov form factor \begin{equation} \Delta^{f_b} (\Phi_n, p_\text{T}) = \prod_{\alpha_r \in \{\alpha_r |f_b \}} \Delta^{f_b}_{\alpha_r} (\Phi_n, p_\text{T}), \end{equation} with \begin{equation} \Delta^{f_b}_{\alpha_r} (\Phi_n, p_\text{T}) = \exp \left\{ - \left[ \int d \Phi_{\text{rad}} \,\frac{R (\Phi_{n+1})}{B^{f_b} (\Phi_n)} \,\theta( k_\text{T} (\Phi_{n+1}) - p_\text{T}) \right]^{\bar{\bf \Phi}_n^{\alpha_r} = \Phi_n}_{\alpha_r} \right\} \end{equation} We expect that an underlying Born flavor structure $f_b$ has been generated with a probability proportional to its contribution to the $\tilde B$ at the given kinematic point. <<[[powheg_matching.f90]]>>= <> module powheg_matching use, intrinsic :: iso_fortran_env <> <> -<> use diagnostics use constants, only: ZERO, ONE, TWO, THREE, FOUR, FIVE, TINY_07, PI, TWOPI - use numeric_utils - use io_units, only: given_output_unit, free_unit - use format_utils, only: write_separator - use format_defs, only: FMT_16, FMT_19 - use string_utils, only: str - use os_interface, only: os_file_exist - use physics_defs, only: CA, BORN, NLO_REAL - use pdg_arrays, only: is_gluon, is_quark use pdf, only: pdf_data_t use lorentz use phs_points, only: assignment(=), operator(*) use sm_qcd, only: qcd_t, alpha_qcd_from_scale_t, alpha_qcd_from_lambda_t - use sf_lhapdf, only: alpha_qcd_lhapdf_t - use sm_physics, only: Li2 - use subevents, only: PRT_INCOMING, PRT_OUTGOING - use colors use particles use grids use solver use rng_base use variables - use process_config, only: COMP_REAL_FIN use phs_fks, only: phs_fks_generator_t, compute_dalitz_bounds, beta_emitter use phs_fks, only: phs_point_set_t, phs_identifier_t, phs_fks_t use phs_fks, only: get_xi_max_isr use phs_fks, only: I_XI, I_Y, I_PLUS, I_MINUS, UBF_FSR_SIMPLE, UBF_FSR_MASSIVE, UBF_FSR_MASSLESS_RECOIL, UBF_ISR use matching_base use instances, only: process_instance_t, process_instance_hook_t use pcm, only: pcm_nlo_t, pcm_nlo_workspace_t -<> <> <> <> <> <> + interface +<> + end interface + contains -<> +<> end module powheg_matching @ %def powheg_matching @ +<<[[powheg_matching_sub.f90]]>>= +<> + +submodule (powheg_matching) powheg_matching_s + +<> +<> + use io_units, only: given_output_unit, free_unit + use format_utils, only: write_separator + use format_defs, only: FMT_16, FMT_19 + use string_utils, only: str + use numeric_utils + use os_interface, only: os_file_exist + use physics_defs, only: CA, BORN, NLO_REAL + use pdg_arrays, only: is_gluon, is_quark + use sf_lhapdf, only: alpha_qcd_lhapdf_t + use sm_physics, only: Li2 + use subevents, only: PRT_INCOMING, PRT_OUTGOING + use colors + use process_config, only: COMP_REAL_FIN + + implicit none + +contains + +<> + +end submodule powheg_matching_s + +@ %def powheg_matching_s +@ \subsection{Base types for settings and data} [[lambda]] enters for now as the lowest scale $2 \Lambda^(5)_{\bar{\text{MS}}}$ where the radiation $\alpha_s^\text{rad}$ is still larger than the true $\alpha_s$. <>= public :: powheg_settings_t <>= type :: powheg_settings_t real(default) :: pt2_min = zero real(default) :: lambda = zero integer :: size_grid_xi = 0 integer :: size_grid_y = 0 integer :: upper_bound_func_type = UBF_FSR_SIMPLE logical :: test_sudakov = .false. logical :: disable_sudakov = .false. logical :: singular_jacobian = .false. contains <> end type powheg_settings_t @ %def powheg_settings_t @ <>= procedure :: init => powheg_settings_init +<>= + module subroutine powheg_settings_init (settings, var_list) + class(powheg_settings_t), intent(out) :: settings + type(var_list_t), intent(in) :: var_list + end subroutine powheg_settings_init <>= - subroutine powheg_settings_init (settings, var_list) + module subroutine powheg_settings_init (settings, var_list) class(powheg_settings_t), intent(out) :: settings type(var_list_t), intent(in) :: var_list settings%pt2_min = & var_list%get_rval (var_str ("powheg_pt_min"))**2 settings%size_grid_xi = & var_list%get_ival (var_str ("powheg_grid_size_xi")) settings%size_grid_y = & var_list%get_ival (var_str ("powheg_grid_size_y")) settings%lambda = var_list%get_rval (var_str ("powheg_lambda")) settings%singular_jacobian = & var_list%get_lval (var_str ("?powheg_use_singular_jacobian")) settings%test_sudakov = & var_list%get_lval (var_str ("?powheg_test_sudakov")) settings%disable_sudakov = & var_list%get_lval (var_str ("?powheg_disable_sudakov")) end subroutine powheg_settings_init @ %def powheg_settings_init @ <>= procedure :: write => powheg_settings_write +<>= + module subroutine powheg_settings_write (powheg_settings, unit) + class(powheg_settings_t), intent(in) :: powheg_settings + integer, intent(in), optional :: unit + end subroutine powheg_settings_write <>= - subroutine powheg_settings_write (powheg_settings, unit) + module subroutine powheg_settings_write (powheg_settings, unit) class(powheg_settings_t), intent(in) :: powheg_settings integer, intent(in), optional :: unit integer :: u u = given_output_unit (unit); if (u < 0) return write (u, "(1X,A)") "POWHEG settings:" write (u, "(3X,A," // FMT_16 //")") "pt2_min = ", powheg_settings%pt2_min write (u, "(3X,A," // FMT_16 //")") "lambda = ", powheg_settings%lambda write (u, "(3X,A,I12)") "size_grid_xi = ", powheg_settings%size_grid_xi write (u, "(3X,A,I12)") "size_grid_y = ", powheg_settings%size_grid_y write (u, "(3X,A,I12)") "upper_bound_func_type = ", powheg_settings%upper_bound_func_type end subroutine powheg_settings_write @ %def powheg_settings_write @ <>= public :: radiation_t <>= type :: radiation_t real(default) :: xi, y, phi, pt2 integer :: alr logical :: valid = .false. contains <> end type radiation_t @ %def radiation_t @ <>= procedure :: write => radiation_write +<>= + module subroutine radiation_write (radiation, unit) + class(radiation_t), intent(in) :: radiation + integer, intent(in), optional :: unit + end subroutine radiation_write <>= - subroutine radiation_write (radiation, unit) + module subroutine radiation_write (radiation, unit) class(radiation_t), intent(in) :: radiation integer, intent(in), optional :: unit integer :: u u = given_output_unit (unit); if (u < 0) return write (u, "(1X, A)") "Radiation:" write (u, "(3X, A," // FMT_16 // ")") "xi = ", radiation%xi write (u, "(3X, A," // FMT_16 // ")") "y = ", radiation%y write (u, "(3X, A," // FMT_16 // ")") "phi = ", radiation%phi write (u, "(3X, A," // FMT_16 // ")") "pt2 = ", radiation%pt2 write (u, "(3X, A, I12)") "alr = ", radiation%alr end subroutine radiation_write @ %def radiation_write @ [[lambda2_gen]] is the scale used in the (log integrated) upper bounding functions. It is not equivalent to [[lambda]] which is the reference scale for the $\alpha_s$ evolution. <>= public :: process_deps_t <>= type :: process_deps_t real(default) :: lambda2_gen, sqrts integer :: n_alr integer :: alpha_power, alphas_power logical :: lab_is_cm = .true. type(pdf_data_t) :: pdf_data type(phs_identifier_t), dimension(:), allocatable :: phs_identifiers integer, dimension(:), allocatable :: alr_to_i_phs integer :: i_term_born integer, dimension(:), allocatable :: i_term_real contains <> end type process_deps_t @ %def process_deps_t @ <>= procedure :: write => process_deps_write +<>= + module subroutine process_deps_write (process_deps, unit) + class(process_deps_t), intent(in) :: process_deps + integer, intent(in), optional :: unit + end subroutine process_deps_write <>= - subroutine process_deps_write (process_deps, unit) + module subroutine process_deps_write (process_deps, unit) class(process_deps_t), intent(in) :: process_deps integer, intent(in), optional :: unit integer :: u, i u = given_output_unit (unit); if (u < 0) return write (u, "(1X,A)") "Process dependencies:" - write (u, "(3X,A," // FMT_19 // ")") "lambda2_gen = ", process_deps%lambda2_gen + write (u, "(3X,A," // FMT_19 // ")") "lambda2_gen = ", & + process_deps%lambda2_gen write (u, "(3X,A, I12)") "n_alr = ", process_deps%n_alr write (u, "(3X,A, L12)") "lab_is_cm = ", process_deps%lab_is_cm write (u, "(3X,A, I10)") "alpha_power = ", process_deps%alpha_power write (u, "(3X,A, I9)") "alphas_power = ", process_deps%alphas_power write (u, "(3X,A, I10)") "i_term_born = ", process_deps%i_term_born do i = 1, size(process_deps%i_term_real) - write (u, "(3X,A,I2,A, I6)") "i_term_real(",i,") = ", process_deps%i_term_real(i) + write (u, "(3X,A,I2,A, I6)") "i_term_real(",i,") = ", & + process_deps%i_term_real(i) end do call process_deps%pdf_data%write(u) end subroutine process_deps_write @ %def process_deps_write @ <>= public :: event_deps_t <>= type :: event_deps_t real(default) :: s_hat real(default), dimension(2) :: x_born real(default) :: s_had type(phs_point_set_t) :: p_born_cms type(phs_point_set_t) :: p_born_lab type(phs_point_set_t) :: p_real_cms type(phs_point_set_t) :: p_real_lab real(default) :: sqme_born contains <> end type event_deps_t @ %def event_deps_t @ <>= procedure :: write => event_deps_write +<>= + module subroutine event_deps_write (event_deps, unit) + class(event_deps_t), intent(in) :: event_deps + integer, intent(in), optional :: unit + end subroutine event_deps_write <>= - subroutine event_deps_write (event_deps, unit) + module subroutine event_deps_write (event_deps, unit) class(event_deps_t), intent(in) :: event_deps integer, intent(in), optional :: unit integer :: u u = given_output_unit (unit); if (u < 0) return write (u, "(1X,A)") "Event dependencies:" write (u, "(3X,A," // FMT_19 // ")") "s_hat = ", event_deps%s_hat write (u, "(3X,A," // FMT_19 // ")") "x(+) = ", event_deps%x_born(I_PLUS) write (u, "(3X,A," // FMT_19 // ")") "x(-) = ", event_deps%x_born(I_MINUS) write (u, "(3X,A," // FMT_19 // ")") "s_had = ", event_deps%s_had write (u, "(3X,A," // FMT_19 // ")") "sqme_born = ", event_deps%sqme_born end subroutine event_deps_write @ %def event_deps_write @ <>= procedure :: update => event_deps_update +<>= + module subroutine event_deps_update (event_deps, sqme_born, & + p_born, x_born, lt_lab_to_cms) + class(event_deps_t), intent(inout) :: event_deps + real(default), intent(in) :: sqme_born + type(vector4_t), dimension(:), intent(in) :: p_born + real(default), dimension(2), intent(in) :: x_born + type(lorentz_transformation_t), intent(in), optional :: lt_lab_to_cms + end subroutine event_deps_update <>= - subroutine event_deps_update (event_deps, sqme_born, p_born, x_born, lt_lab_to_cms) + module subroutine event_deps_update (event_deps, sqme_born, & + p_born, x_born, lt_lab_to_cms) class(event_deps_t), intent(inout) :: event_deps real(default), intent(in) :: sqme_born type(vector4_t), dimension(:), intent(in) :: p_born real(default), dimension(2), intent(in) :: x_born type(lorentz_transformation_t), intent(in), optional :: lt_lab_to_cms integer :: n_born event_deps%sqme_born = sqme_born n_born = size (p_born) if (debug_active (D_MATCHING)) then if (n_born /= event_deps%p_born_lab%get_n_particles (1)) then call msg_fatal & ("event_deps_update: number of born_momenta has changed") end if end if call event_deps%p_born_lab%set_momenta (1, p_born) call event_deps%set_cms (lt_lab_to_cms) event_deps%x_born = x_born event_deps%s_had = event_deps%s_hat / ( x_born(I_PLUS) * x_born(I_MINUS) ) end subroutine event_deps_update @ %def event_deps_update @ Sets the Born momenta in the CMS boosting them if necessary. <>= procedure :: set_cms => event_deps_set_cms +<>= + module subroutine event_deps_set_cms (event_deps, lt_lab_to_cms) + class(event_deps_t), intent(inout) :: event_deps + type(lorentz_transformation_t), intent(in), optional :: lt_lab_to_cms + end subroutine event_deps_set_cms <>= - subroutine event_deps_set_cms (event_deps, lt_lab_to_cms) + module subroutine event_deps_set_cms (event_deps, lt_lab_to_cms) class(event_deps_t), intent(inout) :: event_deps type(lorentz_transformation_t), intent(in), optional :: lt_lab_to_cms associate (pp => event_deps%p_born_lab%phs_point(1)) event_deps%s_hat = pp%get_msq ([1,2]) if (present (lt_lab_to_cms)) then event_deps%p_born_cms%phs_point(1) = lt_lab_to_cms * pp else event_deps%p_born_cms%phs_point(1) = pp end if end associate end subroutine event_deps_set_cms @ %def event_deps_set_cms @ <>= type :: veto_counter_t integer :: n_ubf = 0 integer :: n_first_fail = 0 integer :: n_alpha_s = 0 integer :: n_xi_max = 0 integer :: n_norm = 0 integer :: n_sqme = 0 integer :: n_veto_ubf = 0 integer :: n_veto_alpha_s = 0 integer :: n_veto_xi_max = 0 integer :: n_veto_norm = 0 integer :: n_veto_sqme = 0 integer :: n_fail_ubf = 0 contains <> end type veto_counter_t @ %def veto_counter_t @ <>= procedure :: record_ubf => veto_counter_record_ubf +<>= + pure module subroutine veto_counter_record_ubf (counter, vetoed) + class(veto_counter_t), intent(inout) :: counter + logical, intent(in) :: vetoed + end subroutine veto_counter_record_ubf <>= - pure subroutine veto_counter_record_ubf (counter, vetoed) + pure module subroutine veto_counter_record_ubf (counter, vetoed) class(veto_counter_t), intent(inout) :: counter logical, intent(in) :: vetoed counter%n_ubf = counter%n_ubf + 1 if (vetoed) counter%n_veto_ubf = counter%n_veto_ubf + 1 end subroutine veto_counter_record_ubf @ %def veto_counter_record_ubf @ <>= procedure :: record_first_fail => veto_counter_record_first_fail +<>= + module subroutine veto_counter_record_first_fail (counter) + class(veto_counter_t), intent(inout) :: counter + end subroutine veto_counter_record_first_fail <>= - subroutine veto_counter_record_first_fail (counter) + module subroutine veto_counter_record_first_fail (counter) class(veto_counter_t), intent(inout) :: counter counter%n_first_fail = counter%n_first_fail + 1 end subroutine veto_counter_record_first_fail @ %def veto_counter_record_first_fail @ <>= procedure :: record_alpha_s => veto_counter_record_alpha_s +<>= + module subroutine veto_counter_record_alpha_s (counter, vetoed) + class(veto_counter_t), intent(inout) :: counter + logical, intent(in) :: vetoed + end subroutine veto_counter_record_alpha_s <>= - subroutine veto_counter_record_alpha_s (counter, vetoed) + module subroutine veto_counter_record_alpha_s (counter, vetoed) class(veto_counter_t), intent(inout) :: counter logical, intent(in) :: vetoed counter%n_alpha_s = counter%n_alpha_s + 1 if (vetoed) counter%n_veto_alpha_s = counter%n_veto_alpha_s + 1 end subroutine veto_counter_record_alpha_s @ %def veto_counter_record_alpha_s @ <>= procedure :: record_xi_max => veto_counter_record_xi_max +<>= + module subroutine veto_counter_record_xi_max (counter, vetoed) + class(veto_counter_t), intent(inout) :: counter + logical, intent(in) :: vetoed + end subroutine veto_counter_record_xi_max <>= - subroutine veto_counter_record_xi_max (counter, vetoed) + module subroutine veto_counter_record_xi_max (counter, vetoed) class(veto_counter_t), intent(inout) :: counter logical, intent(in) :: vetoed counter%n_xi_max = counter%n_xi_max + 1 if (vetoed) counter%n_veto_xi_max = counter%n_veto_xi_max + 1 end subroutine veto_counter_record_xi_max @ %def veto_counter_record_xi_max @ <>= procedure :: record_norm => veto_counter_record_norm +<>= + module subroutine veto_counter_record_norm (counter, vetoed) + class(veto_counter_t), intent(inout) :: counter + logical, intent(in) :: vetoed + end subroutine veto_counter_record_norm <>= - subroutine veto_counter_record_norm (counter, vetoed) + module subroutine veto_counter_record_norm (counter, vetoed) class(veto_counter_t), intent(inout) :: counter logical, intent(in) :: vetoed counter%n_norm = counter%n_norm + 1 if (vetoed) counter%n_veto_norm = counter%n_veto_norm + 1 end subroutine veto_counter_record_norm @ %def veto_counter_record_norm @ <>= procedure :: record_sqme => veto_counter_record_sqme +<>= + module subroutine veto_counter_record_sqme (counter, vetoed) + class(veto_counter_t), intent(inout) :: counter + logical, intent(in) :: vetoed + end subroutine veto_counter_record_sqme <>= - subroutine veto_counter_record_sqme (counter, vetoed) + module subroutine veto_counter_record_sqme (counter, vetoed) class(veto_counter_t), intent(inout) :: counter logical, intent(in) :: vetoed counter%n_sqme = counter%n_sqme + 1 if (vetoed) counter%n_veto_sqme = counter%n_veto_sqme + 1 end subroutine veto_counter_record_sqme @ %def veto_counter_record_sqme @ <>= procedure :: record_ubf_fail => veto_counter_record_ubf_fail +<>= + module subroutine veto_counter_record_ubf_fail (counter) + class(veto_counter_t), intent(inout) :: counter + end subroutine veto_counter_record_ubf_fail <>= - subroutine veto_counter_record_ubf_fail (counter) + module subroutine veto_counter_record_ubf_fail (counter) class(veto_counter_t), intent(inout) :: counter counter%n_fail_ubf = counter%n_fail_ubf + 1 end subroutine veto_counter_record_ubf_fail @ %def veto_counter_record_ubf_fail @ This routine fills the POWHEG veto log file with content showing how many events have been vetoed in which step of the sudakov veto algorithm. <>= procedure :: write => veto_counter_write +<>= + module subroutine veto_counter_write (counter, unit) + class(veto_counter_t), intent(in) :: counter + integer, intent(in), optional :: unit + end subroutine veto_counter_write <>= - subroutine veto_counter_write (counter, unit) + module subroutine veto_counter_write (counter, unit) class(veto_counter_t), intent(in) :: counter integer, intent(in), optional :: unit integer :: u u = given_output_unit (unit); if (u < 0) return write (u, "(A29,I10)") "Nr. of ubf-veto calls: ", counter%n_ubf write (u, "(A29,I10)") "Nr. of ubf-vetos: ", counter%n_veto_ubf if (counter%n_ubf > 0) & write (u, "(A29,F9.2,A1)") "Fraction of vetoed points: ", & one * counter%n_veto_ubf / counter%n_ubf * 100, "%" call write_separator (u) write (u, "(A29,I10)") "Nr. of alpha_s-veto calls: ", counter%n_alpha_s write (u, "(A29,I10)") "Nr. of alpha_s-vetos: ", counter%n_veto_alpha_s if (counter%n_alpha_s > 0) & write (u, "(A29,F9.2,A1)") "Fraction of vetoed points: ", & one * counter%n_veto_alpha_s / counter%n_alpha_s * 100, "%" call write_separator (u) write (u, "(A29,I10)") "Nr. of xi_max-veto calls: ", counter%n_xi_max write (u, "(A29,I10)") "Nr. of xi_max-vetos: ", counter%n_veto_xi_max if (counter%n_alpha_s > 0) & write (u, "(A29,F9.2,A1)") "Fraction of vetoed points: ", & one * counter%n_veto_xi_max / counter%n_xi_max * 100, "%" call write_separator (u) write (u, "(A29,I10)") "Nr. of norm-veto calls: ", counter%n_norm write (u, "(A29,I10)") "Nr. of norm-vetos: ", counter%n_veto_norm if (counter%n_norm > 0) & write (u, "(A29,F9.2,A1)") "Fraction of vetoed points: ", & one * counter%n_veto_norm / counter%n_norm * 100, "%" call write_separator (u) write (u, "(A29,I10)") "Nr. of sqme-veto calls: ", counter%n_sqme write (u, "(A29,I10)") "Nr. of sqme-vetos: ", counter%n_veto_sqme if (counter%n_sqme > 0) & write (u, "(A29,F9.2,A1)") "Fraction of vetoed points: ", & one * counter%n_veto_sqme / counter%n_sqme * 100, "%" call write_separator (u) write (u, "(A29,I10)") "Nr. of upper-bound failures: ", & counter%n_fail_ubf if (counter%n_sqme > 0) & write (u, "(A29,F9.2,A1)") "Fraction of failed points: ", & one * counter%n_fail_ubf / counter%n_sqme * 100, "%" end subroutine veto_counter_write @ %def veto_counter_write @ \subsection{Upper bounding functions and [[sudakov]]s} \subsubsection{Abstract version} This contains the pieces that depend on the radiation region $\alpha_r$ <>= public :: sudakov_t <>= type, abstract, extends (solver_function_t) :: sudakov_t type(process_deps_t), pointer :: process_deps => null() type(event_deps_t), pointer :: event_deps => null() type(powheg_settings_t), pointer :: powheg_settings => null() type(phs_fks_generator_t), pointer :: phs_fks_generator => null() type(qcd_t) :: qcd class(rng_t), pointer :: rng => null() real(default) :: xi2_max = zero real(default) :: norm_max = zero real(default) :: current_pt2_max = zero real(default) :: sum_log_rands = zero real(default) :: random = zero type(veto_counter_t) :: veto_counter integer :: i_phs = 0 contains <> end type sudakov_t @ %def sudakov_t @ <>= procedure :: write => sudakov_write +<>= + module subroutine sudakov_write (sudakov, unit) + class(sudakov_t), intent(in) :: sudakov + integer, intent(in), optional :: unit + end subroutine sudakov_write <>= - subroutine sudakov_write (sudakov, unit) + module subroutine sudakov_write (sudakov, unit) class(sudakov_t), intent(in) :: sudakov integer, intent(in), optional :: unit integer :: u u = given_output_unit (unit); if (u < 0) return write (u, "(3X,A," // FMT_19 // ")") "xi2_max = ", sudakov%xi2_max write (u, "(3X,A," // FMT_19 // ")") "norm_max = ", sudakov%norm_max write (u, "(3X,A," // FMT_19 // ")") & "current_pt2_max = ", sudakov%current_pt2_max write (u, "(3X,A," // FMT_19 // ")") "sum_log_rands = ", sudakov%sum_log_rands write (u, "(3X,A," // FMT_19 // ")") "random = ", sudakov%random end subroutine sudakov_write @ %def sudakov_write @ To allow for arrays of this class <>= public :: sudakov_wrapper_t <>= type :: sudakov_wrapper_t class(sudakov_t), allocatable :: s end type sudakov_wrapper_t @ %def sudakov_wrapper_t @ <>= procedure :: init => sudakov_init +<>= + module subroutine sudakov_init (sudakov, process_deps, event_deps, & + powheg_settings, qcd, phs_fks_generator, rng) + class(sudakov_t), intent(out) :: sudakov + type(process_deps_t), target, intent(in) :: process_deps + type(event_deps_t), target, intent(in) :: event_deps + type(powheg_settings_t), target, intent(in) :: powheg_settings + type(qcd_t), intent(in) :: qcd + type(phs_fks_generator_t), target, intent(in) :: phs_fks_generator + class(rng_t), target, intent(in), optional :: rng + end subroutine sudakov_init <>= - subroutine sudakov_init (sudakov, process_deps, event_deps, & + module subroutine sudakov_init (sudakov, process_deps, event_deps, & powheg_settings, qcd, phs_fks_generator, rng) class(sudakov_t), intent(out) :: sudakov type(process_deps_t), target, intent(in) :: process_deps type(event_deps_t), target, intent(in) :: event_deps type(powheg_settings_t), target, intent(in) :: powheg_settings type(qcd_t), intent(in) :: qcd type(phs_fks_generator_t), target, intent(in) :: phs_fks_generator class(rng_t), target, intent(in), optional :: rng sudakov%process_deps => process_deps sudakov%event_deps => event_deps sudakov%powheg_settings => powheg_settings sudakov%qcd = qcd sudakov%phs_fks_generator => phs_fks_generator if (present (rng)) sudakov%rng => rng end subroutine sudakov_init @ %def sudakov_init @ This has to be done after the grids are initialized. <>= procedure :: set_normalization => sudakov_set_normalization +<>= + pure module subroutine sudakov_set_normalization (sudakov, norm_max) + class(sudakov_t), intent(inout) :: sudakov + real(default), intent(in) :: norm_max + end subroutine sudakov_set_normalization <>= - pure subroutine sudakov_set_normalization (sudakov, norm_max) + pure module subroutine sudakov_set_normalization (sudakov, norm_max) class(sudakov_t), intent(inout) :: sudakov real(default), intent(in) :: norm_max sudakov%norm_max = norm_max end subroutine sudakov_set_normalization @ %def sudakov_set_normalization @ <>= procedure :: update_xi2_max => sudakov_update_xi2_max +<>= + pure module subroutine sudakov_update_xi2_max (sudakov, xi2_max) + class(sudakov_t), intent(inout) :: sudakov + real(default), intent(in) :: xi2_max + end subroutine sudakov_update_xi2_max <>= - pure subroutine sudakov_update_xi2_max (sudakov, xi2_max) + pure module subroutine sudakov_update_xi2_max (sudakov, xi2_max) class(sudakov_t), intent(inout) :: sudakov real(default), intent(in) :: xi2_max sudakov%xi2_max = xi2_max end subroutine sudakov_update_xi2_max @ %def sudakov_update_xi2_max @ [[upper_bound_func]] does \emph{not} contain the normalization $N$ which is given by the grids. In the notation of [[1002.2581]], it is thus $\frac 1 N U(\xi,y)$ <>= procedure (sudakov_upper_bound_func), deferred :: upper_bound_func <>= abstract interface pure function sudakov_upper_bound_func (sudakov, xi, y, alpha_s) result (u) import real(default) :: u class(sudakov_t), intent(in) :: sudakov real(default), intent(in) :: xi, y, alpha_s end function sudakov_upper_bound_func end interface @ %def sudakov_upper_bound_func @ Similar to the [[upper_bound_func]], this is $-\frac 1 N \log\Delta(p_T^2)$ where \begin{equation} \Delta^{(U)} (p_T) = \exp\left[-\int U(\xi,y)\theta(k_T - p_T) \; d\xi \, dy \, d\phi \right] \end{equation} <>= procedure (sudakov_log_integrated_ubf), deferred :: log_integrated_ubf <>= abstract interface pure function sudakov_log_integrated_ubf (sudakov, pt2) result (y) import real(default) :: y class(sudakov_t), intent(in) :: sudakov real(default), intent(in) :: pt2 end function sudakov_log_integrated_ubf end interface @ %def sudakov_log_integrated_ubf @ <>= procedure (sudakov_generate_xi_and_y_and_phi), deferred :: generate_xi_and_y_and_phi <>= abstract interface subroutine sudakov_generate_xi_and_y_and_phi (sudakov, r) import class(sudakov_t), intent(inout) :: sudakov type(radiation_t), intent(inout) :: r end subroutine sudakov_generate_xi_and_y_and_phi end interface @ %def sudakov_generate_xi_and_y_and_phi @ Generating $\phi$ can be performed in a generic way for all UBF types. <>= procedure :: generate_phi => sudakov_generate_phi +<>= + module subroutine sudakov_generate_phi (sudakov, r) + class(sudakov_t), intent(inout) :: sudakov + type(radiation_t), intent(inout) :: r + end subroutine sudakov_generate_phi <>= - subroutine sudakov_generate_phi (sudakov, r) + module subroutine sudakov_generate_phi (sudakov, r) class(sudakov_t), intent(inout) :: sudakov type(radiation_t), intent(inout) :: r call sudakov%rng%generate (sudakov%random) r%phi = sudakov%random * twopi end subroutine sudakov_generate_phi @ %def sudakov_generate_phi @ <>= procedure (sudakov_kt2), deferred :: kt2 <>= abstract interface function sudakov_kt2 (sudakov, xi, y) result (kt2) import real(default) :: kt2 class(sudakov_t), intent(in) :: sudakov real(default), intent(in) :: xi, y end function sudakov_kt2 end interface @ %def sudakov_kt2 @ <>= procedure (sudakov_kt2_max), deferred :: kt2_max <>= abstract interface pure function sudakov_kt2_max (sudakov) result (kt2_max) import real(default) :: kt2_max class(sudakov_t), intent(in) :: sudakov end function sudakov_kt2_max end interface @ %def sudakov_kt2_max @ Veto routine for an overestimation of the UBF performed to simplify the log integrated UBF <>= procedure (sudakov_reweight_ubf), deferred :: reweight_ubf <>= abstract interface function sudakov_reweight_ubf (sudakov, pt2) result (accepted) import logical :: accepted class(sudakov_t), intent(inout) :: sudakov real(default), intent(in) :: pt2 end function sudakov_reweight_ubf end interface @ %def sudakov_reweight_ubf @ Veto routine to correct for an overestimation of $\xi_\text{max}$ performed in the case of massive emitters. <>= procedure (sudakov_reweight_xi_max), deferred :: reweight_xi_max <>= abstract interface function sudakov_reweight_xi_max (sudakov, xi) result (accepted) import logical :: accepted class(sudakov_t), intent(in) :: sudakov real(default), intent(in) :: xi end function sudakov_reweight_xi_max end interface @ %def sudakov_reweight_xi_max @ In the generation of $p_T^2$ via [[log_integrated_ubf]], we use the simplified version $\alpha_s^\text{rad}$ while the grids take the improved version. <>= procedure :: alpha_s => sudakov_alpha_s +<>= + module function sudakov_alpha_s (sudakov, kT2, use_correct) result (a) + real(default) :: a + class(sudakov_t), intent(in) :: sudakov + real(default), intent(in) :: kT2 + logical, intent(in), optional :: use_correct + end function sudakov_alpha_s <>= - function sudakov_alpha_s (sudakov, kT2, use_correct) result (a) + module function sudakov_alpha_s (sudakov, kT2, use_correct) result (a) real(default) :: a class(sudakov_t), intent(in) :: sudakov real(default), intent(in) :: kT2 logical, intent(in), optional :: use_correct logical :: yorn yorn = .false.; if (present (use_correct)) yorn = use_correct if (yorn) then a = get_alpha_s (sudakov%qcd, kT2) else a = sudakov%alpha_s_rad (kT2) end if end function sudakov_alpha_s @ %def sudakov_alpha_s @ To generate the transverse momentum, we have to solve the equation \begin{equation*} \frac{\Delta^{(U)}(p_T)}{\Delta^{(U)}(p_T^{\mathrm{max}})} = r_1 \end{equation*} iteratively for $p_T$. In log space and with initially $\Delta^{(U)}(p_T^{\mathrm{max}}) = 1$, this is \begin{equation*} \log\Delta^{(U)}(p_T) = \log r_1 \end{equation*} If no solution is found, we set $p_T = p_T^{\mathrm{min}}$ to exit the event generation loop and generate an emissionless event. If a solution is found but the current emission is not accepted, we set $p_T = p_T^{\mathrm{max}}$ and thus in the next step it is $\Delta^{(U)}(p_T^{\mathrm{max}}) = r_1$, so that we have to solve the equation \begin{equation*} \log\Delta^{(U)}(p_T) = \log r_1 + \log r_2 \end{equation*} using another random number $r_2$. We use [[sum_log_rands]] to remember the sum of the logarithms of all previous random numbers used for this event. <>= procedure :: generate_pt2 => sudakov_generate_pt2 +<>= + module function sudakov_generate_pt2 (sudakov) result (pt2) + real(default) :: pt2 + class(sudakov_t), intent(inout) :: sudakov + end function sudakov_generate_pt2 <>= - function sudakov_generate_pt2 (sudakov) result (pt2) + module function sudakov_generate_pt2 (sudakov) result (pt2) real(default) :: pt2 class(sudakov_t), intent(inout) :: sudakov logical :: success success = .false. if (sudakov%current_pt2_max > sudakov%powheg_settings%pt2_min) then call sudakov%rng%generate (sudakov%random) sudakov%sum_log_rands = sudakov%sum_log_rands + log(sudakov%random) pt2 = solve_interval (sudakov, & sudakov%powheg_settings%pt2_min, & sudakov%current_pt2_max, success, & 0.001_default) end if if (.not. success) then pt2 = sudakov%powheg_settings%pt2_min end if end function sudakov_generate_pt2 @ %def sudakov_generate_pt2 @ This could be activated [[if (debug_active (MATCHING))]]. <>= procedure :: check_solution_interval => sudakov_check_solution_interval +<>= + module subroutine sudakov_check_solution_interval (sudakov) + class(sudakov_t), intent(inout) :: sudakov + end subroutine sudakov_check_solution_interval <>= - subroutine sudakov_check_solution_interval (sudakov) + module subroutine sudakov_check_solution_interval (sudakov) class(sudakov_t), intent(inout) :: sudakov real(default) :: r real(default), parameter :: dr = 0.05 real(default) :: pt2 logical :: success r = 0._default do r = r + dr sudakov%random = r pt2 = solve_interval (sudakov, & sudakov%powheg_settings%pt2_min, & sudakov%current_pt2_max, success, & 0.001_default) if (success) then print *, 'r: ', r, ' zero found' else print *, 'r: ', r, 'no zero found' end if if (r >= 1._default) exit end do end subroutine sudakov_check_solution_interval @ %def sudakov_check_solution_interval @ Generates the FKS variables $(\xi, y,\phi)$ and sets the emission hardness scale. In debug mode, we assert that the [[pt2]] saved in the [[radiation_t]] and the massive hardness scale are equal. This should hold as we compute [[kt2]] with the [[y]] we got from generating [[pt2]] so that this is a full circle. It then initiates the generation of the transverse momentum followed by a sequence of veto steps to get the correct distribution. <>= procedure :: generate_emission => sudakov_generate_emission +<>= + module subroutine sudakov_generate_emission (sudakov, r) + class(sudakov_t), intent(inout) :: sudakov + type(radiation_t), intent(inout) :: r + end subroutine sudakov_generate_emission <>= - subroutine sudakov_generate_emission (sudakov, r) + module subroutine sudakov_generate_emission (sudakov, r) class(sudakov_t), intent(inout) :: sudakov type(radiation_t), intent(inout) :: r logical :: accepted sudakov%current_pt2_max = r%pt2 if (debug_on) call msg_debug2 (D_MATCHING, "sudakov_generate_emission") if (debug_on) call msg_debug2 (D_MATCHING, "sqrt (sudakov%current_pt2_max)", & sqrt (sudakov%current_pt2_max)) if (debug_on) call msg_debug2 (D_MATCHING, "sudakov%sum_log_rands", sudakov%sum_log_rands) LOOP_UNTIL_ACCEPTED: do if (signal_is_pending ()) return r%valid = .false. r%pt2 = sudakov%generate_pt2 () if (debug_on) call msg_debug2 (D_MATCHING, "sudakov_generate_emission: after generate_pt2") if (debug_on) call msg_debug2 (D_MATCHING, "sqrt (r%pt2)", sqrt (r%pt2)) if (debug_on) call msg_debug2 (D_MATCHING, "sudakov%sum_log_rands", sudakov%sum_log_rands) if (r%pt2 <= sudakov%powheg_settings%pt2_min) then exit end if accepted = sudakov%reweight_ubf (r%pt2) call sudakov%veto_counter%record_ubf (.not. accepted) if (.not. accepted) then sudakov%current_pt2_max = r%pt2 cycle end if accepted = sudakov%reweight_alpha_s (r%pt2) call sudakov%veto_counter%record_alpha_s (.not. accepted) if (.not. accepted) then sudakov%current_pt2_max = r%pt2 cycle end if call sudakov%generate_xi_and_y_and_phi (r) accepted = sudakov%reweight_xi_max (r%xi) call sudakov%veto_counter%record_xi_max (.not. accepted) if (.not. accepted) then sudakov%current_pt2_max = r%pt2 cycle end if if (debug_active (D_MATCHING)) then call assert_equal (OUTPUT_UNIT, r%pt2, & sudakov%kt2 (r%xi, r%y), & "sudakov_generate_xi_and_y_and_phi: pt2 inconsistency") end if r%valid = .true. exit end do LOOP_UNTIL_ACCEPTED end subroutine sudakov_generate_emission @ %def sudakov_generate_emission @ Evaluates the Sudakov as needed for [[solve_interval]] to generate a $p_T$ value. <>= procedure :: evaluate => sudakov_evaluate +<>= + module function sudakov_evaluate (solver_f, x) result (f) + complex(default) :: f + class(sudakov_t), intent(in) :: solver_f + real(default), intent(in) :: x + end function sudakov_evaluate <>= - function sudakov_evaluate (solver_f, x) result (f) + module function sudakov_evaluate (solver_f, x) result (f) complex(default) :: f class(sudakov_t), intent(in) :: solver_f real(default), intent(in) :: x - f = solver_f%sum_log_rands + solver_f%norm_max * solver_f%log_integrated_ubf (x) + f = solver_f%sum_log_rands + & + solver_f%norm_max * solver_f%log_integrated_ubf (x) end function sudakov_evaluate @ %def sudakov_evaluate @ <>= procedure :: associated_emitter => sudakov_associated_emitter +<>= + elemental module function sudakov_associated_emitter & + (sudakov) result (emitter) + integer :: emitter + class(sudakov_t), intent(in) :: sudakov + end function sudakov_associated_emitter <>= - elemental function sudakov_associated_emitter (sudakov) result (emitter) + elemental module function sudakov_associated_emitter & + (sudakov) result (emitter) integer :: emitter class(sudakov_t), intent(in) :: sudakov emitter = sudakov%process_deps%phs_identifiers(sudakov%i_phs)%emitter end function sudakov_associated_emitter @ %def sudakov_associated_emitter @ <>= procedure :: set_i_phs => sudakov_set_i_phs +<>= + module subroutine sudakov_set_i_phs (sudakov, alr) + class(sudakov_t), intent(inout) :: sudakov + integer, intent(in) :: alr + end subroutine sudakov_set_i_phs <>= - subroutine sudakov_set_i_phs (sudakov, alr) + module subroutine sudakov_set_i_phs (sudakov, alr) class(sudakov_t), intent(inout) :: sudakov integer, intent(in) :: alr sudakov%i_phs = sudakov%process_deps%alr_to_i_phs(alr) end subroutine sudakov_set_i_phs @ %def sudakov_set_i_phs @ \subsubsection{Simple FSR} This corresponds to Appendix C of [[1002.2581]]. <>= public :: sudakov_simple_fsr_t <>= type, extends (sudakov_t) :: sudakov_simple_fsr_t contains <> end type sudakov_simple_fsr_t @ %def sudakov_simple_fsr_t @ The simplest upper bounding function for final-state radiation is \begin{equation} \mathtt{upper\_bound\_func} = \frac {U(\xi,y)} N = \frac {\alpha_s}{\xi (1-y)} \end{equation} <>= procedure :: upper_bound_func => sudakov_simple_fsr_upper_bound_func +<>= + pure module function sudakov_simple_fsr_upper_bound_func & + (sudakov, xi, y, alpha_s) result (u) + real(default) :: u + class(sudakov_simple_fsr_t), intent(in) :: sudakov + real(default), intent(in) :: xi, y, alpha_s + end function sudakov_simple_fsr_upper_bound_func <>= - pure function sudakov_simple_fsr_upper_bound_func (sudakov, xi, y, alpha_s) result (u) + pure module function sudakov_simple_fsr_upper_bound_func & + (sudakov, xi, y, alpha_s) result (u) real(default) :: u class(sudakov_simple_fsr_t), intent(in) :: sudakov real(default), intent(in) :: xi, y, alpha_s u = alpha_s / (xi * (1 - y)) end function sudakov_simple_fsr_upper_bound_func @ %def sudakov_simple_fsr_upper_bound_func @ The above upper bounding function corresponds to the transverse momentum scale \begin{equation} k_T^2 = \frac{s}{2} \xi^2 (1-y). \end{equation} <>= procedure :: kt2 => sudakov_simple_fsr_kt2 +<>= + module function sudakov_simple_fsr_kt2 (sudakov, xi, y) result (kt2) + real(default) :: kt2 + class(sudakov_simple_fsr_t), intent(in) :: sudakov + real(default), intent(in) :: xi, y + end function sudakov_simple_fsr_kt2 <>= - function sudakov_simple_fsr_kt2 (sudakov, xi, y) result (kt2) + module function sudakov_simple_fsr_kt2 (sudakov, xi, y) result (kt2) real(default) :: kt2 class(sudakov_simple_fsr_t), intent(in) :: sudakov real(default), intent(in) :: xi, y kt2 = sudakov%phs_fks_generator%real_kinematics%kt2 & (sudakov%i_phs, sudakov%associated_emitter (), UBF_FSR_SIMPLE, xi, y) end function sudakov_simple_fsr_kt2 @ %def sudakov_simple_fsr_kt2 @ For massless emitters, the upper bound on the radiated energy is \begin{equation*} t_{\mathrm{max}} = \xi_{\mathrm{max}}^2 \hat{s} \end{equation*} We use this as largest possible scale for the radiation. <>= procedure :: kt2_max => sudakov_simple_fsr_kt2_max +<>= + pure module function sudakov_simple_fsr_kt2_max (sudakov) result (kt2_max) + real(default) :: kt2_max + class(sudakov_simple_fsr_t), intent(in) :: sudakov + end function sudakov_simple_fsr_kt2_max <>= - pure function sudakov_simple_fsr_kt2_max (sudakov) result (kt2_max) - real(default) :: kt2_max - class(sudakov_simple_fsr_t), intent(in) :: sudakov - real(default) :: s_hat - s_hat = sudakov%event_deps%s_hat - kt2_max = sudakov%xi2_max * s_hat + pure module function sudakov_simple_fsr_kt2_max (sudakov) result (kt2_max) + real(default) :: kt2_max + class(sudakov_simple_fsr_t), intent(in) :: sudakov + real(default) :: s_hat + s_hat = sudakov%event_deps%s_hat + kt2_max = sudakov%xi2_max * s_hat end function sudakov_simple_fsr_kt2_max @ %def sudakov_simple_fsr_kt2_max @ This is \begin{equation} - \frac{\log{\Delta^{(U)}}(p_T)}{N} = \frac\pi{b_0} \theta\left(\xi_\text{max}^2-\frac{p_T^2}s\right) \left[\log{\frac{\xi^2_\text{max}s}{\Lambda^2}} \log{\frac{\log{{\xi^2_\text{max}s}/{\Lambda^2}}} {\log{p_T^2/\Lambda^2}}} - \log{\frac{\xi^2_\text{max}s}{p_T^2}}\right] \end{equation} with $p_\text{T,max}^2=\xi_\text{max}^2 s$. <>= procedure :: log_integrated_ubf => sudakov_simple_fsr_log_integrated_ubf +<>= + pure module function sudakov_simple_fsr_log_integrated_ubf & + (sudakov, pt2) result (log_sudakov) + real(default) :: log_sudakov + class(sudakov_simple_fsr_t), intent(in) :: sudakov + real(default), intent(in) :: pt2 + end function sudakov_simple_fsr_log_integrated_ubf <>= - pure function sudakov_simple_fsr_log_integrated_ubf (sudakov, pt2) result (log_sudakov) + pure module function sudakov_simple_fsr_log_integrated_ubf & + (sudakov, pt2) result (log_sudakov) real(default) :: log_sudakov class(sudakov_simple_fsr_t), intent(in) :: sudakov real(default), intent(in) :: pt2 real(default) :: kt2_max, kt2_maxl, pt2l logical :: within_boundaries within_boundaries = pt2 <= sudakov%kt2_max () & .and. pt2 >= sudakov%powheg_settings%pt2_min if (within_boundaries) then kt2_max = sudakov%kt2_max () kt2_maxl = kt2_max / sudakov%process_deps%lambda2_gen pt2l = pt2 / sudakov%process_deps%lambda2_gen log_sudakov = pi / b0rad * (log (kt2_maxl) * & log (log (kt2_maxl) / log (pt2l)) - & log (kt2_max / pt2)) else log_sudakov = 0 end if end function sudakov_simple_fsr_log_integrated_ubf @ %def sudakov_simple_fsr_log_integrated_ubf @ No further veto needed for this upper bounding function. <>= procedure :: reweight_ubf => sudakov_simple_fsr_reweight_ubf +<>= + module function sudakov_simple_fsr_reweight_ubf & + (sudakov, pt2) result (accepted) + logical :: accepted + class(sudakov_simple_fsr_t), intent(inout) :: sudakov + real(default), intent(in) :: pt2 + end function sudakov_simple_fsr_reweight_ubf <>= - function sudakov_simple_fsr_reweight_ubf (sudakov, pt2) result (accepted) + module function sudakov_simple_fsr_reweight_ubf & + (sudakov, pt2) result (accepted) logical :: accepted class(sudakov_simple_fsr_t), intent(inout) :: sudakov real(default), intent(in) :: pt2 accepted = .true. end function sudakov_simple_fsr_reweight_ubf @ %def sudakov_simple_fsr_reweight_ubf @ <>= procedure :: reweight_xi_max => sudakov_simple_fsr_reweight_xi_max +<>= + module function sudakov_simple_fsr_reweight_xi_max & + (sudakov, xi) result (accepted) + logical :: accepted + class(sudakov_simple_fsr_t), intent(in) :: sudakov + real(default), intent(in) :: xi + end function sudakov_simple_fsr_reweight_xi_max <>= - function sudakov_simple_fsr_reweight_xi_max (sudakov, xi) result (accepted) - logical :: accepted - class(sudakov_simple_fsr_t), intent(in) :: sudakov - real(default), intent(in) :: xi - accepted = .true. + module function sudakov_simple_fsr_reweight_xi_max & + (sudakov, xi) result (accepted) + logical :: accepted + class(sudakov_simple_fsr_t), intent(in) :: sudakov + real(default), intent(in) :: xi + accepted = .true. end function sudakov_simple_fsr_reweight_xi_max @ %def sudakov_simple_fsr_reweight_xi_max @ This depends on the choice of $p_T$ and is tested in the assertion. <>= - procedure :: generate_xi_and_y_and_phi => sudakov_simple_fsr_generate_xi_and_y_and_phi + procedure :: generate_xi_and_y_and_phi => & + sudakov_simple_fsr_generate_xi_and_y_and_phi +<>= + module subroutine sudakov_simple_fsr_generate_xi_and_y_and_phi (sudakov, r) + class(sudakov_simple_fsr_t), intent(inout) :: sudakov + type(radiation_t), intent(inout) :: r + end subroutine sudakov_simple_fsr_generate_xi_and_y_and_phi <>= - subroutine sudakov_simple_fsr_generate_xi_and_y_and_phi (sudakov, r) + module subroutine sudakov_simple_fsr_generate_xi_and_y_and_phi (sudakov, r) class(sudakov_simple_fsr_t), intent(inout) :: sudakov type(radiation_t), intent(inout) :: r real(default) :: s s = sudakov%event_deps%s_hat call sudakov%generate_xi (r) r%y = one - (two * r%pt2) / (s * r%xi**2) call sudakov%generate_phi (r) end subroutine sudakov_simple_fsr_generate_xi_and_y_and_phi @ %def sudakov_generate_xi_and_y_and_phi @ We generate $\xi \in [\frac{p_\text{T}}{\sqrt{s}}, \xi_\text{max}]$ with a density $1 / \xi$ i.e. uniformly in $\log(\xi)$. <>= procedure :: generate_xi => sudakov_simple_fsr_generate_xi +<>= + module subroutine sudakov_simple_fsr_generate_xi (sudakov, r) + class(sudakov_simple_fsr_t), intent(inout) :: sudakov + type(radiation_t), intent(inout) :: r + end subroutine sudakov_simple_fsr_generate_xi <>= - subroutine sudakov_simple_fsr_generate_xi (sudakov, r) + module subroutine sudakov_simple_fsr_generate_xi (sudakov, r) class(sudakov_simple_fsr_t), intent(inout) :: sudakov type(radiation_t), intent(inout) :: r real(default) :: s, xi2_max s = sudakov%event_deps%s_hat xi2_max = sudakov%xi2_max call sudakov%rng%generate (sudakov%random) r%xi = exp (((one - sudakov%random) * log (r%pt2 / s) + & sudakov%random * log (xi2_max)) / two) end subroutine sudakov_simple_fsr_generate_xi @ %def sudakov_simple_fsr_generate_xi @ \subsubsection{Dijet production at lepton colliders} In the POWHEG method paper, $e^+e^-\to q\bar{q}$ requires a UBF different from the simple one to account for an additional divergence for $(\xi,y) \to (1,-1)$. <>= public :: sudakov_eeqq_fsr_t <>= type, extends (sudakov_t) :: sudakov_eeqq_fsr_t contains <> end type sudakov_eeqq_fsr_t @ %def sudakov_eeqq_fsr_t @ This $k_T$ measure is the same as in the case for simple FSR up to $\mathcal{O}(\theta^4)$ when $y=\cos(\theta)$. It differs by just a factor of $2$ w.r.t. [[sudakov_simple_fsr_kt2]]. <>= procedure :: kt2 => sudakov_eeqq_fsr_kt2 +<>= + module function sudakov_eeqq_fsr_kt2 (sudakov, xi, y) result (kt2) + real(default) :: kt2 + class(sudakov_eeqq_fsr_t), intent(in) :: sudakov + real(default), intent(in) :: xi, y + end function sudakov_eeqq_fsr_kt2 <>= - function sudakov_eeqq_fsr_kt2 (sudakov, xi, y) result (kt2) + module function sudakov_eeqq_fsr_kt2 (sudakov, xi, y) result (kt2) real(default) :: kt2 class(sudakov_eeqq_fsr_t), intent(in) :: sudakov real(default), intent(in) :: xi, y kt2 = sudakov%phs_fks_generator%real_kinematics%kt2 & - (sudakov%i_phs, sudakov%associated_emitter(), UBF_FSR_MASSLESS_RECOIL, xi, y) + (sudakov%i_phs, sudakov%associated_emitter(), & + UBF_FSR_MASSLESS_RECOIL, xi, y) end function sudakov_eeqq_fsr_kt2 @ %def sudakov_eeqq_fsr_kt2 @ For [[eeqq]], the shower starts at $k_\text{max} = \frac{q^0}{2}$. (c.f. eq. (7.86) in [0709.2092].) <>= procedure :: kt2_max => sudakov_eeqq_fsr_kt2_max +<>= + pure module function sudakov_eeqq_fsr_kt2_max (sudakov) result (kt2_max) + real(default) :: kt2_max + class(sudakov_eeqq_fsr_t), intent(in) :: sudakov + end function sudakov_eeqq_fsr_kt2_max <>= - pure function sudakov_eeqq_fsr_kt2_max (sudakov) result (kt2_max) - real(default) :: kt2_max - class(sudakov_eeqq_fsr_t), intent(in) :: sudakov - kt2_max = 0.25_default * sudakov%event_deps%s_hat + pure module function sudakov_eeqq_fsr_kt2_max (sudakov) result (kt2_max) + real(default) :: kt2_max + class(sudakov_eeqq_fsr_t), intent(in) :: sudakov + kt2_max = 0.25_default * sudakov%event_deps%s_hat end function sudakov_eeqq_fsr_kt2_max @ %def sudakov_eeqq_fsr_kt2_max @ This covers also the singularity at $(\xi,y)\to(1,-1)$ that occurs for a massless recoiling system. <>= procedure :: upper_bound_func => sudakov_eeqq_fsr_upper_bound_func +<>= + pure module function sudakov_eeqq_fsr_upper_bound_func & + (sudakov, xi, y, alpha_s) result (u) + real(default) :: u + class(sudakov_eeqq_fsr_t), intent(in) :: sudakov + real(default), intent(in) :: xi, y, alpha_s + end function sudakov_eeqq_fsr_upper_bound_func <>= - pure function sudakov_eeqq_fsr_upper_bound_func (sudakov, xi, y, alpha_s) result (u) + pure module function sudakov_eeqq_fsr_upper_bound_func & + (sudakov, xi, y, alpha_s) result (u) real(default) :: u class(sudakov_eeqq_fsr_t), intent(in) :: sudakov real(default), intent(in) :: xi, y, alpha_s u = alpha_s / (xi * (1 - y**2)) end function sudakov_eeqq_fsr_upper_bound_func @ %def sudakov_eeqq_fsr_upper_bound_func @ The logarithmic integrated UBF in the [[eeqq]] case as explained in [0709.2092], eq. (7.98), is given by \begin{align} -\frac{1}{N} \log\left[ \Delta^{U}(p_T) \right] &= 2 \pi \int_0^{k_\text{max}} \alpha_s(k_T) \log \left[ \frac{4k_\text{max}^2}{k_T^2} \right] \theta(k_T - p_T) \frac{dk_T}{k_T} \nonumber \\ &= \frac{\pi}{b_0} \left( \log \left[ \frac{4k_\text{max}^2}{k_T^2} \right] \log \left[ \frac{\log(k_\text{max}^2/\Lambda^2)}{\log(p_T^2/\Lambda^2)} \right] - \log \left[ \frac{k_\text{max}^2}{p_T^2} \right] \right) \end{align} with $k_\text{max}=\frac{q^0}{2}$. We see that for these logarithms to be well defined, we need to set a lower bound for the transverse momentum $p_{T,\text{min}} > \Lambda =$ [[lambda2_gen]]. % (PS 2021-05-05) Status: unvalidated <>= procedure :: log_integrated_ubf => sudakov_eeqq_fsr_log_integrated_ubf +<>= + pure module function sudakov_eeqq_fsr_log_integrated_ubf & + (sudakov, pt2) result (log_sudakov) + real(default) :: log_sudakov + class(sudakov_eeqq_fsr_t), intent(in) :: sudakov + real(default), intent(in) :: pt2 + end function sudakov_eeqq_fsr_log_integrated_ubf <>= - pure function sudakov_eeqq_fsr_log_integrated_ubf (sudakov, pt2) result (log_sudakov) + pure module function sudakov_eeqq_fsr_log_integrated_ubf & + (sudakov, pt2) result (log_sudakov) real(default) :: log_sudakov class(sudakov_eeqq_fsr_t), intent(in) :: sudakov real(default), intent(in) :: pt2 real(default) :: k_max2, Lambda2 logical :: within_boundaries k_max2 = sudakov%kt2_max () within_boundaries = pt2 >= sudakov%powheg_settings%pt2_min & .and. pt2 <= k_max2 if (within_boundaries) then Lambda2 = sudakov%process_deps%lambda2_gen log_sudakov = pi / b0rad * ( log (4 * k_max2 / Lambda2) * & log (log (k_max2/Lambda2) / log (pt2/Lambda2)) - log(k_max2/pt2) ) else log_sudakov = 0 end if end function sudakov_eeqq_fsr_log_integrated_ubf @ %def sudakov_eeqq_fsr_log_integrated_ubf @ In order to arrive at the log integrated Sudakov implemented in [[sudakov_eeqq_fsr_log_integrated_ubf]], we overestimated \begin{equation} \log \left[ \frac{1 + \sqrt{1 - (k_T/k_\text{max})^2}} {1 - \sqrt{1 - (k_T/k_\text{max})^2}} \right] \leq \log \left[ \frac{4 k_\text{max}^2}{k_T^2} \right]. \end{equation} We correct for this in this veto step. <>= procedure :: reweight_ubf => sudakov_eeqq_fsr_reweight_ubf +<>= + module function sudakov_eeqq_fsr_reweight_ubf & + (sudakov, pt2) result (accepted) + logical :: accepted + class(sudakov_eeqq_fsr_t), intent(inout) :: sudakov + real(default), intent(in) :: pt2 + end function sudakov_eeqq_fsr_reweight_ubf <>= - function sudakov_eeqq_fsr_reweight_ubf (sudakov, pt2) result (accepted) + module function sudakov_eeqq_fsr_reweight_ubf (sudakov, pt2) result (accepted) logical :: accepted class(sudakov_eeqq_fsr_t), intent(inout) :: sudakov real(default), intent(in) :: pt2 real(default) :: log_bound, k_max2 k_max2 = sudakov%kt2_max () log_bound = log((1 + sqrt(1 - (pt2/k_max2))) & / (1 - sqrt(1 - (pt2/k_max2)))) call sudakov%rng%generate (sudakov%random) accepted = sudakov%random * log( 4 * k_max2 / pt2 ) <= log_bound end function sudakov_eeqq_fsr_reweight_ubf @ %def sudakov_eeqq_fsr_reweight_ubf @ In the eeqq case, we did not overestimate $\xi$ with $\xi_\text{max}$. <>= procedure :: reweight_xi_max => sudakov_eeqq_fsr_reweight_xi_max +<>= + module function sudakov_eeqq_fsr_reweight_xi_max & + (sudakov, xi) result (accepted) + logical :: accepted + class(sudakov_eeqq_fsr_t), intent(in) :: sudakov + real(default), intent(in) :: xi + end function sudakov_eeqq_fsr_reweight_xi_max <>= - function sudakov_eeqq_fsr_reweight_xi_max (sudakov, xi) result (accepted) - logical :: accepted - class(sudakov_eeqq_fsr_t), intent(in) :: sudakov - real(default), intent(in) :: xi - accepted = .true. + module function sudakov_eeqq_fsr_reweight_xi_max & + (sudakov, xi) result (accepted) + logical :: accepted + class(sudakov_eeqq_fsr_t), intent(in) :: sudakov + real(default), intent(in) :: xi + accepted = .true. end function sudakov_eeqq_fsr_reweight_xi_max @ %def sudakov_eeqq_fsr_reweight_xi_max @ In the eeqq case, we generate $y$ from a random number and afterwards compute $\xi(p_T,y)$ so the sequence is reversed in comparison to the simple case. <>= - procedure :: generate_xi_and_y_and_phi => sudakov_eeqq_fsr_generate_xi_and_y_and_phi + procedure :: generate_xi_and_y_and_phi => & + sudakov_eeqq_fsr_generate_xi_and_y_and_phi +<>= + module subroutine sudakov_eeqq_fsr_generate_xi_and_y_and_phi (sudakov, r) + class(sudakov_eeqq_fsr_t), intent(inout) :: sudakov + type(radiation_t), intent(inout) :: r + end subroutine sudakov_eeqq_fsr_generate_xi_and_y_and_phi <>= - subroutine sudakov_eeqq_fsr_generate_xi_and_y_and_phi (sudakov, r) + module subroutine sudakov_eeqq_fsr_generate_xi_and_y_and_phi (sudakov, r) class(sudakov_eeqq_fsr_t), intent(inout) :: sudakov type(radiation_t), intent(inout) :: r call sudakov%generate_y (r) call sudakov%generate_xi (r) call sudakov%generate_phi (r) end subroutine sudakov_eeqq_fsr_generate_xi_and_y_and_phi @ %def sudakov_eeqq_fsr_generate_xi_and_y_and_phi @ To generate $y$ for the [[eeqq]] UBF, we take a random number \begin{equation} -\log \left[ \frac{1 + \sqrt{1 - (p_T / k_\text{max})^2 }} {1 - \sqrt{1 - (p_T / k_\text{max})^2 }} \right] < r_y < \log \left[ \frac{1 + \sqrt{1 - (p_T / k_\text{max})^2 }} {1 - \sqrt{1 - (p_T / k_\text{max})^2 }} \right] \end{equation} and then compute \begin{equation} y(r_y) = \frac{e^{r_y} - 1}{e^{r_y} + 1} \end{equation} <>= procedure :: generate_y => sudakov_eeqq_fsr_generate_y +<>= + module subroutine sudakov_eeqq_fsr_generate_y (sudakov, r) + class(sudakov_eeqq_fsr_t), intent(inout) :: sudakov + type(radiation_t), intent(inout) :: r + end subroutine sudakov_eeqq_fsr_generate_y <>= - subroutine sudakov_eeqq_fsr_generate_y (sudakov, r) + module subroutine sudakov_eeqq_fsr_generate_y (sudakov, r) class(sudakov_eeqq_fsr_t), intent(inout) :: sudakov type(radiation_t), intent(inout) :: r real(default) :: k_max2, log_bound, r_y k_max2 = sudakov%kt2_max () call sudakov%rng%generate (sudakov%random) log_bound = log((1 + sqrt(1 - (r%pt2/k_max2))) & / (1 - sqrt(1 - (r%pt2/k_max2)))) r_y = -log_bound + sudakov%random * (two * log_bound) r%y = (exp(r_y) - 1) / (exp(r_y) + 1) end subroutine sudakov_eeqq_fsr_generate_y @ %def sudakov_eeqq_fsr_generate_y @ We generate $\xi$ for the [[eeqq]] UBF according to \begin{equation} \xi(p_T, y) = \frac{2 p_T}{\sqrt{s} \sqrt{1-y^2}} \end{equation} <>= procedure :: generate_xi => sudakov_eeqq_fsr_generate_xi +<>= + module subroutine sudakov_eeqq_fsr_generate_xi (sudakov, r) + class(sudakov_eeqq_fsr_t), intent(inout) :: sudakov + type(radiation_t), intent(inout) :: r + end subroutine sudakov_eeqq_fsr_generate_xi <>= - subroutine sudakov_eeqq_fsr_generate_xi (sudakov, r) + module subroutine sudakov_eeqq_fsr_generate_xi (sudakov, r) class(sudakov_eeqq_fsr_t), intent(inout) :: sudakov type(radiation_t), intent(inout) :: r real(default) :: q0 q0 = sqrt(sudakov%event_deps%s_hat) r%xi = (2 * sqrt(r%pt2)) / (q0 * sqrt(1 - r%y**2)) end subroutine sudakov_eeqq_fsr_generate_xi @ %def sudakov_eeqq_fsr_generate_xi @ \subsubsection{Massive FSR} <>= public :: sudakov_massive_fsr_t <>= type, extends (sudakov_t) :: sudakov_massive_fsr_t real(default) :: z, z1, z2 = 0._default real(default) :: xi_1, xi_min, xi_m = 0._default real(default) :: xi_max_extended = 1._default contains <> end type sudakov_massive_fsr_t @ %def sudakov_massive_fsr_t @ According to eq. A.42 [[1202.0465]], during the radiation generation, an alternative expression for $\xi_{\mathrm{max}}$, \begin{equation*} \xi_{\mathrm{max}} = 1 - \frac{(m+m_{\mathrm{rec}})^2}{q^2}, \end{equation*} is used, which corresponds to an extended Dalitz region. Phase space points outside of the original Dalitz region will be vetoed afterwards in [[sudakov_massive_fsr_reweight_xi_max]]. <>= procedure :: compute_xi_max_extended & => sudakov_massive_fsr_compute_xi_max_extended +<>= + module subroutine sudakov_massive_fsr_compute_xi_max_extended & + (sudakov, i_phs) + class(sudakov_massive_fsr_t), intent(inout) :: sudakov + integer, intent(in) :: i_phs + end subroutine sudakov_massive_fsr_compute_xi_max_extended <>= - subroutine sudakov_massive_fsr_compute_xi_max_extended (sudakov, i_phs) - class(sudakov_massive_fsr_t), intent(inout) :: sudakov - integer, intent(in) :: i_phs - real(default) :: m, mrec - real(default) :: q0 - type(vector4_t) :: p - q0 = sqrt(sudakov%event_deps%s_hat) - p = sudakov%event_deps%p_born_lab%get_momentum (1, sudakov%associated_emitter()) - m = p**1 - mrec = sqrt ((q0 - p%p(0))**2 - p%p(1)**2 - p%p(2)**2 - p%p(3)**2) - sudakov%xi_max_extended = one - (m + mrec)**2 / q0**2 + module subroutine sudakov_massive_fsr_compute_xi_max_extended (sudakov, i_phs) + class(sudakov_massive_fsr_t), intent(inout) :: sudakov + integer, intent(in) :: i_phs + real(default) :: m, mrec + real(default) :: q0 + type(vector4_t) :: p + q0 = sqrt(sudakov%event_deps%s_hat) + p = sudakov%event_deps%p_born_lab%get_momentum & + (1, sudakov%associated_emitter()) + m = p**1 + mrec = sqrt ((q0 - p%p(0))**2 - p%p(1)**2 - p%p(2)**2 - p%p(3)**2) + sudakov%xi_max_extended = one - (m + mrec)**2 / q0**2 end subroutine sudakov_massive_fsr_compute_xi_max_extended @ %def sudakov_massive_fsr_compute_xi_max_extended @ As described in [[1202.0465]], App. A.5., for massive emitters, the radiation variable $\xi$ is constructed as follows. First, \begin{equation} \xi_{\mathrm{min}}(K_T^2) = \frac{\sqrt{K_T^2 \left(K_T^2z_2^2 + 8\bar{p}^0q(1-z_2)\right)} - K_T^2z_2}{2q^2(1-z_2)} \end{equation} is computed where $K_T$ is a hardness scale for the radiation that determines the upper limit of the integral for the Sudakov form factor and coincides with the usual one, i.e. the transverse momentum of the radiated particle, in the massless limit. Then $\xi_1$ is computed according to the same equation with $z_2 \leftrightarrow z_1$. Finally, $\xi$ is generated according to \begin{equation} \xi = \frac{1}{q^2}\left(\exp\left[\log\left(\xi_{\rm{min}}q^2-K_T^2\right) + r\log\frac{\xi_m q^2 - K_T^2}{\xi_{\rm{min}}q^2-K_T^2}\right] + K_T^2\right), \end{equation} where $\xi_m = \rm{min}\left(\xi_{\rm{max}}, \xi_1\right)$. <>= procedure :: generate_xi => sudakov_massive_fsr_generate_xi +<>= + module subroutine sudakov_massive_fsr_generate_xi (sudakov, r) + class(sudakov_massive_fsr_t), intent(inout) :: sudakov + type(radiation_t), intent(inout) :: r + end subroutine sudakov_massive_fsr_generate_xi <>= - subroutine sudakov_massive_fsr_generate_xi (sudakov, r) + module subroutine sudakov_massive_fsr_generate_xi (sudakov, r) class(sudakov_massive_fsr_t), intent(inout) :: sudakov type(radiation_t), intent(inout) :: r real(default) :: pt2, q0, q02 real(default) :: E_em, xi_max_ext real(default) :: xi_1, xi_min, xi_m pt2 = r%pt2 E_em = sudakov%event_deps%p_born_lab%get_energy & (1, sudakov%associated_emitter()) q02 = sudakov%event_deps%s_hat; q0 = sqrt(q02) xi_max_ext = sudakov%xi_max_extended associate (z1 => sudakov%z1, z2 => sudakov%z2) xi_1 = (sqrt(pt2 * (pt2 * z1**2 + 8 * E_em * q0 * (one - z1))) - pt2 * z1) / & (two * q02 * (one - z1)) xi_min = (sqrt(pt2 * (pt2 * z2**2 + 8 * E_em * q0 * (one - z2))) - pt2 * z2) / & (two * q02 * (one - z2)) end associate xi_m = min (xi_max_ext, xi_1) call sudakov%rng%generate (sudakov%random) r%xi = (exp (log(xi_min * q02 - pt2) + sudakov%random * & log((xi_m * q02 - pt2) / (xi_min * q02 - pt2))) + pt2) / q02 end subroutine sudakov_massive_fsr_generate_xi @ %def sudakov_massive_fsr_generate_xi @ Generates the FKS variables in the case of a massive emitter and additionally computes $z$ and $z_{1/2}$ as well as $\xi_\text{max}$ to be stored in the Sudakov for later use when reweighting [[xi_max]] and retrieving the norm from the POWHEG grid. The corresponding derivations can be found in [[1202.0465]], App. A. <>= - procedure :: generate_xi_and_y_and_phi => sudakov_massive_fsr_generate_xi_and_y_and_phi + procedure :: generate_xi_and_y_and_phi => & + sudakov_massive_fsr_generate_xi_and_y_and_phi +<>= + module subroutine sudakov_massive_fsr_generate_xi_and_y_and_phi (sudakov, r) + class(sudakov_massive_fsr_t), intent(inout) :: sudakov + type(radiation_t), intent(inout) :: r + end subroutine sudakov_massive_fsr_generate_xi_and_y_and_phi <>= - subroutine sudakov_massive_fsr_generate_xi_and_y_and_phi (sudakov, r) + module subroutine sudakov_massive_fsr_generate_xi_and_y_and_phi (sudakov, r) class(sudakov_massive_fsr_t), intent(inout) :: sudakov type(radiation_t), intent(inout) :: r real(default) :: q0 real(default) :: m2, mrec2, k0_rec_max real(default) :: E_em type(vector4_t) :: p_emitter q0 = sqrt (sudakov%event_deps%s_hat) p_emitter = sudakov%event_deps%p_born_lab%get_momentum & (1, sudakov%associated_emitter()) associate (p => p_emitter%p) mrec2 = (q0 - p(0))**2 - p(1)**2 - p(2)**2 - p(3)**2 E_em = p(0) end associate m2 = p_emitter**2 call compute_dalitz_bounds (q0, m2, mrec2, sudakov%z1, sudakov%z2, k0_rec_max) call sudakov%generate_xi (r) sudakov%z = (two * r%pt2 * E_em - r%xi**2 * q0**3) / & (r%pt2 * r%xi * q0 - r%xi**2 * q0**3) sudakov%xi2_max = - (q0**2 * sudakov%z**2 - two * q0 * k0_rec_max * sudakov%z + mrec2) / & (q0**2 * sudakov%z * (one - sudakov%z)) sudakov%xi2_max = sudakov%xi2_max**2 r%y = two * (sudakov%z2 - sudakov%z) / (sudakov%z2 - sudakov%z1) - one call sudakov%generate_phi (r) end subroutine sudakov_massive_fsr_generate_xi_and_y_and_phi @ %def sudakov_massive_fsr_generate_xi_and_y_and_phi @ Computes the hardness scale as discussed in [[1202.0465]], App. A.2. \begin{equation} K_T^2 = \frac{\xi^2q^3 (1-z)}{2\bar{p}_{\rm{em}}^0 - z\xi q} \label{eq:HardnessMassiveFSR} \end{equation} <>= procedure :: kt2 => sudakov_massive_fsr_kt2 +<>= + module function sudakov_massive_fsr_kt2 (sudakov, xi, y) result (kt2) + real(default) :: kt2 + class(sudakov_massive_fsr_t), intent(in) :: sudakov + real(default), intent(in) :: xi, y + end function sudakov_massive_fsr_kt2 <>= - function sudakov_massive_fsr_kt2 (sudakov, xi, y) result (kt2) + module function sudakov_massive_fsr_kt2 (sudakov, xi, y) result (kt2) real(default) :: kt2 class(sudakov_massive_fsr_t), intent(in) :: sudakov real(default), intent(in) :: xi, y kt2 = sudakov%phs_fks_generator%real_kinematics%kt2 & (sudakov%i_phs, sudakov%associated_emitter(), UBF_FSR_MASSIVE, xi, y) end function sudakov_massive_fsr_kt2 @ %def sudakov_massive_fsr_kt2 @ For massive emitters, the upper bound on the radiated $p_T$ is \begin{equation*} t_{\mathrm{max}} = \frac{\xi_{\mathrm{max}}^2q^3(1-z_2)}{2 \cdot \bar{p}^0 - z_2\xi_{\mathrm{max}}q} \end{equation*} For this, we have to compute $z_2$ from the Dalitz bounds as [[sudakov%z2]] has not yet been set for the current event. <>= procedure :: kt2_max => sudakov_massive_fsr_kt2_max +<>= + pure module function sudakov_massive_fsr_kt2_max (sudakov) result (kt2_max) + real(default) :: kt2_max + class(sudakov_massive_fsr_t), intent(in) :: sudakov + end function sudakov_massive_fsr_kt2_max <>= - pure function sudakov_massive_fsr_kt2_max (sudakov) result (kt2_max) - real(default) :: kt2_max - class(sudakov_massive_fsr_t), intent(in) :: sudakov - real(default) :: q0, E_em, xi_max, z1, z2, m2, mrec2, k0_rec_max - type(vector4_t) :: p_emitter - q0 = sqrt(sudakov%event_deps%s_hat) - p_emitter = sudakov%event_deps%p_born_lab%get_momentum & - (1, sudakov%associated_emitter()) - associate (p => p_emitter%p) - mrec2 = (q0 - p(0))**2 - p(1)**2 - p(2)**2 - p(3)**2 - E_em = p(0) - end associate - m2 = p_emitter**2 - call compute_dalitz_bounds (q0, m2, mrec2, z1, z2, k0_rec_max) - xi_max = sudakov%xi_max_extended - kt2_max = (xi_max**2 * q0**3 * (one - z2)) / (two * E_em - z2 * xi_max * q0) + pure module function sudakov_massive_fsr_kt2_max (sudakov) result (kt2_max) + real(default) :: kt2_max + class(sudakov_massive_fsr_t), intent(in) :: sudakov + real(default) :: q0, E_em, xi_max, z1, z2, m2, mrec2, k0_rec_max + type(vector4_t) :: p_emitter + q0 = sqrt(sudakov%event_deps%s_hat) + p_emitter = sudakov%event_deps%p_born_lab%get_momentum & + (1, sudakov%associated_emitter()) + associate (p => p_emitter%p) + mrec2 = (q0 - p(0))**2 - p(1)**2 - p(2)**2 - p(3)**2 + E_em = p(0) + end associate + m2 = p_emitter**2 + call compute_dalitz_bounds (q0, m2, mrec2, z1, z2, k0_rec_max) + xi_max = sudakov%xi_max_extended + kt2_max = (xi_max**2 * q0**3 * (one - z2)) / (two * E_em - z2 * xi_max * q0) end function sudakov_massive_fsr_kt2_max @ %def sudakov_massive_fsr_kt2_max @ The upper bounding function for massive emitters according to Eq. A.56 in [1202.0465] (disregarding a possible factor of $\alpha_s$) is \begin{equation} U(\xi, y) \sim \frac{\sqrt{s}}{\bar{p}_{\rm{em}}} \frac{1}{\xi(1-z)} \label{eq:UBFMassiveFSR} \end{equation} <>= procedure :: upper_bound_func => sudakov_massive_fsr_upper_bound_func +<>= + pure module function sudakov_massive_fsr_upper_bound_func & + (sudakov, xi, y, alpha_s) result (u) + real(default) :: u + class(sudakov_massive_fsr_t), intent(in) :: sudakov + real(default), intent(in) :: xi, y, alpha_s + end function sudakov_massive_fsr_upper_bound_func <>= - pure function sudakov_massive_fsr_upper_bound_func (sudakov, xi, y, alpha_s) result (u) + pure module function sudakov_massive_fsr_upper_bound_func & + (sudakov, xi, y, alpha_s) result (u) real(default) :: u class(sudakov_massive_fsr_t), intent(in) :: sudakov real(default), intent(in) :: xi, y, alpha_s real(default) :: q, p_em q = sqrt (sudakov%event_deps%s_hat) p_em = space_part_norm & (sudakov%event_deps%p_born_lab%get_momentum (1, & sudakov%associated_emitter())) u = alpha_s * q / p_em * one / (xi * (one - sudakov%z)) end function sudakov_massive_fsr_upper_bound_func @ %def sudakov_massive_fsr_upper_bound_func @ The integrated upper-bounding function for massive final-state emitters is given by \begin{align*} I(t) &= \frac{q}{\bar{p}_{\rm{em}}}\left[\log\xi\log\left[(1-z_2)\frac{q}{k_T^2}\right] + \frac{1}{2} \log^2\xi + G(-k_T^2,q^2,\xi) - G(2\bar{p}_{\rm{em}},-q,\xi)\right] ^{\rm{min}(\xi_1(k_T^2),\xi_{\rm{max}})}_{\xi_{\rm{min}}} \\ &+ \frac{q}{\bar{p}_{\rm{em}}} \theta (\xi_{\rm{max}} - \xi_1(k_T^2)) \log\frac{\xi_{\rm{max}}}{\xi_1(k_T^2)}\log\frac{1-z_2}{1-z_1}, \end{align*} where the function $G(a,b,\xi)$ is given by \begin{equation} G(a,b,\xi) = \log(a+b\xi)\log\left(1-\frac{a+b\xi}{a}\right) + Li_2\left(\frac{a+b\xi}{a}\right), \label{GMinusMassiveFSR} \end{equation} for $a < 0$ and by \begin{equation} G(a,b,\xi) = \log\left|\frac{b\xi}{a}\right|\log a - Li_2\left(-\frac{b\xi}{a}\right) + \frac{\pi^2}{6}, \label{GPlusMassiveFSR} \end{equation} for $a > 0$. <>= procedure :: log_integrated_ubf => sudakov_massive_fsr_log_integrated_ubf +<>= + pure module function sudakov_massive_fsr_log_integrated_ubf & + (sudakov, pt2) result (log_sudakov) + real(default) :: log_sudakov + class(sudakov_massive_fsr_t), intent(in) :: sudakov + real(default), intent(in) :: pt2 + end function sudakov_massive_fsr_log_integrated_ubf <>= - pure function sudakov_massive_fsr_log_integrated_ubf (sudakov, pt2) result (log_sudakov) + pure module function sudakov_massive_fsr_log_integrated_ubf & + (sudakov, pt2) result (log_sudakov) real(default) :: log_sudakov class(sudakov_massive_fsr_t), intent(in) :: sudakov real(default), intent(in) :: pt2 real(default) :: xi, xi_max, xi_1, xi_min real(default) :: q0, p_em, E_em, m2, mrec2, k0_rec_max type(vector4_t) :: p_emitter real(default) :: log_sudakov_up, log_sudakov_down, z1, z2 logical :: within_boundaries within_boundaries = pt2 >= sudakov%powheg_settings%pt2_min & .and. pt2 <= sudakov%kt2_max () if (within_boundaries) then q0 = sqrt (sudakov%event_deps%s_hat) p_emitter = sudakov%event_deps%p_born_lab%get_momentum & (1, sudakov%associated_emitter()) associate (p => p_emitter%p) mrec2 = (q0 - p(0))**2 - p(1)**2 - p(2)**2 - p(3)**2 E_em = p(0) end associate m2 = p_emitter**2 call compute_dalitz_bounds (q0, m2, mrec2, z1, z2, k0_rec_max) p_em = space_part_norm (p_emitter) xi_max = sudakov%xi_max_extended xi_1 = (sqrt (pt2 * (pt2 * z1**2 + 8 * E_em * q0 * (one - z1))) - pt2 * z1) / & (two * q0**2 * (one - z1)) xi_min = (sqrt (pt2 * (pt2 * z2**2 + 8 * E_em * q0 * (one - z2))) - pt2 * z2) / & (two * q0**2 * (one - z2)) xi = min (xi_1, xi_max) log_sudakov_up = log(xi) * log((one - z2) * q0 / pt2) + log(xi)**2 / two + & G_FSR(-pt2, q0**2, xi) - G_FSR(two * E_em, -q0, xi) xi = xi_min log_sudakov_down = log(xi) * log((one - z2) * q0 / pt2) + log(xi)**2/two + & G_FSR(-pt2, q0**2, xi) - G_FSR(two * E_em, -q0, xi) log_sudakov = log_sudakov_up - log_sudakov_down if (xi_max > xi_1) & log_sudakov = log_sudakov + log(xi_max / xi_1) * log((one - z2) / (one - z1)) log_sudakov = twopi * q0 / p_em * log_sudakov else log_sudakov = 0 end if end function sudakov_massive_fsr_log_integrated_ubf @ %def sudakov_massive_fsr_log_integrated_ubf @ No further ubf veto needed for now. <>= procedure :: reweight_ubf => sudakov_massive_fsr_reweight_ubf +<>= + module function sudakov_massive_fsr_reweight_ubf & + (sudakov, pt2) result (accepted) + logical :: accepted + class(sudakov_massive_fsr_t), intent(inout) :: sudakov + real(default), intent(in) :: pt2 + end function sudakov_massive_fsr_reweight_ubf <>= - function sudakov_massive_fsr_reweight_ubf (sudakov, pt2) result (accepted) + module function sudakov_massive_fsr_reweight_ubf & + (sudakov, pt2) result (accepted) logical :: accepted class(sudakov_massive_fsr_t), intent(inout) :: sudakov real(default), intent(in) :: pt2 accepted = .true. end function sudakov_massive_fsr_reweight_ubf @ %def sudakov_massive_fsr_reweight_ubf @ In the massive case, we generated $0 < \xi < $[[xi_max_extended]] and thus we need to veto values of $\xi$ larger than $\xi_{\mathrm{max}}$ in this extra veto step. <>= procedure :: reweight_xi_max => sudakov_massive_fsr_reweight_xi_max +<>= + module function sudakov_massive_fsr_reweight_xi_max & + (sudakov, xi) result (accepted) + logical :: accepted + class(sudakov_massive_fsr_t), intent(in) :: sudakov + real(default), intent(in) :: xi + end function sudakov_massive_fsr_reweight_xi_max <>= - function sudakov_massive_fsr_reweight_xi_max (sudakov, xi) result (accepted) + module function sudakov_massive_fsr_reweight_xi_max & + (sudakov, xi) result (accepted) logical :: accepted class(sudakov_massive_fsr_t), intent(in) :: sudakov real(default), intent(in) :: xi accepted = xi < sqrt (sudakov%xi2_max) end function sudakov_massive_fsr_reweight_xi_max @ %def sudakov_massive_fsr_reweight_xi_max @ \subsubsection{Auxiliary functions} Implements the function $G(a,b,\xi)$ given in eq. (\ref{GPlusMassiveFSR}) and eq. (\ref{GMinusMassiveFSR}). Cannot be evaluated for $a = 0$. <>= elemental function G_FSR (a,b,xi) real(default) :: G_FSR real(default), intent(in) :: a, b, xi if (a > 0) then G_FSR = G_FSR_Plus (a,b,xi) else if (a < 0) then G_FSR = G_FSR_Minus (a,b,xi) end if end function G_FSR elemental function G_FSR_Minus (a,b,xi) real(default) :: G_FSR_Minus real(default), intent(in) :: a, b, xi G_FSR_Minus = log(a+b*xi)*log(one - (a+b*xi)/a) + Li2((a+b*xi)/a) end function G_FSR_Minus elemental function G_FSR_Plus (a,b,xi) real(default) :: G_FSR_Plus real(default), intent(in) :: a, b, xi G_FSR_Plus = log(abs(b*xi/a))*log(a) - Li2(-b*xi/a) + pi**2/6 end function G_FSR_Plus @ %def G_FSR, G_FSR_Minus, G_FSR_Plus @ \subsubsection{ISR} This is the Sudakov for ISR. This corresponds to [[1002.2581]], App. D. <>= public :: sudakov_isr_t <>= type, extends (sudakov_t) :: sudakov_isr_t real(default) :: xi_max_extended = 1._default contains <> end type sudakov_isr_t @ %def sudakov_isr_t @ This $k_T$ measure differs by a factor of $\frac{1}{1-\xi}$ w.r.t. [[sudakov_eeqq_fsr_kt2]] to account for $s_\text{born} = s_\text{real} * (1 - \xi)$. <>= procedure :: kt2 => sudakov_isr_kt2 +<>= + module function sudakov_isr_kt2 (sudakov, xi, y) result (kt2) + real(default) :: kt2 + class(sudakov_isr_t), intent(in) :: sudakov + real(default), intent(in) :: xi, y + end function sudakov_isr_kt2 <>= - function sudakov_isr_kt2 (sudakov, xi, y) result (kt2) + module function sudakov_isr_kt2 (sudakov, xi, y) result (kt2) real(default) :: kt2 class(sudakov_isr_t), intent(in) :: sudakov real(default), intent(in) :: xi, y kt2 = sudakov%phs_fks_generator%real_kinematics%kt2 & (sudakov%i_phs, sudakov%associated_emitter(), UBF_ISR, xi, y) end function sudakov_isr_kt2 @ %def sudakov_isr_kt2 @ The maximal transverse momentum is given by (D.4) in [1002.2581]. <>= procedure :: kt2_max => sudakov_isr_kt2_max +<>= + pure module function sudakov_isr_kt2_max (sudakov) result (kt2_max) + real(default) :: kt2_max + class(sudakov_isr_t), intent(in) :: sudakov + end function sudakov_isr_kt2_max <>= - pure function sudakov_isr_kt2_max (sudakov) result (kt2_max) - real(default) :: kt2_max - class(sudakov_isr_t), intent(in) :: sudakov - real(default) :: s_hat, xb_plus, xb_minus - s_hat = sudakov%event_deps%s_hat - xb_plus = sudakov%event_deps%x_born(I_PLUS) - xb_minus = sudakov%event_deps%x_born(I_MINUS) - kt2_max = s_hat * (1 - xb_plus**2) * (1 - xb_minus**2) & - / (xb_plus + xb_minus)**2 + pure module function sudakov_isr_kt2_max (sudakov) result (kt2_max) + real(default) :: kt2_max + class(sudakov_isr_t), intent(in) :: sudakov + real(default) :: s_hat, xb_plus, xb_minus + s_hat = sudakov%event_deps%s_hat + xb_plus = sudakov%event_deps%x_born(I_PLUS) + xb_minus = sudakov%event_deps%x_born(I_MINUS) + kt2_max = s_hat * (1 - xb_plus**2) * (1 - xb_minus**2) & + / (xb_plus + xb_minus)**2 end function sudakov_isr_kt2_max @ %def sudakov_isr_kt2_max -@ This covers the soft singularity for $\xi \to 0$ and the collinear singularities to both beam axes $y \to \pm 1$. +@ This covers the soft singularity for $\xi \to 0$ and the collinear +singularities to both beam axes $y \to \pm 1$. <>= procedure :: upper_bound_func => sudakov_isr_upper_bound_func +<>= + pure module function sudakov_isr_upper_bound_func & + (sudakov, xi, y, alpha_s) result (u) + real(default) :: u + class(sudakov_isr_t), intent(in) :: sudakov + real(default), intent(in) :: xi, y, alpha_s + end function sudakov_isr_upper_bound_func <>= - pure function sudakov_isr_upper_bound_func (sudakov, xi, y, alpha_s) result (u) + pure module function sudakov_isr_upper_bound_func & + (sudakov, xi, y, alpha_s) result (u) real(default) :: u class(sudakov_isr_t), intent(in) :: sudakov real(default), intent(in) :: xi, y, alpha_s u = alpha_s / (xi * (1 - y**2)) end function sudakov_isr_upper_bound_func @ %def sudakov_isr_upper_bound_func @ This is eq. (D.14) in [1002.2581]. It is similar to the eeqq case with the only difference beeing that the first log here includes both, $k^2_\text{T, max}$ and $s_b$ via $q^2$. <>= procedure :: log_integrated_ubf => sudakov_isr_log_integrated_ubf +<>= + pure module function sudakov_isr_log_integrated_ubf & + (sudakov, pt2) result (log_sudakov) + real(default) :: log_sudakov + class(sudakov_isr_t), intent(in) :: sudakov + real(default), intent(in) :: pt2 + end function sudakov_isr_log_integrated_ubf <>= - pure function sudakov_isr_log_integrated_ubf (sudakov, pt2) result (log_sudakov) + pure module function sudakov_isr_log_integrated_ubf & + (sudakov, pt2) result (log_sudakov) real(default) :: log_sudakov class(sudakov_isr_t), intent(in) :: sudakov real(default), intent(in) :: pt2 real(default) :: kT_max2, Lambda2, s_hat, q2 logical :: within_boundaries kT_max2 = sudakov%kt2_max () s_hat = sudakov%event_deps%s_hat q2 = kT_max2 + s_hat within_boundaries = pt2 >= sudakov%powheg_settings%pt2_min & .and. pt2 <= kT_max2 if (within_boundaries) then Lambda2 = sudakov%process_deps%lambda2_gen log_sudakov = pi / b0rad * ( log (q2 / Lambda2) * & log (log (kT_max2/Lambda2) / log (pt2/Lambda2)) - log(kT_max2/pt2) ) else log_sudakov = 0 end if end function sudakov_isr_log_integrated_ubf @ %def sudakov_isr_log_integrated_ubf @ In order to arrive at the log integrated Sudakov implemented in [[sudakov_isr_log_integrated_ubf]], we overestimated \begin{equation} \log \left[ \frac{\sqrt{x_+ - \rho} + \sqrt{x_- - \rho}} {\sqrt{x_+ - \rho} - \sqrt{x_- - \rho}} \right] \leq \log \left[ \frac{\sqrt{x_+} + \sqrt{x_-}} {\sqrt{x_+} - \sqrt{x_-}} \right] = \frac12 \log \left[ \frac{k_T^2 + s_b}{k_T^2} \right] \leq \frac12 \log \left[ \frac{q^2}{k_T^2} \right] \end{equation} with \begin{equation} x_\pm = \left( \sqrt{1 + \frac{k_T^2}{s_b}} \pm \frac{k_T}{\sqrt{s_b}} \right)^2 \quad , \quad \rho = \frac{s_b}{S} \quad \text{and} \quad q^2 = k^2_{T, \text{max}} + s_b. \end{equation} We correct for this in this veto step. <>= procedure :: reweight_ubf => sudakov_isr_reweight_ubf +<>= + module function sudakov_isr_reweight_ubf (sudakov, pt2) result (accepted) + logical :: accepted + class(sudakov_isr_t), intent(inout) :: sudakov + real(default), intent(in) :: pt2 + end function sudakov_isr_reweight_ubf <>= - function sudakov_isr_reweight_ubf (sudakov, pt2) result (accepted) + module function sudakov_isr_reweight_ubf (sudakov, pt2) result (accepted) logical :: accepted class(sudakov_isr_t), intent(inout) :: sudakov real(default), intent(in) :: pt2 real(default) :: log_bound, kT_max2, s_hat, s_had, rho, x_plus, x_minus, q2 s_hat = sudakov%event_deps%s_hat s_had = sudakov%event_deps%s_had rho = s_hat / s_had kT_max2 = sudakov%kt2_max () q2 = kT_max2 + s_hat x_plus = ( sqrt( 1 + pt2 / s_hat ) + sqrt( pt2 / s_hat ) )**2 x_minus = ( sqrt( 1 + pt2 / s_hat ) - sqrt( pt2 / s_hat ) )**2 !!! These are not the ISR x values! log_bound = log ((sqrt( x_plus - rho ) + sqrt( x_minus - rho )) & / (sqrt( x_plus - rho ) - sqrt( x_minus - rho ))) call sudakov%rng%generate (sudakov%random) accepted = sudakov%random * log( q2 / pt2 ) <= 2 * log_bound end function sudakov_isr_reweight_ubf @ %def sudakov_isr_reweight_ubf -@ In the ISR case, we overestimated $\xi_\text{max}$ with $\xi_\text{max}^{\text{ext}}$. +@ In the ISR case, we overestimated $\xi_\text{max}$ with +$\xi_\text{max}^{\text{ext}}$. We veto too large $\xi$ in this veto step. <>= procedure :: reweight_xi_max => sudakov_isr_reweight_xi_max +<>= + module function sudakov_isr_reweight_xi_max (sudakov, xi) result (accepted) + logical :: accepted + class(sudakov_isr_t), intent(in) :: sudakov + real(default), intent(in) :: xi + end function sudakov_isr_reweight_xi_max <>= - function sudakov_isr_reweight_xi_max (sudakov, xi) result (accepted) - logical :: accepted - class(sudakov_isr_t), intent(in) :: sudakov - real(default), intent(in) :: xi - accepted = xi < sqrt (sudakov%xi2_max) + module function sudakov_isr_reweight_xi_max (sudakov, xi) result (accepted) + logical :: accepted + class(sudakov_isr_t), intent(in) :: sudakov + real(default), intent(in) :: xi + accepted = xi < sqrt (sudakov%xi2_max) end function sudakov_isr_reweight_xi_max @ %def sudakov_isr_reweight_xi_max @ The actual $\xi_\text{max}$ for ISR is given by eq. (\ref{eqn:xi_max_isr}) and depends on $y$ which will be known only after we generated $\xi$. We thus use an overestimated \begin{equation*} \xi_{\mathrm{max}}^{\text{ext}} = 1 - \frac{s_b}{S_\text{had}} \end{equation*} to generate both and veto $\xi > \xi_\text{max}$ later on in [[sudakov_isr_reweight_xi_max]]. <>= procedure :: compute_xi_max_extended & => sudakov_isr_compute_xi_max_extended +<>= + module subroutine sudakov_isr_compute_xi_max_extended (sudakov) + class(sudakov_isr_t), intent(inout) :: sudakov + real(default) :: s_hat, s_had + end subroutine sudakov_isr_compute_xi_max_extended <>= - subroutine sudakov_isr_compute_xi_max_extended (sudakov) - class(sudakov_isr_t), intent(inout) :: sudakov - real(default) :: s_hat, s_had - s_hat = sudakov%event_deps%s_hat - s_had = sudakov%event_deps%s_had + module subroutine sudakov_isr_compute_xi_max_extended (sudakov) + class(sudakov_isr_t), intent(inout) :: sudakov + real(default) :: s_hat, s_had + s_hat = sudakov%event_deps%s_hat + s_had = sudakov%event_deps%s_had sudakov%xi_max_extended = one - s_hat / s_had - end subroutine sudakov_isr_compute_xi_max_extended + end subroutine sudakov_isr_compute_xi_max_extended @ %def sudakov_isr_compute_xi_max_extended -@ In the ISR case, we generate $\xi$ from a random number and afterwards compute -$y(p_T,\xi)$ as in the simple case. +@ In the ISR case, we generate $\xi$ from a random number and +afterwards compute $y(p_T,\xi)$ as in the simple case. <>= - procedure :: generate_xi_and_y_and_phi => sudakov_isr_generate_xi_and_y_and_phi + procedure :: generate_xi_and_y_and_phi => & + sudakov_isr_generate_xi_and_y_and_phi +<>= + module subroutine sudakov_isr_generate_xi_and_y_and_phi (sudakov, r) + class(sudakov_isr_t), intent(inout) :: sudakov + type(radiation_t), intent(inout) :: r + end subroutine sudakov_isr_generate_xi_and_y_and_phi <>= - subroutine sudakov_isr_generate_xi_and_y_and_phi (sudakov, r) + module subroutine sudakov_isr_generate_xi_and_y_and_phi (sudakov, r) class(sudakov_isr_t), intent(inout) :: sudakov type(radiation_t), intent(inout) :: r call sudakov%generate_xi (r) call sudakov%generate_y (r) call sudakov%generate_phi (r) call sudakov%compute_xi2_max (r) end subroutine sudakov_isr_generate_xi_and_y_and_phi @ %def sudakov_isr_generate_xi_and_y_and_phi @ Inverting [[real_kinematics_kt2]], we find \begin{equation} |y(k_T^2, \xi)| = \sqrt{1 - \frac{4 (1 - \xi)}{\xi^2}\frac{k_T^2}{s_b}}. \end{equation} We then determine the sign randomly as there is no preferred direction. <>= procedure :: generate_y => sudakov_isr_generate_y +<>= + module subroutine sudakov_isr_generate_y (sudakov, r) + class(sudakov_isr_t), intent(inout) :: sudakov + type(radiation_t), intent(inout) :: r + end subroutine sudakov_isr_generate_y <>= - subroutine sudakov_isr_generate_y (sudakov, r) + module subroutine sudakov_isr_generate_y (sudakov, r) class(sudakov_isr_t), intent(inout) :: sudakov type(radiation_t), intent(inout) :: r real(default) :: s_hat s_hat = sudakov%event_deps%s_hat r%y = sqrt (1 - (4 * (1 - r%xi) * r%pt2 / (s_hat * r%xi**2))) call sudakov%rng%generate (sudakov%random) if (sudakov%random > 0.5_default) then r%y = - r%y end if end subroutine sudakov_isr_generate_y @ %def sudakov_isr_generate_y @ We want to generate $x := 1 - \xi$ for the ISR UBF with probability density \begin{equation} f(x) = \frac{1}{\sqrt{(x_+ - x)(x_- - x)}}. \end{equation} within the limits \begin{equation} \frac{s_b}{S} = x_{b,\oplus} x_{b,\ominus} = \rho \leq x \leq x_- = x_- = \left( \sqrt{1 + \frac{k_T^2}{s_b}} - \frac{k_T}{\sqrt{s_b}} \right)^2. \end{equation} To generate $x$ according to the distribution \begin{equation} F(x) = \int f(x) \, dx = = 2 \sinh^{-1} \left( \sqrt{ \frac{x_- - x}{x_+ - x_-} } \right) \end{equation} we use \begin{equation} x = x_- + (x_- - x_+) \sinh^2(\frac{r_x}{2}) \end{equation} with a random number $r_x \in \text{Unif}[F(x_-),F(\rho)] = \text{Unif}\left[0,2 \sinh^{-1} \left( \sqrt{ \frac{x_- - \rho}{x_+ - x_-} } \right) \right]$.\\ <>= procedure :: generate_xi => sudakov_isr_generate_xi +<>= + module subroutine sudakov_isr_generate_xi (sudakov, r) + class(sudakov_isr_t), intent(inout) :: sudakov + type(radiation_t), intent(inout) :: r + end subroutine sudakov_isr_generate_xi <>= - subroutine sudakov_isr_generate_xi (sudakov, r) + module subroutine sudakov_isr_generate_xi (sudakov, r) class(sudakov_isr_t), intent(inout) :: sudakov type(radiation_t), intent(inout) :: r real(default) :: s_hat, s_had, r_x, x, x_minus, x_plus, rho s_hat = sudakov%event_deps%s_hat s_had = sudakov%event_deps%s_had rho = 1 - sudakov%xi_max_extended x_plus = ( sqrt( 1 + r%pt2 / s_hat ) + sqrt( r%pt2 / s_hat ) )**2 x_minus = ( sqrt( 1 + r%pt2 / s_hat ) - sqrt( r%pt2 / s_hat ) )**2 call sudakov%rng%generate (sudakov%random) r_x = sudakov%random r_x = r_x * 2 * asinh (sqrt ( (x_minus - rho) / (x_plus - x_minus) ) ) x = x_minus + (x_minus - x_plus) * sinh (r_x / 2)**2 r%xi = 1 - x end subroutine sudakov_isr_generate_xi @ %def sudakov_isr_generate_xi @ Computes the actual [[xi2_max]] for ISR. Can only be computed after y has been generated. Needs to be called before [[sudakov_isr_reweight_xi_max]]. <>= procedure :: compute_xi2_max => sudakov_isr_compute_xi2_max +<>= + module subroutine sudakov_isr_compute_xi2_max (sudakov, r) + class(sudakov_isr_t), intent(inout) :: sudakov + type(radiation_t), intent(inout) :: r + end subroutine sudakov_isr_compute_xi2_max <>= - subroutine sudakov_isr_compute_xi2_max (sudakov, r) + module subroutine sudakov_isr_compute_xi2_max (sudakov, r) class(sudakov_isr_t), intent(inout) :: sudakov type(radiation_t), intent(inout) :: r real(default) :: xi_max xi_max = get_xi_max_isr (sudakov%event_deps%x_born, r%y) sudakov%xi2_max = xi_max**2 end subroutine sudakov_isr_compute_xi2_max @ %def sudakov_isr_compute_xi2_max @ \subsection{Main POWHEG class} <>= public :: powheg_matching_t <>= type, extends(matching_t) :: powheg_matching_t type(grid_t) :: grid type(phs_fks_generator_t) :: phs_fks_generator type(powheg_settings_t) :: settings type(event_deps_t) :: event_deps type(process_deps_t) :: process_deps type(sudakov_wrapper_t), dimension(:), allocatable :: sudakov integer :: n_emissions = 0 logical :: active = .true. contains <> end type powheg_matching_t @ %def powheg_matching_t @ <>= procedure :: get_method => powheg_matching_get_method +<>= + module function powheg_matching_get_method (matching) result (method) + type(string_t) :: method + class(powheg_matching_t), intent(in) :: matching + end function powheg_matching_get_method <>= - function powheg_matching_get_method (matching) result (method) - type(string_t) :: method - class(powheg_matching_t), intent(in) :: matching - method = matching_method (MATCH_POWHEG) + module function powheg_matching_get_method (matching) result (method) + type(string_t) :: method + class(powheg_matching_t), intent(in) :: matching + method = matching_method (MATCH_POWHEG) end function powheg_matching_get_method @ %def powheg_matching_get_method @ <>= procedure :: before_shower => powheg_matching_before_shower +<>= + module subroutine powheg_matching_before_shower & + (matching, particle_set, vetoed) + class(powheg_matching_t), intent(inout) :: matching + type(particle_set_t), intent(inout) :: particle_set + logical, intent(out) :: vetoed + end subroutine powheg_matching_before_shower <>= - subroutine powheg_matching_before_shower & + module subroutine powheg_matching_before_shower & (matching, particle_set, vetoed) class(powheg_matching_t), intent(inout) :: matching type(particle_set_t), intent(inout) :: particle_set logical, intent(out) :: vetoed if (debug_on) call msg_debug2 (D_MATCHING, "powheg_matching_before_shower") if (signal_is_pending ()) return if (.not. matching%active) return call matching%update (particle_set) if (matching%settings%test_sudakov) then call matching%test_sudakov () stop end if if (.not. matching%settings%disable_sudakov) & call matching%generate_emission (particle_set = particle_set) vetoed = .false. end subroutine powheg_matching_before_shower @ %def powheg_matching_before_shower @ <>= procedure :: first_event => powheg_matching_first_event +<>= + module subroutine powheg_matching_first_event (matching) + class(powheg_matching_t), intent(inout), target :: matching + end subroutine powheg_matching_first_event <>= - subroutine powheg_matching_first_event (matching) + module subroutine powheg_matching_first_event (matching) class(powheg_matching_t), intent(inout), target :: matching associate (instance => matching%process_instance) matching%process_deps%lab_is_cm = instance%lab_is_cm (1) end associate call matching%setup_grids () end subroutine powheg_matching_first_event @ %def powheg_matching_first_event @ <>= procedure :: after_shower => powheg_matching_after_shower +<>= + module subroutine powheg_matching_after_shower & + (matching, particle_set, vetoed) + class(powheg_matching_t), intent(inout) :: matching + type(particle_set_t), intent(inout) :: particle_set + logical, intent(out) :: vetoed + end subroutine powheg_matching_after_shower <>= - subroutine powheg_matching_after_shower (matching, particle_set, vetoed) + module subroutine powheg_matching_after_shower & + (matching, particle_set, vetoed) class(powheg_matching_t), intent(inout) :: matching type(particle_set_t), intent(inout) :: particle_set logical, intent(out) :: vetoed vetoed = .false. end subroutine powheg_matching_after_shower @ %def powheg_matching_after_shower @ \subsubsection{Output} <>= procedure :: write => powheg_write +<>= + module subroutine powheg_write (matching, unit) + class(powheg_matching_t), intent(in) :: matching + integer, intent(in), optional :: unit + end subroutine powheg_write <>= - subroutine powheg_write (matching, unit) + module subroutine powheg_write (matching, unit) class(powheg_matching_t), intent(in) :: matching integer, intent(in), optional :: unit integer :: u, alr u = given_output_unit (unit); if (u < 0) return call write_separator (u, 2) write (u, "(1X,A)") "POWHEG Emission Generator" write (u, "(1X,A)") "Process name: " // char (matching%process_name) if (allocated (matching%rng)) then call matching%rng%write (u) else write (u, "(1X,A)") "RNG not allocated" end if call matching%qcd%write (u) call matching%settings%write (u) call matching%event_deps%write (u) call matching%process_deps%write (u) do alr = 1, size (matching%sudakov) call write_separator (u) write (u, "(1X,A,I12,A)") "sudakov (alr = ", alr, ")" call matching%sudakov(alr)%s%write (u) end do call write_separator (u, 2) end subroutine powheg_write @ %def powheg_write @ Finalization of the POWHEG matching. It creates a log file showing how many events failed at which stage of the -POWHEG veto algorithm. -It also issues a warning if too many POWHEG excess events occured. If this is the -case, the chosen upper-bounding function is unsuitable for the considered process. +POWHEG veto algorithm. It also issues a warning if too many POWHEG +excess events occured. If this is the case, the chosen upper-bounding +function is unsuitable for the considered process. <>= procedure :: final => powheg_matching_final +<>= + module subroutine powheg_matching_final (matching) + class(powheg_matching_t), intent(in) :: matching + end subroutine powheg_matching_final <>= - subroutine powheg_matching_final (matching) + module subroutine powheg_matching_final (matching) class(powheg_matching_t), intent(in) :: matching integer :: u, alr, n_fail_ubf_total, n_sqme_total real :: frac_ubf_fail type(string_t) :: filename n_fail_ubf_total = 0 n_sqme_total = 0 u = free_unit () filename = matching%process_name // "_veto.log" open (file=char(filename), unit=u, action='write') write (u, '(A)') "Summary of POWHEG veto procedure" do alr = 1, matching%process_deps%n_alr write(u,'(A,I0)') 'alr: ', alr call matching%sudakov(alr)%s%veto_counter%write (u) call write_separator (u) n_fail_ubf_total = n_fail_ubf_total & + matching%sudakov(alr)%s%veto_counter%n_fail_ubf n_sqme_total = n_sqme_total & + matching%sudakov(alr)%s%veto_counter%n_sqme end do write (u,'(A,I0)') "Total number of events which radiate a gluon: ", & matching%n_emissions close (u) if (n_sqme_total > 0) then frac_ubf_fail = one * n_fail_ubf_total / n_sqme_total * 100 if (frac_ubf_fail > 1) then - write (msg_buffer, "(A16,I4,A2,F6.2,A24)") "There have been ", n_fail_ubf_total,& + write (msg_buffer, "(A16,I4,A2,F6.2,A24)") "There have been ", & + n_fail_ubf_total,& " (", frac_ubf_fail, "%) POWHEG grid excesses." call msg_warning end if end if end subroutine powheg_matching_final @ %def powheg_matching_final @ \subsubsection{Initialization and Finalization} -Setup the POWHEG grids. As they are filled during integration using the POWHEG hook, -we can load them from file here. +Setup the POWHEG grids. As they are filled during integration using +the POWHEG hook, we can load them from file here. <>= procedure :: setup_grids => powheg_matching_setup_grids +<>= + module subroutine powheg_matching_setup_grids (matching) + class(powheg_matching_t), intent(inout), target :: matching + end subroutine powheg_matching_setup_grids <>= - subroutine powheg_matching_setup_grids (matching) + module subroutine powheg_matching_setup_grids (matching) class(powheg_matching_t), intent(inout), target :: matching call matching%prepare_for_events () if (matching%check_grids ()) then call matching%load_grids () else call msg_fatal ("POWHEG: POWHEG grids are invalid.") end if call matching%grid%compute_and_write_mean_and_max () call matching%import_norms_from_grid () end subroutine powheg_matching_setup_grids @ %def powheg_matching_setup_grids @ Sets up the Sudakovs and chooses the correct UBF type. To determine the UBF type, we check for massive or initial state emitters for each ALR. For consistency, we thus also check the eeqq case for each ALR. +Gfortran 7/8/9 bug, has to remain in the main module: <>= procedure :: setup_sudakovs => powheg_matching_setup_sudakovs -<>= +<>= subroutine powheg_matching_setup_sudakovs (powheg) class(powheg_matching_t), intent(inout), target :: powheg integer :: alr, ubf_type, n_in, emitter logical :: is_fsr, is_massive, is_eeqq allocate (powheg%sudakov (powheg%process_deps%n_alr)) do alr = 1, powheg%process_deps%n_alr select type (pcm => powheg%process_instance%pcm) type is (pcm_nlo_t) associate(reg_data => pcm%region_data, & phs => powheg%phs_fks_generator) n_in = reg_data%get_n_in() emitter = reg_data%get_emitter (alr) is_fsr = emitter > n_in if (is_fsr) then is_massive = phs%is_massive (emitter) else if (emitter == 0) then is_massive = phs%is_massive (1) .or. phs%is_massive (2) else is_massive = phs%is_massive (emitter) end if - if (is_massive) call msg_bug ("setup_sudakovs: ISR off massive emitters not implemented.") + if (is_massive) call msg_bug ("setup_sudakovs: ISR " // & + "off massive emitters not implemented.") end if is_eeqq = n_in == 2 .and. & reg_data%n_legs_born == 4 .and. & .not. phs%is_massive(3) .and. & .not. phs%is_massive(4) ! (PS 2021-05-28) This includes FSR regions of pp -> jj. end associate end select if (is_fsr) then if (is_eeqq) then ubf_type = UBF_FSR_MASSLESS_RECOIL else if (is_massive) then ubf_type = UBF_FSR_MASSIVE else ubf_type = powheg%settings%upper_bound_func_type end if select case (ubf_type) case (UBF_FSR_SIMPLE) allocate (sudakov_simple_fsr_t :: powheg%sudakov(alr)%s) case (UBF_FSR_MASSLESS_RECOIL) allocate (sudakov_eeqq_fsr_t :: powheg%sudakov(alr)%s) case (UBF_FSR_MASSIVE) allocate (sudakov_massive_fsr_t :: powheg%sudakov(alr)%s) case default - call msg_fatal ("powheg_setup_sudakovs: Please choose upper bounding function!") + call msg_fatal ("powheg_setup_sudakovs: Please choose " // & + "upper bounding function!") end select else allocate (sudakov_isr_t :: powheg%sudakov(alr)%s) end if if (allocated (powheg%rng)) then !!! generator mode call powheg%sudakov(alr)%s%init (powheg%process_deps, & powheg%event_deps, powheg%settings, & powheg%qcd, powheg%phs_fks_generator, powheg%rng) else !!! lookup mode call powheg%sudakov(alr)%s%init (powheg%process_deps, & powheg%event_deps, powheg%settings, & powheg%qcd, powheg%phs_fks_generator) end if call powheg%sudakov(alr)%s%set_i_phs (alr) end do end subroutine powheg_matching_setup_sudakovs @ %def powheg_matching_setup_sudakovs @ <>= procedure :: init => powheg_matching_init +<>= + module subroutine powheg_matching_init (matching, var_list, process_name) + class(powheg_matching_t), intent(out) :: matching + type(var_list_t), intent(in) :: var_list + type(string_t), intent(in) :: process_name + end subroutine powheg_matching_init <>= - subroutine powheg_matching_init (matching, var_list, process_name) + module subroutine powheg_matching_init (matching, var_list, process_name) class(powheg_matching_t), intent(out) :: matching <> end subroutine powheg_matching_init @ %def powheg_matching_init @ Updates the Born momenta and the Born matrix element stored in the event dependencies. <>= generic :: update => update_momenta, & update_particle_set procedure :: update_momenta => powheg_matching_update_momenta procedure :: update_particle_set => powheg_matching_update_particle_set +<>= + module subroutine powheg_matching_update_momenta (powheg, p_born) + class(powheg_matching_t), intent(inout) :: powheg + type(vector4_t), dimension(:), intent(in) :: p_born + end subroutine powheg_matching_update_momenta + module subroutine powheg_matching_update_particle_set (powheg, particle_set) + class(powheg_matching_t), intent(inout) :: powheg + type(particle_set_t), intent(in) :: particle_set + end subroutine powheg_matching_update_particle_set <>= - subroutine powheg_matching_update_momenta (powheg, p_born) + module subroutine powheg_matching_update_momenta (powheg, p_born) class(powheg_matching_t), intent(inout) :: powheg type(vector4_t), dimension(:), intent(in) :: p_born type(lorentz_transformation_t) :: lt_lab_to_cms real(default) :: sqme_born real(default), dimension(2) :: x_born - sqme_born = powheg%process_instance%get_sqme (powheg%process_deps%i_term_born) + sqme_born = powheg%process_instance%get_sqme & + (powheg%process_deps%i_term_born) x_born = powheg%phs_fks_generator%isr_kinematics%x if (.not. powheg%process_deps%lab_is_cm) then lt_lab_to_cms = powheg%process_instance%get_boost_to_cms (1) call powheg%update_event_deps (sqme_born, p_born, x_born, lt_lab_to_cms) else call powheg%update_event_deps (sqme_born, p_born, x_born) end if end subroutine powheg_matching_update_momenta - subroutine powheg_matching_update_particle_set (powheg, particle_set) + module subroutine powheg_matching_update_particle_set (powheg, particle_set) class(powheg_matching_t), intent(inout) :: powheg type(particle_set_t), intent(in) :: particle_set integer, dimension(:), allocatable :: indices logical, dimension(:), allocatable :: in_out_mask integer :: i allocate (in_out_mask (particle_set%get_n_tot())) do i = 1, particle_set%get_n_tot() in_out_mask(i) = particle_set%prt(i)%get_status () == PRT_INCOMING & .or. particle_set%prt(i)%get_status () == PRT_OUTGOING end do allocate (indices (size (particle_set%get_indices (in_out_mask)))) indices = particle_set%get_indices (in_out_mask) call powheg%update_momenta (particle_set%get_momenta (indices)) end subroutine powheg_matching_update_particle_set @ %def powheg_matching_update @ <>= procedure :: update_event_deps => powheg_matching_update_event_deps +<>= + module subroutine powheg_matching_update_event_deps & + (powheg, sqme_born, p_born, x_born, lt_lab_to_cms) + class(powheg_matching_t), intent(inout) :: powheg + real(default), intent(in) :: sqme_born + type(vector4_t), dimension(:), intent(in) :: p_born + real(default), dimension(2), intent(in) :: x_born + type(lorentz_transformation_t), intent(in), optional :: lt_lab_to_cms + end subroutine powheg_matching_update_event_deps <>= - subroutine powheg_matching_update_event_deps (powheg, sqme_born, p_born, x_born, lt_lab_to_cms) + module subroutine powheg_matching_update_event_deps & + (powheg, sqme_born, p_born, x_born, lt_lab_to_cms) class(powheg_matching_t), intent(inout) :: powheg real(default), intent(in) :: sqme_born type(vector4_t), dimension(:), intent(in) :: p_born real(default), dimension(2), intent(in) :: x_born type(lorentz_transformation_t), intent(in), optional :: lt_lab_to_cms call powheg%event_deps%update (sqme_born, p_born, x_born, lt_lab_to_cms) end subroutine powheg_matching_update_event_deps @ %def powheg_matching_update_event_deps @ <>= - procedure :: boost_preal_to_lab_frame => powheg_matching_boost_preal_to_lab_frame + procedure :: boost_preal_to_lab_frame => & + powheg_matching_boost_preal_to_lab_frame +<>= + module subroutine powheg_matching_boost_preal_to_lab_frame (powheg, i_phs) + class(powheg_matching_t), intent(inout) :: powheg + integer, intent(in) :: i_phs + end subroutine powheg_matching_boost_preal_to_lab_frame <>= - subroutine powheg_matching_boost_preal_to_lab_frame (powheg, i_phs) + module subroutine powheg_matching_boost_preal_to_lab_frame (powheg, i_phs) class(powheg_matching_t), intent(inout) :: powheg type(lorentz_transformation_t) :: lt_cms_to_lab integer, intent(in) :: i_phs associate (event_deps => powheg%event_deps) if (powheg%process_deps%lab_is_cm) then - event_deps%p_real_lab%phs_point(i_phs) = event_deps%p_real_cms%phs_point(i_phs) + event_deps%p_real_lab%phs_point(i_phs) = & + event_deps%p_real_cms%phs_point(i_phs) else lt_cms_to_lab = powheg%process_instance%get_boost_to_lab (1) event_deps%p_real_lab%phs_point(i_phs) = & lt_cms_to_lab * event_deps%p_real_cms%phs_point(i_phs) end if end associate end subroutine powheg_matching_boost_preal_to_lab_frame @ %def powheg_matching_boost_preal_to_lab_frame @ <>= procedure :: boost_preal_to_cms => powheg_matching_boost_preal_to_cms +<>= + module subroutine powheg_matching_boost_preal_to_cms (powheg, i_phs) + class(powheg_matching_t), intent(inout) :: powheg + integer, intent(in) :: i_phs + end subroutine powheg_matching_boost_preal_to_cms <>= - subroutine powheg_matching_boost_preal_to_cms (powheg, i_phs) + module subroutine powheg_matching_boost_preal_to_cms (powheg, i_phs) class(powheg_matching_t), intent(inout) :: powheg type(lorentz_transformation_t) :: lt type(vector4_t), dimension(:), allocatable :: p_real integer, intent(in) :: i_phs real(default) :: sqrts_real type(vector4_t) :: p0, p1 associate (event_deps => powheg%event_deps) if (powheg%process_deps%lab_is_cm) then event_deps%p_real_cms%phs_point(i_phs) = event_deps%p_real_lab%phs_point(i_phs) else p_real = event_deps%p_real_lab%phs_point(i_phs) p0 = p_real(1) + p_real(2) sqrts_real = (p0)**1 lt = boost (p0, sqrts_real) p1 = inverse(lt) * p_real(1) lt = lt * rotation_to_2nd (3, space_part (p1)) p_real = inverse (lt) * p_real event_deps%p_real_cms%phs_point(i_phs) = p_real end if end associate end subroutine powheg_matching_boost_preal_to_cms @ %def powheg_matching_boost_preal_to_cms @ This is the last step of the Sudakov veto algorithm removing the dependence on the normalization factor $N(\xi_i,y_i)$ and the upper bounding function $U(\xi,y,\alpha_s)$ from the final result. In cases where the Born kinematic failed the cuts, $\mathcal{B}$ was artificially set to zero, so we cannot use it in a denominator. We assume that the generator cuts applied for parton showering are weak cuts to avoid divergent phase space regions, thus $\mathcal{B}$ would have been large and the reweighting would have failed. Here, we can compute [[sqme_real]] directly with scale $p_T$ but we still need to correct [[sqme_born]]. <>= - procedure :: reweight_matrix_elements => powheg_matching_reweight_matrix_elements + procedure :: reweight_matrix_elements => & + powheg_matching_reweight_matrix_elements +<>= + module function powheg_matching_reweight_matrix_elements & + (powheg, r) result (accepted) + logical :: accepted + class(powheg_matching_t), intent(inout) :: powheg + type(radiation_t), intent(in) :: r + end function powheg_matching_reweight_matrix_elements <>= - function powheg_matching_reweight_matrix_elements (powheg, r) result (accepted) + module function powheg_matching_reweight_matrix_elements & + (powheg, r) result (accepted) logical :: accepted class(powheg_matching_t), intent(inout) :: powheg type(radiation_t), intent(in) :: r integer :: emitter, n_in, i_phs, i_term_born, alphas_power, em real(default) :: sqme_real_x_jacobian, sqme_born real(default) :: norm, ubf, ubound, random, weight real(default) :: alpha_s_kt, alpha_s_muren, muren, mufac, pt2 integer, dimension(2) :: flv_born real(double), dimension(-6:6) :: pdf_dbl real(double) :: x_dbl, q_dbl type(vector4_t), dimension(:), allocatable :: p_real_cms, p_real_lab real(default), dimension(2) :: pdf_born_mufac, pdf_born_kt if (debug_on) call msg_debug (D_MATCHING, "reweight_matrix_elements") sqme_born = powheg%event_deps%sqme_born if (sqme_born == 0) then accepted = .false. return end if call powheg%rng%generate (random) i_phs = powheg%process_deps%alr_to_i_phs (r%alr) select type (pcm => powheg%process_instance%pcm) type is (pcm_nlo_t) emitter = pcm%region_data%get_emitter (r%alr) n_in = pcm%region_data%get_n_in() if (emitter <= n_in) then allocate(p_real_lab(pcm%region_data%get_n_legs_real())) call powheg%phs_fks_generator%generate_isr_from_xi_and_y (& r%xi, sqrt (powheg%sudakov(r%alr)%s%xi2_max), r%y, & r%phi, i_phs, powheg%event_deps%p_born_lab%get_momenta(1), & p_real_lab) powheg%event_deps%p_real_lab%phs_point(i_phs) = p_real_lab call powheg%boost_preal_to_cms (i_phs) else allocate(p_real_cms(pcm%region_data%get_n_legs_real())) call powheg%phs_fks_generator%generate_fsr_from_xi_and_y (r%xi, r%y, & r%phi, emitter, i_phs, & powheg%event_deps%p_born_cms%get_momenta(1), & p_real_cms) powheg%event_deps%p_real_cms%phs_point(i_phs) = p_real_cms call powheg%boost_preal_to_lab_frame (i_phs) end if if (debug_active (D_MATCHING)) then if (emitter <= n_in) then pt2 = (transverse_part (p_real_lab (size(p_real_lab))))**2 else pt2 = (transverse_part (p_real_cms (size(p_real_cms))))**2 end if call assert_equal (OUTPUT_UNIT, r%pt2, pt2, & "reweight_matrix_elements: generated p_real does not fit to sudakovs pt2") end if call powheg%copy_momenta (i_phs) norm = powheg%norm_from_xi_and_y (r) associate (s => powheg%sudakov(r%alr)%s) alpha_s_kt = s%alpha_s (r%pt2, use_correct=.true.) i_term_born = powheg%process_deps%i_term_born muren = powheg%process_instance%term(i_term_born)%get_ren_scale () mufac = powheg%process_instance%term(i_term_born)%get_fac_scale () alpha_s_muren = s%alpha_s (muren**2, use_correct=.true.) ubf = s%upper_bound_func (r%xi, r%y, alpha_s_kt) sqme_real_x_jacobian = powheg%compute_sqme_real (r%alr, sqrt(r%pt2)) !!! Correct all factors of alphas(muren) to alphas(kt). alphas_power = powheg%process_deps%alphas_power sqme_born = (alpha_s_kt / alpha_s_muren)**alphas_power * sqme_born !!! Also correct the PDFs previously computed at mufac instead of kt select type (pcm => powheg%process_instance%term(i_term_born)%pcm) type is (pcm_nlo_t) if (pcm%has_pdfs) then associate (reg_data => pcm%region_data) flv_born = reg_data%flv_born(i_term_born)%flst(1:2) where (flv_born == 21) flv_born = 0 end associate associate (pdf_data => powheg%process_deps%pdf_data, & x_born => powheg%event_deps%x_born) do em = 1, 2 x_dbl = x_born(em); q_dbl = mufac call pdf_data%evolve(x_dbl, q_dbl, pdf_dbl) pdf_born_mufac(em) = pdf_dbl(flv_born(em)) / x_born(em) x_dbl = x_born(em); q_dbl = sqrt(r%pt2) call pdf_data%evolve(x_dbl, q_dbl, pdf_dbl) pdf_born_kt(em) = pdf_dbl(flv_born(em)) / x_born(em) end do sqme_born = pdf_born_kt(1) * pdf_born_kt(2) / & (pdf_born_mufac(1) * pdf_born_mufac(2)) * sqme_born end associate end if end select ubound = sqme_born * ubf * norm weight = sqme_real_x_jacobian / ubound if (weight > 1) call s%veto_counter%record_ubf_fail() if (debug_active (D_MATCHING)) then if (weight < 0) call msg_warning ("R/B < 0!") end if accepted = random < weight end associate if (debug_active (D_MATCHING)) then print *, ' r%alr = ', r%alr print *, ' r%xi = ', r%xi print *, ' r%y = ', r%y print *, ' emitter = ', emitter print *, ' random = ', random print *, ' sqme_real_x_jacobian = ', sqme_real_x_jacobian print *, ' sqme_born = ', sqme_born print *, ' ubf = ', ubf print *, ' norm = ', norm print *, ' ubound = ', ubound print *, ' matrix element accepted = ', accepted end if end select end function powheg_matching_reweight_matrix_elements @ %def powheg_matching_reweight_matrix_element @ \subsubsection{Generation algorithm and grid initialization} [[compute_sqme_real]] is the projected real matrix element $R_{\alpha_r} = S_{\alpha_r} R$ whereby the current $\alpha_r$ is implied by the [[emitter]]. Furthermore, it is multiplied by the Jacobian. <>= procedure :: compute_sqme_real => powheg_matching_compute_sqme_real +<>= + module function powheg_matching_compute_sqme_real & + (powheg, alr, scale) result (sqme) + real(default) :: sqme + class(powheg_matching_t), intent(inout) :: powheg + integer, intent(in) :: alr + real(default), intent(in), optional :: scale + end function powheg_matching_compute_sqme_real <>= - function powheg_matching_compute_sqme_real (powheg, alr, scale) result (sqme) + module function powheg_matching_compute_sqme_real & + (powheg, alr, scale) result (sqme) real(default) :: sqme class(powheg_matching_t), intent(inout) :: powheg integer, intent(in) :: alr real(default), intent(in), optional :: scale real(default), allocatable :: q integer :: i_phs, i_term logical :: is_subtraction is_subtraction = .false. select type (pcm_work => powheg%process_instance%pcm_work) class is (pcm_nlo_workspace_t) i_phs = powheg%process_deps%alr_to_i_phs (alr) i_term = powheg%process_deps%i_term_real(i_phs) associate (instance => powheg%process_instance) select type (pcm_work => instance%term(i_term)%pcm_work) type is (pcm_nlo_workspace_t) call pcm_work%set_powheg_mode () end select if (present(scale)) then if (.not. allocated (q)) then allocate (q, source = scale) else q = scale end if - call instance%compute_sqme_rad (i_term, i_phs, is_subtraction, scale_forced=q) + call instance%compute_sqme_rad & + (i_term, i_phs, is_subtraction, scale_forced=q) else call instance%compute_sqme_rad (i_term, i_phs, is_subtraction) end if sqme = instance%get_sqme (i_term) end associate end select end function powheg_matching_compute_sqme_real @ %def powheg_matching_compute_sqme_real @ <>= procedure :: set_scale => powheg_matching_set_scale +<>= + module subroutine powheg_matching_set_scale (powheg, pT2) + class(powheg_matching_t), intent(inout) :: powheg + real(default), intent(in) :: pT2 + end subroutine powheg_matching_set_scale <>= - subroutine powheg_matching_set_scale (powheg, pT2) + module subroutine powheg_matching_set_scale (powheg, pT2) class(powheg_matching_t), intent(inout) :: powheg real(default), intent(in) :: pT2 call powheg%process_instance%set_fac_scale (sqrt(pT2)) end subroutine powheg_matching_set_scale @ %def powheg_matching_set_scale @ <>= procedure :: update_sudakovs => powheg_matching_update_sudakovs +<>= + module subroutine powheg_matching_update_sudakovs (powheg, alr, i_phs, y) + class(powheg_matching_t), intent(inout) :: powheg + integer, intent(in) :: alr, i_phs + real(default), intent(in) :: y + end subroutine powheg_matching_update_sudakovs <>= - subroutine powheg_matching_update_sudakovs (powheg, alr, i_phs, y) + module subroutine powheg_matching_update_sudakovs (powheg, alr, i_phs, y) class(powheg_matching_t), intent(inout) :: powheg integer, intent(in) :: alr, i_phs real(default), intent(in) :: y real(default) :: q0, m2, mrec2, k0_rec_max type(vector4_t) :: p_emitter select type (s => powheg%sudakov(alr)%s) type is (sudakov_massive_fsr_t) q0 = sqrt (s%event_deps%s_hat) p_emitter = s%event_deps%p_born_lab%get_momentum (1, & s%associated_emitter ()) associate (p => p_emitter%p) mrec2 = (q0 - p(0))**2 - p(1)**2 - p(2)**2 - p(3)**2 end associate m2 = p_emitter**2 call compute_dalitz_bounds (q0, m2, mrec2, s%z1, s%z2, k0_rec_max) s%z = s%z2 - (s%z2 - s%z1) * (one + y) / two end select end subroutine powheg_matching_update_sudakovs @ %def powheg_matching_update_sudakovs @ <>= procedure :: import_norms_from_grid => powheg_matching_import_norms_from_grid +<>= + module subroutine powheg_matching_import_norms_from_grid (powheg) + class(powheg_matching_t), intent(inout) :: powheg + end subroutine powheg_matching_import_norms_from_grid <>= - subroutine powheg_matching_import_norms_from_grid (powheg) + module subroutine powheg_matching_import_norms_from_grid (powheg) class(powheg_matching_t), intent(inout) :: powheg integer :: alr real(default) :: norm_max do alr = 1, powheg%process_deps%n_alr norm_max = powheg%grid%get_maximum_in_3d (alr) call powheg%sudakov(alr)%s%set_normalization (norm_max) end do end subroutine powheg_matching_import_norms_from_grid @ %def powheg_matching_import_norms_from_grid @ We save the POWHEG grid to a file to be used for the event generation. If it has no non-zero entries, we assume that the integration was skipped because there were existing VAMP(2) and POWHEG grid files. <>= procedure :: save_grids => powheg_matching_save_grids +<>= + module subroutine powheg_matching_save_grids (powheg) + class(powheg_matching_t), intent(inout) :: powheg + end subroutine powheg_matching_save_grids <>= - subroutine powheg_matching_save_grids (powheg) + module subroutine powheg_matching_save_grids (powheg) class(powheg_matching_t), intent(inout) :: powheg type(string_t) :: filename, n_points filename = powheg%process_name // ".pg" if (powheg%grid%has_non_zero_entries()) then call powheg%grid%save_to_file (char (filename)) else call msg_warning("POWHEG grid not saved to file. An existing " // & & char(filename) // " will be used for the events.") end if end subroutine powheg_matching_save_grids @ %def powheg_matching_save_grids @ <>= procedure :: load_grids => powheg_matching_load_grids +<>= + module subroutine powheg_matching_load_grids (powheg) + class(powheg_matching_t), intent(inout) :: powheg + end subroutine powheg_matching_load_grids <>= - subroutine powheg_matching_load_grids (powheg) + module subroutine powheg_matching_load_grids (powheg) class(powheg_matching_t), intent(inout) :: powheg type(string_t) :: filename, n_points filename = powheg%process_name // ".pg" call powheg%grid%load_from_file (char (filename)) if (powheg%grid%has_non_zero_entries()) then write (msg_buffer, "(A,A,A)") "POWHEG: using grids from file '", & char (filename), "'" else call msg_fatal("POWHEG grid in " // char(filename) // & & " is zero and cannot be used for event generation.") end if call msg_message () end subroutine powheg_matching_load_grids @ %def powheg_matching_load_grids @ <>= procedure :: check_grids => powheg_matching_check_grids +<>= + module function powheg_matching_check_grids (powheg) result (ok) + logical :: ok + class(powheg_matching_t), intent(in) :: powheg + end function powheg_matching_check_grids <>= - function powheg_matching_check_grids (powheg) result (ok) + module function powheg_matching_check_grids (powheg) result (ok) logical :: ok class(powheg_matching_t), intent(in) :: powheg type(string_t) :: filename, n_points filename = powheg%process_name // ".pg" ok = os_file_exist (filename) .and. & verify_points_for_grid (char (filename), & [powheg%settings%size_grid_xi, & powheg%settings%size_grid_y, & powheg%process_deps%n_alr]) end function powheg_matching_check_grids @ %def powheg_matching_check_grids @ This routine implements the Sudakov veto algorithm as explained in [[[1002.2581]]], Sec. 7.1 and Bijans thesis, Sec. 3.3 and B.2. Here, we veto in order to correct for \begin{enumerate} \item an overestimation of parts of the upper bounding function for complicated UBFs in [[reweight_ubf]] \item the usage of $\alpha_s^\text{rad}$ instead of $\alpha_s^\text{true}$ in [[reweight_alpha_s]] \item the usage of [[xi_max_extended]] instead of [[xi_max]] in [[reweight_xi_max]] which is only non-trivial for massive FSR and ISR. \item using the maximum of the entire POWHEG normalization grid $N_\text{max}$ when generating $p_T$ instead of the correct bin value $N(\xi_i,y_i)$ in [[reweight_norm]]. \item the normalization grid and the upper bounding function in the first place in [[reweight_matrix_elements]] to retrieve the correct fraction $\frac{\mathcal{R}(\xi,y) \mathcal{J}(\xi,y)}{\mathcal{B}}$. \end{enumerate} By keeping the radiation with the largest [[pt2]], we are effectively implementing the highest bid procedure. The loop over ALRs implemented here should instead loop over all radiation regions [[i_phs]] for each underlying Born flavor structure $f_b$. We should then pick the ALR with probability proportional to the corresponding $\mathcal{R}_{\alpha_r}$. The version here works only for $ee \to t \overline{t}$ as there is just a single $f_b$ and each ALR corresponds to its own real phase space. We need to compute [[xi2_max]] before [[kt2_max]] here. In the simple case, it is given by $\xi_\text{max} = \frac{s-M_\text{rec}^2}{s}$ but we can also take the $\xi_\text{max}$ from the phase space computed before; they coincide. This is possible already at this point because in the simple case, $\xi_\text{max}$ only depends on the Born momentum of the emitter and not on the real kinematics. In the eeqq case, $M_\text{rec} = 0$, so $\xi_\text{max} = 1$ and in the massive case as well as in the ISR case, we compute an extended $\xi_\text{max}^\text{ext}$ because the actual $\xi_\text{max}$ depends on the real kinematics. <>= procedure :: generate_emission => powheg_matching_generate_emission +<>= + module subroutine powheg_matching_generate_emission & + (powheg, particle_set, pt2_generated) + class(powheg_matching_t), intent(inout) :: powheg + type(particle_set_t), intent(inout), optional :: particle_set + real(default), intent(out), optional :: pt2_generated + end subroutine powheg_matching_generate_emission <>= - subroutine powheg_matching_generate_emission (powheg, particle_set, pt2_generated) + module subroutine powheg_matching_generate_emission & + (powheg, particle_set, pt2_generated) class(powheg_matching_t), intent(inout) :: powheg type(particle_set_t), intent(inout), optional :: particle_set real(default), intent(out), optional :: pt2_generated type(radiation_t) :: r, r_max real(default) :: xi2_max integer :: alr logical :: accepted type(vector4_t), dimension(:), allocatable :: p_real_max if (signal_is_pending ()) return r_max%pt2 = zero r_max%alr = 0 if (debug_on) call msg_debug (D_MATCHING, "powheg_matching_generate_emission") select type (pcm => powheg%process_instance%pcm) type is (pcm_nlo_t) allocate (p_real_max (pcm%region_data%n_legs_real)) end select select type (pcm_work => powheg%process_instance%pcm_work) class is (pcm_nlo_workspace_t) do alr = 1, powheg%process_deps%n_alr if (signal_is_pending ()) return associate (sudakov => powheg%sudakov(alr)%s) select type (sudakov) type is (sudakov_simple_fsr_t) xi2_max = pcm_work%get_xi_max (alr)**2 call sudakov%update_xi2_max (xi2_max) if (debug2_active (D_MATCHING)) then call check_xi2_max (sudakov) end if type is (sudakov_eeqq_fsr_t) call sudakov%update_xi2_max (one) type is (sudakov_massive_fsr_t) call sudakov%compute_xi_max_extended (sudakov%i_phs) type is (sudakov_isr_t) call sudakov%compute_xi_max_extended () class default call msg_fatal ("powheg_matching_generate_emission: unknown Sudakov!") end select r%alr = alr r%pt2 = sudakov%kt2_max () sudakov%sum_log_rands = 0 if (debug_on) call msg_debug (D_MATCHING, "Starting evolution at r%pt2", r%pt2) PT_EVOLUTION: do if (signal_is_pending ()) return call sudakov%generate_emission (r) if (signal_is_pending ()) return if (r%valid) then accepted = powheg%reweight_norm (r) call sudakov%veto_counter%record_norm (.not. accepted) if (.not. accepted) cycle PT_EVOLUTION accepted = powheg%reweight_matrix_elements (r) call sudakov%veto_counter%record_sqme (.not. accepted) if (.not. accepted) cycle PT_EVOLUTION end if exit end do PT_EVOLUTION if (r%pt2 > r_max%pt2 .and. r%valid) then r_max = r p_real_max = powheg%event_deps%p_real_lab%get_momenta (sudakov%i_phs) end if end associate end do if (r_max%pt2 > powheg%settings%pt2_min) then powheg%n_emissions = powheg%n_emissions + 1 call powheg%set_scale (r_max%pt2) if (present (particle_set)) then call powheg%build_particle_set (particle_set, & powheg%event_deps%p_born_lab%get_momenta (1), & p_real_max, r_max%alr, r_max%y) end if if (present (pt2_generated)) pt2_generated = r_max%pt2 else call powheg%set_scale (powheg%settings%pt2_min) if (present (pt2_generated)) pt2_generated = powheg%settings%pt2_min end if end select contains subroutine check_xi2_max (sudakov) class(sudakov_t), intent(in) :: sudakov real(default) :: s_hat, mrec2, xi2_max type(vector4_t) :: p_emitter s_hat = sudakov%event_deps%s_hat p_emitter = sudakov%event_deps%p_born_lab%get_momentum & (1, sudakov%associated_emitter()) associate (p => p_emitter%p) mrec2 = (sqrt(s_hat) - p(0))**2 - p(1)**2 - p(2)**2 - p(3)**2 end associate xi2_max = (s_hat - mrec2) / s_hat xi2_max = xi2_max**2 call assert_equal (OUTPUT_UNIT, sudakov%xi2_max, xi2_max, & "powheg_matching_generate_emission: xi_max inconsistent") end subroutine check_xi2_max end subroutine powheg_matching_generate_emission @ %def generate_emission @ <>= procedure :: build_particle_set => powheg_matching_build_particle_set +<>= + module subroutine powheg_matching_build_particle_set & + (powheg, particle_set, p_born, p_real, alr, y) + class(powheg_matching_t), intent(inout) :: powheg + type(particle_set_t), intent(inout) :: particle_set + type(vector4_t), dimension(:), intent(in) :: p_born, p_real + integer, intent(in) :: alr + real(default), intent(in) :: y + end subroutine powheg_matching_build_particle_set <>= - subroutine powheg_matching_build_particle_set & + module subroutine powheg_matching_build_particle_set & (powheg, particle_set, p_born, p_real, alr, y) class(powheg_matching_t), intent(inout) :: powheg type(particle_set_t), intent(inout) :: particle_set type(vector4_t), dimension(:), intent(in) :: p_born, p_real integer, intent(in) :: alr real(default), intent(in) :: y integer, dimension(:), allocatable :: i_flvs_real, flv_radiated real(default) :: r_col integer :: emitter, i_flv_real select type (pcm => powheg%process_instance%pcm) type is (pcm_nlo_t) i_flvs_real = pcm%region_data%get_flavor_indices (.false.) i_flv_real = i_flvs_real (alr) allocate (flv_radiated (size (pcm%region_data%get_flv_states_real (i_flv_real)))) flv_radiated = pcm%region_data%get_flv_states_real (i_flv_real) emitter = pcm%region_data%get_emitter(alr) if (emitter == 0) then if (y > 0) then emitter = 1 else emitter = 2 end if end if end select call powheg%rng%generate (r_col) call particle_set%build_radiation (p_real, emitter, flv_radiated, & powheg%process_instance%process%get_model_ptr (), r_col) end subroutine powheg_matching_build_particle_set @ %def powheg_matching_build_particle_set @ When generating the transverse momentum in [[sudakov_generate_pt2]], we used $N^\mathrm{max}$ as normalization. In this veto step, we correct for this by keeping events with probability \begin{equation*} \frac{N(\xi_i,y_i)}{N^\mathrm{max}}. \end{equation*} <>= procedure :: reweight_norm => powheg_matching_reweight_norm +<>= + module function powheg_matching_reweight_norm (powheg, r) result (accepted) + logical :: accepted + class(powheg_matching_t), intent(inout) :: powheg + type(radiation_t), intent(in) :: r + end function powheg_matching_reweight_norm <>= - function powheg_matching_reweight_norm (powheg, r) result (accepted) + module function powheg_matching_reweight_norm (powheg, r) result (accepted) logical :: accepted class(powheg_matching_t), intent(inout) :: powheg type(radiation_t), intent(in) :: r real(default) :: random, norm_max, norm_true if (debug_on) call msg_debug2 (D_MATCHING, "reweight_norm") call powheg%rng%generate (random) norm_true = powheg%norm_from_xi_and_y (r) norm_max = powheg%sudakov(r%alr)%s%norm_max accepted = random < norm_true / norm_max if (debug2_active (D_MATCHING)) then print *, ' r%alr = ', r%alr print *, ' random = ', random print *, ' norm_true = ', norm_true print *, ' norm_max = ', norm_max print *, ' norm accepted = ', accepted end if if (debug_active (D_MATCHING)) then if (.not. (zero < r%xi .and. & r%xi < sqrt(powheg%sudakov(r%alr)%s%xi2_max))) then call msg_bug ("powheg_matching_reweight_norm: xi is out of bounds") end if if (norm_true > norm_max) then call msg_bug ("powheg_matching_reweight_norm: norm shouldnt be larger than norm_max") end if end if end function powheg_matching_reweight_norm @ %def powheg_matching_reweight_norm @ Retrieves the value of the norm for given $\xi$ and $y$ from the pre-computed POWHEG grid. To find the correct bin in the grid, we need to invert the generation of $\xi$ and $y$ from their corresponding random numbers. As the grid is filled during integration, the formulas here are the inverse of those in [[compute_y_from_emitter]] instead of those specific for POWHEG matching. There, we coompute $y$ with [[y = (one - two * r_y)]] and [[y = 1.5_default * (y - y**3 / 3)]]. The inverse of this operation is given by \begin{align} r_y &= \begin{cases} \frac14 \left[ (i\sqrt{3}-1)A^{\frac13} - (i\sqrt{3}+1)A^{-\frac13} + 2 \right] \quad \text{for} \quad y < 0 \\ \frac14 \left[-(i\sqrt{3}+1)A^{\frac13} + (i\sqrt{3}-1)A^{-\frac13} + 2 \right] \quad \text{for} \quad y > 0 \end{cases}\nonumber\\ &\text{with} \quad A = \sqrt{y^2-1} + y \end{align} <>= procedure :: norm_from_xi_and_y => powheg_matching_norm_from_xi_and_y +<>= + module function powheg_matching_norm_from_xi_and_y & + (powheg, r) result (norm_true) + real(default) :: norm_true + class(powheg_matching_t), intent(inout) :: powheg + type(radiation_t), intent(in) :: r + end function powheg_matching_norm_from_xi_and_y <>= - function powheg_matching_norm_from_xi_and_y (powheg, r) result (norm_true) + module function powheg_matching_norm_from_xi_and_y & + (powheg, r) result (norm_true) real(default) :: norm_true class(powheg_matching_t), intent(inout) :: powheg type(radiation_t), intent(in) :: r real(default) :: f_alr real(default), dimension(2) :: rands real(default) :: beta f_alr = (one * r%alr) / powheg%process_deps%n_alr - tiny_07 rands(I_XI) = r%xi / sqrt (powheg%sudakov(r%alr)%s%xi2_max) select type (s => powheg%sudakov(r%alr)%s) type is (sudakov_simple_fsr_t) rands(I_Y) = map_y_to_rands_massless(r%y) type is (sudakov_massive_fsr_t) beta = beta_emitter (sqrt (powheg%event_deps%s_hat), & powheg%event_deps%p_born_lab%get_momentum (1, s%associated_emitter())) rands(I_Y) = - log((one - r%y * beta) / (one + beta)) / & log((one + beta) / (one - beta)) type is (sudakov_eeqq_fsr_t) rands(I_Y) = map_y_to_rands_massless(r%y) type is (sudakov_isr_t) rands(I_Y) = map_y_to_rands_massless(r%y) class default call msg_fatal ("powheg_matching_norm_from_xi_and_y: unknown Sudakov!") end select norm_true = powheg%grid%get_value ([rands, f_alr]) contains function map_y_to_rands_massless (y) result (r_y) real(default), intent(in) :: y real(default) :: r_y complex(default) :: y_cmplx, A, r_y_cmplx y_cmplx = cmplx(y, 0, default) A = sqrt(y_cmplx**2 - 1) + y_cmplx if (y < 0) then r_y_cmplx = 0.25_default * (cmplx(-one,sqrt(three)) * A**(one/three) & - cmplx(one,sqrt(three)) * A**(-one/three) + two) else if (y == 0) then r_y_cmplx = cmplx(0.5_default,0) else if (y > 0) then r_y_cmplx = 0.25_default * (-cmplx(one,sqrt(three)) * A**(one/three) & + cmplx(-one,sqrt(three)) * A**(-one/three) + two) end if r_y = real(r_y_cmplx) end function map_y_to_rands_massless end function powheg_matching_norm_from_xi_and_y @ %def powheg_matching_norm_from_xi_and_y @ \subsection{$\alpha_s$ and its reweighting} The main point to ensure here is that the simple fixed-flavor-1-loop expression $\alpha_s^\text{rad}$ is larger than the more accurate $\alpha_s$ such that we can use a reweighting veto and use $\alpha_s^\text{rad}$ for the generation of the emission. This can be done by setting \begin{equation} \alpha_s^\text{rad}(\mu_0) = \alpha_s (\mu_0) \end{equation} whereby $\mu_0^2$ is the [[scale_to_relate2]] that is taken to be $p_{T,\text{min}}^2$. <>= procedure :: prepare_for_events => powheg_matching_prepare_for_events +<>= + module subroutine powheg_matching_prepare_for_events (matching) + class(powheg_matching_t), intent(inout), target :: matching + end subroutine powheg_matching_prepare_for_events <>= - subroutine powheg_matching_prepare_for_events (matching) + module subroutine powheg_matching_prepare_for_events (matching) class(powheg_matching_t), intent(inout), target :: matching - if (debug_on) call msg_debug (D_MATCHING, "powheg_matching_prepare_for_events") + if (debug_on) call msg_debug & + (D_MATCHING, "powheg_matching_prepare_for_events") call matching%setup_nlo_environment () call matching%grid%init ([matching%settings%size_grid_xi, & matching%settings%size_grid_y, matching%process_deps%n_alr]) call matching%compute_lambda2_gen () call matching%setup_sudakovs () end subroutine powheg_matching_prepare_for_events @ %def powheg_matching_prepare_for_events @ Computes the scale $\Lambda$ used for the (log integrated) UBFs. By construction, it is always $p_{T,\text{min}} > \Lambda$. <>= procedure :: compute_lambda2_gen => powheg_matching_compute_lambda2_gen +<>= + module subroutine powheg_matching_compute_lambda2_gen (matching) + class(powheg_matching_t), intent(inout) :: matching + end subroutine powheg_matching_compute_lambda2_gen <>= - subroutine powheg_matching_compute_lambda2_gen (matching) + module subroutine powheg_matching_compute_lambda2_gen (matching) class(powheg_matching_t), intent(inout) :: matching real(default) :: scale_to_relate2, alpha_s scale_to_relate2 = matching%settings%pt2_min alpha_s = get_alpha_s (matching%qcd, scale_to_relate2) matching%process_deps%lambda2_gen = exp (- one / (b0rad * alpha_s)) * & scale_to_relate2 end subroutine powheg_matching_compute_lambda2_gen @ %def powheg_matching_compute_lambda2_gen @ This is the setup of the process dependencies stored for the matching. <>= procedure :: setup_nlo_environment => powheg_matching_setup_nlo_environment +<>= + module subroutine powheg_matching_setup_nlo_environment (matching) + class(powheg_matching_t), intent(inout) :: matching + end subroutine powheg_matching_setup_nlo_environment <>= - subroutine powheg_matching_setup_nlo_environment (matching) + module subroutine powheg_matching_setup_nlo_environment (matching) class(powheg_matching_t), intent(inout) :: matching integer :: n_tot_born, n_tot_real integer :: i_phs, i_term_real, i_term integer :: n_phs, nlo_type - if (debug_on) call msg_debug (D_MATCHING, "powheg_matching_setup_nlo_environment") + if (debug_on) call msg_debug & + (D_MATCHING, "powheg_matching_setup_nlo_environment") select type (pcm_work => matching%process_instance%pcm_work) class is (pcm_nlo_workspace_t) matching%process_deps%sqrts = matching%process_instance%get_sqrts () select type (pcm => matching%process_instance%pcm) type is (pcm_nlo_t) matching%process_deps%n_alr = pcm%region_data%n_regions n_tot_born = pcm%region_data%n_legs_born n_tot_real = pcm%region_data%n_legs_real call pcm%setup_phs_generator (pcm_work, & matching%phs_fks_generator, matching%process_deps%sqrts, & singular_jacobian = matching%settings%singular_jacobian) end select associate (instance => matching%process_instance) i_term_real = instance%process%get_first_real_component () associate (process_deps => matching%process_deps) select type (phs => instance%kin(i_term_real)%phs) type is (phs_fks_t) n_phs = size (phs%phs_identifiers) allocate (process_deps%phs_identifiers (n_phs)) process_deps%phs_identifiers = phs%phs_identifiers end select call instance%process%get_coupling_powers(process_deps%alpha_power, & process_deps%alphas_power) allocate (matching%process_deps%alr_to_i_phs & (size (pcm_work%real_kinematics%alr_to_i_phs))) process_deps%alr_to_i_phs = pcm_work%real_kinematics%alr_to_i_phs allocate (process_deps%i_term_real (n_phs)) i_phs = 1 do i_term = 1, size (instance%term) nlo_type = instance%term(i_term)%nlo_type if (nlo_type == BORN) then process_deps%i_term_born = i_term else if (nlo_type == NLO_REAL) then if (instance%kin(i_term)%emitter >= 0) then process_deps%i_term_real(i_phs) = i_term i_phs = i_phs + 1 end if end if end do end associate end associate call matching%event_deps%p_born_lab%init (n_tot_born, 1) call matching%event_deps%p_born_cms%init (n_tot_born, 1) call matching%event_deps%p_real_lab%init (n_tot_real, n_phs) call matching%event_deps%p_real_cms%init (n_tot_real, n_phs) end select end subroutine powheg_matching_setup_nlo_environment @ %def powheg_matching_setup_nlo_environment @ Copy momenta from [[event_deps]] to [[real_kinematics]]. <>= procedure :: copy_momenta => powheg_matching_copy_momenta +<>= + module subroutine powheg_matching_copy_momenta (matching, i_phs) + class(powheg_matching_t), intent(inout) :: matching + integer, intent(in) :: i_phs + end subroutine powheg_matching_copy_momenta <>= - subroutine powheg_matching_copy_momenta (matching, i_phs) - class(powheg_matching_t), intent(inout) :: matching - integer, intent(in) :: i_phs - select type (pcm_work => matching%process_instance%pcm_work) - class is (pcm_nlo_workspace_t) - call pcm_work%real_kinematics%p_real_cms%set_momenta & - (i_phs, matching%event_deps%p_real_cms%get_momenta (i_phs)) - call pcm_work%real_kinematics%p_real_lab%set_momenta & - (i_phs, matching%event_deps%p_real_lab%get_momenta (i_phs)) - end select + module subroutine powheg_matching_copy_momenta (matching, i_phs) + class(powheg_matching_t), intent(inout) :: matching + integer, intent(in) :: i_phs + select type (pcm_work => matching%process_instance%pcm_work) + class is (pcm_nlo_workspace_t) + call pcm_work%real_kinematics%p_real_cms%set_momenta & + (i_phs, matching%event_deps%p_real_cms%get_momenta (i_phs)) + call pcm_work%real_kinematics%p_real_lab%set_momenta & + (i_phs, matching%event_deps%p_real_lab%get_momenta (i_phs)) + end select end subroutine powheg_matching_copy_momenta @ %def powheg_matching_copy_momenta @ [[qcd%alpha%get]] should implement a variable-flavor result and optionally return [[n_flavors]] that are active at the scale... <>= function get_alpha_s (qcd, scale2) result (alpha_s) real(default) :: alpha_s class(qcd_t), intent(in) :: qcd real(default), intent(in) :: scale2 integer :: nf, order ! TODO: (bcn 2015-01-30) implement variable flavor alpha_s alpha_s = qcd%alpha%get (sqrt(scale2)) select type (alpha => qcd%alpha) type is (alpha_qcd_from_scale_t) nf = alpha%nf order = alpha%order type is (alpha_qcd_from_lambda_t) nf = alpha%nf order = alpha%order type is (alpha_qcd_lhapdf_t) call msg_warning ("get_alpha_s for LHADPF not supported!" // & " Assuming 5-flavors and LO (1-loop) running!") !!! TODO (PS-2021-09-16) We'd need to get nf from LHAPDF to continue here. nf = 5 order = 0 class default call msg_warning ("get_alpha_s: QCD type is not running!" // & " Assuming 5-flavors and LO (1-loop) running!") nf = 5 order = 0 end select if (order > 0) alpha_s = improve_nll_accuracy (alpha_s, nf) end function get_alpha_s @ %def get_alpha_s @ See Eq. (4.31) in [[0709.2092]]. Should be used everywhere in the Sudakov exponent. <>= pure function improve_nll_accuracy (alpha_s, n_flavors) result (alpha_s_imp) real(default) :: alpha_s_imp real(default), intent(in) :: alpha_s integer, intent(in) :: n_flavors alpha_s_imp = alpha_s * (one + alpha_s / (two*pi) * & ((67.0_default/18 - pi**2/6) * CA - five/9 * n_flavors)) end function improve_nll_accuracy @ %def improve_nll_accuracy @ This is fixed to $n_f=5$ for radiation generation. It will be reweighted to the more precise $\alpha_s$. <>= real(default), parameter :: b0rad = (33 - 2 * 5) / (12 * pi) @ %def b0rad @ <>= procedure :: alpha_s_rad => sudakov_alpha_s_rad +<>= + elemental module function sudakov_alpha_s_rad (sudakov, scale2) result (y) + real(default) :: y + class(sudakov_t), intent(in) :: sudakov + real(default), intent(in) :: scale2 + end function sudakov_alpha_s_rad <>= - elemental function sudakov_alpha_s_rad (sudakov, scale2) result (y) + elemental module function sudakov_alpha_s_rad (sudakov, scale2) result (y) real(default) :: y class(sudakov_t), intent(in) :: sudakov real(default), intent(in) :: scale2 y = one / (b0rad * log (scale2 / sudakov%process_deps%lambda2_gen)) end function sudakov_alpha_s_rad @ %def sudakov_alpha_s_rad @ <>= procedure :: reweight_alpha_s => sudakov_reweight_alpha_s +<>= + module function sudakov_reweight_alpha_s (sudakov, pt2) result (accepted) + logical :: accepted + class(sudakov_t), intent(inout) :: sudakov + real(default), intent(in) :: pt2 + end function sudakov_reweight_alpha_s <>= - function sudakov_reweight_alpha_s (sudakov, pt2) result (accepted) + module function sudakov_reweight_alpha_s (sudakov, pt2) result (accepted) logical :: accepted class(sudakov_t), intent(inout) :: sudakov real(default), intent(in) :: pt2 real(default) :: alpha_s_true, alpha_s_rad logical :: alpha_s_equal if (debug_on) call msg_debug2 (D_MATCHING, "reweight_alpha_s") alpha_s_true = get_alpha_s (sudakov%qcd, pt2) alpha_s_rad = sudakov%alpha_s_rad (pt2) call sudakov%rng%generate (sudakov%random) alpha_s_equal = nearly_equal (alpha_s_true, alpha_s_rad) accepted = alpha_s_equal .or. sudakov%random < alpha_s_true / alpha_s_rad if (debug2_active (D_MATCHING)) then print *, ' sudakov%random = ', sudakov%random print *, ' alpha_s_true = ', alpha_s_true print *, ' alpha_s_rad = ', alpha_s_rad print *, ' alpha_s accepted = ', accepted if (alpha_s_rad < alpha_s_true .and. .not. alpha_s_equal) then print *, 'pt2 = ', pt2 print *, 'sudakov%process_deps%lambda2_gen = ', & sudakov%process_deps%lambda2_gen call msg_fatal ("sudakov_reweight_alpha_s: This should never happen. & &Have you chosen a running alpha_s?") end if end if end function sudakov_reweight_alpha_s @ %def sudakov_reweight_alpha_s @ \subsection{POWHEG hook} We provide a POWHEG hook to be called by [[process_instance_evaluate]] to prefill the adaptation grid. We store the actual [[powheg]] object, which does the computations. <>= public :: powheg_matching_hook_t <>= type, extends(process_instance_hook_t) :: powheg_matching_hook_t type(string_t) :: process_name type(powheg_matching_t) :: powheg contains <> end type powheg_matching_hook_t @ %def powheg_matching_t @ Init the hook. The init procedure will be called in [[setup_process]], after everything is set up. Additionally, we have to include [[var_list]] in order to retrieve the grid size in [[xi]] and [[y]]. <>= procedure :: init => powheg_matching_hook_init -<>= - subroutine powheg_matching_hook_init (hook, var_list, instance, pdf_data) - class(powheg_matching_hook_t), intent(inout), target :: hook - type(var_list_t), intent(in) :: var_list - class(process_instance_t), intent(in), target :: instance - type(pdf_data_t), intent(in), optional :: pdf_data - if (debug_on) call msg_debug (D_MATCHING, "powheg_matching_hook_init") - hook%process_name = instance%get_process_name () - call hook%powheg%init (var_list, hook%process_name) - hook%powheg%qcd = instance%get_qcd () - call hook%powheg%connect (instance) - hook%powheg%process_deps%lab_is_cm = & - hook%powheg%process_instance%lab_is_cm (1) - if (present(pdf_data)) then - hook%powheg%process_deps%pdf_data = pdf_data - end if - call hook%powheg%prepare_for_events () - end subroutine powheg_matching_hook_init +<>= + module subroutine powheg_matching_hook_init (hook, var_list, & + instance, pdf_data) + class(powheg_matching_hook_t), intent(inout), target :: hook + type(var_list_t), intent(in) :: var_list + class(process_instance_t), intent(in), target :: instance + type(pdf_data_t), intent(in), optional :: pdf_data + end subroutine powheg_matching_hook_init +<>= + module subroutine powheg_matching_hook_init (hook, var_list, & + instance, pdf_data) + class(powheg_matching_hook_t), intent(inout), target :: hook + type(var_list_t), intent(in) :: var_list + class(process_instance_t), intent(in), target :: instance + type(pdf_data_t), intent(in), optional :: pdf_data + if (debug_on) call msg_debug (D_MATCHING, "powheg_matching_hook_init") + hook%process_name = instance%get_process_name () + call hook%powheg%init (var_list, hook%process_name) + hook%powheg%qcd = instance%get_qcd () + call hook%powheg%connect (instance) + hook%powheg%process_deps%lab_is_cm = & + hook%powheg%process_instance%lab_is_cm (1) + if (present(pdf_data)) then + hook%powheg%process_deps%pdf_data = pdf_data + end if + call hook%powheg%prepare_for_events () + end subroutine powheg_matching_hook_init @ %def powheg_matching_hook_init @ We save the filled grid to file, such that it can be retrieved later on. The hook object will be deallocated when the instance gets deallocated. <>= procedure :: final => powheg_matching_hook_final +<>= + module subroutine powheg_matching_hook_final (hook) + class(powheg_matching_hook_t), intent(inout) :: hook + end subroutine powheg_matching_hook_final <>= - subroutine powheg_matching_hook_final (hook) - class(powheg_matching_hook_t), intent(inout) :: hook - type(string_t) :: filename - if (debug_on) call msg_debug (D_MATCHING, "powheg_matching_hook_final") + module subroutine powheg_matching_hook_final (hook) + class(powheg_matching_hook_t), intent(inout) :: hook + type(string_t) :: filename + if (debug_on) call msg_debug (D_MATCHING, "powheg_matching_hook_final") <> - call hook%powheg%save_grids () + call hook%powheg%save_grids () end subroutine powheg_matching_hook_final @ %def powheg_matching_hook_final @ <>= @ Reduce all grids to a single grid by using [[MPI_MAX]] on each element. <>= call hook%powheg%grid%mpi_reduce (MPI_MAX) @ This routine fills each bin of the POWHEG normalization grid $\{\tilde\xi, \tilde y\}$ @ such that \begin{equation} N(\{\tilde\xi, \tilde y\}) = \max_{\forall \xi,y \in \{\tilde\xi, \tilde y\}} \frac{\mathcal{R}(\xi,y)\mathcal{J}(\xi,y)}{\mathcal{B}\,U(\xi,y,\alpha_s^\text{true})} \end{equation} where $\mathcal{R}(\xi,y)$ is the real cross section [[sqme_real]], $\mathcal{B}$ is the Born cross section [[sqme_born]], $\mathcal{J}$ is the Jacobian and $U$ the upper bounding function [[ubf]]. For each underlying Born $f_b$, there is a number of radiation regions. A radiation region rr may correspond to multiple $\alpha_r$s. The phase space only depends upon the radiation region kinematics rr and not on the specific $\alpha_r$. $\alpha_r$ can be picked in the set $\{\alpha_r|f_b,\text{rr}\}$ proportional to their $R_{\alpha_r}$. For now, we simplify things though and just work with the $\alpha_r$. References: \begin{itemize} \item \texttt{[1002.2581]}, Sec. 7.1 \item Bijans Thesis, Sec. 3.3 and B.2 \end{itemize} In cases where just the Born kinematic failed the cuts, $\mathcal{B}$ was artificially set to zero. Unfortunately, we lost information about the value of [[sqme_born]] in these cases but we assume that the generator cuts applied for parton showering are weak cuts to avoid divergent phase space regions, thus $\mathcal{B}$ is large and does not influence $N$ so we can skip these phase space points. A more appropriate way to solve this problem of Born zeroes is given by the real partition. With real partition, there are two possible cases if $\mathcal{B} = 0$: if the real phase space point belongs to the real singular, the real kinematic is Born-like and very likely also fails the cuts. If the real kinematic is not Born-like, it belongs to the real finite. In both cases, it is $\mathcal{R} = \mathcal{B} = 0$ and we can savely skip the phase space point. In cases where just the real kinematic failed the cuts, $\mathcal{R}$ was artificially set to zero and we have no way to recover it. We skip these phase space points too. If we apply the same cuts during integration and event generation, these points are not relevant. We also skip points where either matrix element is negative. If just one of both is negative, their ratio is negative and thus cannot possibly be relevant for the grid. If both are negative, we are in an unphysical region for the PDFs. We need to divide the real matrix element by [[xi_max]] to correct for the difference in [[apply_kinematic_factors_radiation]] between the real matrix element when integrating and when POWHEG matching. We also need to correct all scale dependences of [[sqme_born]] and [[sqme_real]] from the process' scale to $p_T$ as this is the scale of the entire Sudakov. <>= procedure :: evaluate => powheg_matching_hook_evaluate +<>= + module subroutine powheg_matching_hook_evaluate (hook, instance) + class(powheg_matching_hook_t), intent(inout) :: hook + class(process_instance_t), intent(in), target :: instance + end subroutine powheg_matching_hook_evaluate <>= - subroutine powheg_matching_hook_evaluate (hook, instance) - class(powheg_matching_hook_t), intent(inout) :: hook - class(process_instance_t), intent(in), target :: instance - type(vector4_t), dimension(:), allocatable :: p_hard_born - real(default), dimension(:), allocatable :: x - real(default) :: kt2, xi, y, xi_max, onepy, onemy - real(default), dimension(2) :: x_real - real(default) :: sqme_real, sqme_born - real(default) :: muren_born, mufac_born, muren_real, mufac_real - real(default) :: alpha_s_kt, alpha_s_muren_born, alpha_s_muren_real - real(default) :: f_alr, norm, ubf - real(double), dimension(-6:6) :: pdf_dbl - real(default), dimension(2) :: pdf_born_mufac, pdf_born_kt, pdf_real_mufac, pdf_real_kt - real(double) :: x_dbl, q_dbl - integer, dimension(2) :: flv_born, flv_real - integer :: alr, i_phs, n_x, i_term_born, i_term_real, i_real, emitter, em, n_in - integer :: alphas_power - logical :: is_isr - if (instance%get_active_component_type () == COMP_REAL_FIN) return - associate (powheg => hook%powheg) - allocate (p_hard_born (size (instance%get_p_hard (1)))) - p_hard_born = instance%get_p_hard (1) - call powheg%update (p_hard_born) - do alr = 1, powheg%process_deps%n_alr - i_phs = powheg%process_deps%alr_to_i_phs (alr) - i_term_real = powheg%process_deps%i_term_real(i_phs) - i_term_born = powheg%process_deps%i_term_born - select type (phs => instance%kin(i_term_real)%phs) - type is (phs_fks_t) - call phs%generator%get_radiation_variables (i_phs, xi, y) - end select - call powheg%update_sudakovs (alr, i_phs, y) - if (powheg%sudakov(alr)%s%kt2 (xi, y) >= powheg%settings%pt2_min) then - select type (pcm_work => instance%pcm_work) - class is (pcm_nlo_workspace_t) - xi_max = pcm_work%get_xi_max (alr) - end select - sqme_born = instance%get_sqme (powheg%process_deps%i_term_born) - if (debug_active (D_MATCHING)) then - call assert_equal (OUTPUT_UNIT, sqme_born, & - powheg%event_deps%sqme_born, & - "sqme_born stored in instance and event_deps do not coincide!") - end if - sqme_real = instance%get_sqme (i_term_real) - if (sqme_born <= 0 .or. sqme_real <= 0) then - cycle - end if - n_in = instance%kin(i_term_real)%n_in - emitter = instance%kin(i_term_real)%emitter - is_isr = emitter <= n_in - sqme_real = sqme_real / xi_max - associate (s => powheg%sudakov(alr)%s) - kt2 = s%kt2(xi, y) - alpha_s_kt = s%alpha_s (kt2, use_correct=.true.) - muren_born = instance%term(i_term_born)%get_ren_scale () - mufac_born = instance%term(i_term_born)%get_fac_scale () - muren_real = instance%term(i_term_real)%get_ren_scale () - mufac_real = instance%term(i_term_real)%get_fac_scale () - alpha_s_muren_born = s%alpha_s (muren_born**2, use_correct=.true.) - alpha_s_muren_real = s%alpha_s (muren_real**2, use_correct=.true.) - ubf = s%upper_bound_func (xi, y, alpha_s_kt) - end associate - - !!! Correct all factors of alphas(muren) to alphas(kt). - alphas_power = powheg%process_deps%alphas_power - sqme_born = (alpha_s_kt / alpha_s_muren_born)**alphas_power * sqme_born - sqme_real = (alpha_s_kt / alpha_s_muren_real)**(alphas_power+1) * sqme_real + module subroutine powheg_matching_hook_evaluate (hook, instance) + class(powheg_matching_hook_t), intent(inout) :: hook + class(process_instance_t), intent(in), target :: instance + type(vector4_t), dimension(:), allocatable :: p_hard_born + real(default), dimension(:), allocatable :: x + real(default) :: kt2, xi, y, xi_max, onepy, onemy + real(default), dimension(2) :: x_real + real(default) :: sqme_real, sqme_born + real(default) :: muren_born, mufac_born, muren_real, mufac_real + real(default) :: alpha_s_kt, alpha_s_muren_born, alpha_s_muren_real + real(default) :: f_alr, norm, ubf + real(double), dimension(-6:6) :: pdf_dbl + real(default), dimension(2) :: pdf_born_mufac, pdf_born_kt, & + pdf_real_mufac, pdf_real_kt + real(double) :: x_dbl, q_dbl + integer, dimension(2) :: flv_born, flv_real + integer :: alr, i_phs, n_x, i_term_born, i_term_real, i_real, & + emitter, em, n_in + integer :: alphas_power + logical :: is_isr + if (instance%get_active_component_type () == COMP_REAL_FIN) return + associate (powheg => hook%powheg) + allocate (p_hard_born (size (instance%get_p_hard (1)))) + p_hard_born = instance%get_p_hard (1) + call powheg%update (p_hard_born) + do alr = 1, powheg%process_deps%n_alr + i_phs = powheg%process_deps%alr_to_i_phs (alr) + i_term_real = powheg%process_deps%i_term_real(i_phs) + i_term_born = powheg%process_deps%i_term_born + select type (phs => instance%kin(i_term_real)%phs) + type is (phs_fks_t) + call phs%generator%get_radiation_variables (i_phs, xi, y) + end select + call powheg%update_sudakovs (alr, i_phs, y) + if (powheg%sudakov(alr)%s%kt2 (xi, y) >= powheg%settings%pt2_min) then + select type (pcm_work => instance%pcm_work) + class is (pcm_nlo_workspace_t) + xi_max = pcm_work%get_xi_max (alr) + end select + sqme_born = instance%get_sqme (powheg%process_deps%i_term_born) + if (debug_active (D_MATCHING)) then + call assert_equal (OUTPUT_UNIT, sqme_born, & + powheg%event_deps%sqme_born, "sqme_born " // & + "stored in instance and event_deps do not coincide!") + end if + sqme_real = instance%get_sqme (i_term_real) + if (sqme_born <= 0 .or. sqme_real <= 0) then + cycle + end if + n_in = instance%kin(i_term_real)%n_in + emitter = instance%kin(i_term_real)%emitter + is_isr = emitter <= n_in + sqme_real = sqme_real / xi_max + associate (s => powheg%sudakov(alr)%s) + kt2 = s%kt2(xi, y) + alpha_s_kt = s%alpha_s (kt2, use_correct=.true.) + muren_born = instance%term(i_term_born)%get_ren_scale () + mufac_born = instance%term(i_term_born)%get_fac_scale () + muren_real = instance%term(i_term_real)%get_ren_scale () + mufac_real = instance%term(i_term_real)%get_fac_scale () + alpha_s_muren_born = s%alpha_s (muren_born**2, use_correct=.true.) + alpha_s_muren_real = s%alpha_s (muren_real**2, use_correct=.true.) + ubf = s%upper_bound_func (xi, y, alpha_s_kt) + end associate + + !!! Correct all factors of alphas(muren) to alphas(kt). + alphas_power = powheg%process_deps%alphas_power + sqme_born = (alpha_s_kt / alpha_s_muren_born)**alphas_power * sqme_born + sqme_real = (alpha_s_kt / alpha_s_muren_real)**(alphas_power+1) * sqme_real !!! Also correct the PDFs previously computed at mufac instead of kt - select type (pcm => instance%term(i_term_real)%pcm) - type is (pcm_nlo_t) - if (pcm%has_pdfs) then - associate (reg_data => pcm%region_data) - flv_born = reg_data%flv_born(i_term_born)%flst(1:2) - where (flv_born == 21) flv_born = 0 - i_real = reg_data%regions(alr)%real_index - flv_real = reg_data%flv_real(i_real)%flst(1:2) - where (flv_real == 21) flv_real = 0 - end associate - associate (pdf_data => powheg%process_deps%pdf_data, & - x_born => powheg%event_deps%x_born) - if (.not. pdfs_valid(instance, pdf_data, i_term_born, i_term_real, & - sqrt(kt2), x_born, is_isr)) cycle - do em = 1, 2 - x_dbl = x_born(em) ; q_dbl = mufac_born - call pdf_data%evolve(x_dbl, q_dbl, pdf_dbl) - pdf_born_mufac(em) = pdf_dbl(flv_born(em)) / x_born(em) - x_dbl = x_born(em) ; q_dbl = sqrt(kt2) - call pdf_data%evolve(x_dbl, q_dbl, pdf_dbl) - pdf_born_kt(em) = pdf_dbl(flv_born(em)) / x_born(em) - end do - onepy = one + y; onemy = one - y - x_real(1) = x_born(1) * sqrt ((two - xi * onemy) / (two - xi * onepy)) - x_real(2) = x_born(2) * sqrt ((two - xi * onepy) / (two - xi * onemy)) - x_real = x_real / sqrt (one - xi) - do em = 1, 2 - x_dbl = x_real(em) ; q_dbl = mufac_real - call pdf_data%evolve(x_dbl, q_dbl, pdf_dbl) - pdf_real_mufac(em) = pdf_dbl(flv_real(em)) / x_real(em) - x_dbl = x_real(em) ; q_dbl = sqrt(kt2) - call pdf_data%evolve(x_dbl, q_dbl, pdf_dbl) - pdf_real_kt(em) = pdf_dbl(flv_real(em)) / x_real(em) - end do - sqme_born = pdf_born_kt(1) * pdf_born_kt(2) / & - (pdf_born_mufac(1) * pdf_born_mufac(2)) * sqme_born - sqme_real = pdf_real_kt(1) * pdf_real_kt(2) / & - (pdf_real_mufac(1) * pdf_real_mufac(2)) * sqme_real - end associate - end if - end select - - - norm = sqme_real / (sqme_born * ubf) - f_alr = (one * alr) / powheg%process_deps%n_alr - tiny_07 - if (.not. allocated (x)) & - allocate (x(size (instance%get_x_process ()))) - x = instance%get_x_process () - n_x = size (x) - call powheg%grid%update_maxima & - ([x(n_x - 2 : n_x - 1), f_alr], norm) - end if - end do - end associate + select type (pcm => instance%term(i_term_real)%pcm) + type is (pcm_nlo_t) + if (pcm%has_pdfs) then + associate (reg_data => pcm%region_data) + flv_born = reg_data%flv_born(i_term_born)%flst(1:2) + where (flv_born == 21) flv_born = 0 + i_real = reg_data%regions(alr)%real_index + flv_real = reg_data%flv_real(i_real)%flst(1:2) + where (flv_real == 21) flv_real = 0 + end associate + associate (pdf_data => powheg%process_deps%pdf_data, & + x_born => powheg%event_deps%x_born) + if (.not. pdfs_valid(instance, pdf_data, i_term_born, i_term_real, & + sqrt(kt2), x_born, is_isr)) cycle + do em = 1, 2 + x_dbl = x_born(em) ; q_dbl = mufac_born + call pdf_data%evolve(x_dbl, q_dbl, pdf_dbl) + pdf_born_mufac(em) = pdf_dbl(flv_born(em)) / x_born(em) + x_dbl = x_born(em) ; q_dbl = sqrt(kt2) + call pdf_data%evolve(x_dbl, q_dbl, pdf_dbl) + pdf_born_kt(em) = pdf_dbl(flv_born(em)) / x_born(em) + end do + onepy = one + y; onemy = one - y + x_real(1) = x_born(1) * sqrt ((two - xi * onemy) / (two - xi * onepy)) + x_real(2) = x_born(2) * sqrt ((two - xi * onepy) / (two - xi * onemy)) + x_real = x_real / sqrt (one - xi) + do em = 1, 2 + x_dbl = x_real(em) ; q_dbl = mufac_real + call pdf_data%evolve(x_dbl, q_dbl, pdf_dbl) + pdf_real_mufac(em) = pdf_dbl(flv_real(em)) / x_real(em) + x_dbl = x_real(em) ; q_dbl = sqrt(kt2) + call pdf_data%evolve(x_dbl, q_dbl, pdf_dbl) + pdf_real_kt(em) = pdf_dbl(flv_real(em)) / x_real(em) + end do + sqme_born = pdf_born_kt(1) * pdf_born_kt(2) / & + (pdf_born_mufac(1) * pdf_born_mufac(2)) * sqme_born + sqme_real = pdf_real_kt(1) * pdf_real_kt(2) / & + (pdf_real_mufac(1) * pdf_real_mufac(2)) * sqme_real + end associate + end if + end select + + norm = sqme_real / (sqme_born * ubf) + f_alr = (one * alr) / powheg%process_deps%n_alr - tiny_07 + if (.not. allocated (x)) & + allocate (x(size (instance%get_x_process ()))) + x = instance%get_x_process () + n_x = size (x) + call powheg%grid%update_maxima & + ([x(n_x - 2 : n_x - 1), f_alr], norm) + end if + end do + end associate contains <> end subroutine powheg_matching_hook_evaluate @ %def powheg_matching_hook_evaluate <>= @ In some phase space point, the quark PDF may vanish although the gluon PDF does not. This may lead to very large ratios $\mathcal{R} / \mathcal{B}$ which may spoil the reweighting performance of the POWHEG grid. Excluding these values when filling the grid is a minor approximation with a huge performance gain. The PDF veto we perform is independent on the actual emitter to not artificially favor specific ALRs. <>= function pdfs_valid (instance, pdf_data, i_term_born, i_term_real, q, x, is_isr) result (valid) logical :: valid class(process_instance_t), intent(in), target :: instance type(pdf_data_t), intent(inout) :: pdf_data integer, intent(in) :: i_term_born, i_term_real real(default), intent(in) :: q real(default), dimension(2), intent(in) :: x real(double) :: q_dbl real(double), dimension(2) :: x_dbl logical, intent(in) :: is_isr real(default) :: sum_pdf_q real(default), dimension(-6:6) :: pdf real(double), dimension(-6:6) :: pdf_dbl integer :: flv_born_em, em, i_q logical :: warned_once = .false. valid = .true. if (is_isr) then if (q**2 < two .or. any(x > 0.9_default)) then valid = .false. return end if x_dbl = x q_dbl = q do em = 1, 2 call pdf_data%evolve(x_dbl(em), q_dbl, pdf_dbl) pdf = pdf_dbl / x(em) select type (pcm => instance%term(i_term_real)%pcm) type is (pcm_nlo_t) flv_born_em = pcm%region_data%flv_born(i_term_born)%flst(em) end select if (is_gluon(flv_born_em)) then sum_pdf_q = 0 do i_q = 1, 6 sum_pdf_q = sum_pdf_q + pdf(-i_q) + pdf(i_q) end do if (sum_pdf_q * x(em) * (1-x(em)) > 30._default * pdf(0)) then valid = .false. return end if elseif (is_quark(flv_born_em)) then if (pdf(0) * x(em) * (1-x(em)) > 30._default * pdf(flv_born_em)) then valid = .false. return end if else if (.not. warned_once) then call msg_warning ("powheg_matching_hook_evaluate: unexpected emitter flavor") warned_once = .true. end if end if end do end if end function pdfs_valid @ %def pdfs_valid @ \subsection{Unit tests} Test module, followed by the corresponding implementation module. <<[[powheg_matching_ut.f90]]>>= <> module powheg_matching_ut use unit_tests use powheg_matching_uti <> <> contains <> end module powheg_matching_ut @ %def powheg_matching_ut @ <<[[powheg_matching_uti.f90]]>>= <> module powheg_matching_uti <> <> use constants, only: zero, one use lorentz use physics_defs, only: LAMBDA_QCD_REF use sm_qcd use subevents, only: PRT_INCOMING, PRT_OUTGOING use model_data use particles use rng_base use variables use shower_base use shower_core use powheg_matching use rng_base_ut, only: rng_test_factory_t <> <> contains <> end module powheg_matching_uti @ %def powheg_matching_ut @ API: driver for the unit tests below. <>= public :: powheg_test <>= subroutine powheg_test (u, results) integer, intent(in) :: u type(test_results_t), intent(inout) :: results <> end subroutine powheg_test @ %def powheg_test @ \subsubsection{Initialization} There are no Powheg unit tests so far. \subsubsection{Compare generated emission with Sudakov form factor} This is a nontrivial test of the generation algorithm and should be independent of the used upper bounding function (as long as all singularities are included). It runs for an extensive amount of time generating a file [[sudakov.dat]] containing a histogram for the Sudakov factor in bins of the radiation's transverse momentum $p_T^2$. So far, it is only designed to work for lepton collisions without beam spectra and cuts. In its current state however, it fills the histogram with either $1$ or [[NaN]] as the [[p_hard]] received in [[term_instance_evaluate_interaction]] is identically zero. <>= procedure :: test_sudakov => powheg_test_sudakov +<>= + module subroutine powheg_test_sudakov (powheg) + class(powheg_matching_t), intent(inout) :: powheg + end subroutine powheg_test_sudakov <>= - subroutine powheg_test_sudakov (powheg) + module subroutine powheg_test_sudakov (powheg) class(powheg_matching_t), intent(inout) :: powheg integer :: n_calls1, n_calls2 integer, parameter :: n_bins = 20 real(default) :: sqme_real_x_jacobian, sqme_born type(vector4_t), dimension(:), allocatable :: p_born real(default), dimension(3) :: random real(default) :: xi, y integer :: i_call, i_bin, alr, emitter, i_phs real(default) :: alpha_s, kT2, weight real(default) :: pt2_min, s, random_jacobian real(default), dimension(n_bins) :: histo1, histo2, histo1sq, histo2sq real(default), dimension(n_bins) :: tmp integer :: i_strip, n_in_strip, n_strips real(default), dimension(n_bins) :: average, average_sq, error real(default), dimension(n_bins) :: & sudakov_0, sudakov_p, sudakov_m, rel_error integer :: u p_born = powheg%event_deps%p_born_lab%get_momenta (1) sqme_born = powheg%event_deps%sqme_born s = powheg%event_deps%s_hat pt2_min = powheg%settings%pt2_min n_calls1 = 100000; n_calls2 = 1000000 histo1 = zero; histo2 = zero; histo1sq = zero; histo2sq = zero n_strips = 10 call compute_integrals () call generate_emissions () call write_to_screen_and_file () contains <> end subroutine powheg_test_sudakov @ %def powheg_test_sudakov @ This determines the binning of the Sudakov histogram. Linear and logarithmic binning are available. <>= pure function binning (i) result (pt2) real(default) :: pt2 integer, intent(in) :: i !pt2 = pt2_min + (s-pt2_min) * (i-1) / (n_bins-1) pt2 = pt2_min * exp (log (s / pt2_min) * (i-1) / (n_bins-1)) end function @ %def binning @ <>= subroutine compute_integrals () write (msg_buffer, "(A)") "POWHEG: test_sudakov: Computing integrals" call msg_message () select type (pcm_work => powheg%process_instance%pcm_work) class is (pcm_nlo_workspace_t) associate (fks => powheg%phs_fks_generator) do i_call = 1, n_calls1 do alr = 1, powheg%process_deps%n_alr call powheg%rng%generate (random) select type (pcm => powheg%process_instance%pcm) type is (pcm_nlo_t) emitter = pcm%region_data%get_emitter (alr) i_phs = powheg%process_deps%alr_to_i_phs(alr) end select !!! The sudakov test works only with lepton collisions without beam spectra !!! so we can identify the cms and lab momenta. powheg%event_deps%p_real_lab = powheg%event_deps%p_real_cms call powheg%copy_momenta (i_phs) call fks%get_radiation_variables (i_phs, xi, y) kT2 = powheg%sudakov(alr)%s%kt2(xi, y) if (kT2 >= pt2_min .and. xi < one - tiny_07) then alpha_s = get_alpha_s (powheg%qcd, kT2) sqme_real_x_jacobian = powheg%compute_sqme_real (alr, sqrt(kT2)) random_jacobian = pcm_work%real_kinematics%jac_rand (emitter) weight = sqme_real_x_jacobian * random_jacobian / sqme_born do i_bin = 1, n_bins if (kT2 > binning(i_bin)) then histo1(i_bin) = histo1(i_bin) + weight histo1sq(i_bin) = histo1sq(i_bin) + weight**2 end if end do end if ! Do not cycle since there is a Heaviside in the exponent end do call msg_show_progress (i_call, n_calls1) end do end associate end select average = histo1 / n_calls1 average_sq = histo1sq / n_calls1 error = sqrt ((average_sq - average**2) / n_calls1) sudakov_0 = exp(-average) sudakov_p = exp(-(average + error)) sudakov_m = exp(-(average - error)) rel_error = (sudakov_m - sudakov_p) / (2 * sudakov_0) * 100 end subroutine compute_integrals @ %def compute_integrals @ <>= subroutine generate_emissions () write (msg_buffer, "(A)") "POWHEG: test_sudakov: Generating emissions" call msg_message () do i_strip = 1, n_strips tmp = 0 n_in_strip = n_calls2 / n_strips do i_call = 1, n_in_strip if (signal_is_pending ()) return call powheg%generate_emission (pt2_generated = kT2) do i_bin = 1, n_bins if (kT2 > binning(i_bin)) then tmp(i_bin) = tmp(i_bin) + 1 end if end do end do tmp = one - (one * tmp) / n_in_strip histo2 = histo2 + tmp histo2sq = histo2sq + tmp**2 call msg_show_progress (i_strip, n_strips) end do average = histo2 / n_strips average_sq = histo2sq / n_strips error = sqrt ((average_sq - average**2) / n_strips) end subroutine generate_emissions @ %def generate_emissions @ <>= subroutine write_to_screen_and_file () u = free_unit () open (file='sudakov.dat', unit=u, action='write') print *, 'exp(-Integrated R/B)-distribution: ' print *, 'pT2 sudakov_+ sudakov_0 sudakov_- rel_err[%]: ' do i_bin = 1, n_bins print *, binning(i_bin), & sudakov_p(i_bin), sudakov_0(i_bin), sudakov_m(i_bin), & rel_error(i_bin) write (u, "(6(" // FMT_16 // ",2X))") binning(i_bin), & sudakov_p(i_bin), sudakov_0(i_bin), sudakov_m(i_bin), & average(i_bin), error(i_bin) end do close (u) print *, '*******************************' print *, 'Noemission probability: ' do i_bin = 1, n_bins print *, binning (i_bin), average (i_bin), error(i_bin) end do end subroutine write_to_screen_and_file @ %def write_to_screen_and_file @ Index: trunk/share/debug/Makefile_full =================================================================== --- trunk/share/debug/Makefile_full (revision 8803) +++ trunk/share/debug/Makefile_full (revision 8804) @@ -1,687 +1,690 @@ FC=pgfortran_2019 FCFLAGS=-Mbackslash CC=gcc CCFLAGS= MODELS = \ SM.mdl \ SM_hadrons.mdl \ Test.mdl CC_SRC = \ sprintf_interface.c \ signal_interface.c F77_SRC = \ pythia.F \ pythia_pdf.f \ pythia6_up.f \ toppik.f \ toppik_axial.f FC0_SRC = FC_SRC = \ format_defs.f90 \ io_units.f90 \ kinds.f90 \ constants.f90 \ iso_varying_string.f90 \ unit_tests.f90 \ unit_tests_sub.f90 \ numeric_utils.f90 \ numeric_utils_sub.f90 \ system_dependencies.f90 \ string_utils.f90 \ string_utils_sub.f90 \ system_defs.f90 \ system_defs_sub.f90 \ debug_master.f90 \ diagnostics.f90 \ diagnostics_sub.f90 \ sorting.f90 \ physics_defs.f90 \ physics_defs_sub.f90 \ pdg_arrays.f90 \ bytes.f90 \ hashes.f90 \ md5.f90 \ model_data.f90 \ model_data_sub.f90 \ auto_components.f90 \ auto_components_sub.f90 \ var_base.f90 \ model_testbed.f90 \ auto_components_uti.f90 \ auto_components_ut.f90 \ os_interface.f90 \ os_interface_sub.f90 \ c_particles.f90 \ c_particles_sub.f90 \ format_utils.f90 \ lorentz.f90 \ lorentz_sub.f90 \ phs_points.f90 \ phs_points_sub.f90 \ colors.f90 \ colors_sub.f90 \ flavors.f90 \ flavors_sub.f90 \ helicities.f90 \ helicities_sub.f90 \ quantum_numbers.f90 \ quantum_numbers_sub.f90 \ state_matrices.f90 \ state_matrices_sub.f90 \ interactions.f90 \ interactions_sub.f90 \ CppStringsWrap_dummy.f90 \ FastjetWrap_dummy.f90 \ cpp_strings.f90 \ cpp_strings_sub.f90 \ fastjet.f90 \ fastjet_sub.f90 \ jets.f90 \ subevents.f90 \ su_algebra.f90 \ su_algebra_sub.f90 \ bloch_vectors.f90 \ bloch_vectors_sub.f90 \ polarizations.f90 \ polarizations_sub.f90 \ particles.f90 \ particles_sub.f90 \ event_base.f90 \ event_base_sub.f90 \ eio_data.f90 \ eio_data_sub.f90 \ event_handles.f90 \ eio_base.f90 \ eio_base_sub.f90 \ eio_base_uti.f90 \ eio_base_ut.f90 \ variables.f90 \ variables_sub.f90 \ rng_base.f90 \ rng_base_sub.f90 \ tao_random_numbers.f90 \ rng_tao.f90 \ rng_tao_sub.f90 \ rng_stream.f90 \ rng_stream_sub.f90 \ rng_base_uti.f90 \ rng_base_ut.f90 \ dispatch_rng.f90 \ dispatch_rng_sub.f90 \ dispatch_rng_uti.f90 \ dispatch_rng_ut.f90 \ beam_structures.f90 \ beam_structures_sub.f90 \ evaluators.f90 \ evaluators_sub.f90 \ beams.f90 \ beams_sub.f90 \ sm_physics.f90 \ sm_physics_sub.f90 \ file_registries.f90 \ file_registries_sub.f90 \ sf_aux.f90 \ sf_aux_sub.f90 \ sf_mappings.f90 \ sf_mappings_sub.f90 \ sf_base.f90 \ sf_base_sub.f90 \ electron_pdfs.f90 \ electron_pdfs_sub.f90 \ sf_isr.f90 \ sf_isr_sub.f90 \ sf_epa.f90 \ sf_epa_sub.f90 \ sf_ewa.f90 \ sf_ewa_sub.f90 \ sf_escan.f90 \ sf_escan_sub.f90 \ sf_gaussian.f90 \ sf_gaussian_sub.f90 \ sf_beam_events.f90 \ sf_beam_events_sub.f90 \ circe1.f90 \ sf_circe1.f90 \ sf_circe1_sub.f90 \ circe2.f90 \ selectors.f90 \ selectors_sub.f90 \ sf_circe2.f90 \ sf_circe2_sub.f90 \ sm_qcd.f90 \ sm_qcd_sub.f90 \ sm_qed.f90 \ sm_qed_sub.f90 \ mrst2004qed.f90 \ cteq6pdf.f90 \ mstwpdf.f90 \ ct10pdf.f90 \ CJpdf.f90 \ ct14pdf.f90 \ pdf_builtin.f90 \ pdf_builtin_sub.f90 \ LHAPDFWrap_dummy.f90 \ lhapdf5_full_dummy.f90 \ lhapdf5_has_photon_dummy.f90 \ lhapdf.f90 \ hoppet_dummy.f90 \ hoppet_interface.f90 \ sf_pdf_builtin.f90 \ sf_pdf_builtin_sub.f90 \ sf_lhapdf.f90 \ sf_lhapdf_sub.f90 \ dispatch_beams.f90 \ dispatch_beams_sub.f90 \ process_constants.f90 \ process_constants_sub.f90 \ prclib_interfaces.f90 \ prc_core_def.f90 \ prc_core_def_sub.f90 \ particle_specifiers.f90 \ particle_specifiers_sub.f90 \ process_libraries.f90 \ process_libraries_sub.f90 \ prc_test.f90 \ prc_test_sub.f90 \ prc_core.f90 \ prc_core_sub.f90 \ prc_test_core.f90 \ prc_test_core_sub.f90 \ sm_qed.f90 \ sm_qed_sub.f90 \ prc_omega.f90 \ prc_omega_sub.f90 \ phs_base.f90 \ phs_base_sub.f90 \ ifiles.f90 \ lexers.f90 \ syntax_rules.f90 \ parser.f90 \ expr_base.f90 \ formats.f90 \ formats_sub.f90 \ analysis.f90 \ user_code_interface.f90 \ observables.f90 \ observables_sub.f90 \ eval_trees.f90 \ eval_trees_sub.f90 \ interpolation.f90 \ interpolation_sub.f90 \ nr_tools.f90 \ ttv_formfactors.f90 \ ttv_formfactors_use.f90 \ ttv_formfactors_uti.f90 \ ttv_formfactors_ut.f90 \ models.f90 \ prclib_stacks.f90 \ prclib_stacks_sub.f90 \ user_files.f90 \ cputime.f90 \ cputime_sub.f90 \ mci_base.f90 \ mci_base_sub.f90 \ integration_results.f90 \ integration_results_sub.f90 \ integration_results_uti.f90 \ integration_results_ut.f90 \ mappings.f90 \ mappings_sub.f90 \ permutations.f90 \ permutations_sub.f90 \ resonances.f90 \ resonances_sub.f90 \ phs_trees.f90 \ phs_trees_sub.f90 \ phs_forests.f90 \ phs_forests_sub.f90 \ prc_external.f90 \ blha_config.f90 \ blha_config_sub.f90 \ blha_olp_interfaces.f90 \ blha_olp_interfaces_sub.f90 \ prc_openloops.f90 \ prc_openloops_sub.f90 \ prc_threshold.f90 \ prc_threshold_sub.f90 \ process_config.f90 \ process_config_sub.f90 \ process_counter.f90 \ process_counter_sub.f90 \ process_mci.f90 \ pcm_base.f90 \ pcm_base_sub.f90 \ nlo_data.f90 \ nlo_data_sub.f90 \ cascades.f90 \ cascades_sub.f90 \ cascades2_lexer.f90 \ cascades2_lexer_sub.f90 \ cascades2_lexer_uti.f90 \ cascades2_lexer_ut.f90 \ cascades2.f90 \ cascades2_sub.f90 \ cascades2_uti.f90 \ cascades2_ut.f90 \ phs_none.f90 \ phs_none_sub.f90 \ phs_rambo.f90 \ phs_rambo_sub.f90 \ phs_wood.f90 \ phs_wood_sub.f90 \ phs_fks.f90 \ phs_fks_sub.f90 \ phs_single.f90 \ phs_single_sub.f90 \ fks_regions.f90 \ fks_regions_sub.f90 \ virtual.f90 \ virtual_sub.f90 \ pdf.f90 \ pdf_sub.f90 \ real_subtraction.f90 \ real_subtraction_sub.f90 \ dglap_remnant.f90 \ dispatch_fks.f90 \ dispatch_fks_sub.f90 \ dispatch_phase_space.f90 \ dispatch_phase_space_sub.f90 \ pcm.f90 \ pcm_sub.f90 \ recola_wrapper_dummy.f90 \ prc_recola.f90 \ subevt_expr.f90 \ subevt_expr_sub.f90 \ parton_states.f90 \ parton_states_sub.f90 \ prc_template_me.f90 \ prc_template_me_sub.f90 \ process.f90 \ process_sub.f90 \ process_stacks.f90 \ process_stacks_sub.f90 \ iterations.f90 \ rt_data.f90 \ file_utils.f90 \ file_utils_sub.f90 \ prc_gosam.f90 \ prc_gosam_sub.f90 \ dispatch_me_methods.f90 \ sf_base_uti.f90 \ sf_base_ut.f90 \ dispatch_uti.f90 \ dispatch_ut.f90 \ formats_uti.f90 \ formats_ut.f90 \ md5_uti.f90 \ md5_ut.f90 \ os_interface_uti.f90 \ os_interface_ut.f90 \ sorting_uti.f90 \ sorting_ut.f90 \ grids.f90 \ grids_uti.f90 \ grids_ut.f90 \ solver.f90 \ solver_uti.f90 \ solver_ut.f90 \ cputime_uti.f90 \ cputime_ut.f90 \ sm_qcd_uti.f90 \ sm_qcd_ut.f90 \ sm_physics_uti.f90 \ sm_physics_ut.f90 \ lexers_uti.f90 \ lexers_ut.f90 \ parser_uti.f90 \ parser_ut.f90 \ xml.f90 \ xml_uti.f90 \ xml_ut.f90 \ colors_uti.f90 \ colors_ut.f90 \ state_matrices_uti.f90 \ state_matrices_ut.f90 \ analysis_uti.f90 \ analysis_ut.f90 \ particles_uti.f90 \ particles_ut.f90 \ radiation_generator.f90 \ radiation_generator_sub.f90 \ radiation_generator_uti.f90 \ radiation_generator_ut.f90 \ blha_uti.f90 \ blha_ut.f90 \ evaluators_uti.f90 \ evaluators_ut.f90 \ models_uti.f90 \ models_ut.f90 \ eval_trees_uti.f90 \ eval_trees_ut.f90 \ resonances_uti.f90 \ resonances_ut.f90 \ phs_trees_uti.f90 \ phs_trees_ut.f90 \ phs_forests_uti.f90 \ phs_forests_ut.f90 \ beams_uti.f90 \ beams_ut.f90 \ su_algebra_uti.f90 \ su_algebra_ut.f90 \ bloch_vectors_uti.f90 \ bloch_vectors_ut.f90 \ polarizations_uti.f90 \ polarizations_ut.f90 \ sf_aux_uti.f90 \ sf_aux_ut.f90 \ sf_mappings_uti.f90 \ sf_mappings_ut.f90 \ sf_pdf_builtin_uti.f90 \ sf_pdf_builtin_ut.f90 \ sf_lhapdf_uti.f90 \ sf_lhapdf_ut.f90 \ sf_isr_uti.f90 \ sf_isr_ut.f90 \ sf_epa_uti.f90 \ sf_epa_ut.f90 \ sf_ewa_uti.f90 \ sf_ewa_ut.f90 \ sf_circe1_uti.f90 \ sf_circe1_ut.f90 \ sf_circe2_uti.f90 \ sf_circe2_ut.f90 \ sf_gaussian_uti.f90 \ sf_gaussian_ut.f90 \ sf_beam_events_uti.f90 \ sf_beam_events_ut.f90 \ sf_escan_uti.f90 \ sf_escan_ut.f90 \ phs_base_uti.f90 \ phs_base_ut.f90 \ phs_none_uti.f90 \ phs_none_ut.f90 \ phs_single_uti.f90 \ phs_single_ut.f90 \ phs_rambo_uti.f90 \ phs_rambo_ut.f90 \ phs_wood_uti.f90 \ phs_wood_ut.f90 \ phs_fks_uti.f90 \ phs_fks_ut.f90 \ fks_regions_uti.f90 \ fks_regions_ut.f90 \ mci_midpoint.f90 \ mci_midpoint_sub.f90 \ mci_base_uti.f90 \ mci_base_ut.f90 \ mci_midpoint_uti.f90 \ mci_midpoint_ut.f90 \ kinematics.f90 \ kinematics_sub.f90 \ instances.f90 \ instances_sub.f90 \ mci_none.f90 \ mci_none_sub.f90 \ mci_none_uti.f90 \ mci_none_ut.f90 \ processes_uti.f90 \ processes_ut.f90 \ process_stacks_uti.f90 \ process_stacks_ut.f90 \ prc_recola_uti.f90 \ prc_recola_ut.f90 \ rng_tao_uti.f90 \ rng_tao_ut.f90 \ rng_stream_uti.f90 \ rng_stream_ut.f90 \ selectors_uti.f90 \ selectors_ut.f90 \ vegas.f90 \ vegas_sub.f90 \ vegas_uti.f90 \ vegas_ut.f90 \ vamp2.f90 \ vamp2_sub.f90 \ vamp2_uti.f90 \ vamp2_ut.f90 \ exceptions.f90 \ vamp_stat.f90 \ utils.f90 \ divisions.f90 \ linalg.f90 \ vamp.f90 \ mci_vamp.f90 \ mci_vamp_sub.f90 \ mci_vamp_uti.f90 \ mci_vamp_ut.f90 \ mci_vamp2.f90 \ mci_vamp2_sub.f90 \ mci_vamp2_uti.f90 \ mci_vamp2_ut.f90 \ prclib_interfaces_uti.f90 \ prclib_interfaces_ut.f90 \ particle_specifiers_uti.f90 \ particle_specifiers_ut.f90 \ process_libraries_uti.f90 \ process_libraries_ut.f90 \ prclib_stacks_uti.f90 \ prclib_stacks_ut.f90 \ slha_interface.f90 \ slha_interface_sub.f90 \ slha_interface_uti.f90 \ slha_interface_ut.f90 \ cascades_uti.f90 \ cascades_ut.f90 \ prc_test_uti.f90 \ prc_test_ut.f90 \ prc_template_me_uti.f90 \ prc_template_me_ut.f90 \ prc_omega_uti.f90 \ prc_omega_ut.f90 \ event_transforms.f90 \ event_transforms_uti.f90 \ event_transforms_ut.f90 \ hep_common.f90 \ hep_common_sub.f90 \ hepev4_aux.f90 \ tauola_dummy.f90 \ tauola_interface.f90 \ tauola_interface_sub.f90 \ shower_base.f90 \ shower_base_sub.f90 \ shower_partons.f90 \ shower_partons_sub.f90 \ muli.f90 \ matching_base.f90 \ + matching_base_sub.f90 \ powheg_matching.f90 \ shower_core.f90 \ shower_core_sub.f90 \ shower_base_uti.f90 \ shower_base_ut.f90 \ shower.f90 \ shower_uti.f90 \ shower_ut.f90 \ shower_pythia6.f90 \ shower_pythia6_sub.f90 \ whizard_lha.f90 \ whizard_lha_uti.f90 \ whizard_lha_ut.f90 \ LHAWhizard_dummy.f90 \ Pythia8Wrap_dummy.f90 \ pythia8.f90 \ pythia8_uti.f90 \ pythia8_ut.f90 \ shower_pythia8.f90 \ shower_pythia8_sub.f90 \ hadrons.f90 \ ktclus.f90 \ mlm_matching.f90 \ + mlm_matching_sub.f90 \ ckkw_matching.f90 \ + ckkw_matching_sub.f90 \ jets_uti.f90 \ jets_ut.f90 \ pdg_arrays_uti.f90 \ pdg_arrays_ut.f90 \ interactions_uti.f90 \ interactions_ut.f90 \ decays.f90 \ decays_uti.f90 \ decays_ut.f90 \ evt_nlo.f90 \ events.f90 \ events_uti.f90 \ events_ut.f90 \ HepMCWrap_dummy.f90 \ hepmc_interface.f90 \ hepmc_interface_sub.f90 \ hepmc_interface_uti.f90 \ hepmc_interface_ut.f90 \ LCIOWrap_dummy.f90 \ lcio_interface.f90 \ lcio_interface_sub.f90 \ lcio_interface_uti.f90 \ lcio_interface_ut.f90 \ hep_events.f90 \ hep_events_sub.f90 \ hep_events_uti.f90 \ hep_events_ut.f90 \ expr_tests_uti.f90 \ expr_tests_ut.f90 \ parton_states_uti.f90 \ parton_states_ut.f90 \ eio_data_uti.f90 \ eio_data_ut.f90 \ eio_raw.f90 \ eio_raw_uti.f90 \ eio_raw_ut.f90 \ eio_checkpoints.f90 \ eio_checkpoints_sub.f90 \ eio_checkpoints_uti.f90 \ eio_checkpoints_ut.f90 \ eio_lhef.f90 \ eio_lhef_sub.f90 \ eio_lhef_uti.f90 \ eio_lhef_ut.f90 \ eio_hepmc.f90 \ eio_hepmc_sub.f90 \ eio_hepmc_uti.f90 \ eio_hepmc_ut.f90 \ eio_lcio.f90 \ eio_lcio_sub.f90 \ eio_lcio_uti.f90 \ eio_lcio_ut.f90 \ stdhep_dummy.f90 \ xdr_wo_stdhep.f90 \ eio_stdhep.f90 \ eio_stdhep_sub.f90 \ eio_stdhep_uti.f90 \ eio_stdhep_ut.f90 \ eio_ascii.f90 \ eio_ascii_sub.f90 \ eio_ascii_uti.f90 \ eio_ascii_ut.f90 \ eio_weights.f90 \ eio_weights_sub.f90 \ eio_weights_uti.f90 \ eio_weights_ut.f90 \ eio_dump.f90 \ eio_dump_sub.f90 \ eio_dump_uti.f90 \ eio_dump_ut.f90 \ eio_callback.f90 \ eio_callback_sub.f90 \ real_subtraction_uti.f90 \ real_subtraction_ut.f90 \ iterations_uti.f90 \ iterations_ut.f90 \ rt_data_uti.f90 \ rt_data_ut.f90 \ dispatch_mci.f90 \ dispatch_mci_sub.f90 \ dispatch_mci_uti.f90 \ dispatch_mci_ut.f90 \ dispatch_phs_uti.f90 \ dispatch_phs_ut.f90 \ resonance_insertion.f90 \ resonance_insertion_uti.f90 \ resonance_insertion_ut.f90 \ recoil_kinematics.f90 \ recoil_kinematics_uti.f90 \ recoil_kinematics_ut.f90 \ isr_epa_handler.f90 \ isr_epa_handler_uti.f90 \ isr_epa_handler_ut.f90 \ dispatch_transforms.f90 \ dispatch_transforms_uti.f90 \ dispatch_transforms_ut.f90 \ beam_structures_uti.f90 \ beam_structures_ut.f90 \ process_configurations.f90 \ process_configurations_uti.f90 \ process_configurations_ut.f90 \ compilations.f90 \ compilations_uti.f90 \ compilations_ut.f90 \ integrations.f90 \ integrations_uti.f90 \ integrations_ut.f90 \ event_streams.f90 \ event_streams_uti.f90 \ event_streams_ut.f90 \ restricted_subprocesses.f90 \ eio_direct.f90 \ eio_direct_sub.f90 \ eio_direct_uti.f90 \ eio_direct_ut.f90 \ simulations.f90 \ restricted_subprocesses_uti.f90 \ restricted_subprocesses_ut.f90 \ simulations_uti.f90 \ simulations_ut.f90 \ commands.f90 \ commands_uti.f90 \ commands_ut.f90 \ cmdline_options.f90 \ libmanager.f90 \ features.f90 \ whizard.f90 \ api.f90 \ api_hepmc_uti.f90 \ api_hepmc_ut.f90 \ api_lcio_uti.f90 \ api_lcio_ut.f90 \ api_uti.f90 \ api_ut.f90 FC_OBJ = $(FC0_SRC:.f90=.o) $(F77_SRC:.f=.o) $(FC_SRC:.f90=.o) CC_OBJ = $(CC_SRC:.c=.o) all: whizard_test check: whizard_test ./whizard_test --check resonances whizard_test: $(FC_OBJ) $(CC_OBJ) main_ut.f90 $(FC) $(FC_OBJ) $(CC_OBJ) -ldl -o $@ main_ut.f90 whizard: $(FC_OBJ) $(CC_OBJ) main.f90 $(FC) $(FC_OBJ) $(CC_OBJ) -ldl -o $@ main.f90 %.o: %.f90 $(FC) $(FCFLAGS) -c $< %.o: %.f $(FC) $(FCFLAGS) -c $< %.o: %.c $(CC) $(CCFLAGS) -c $< tar: $(FC_SRC) $(F77_SRC) $(FC0_SRC) $(CC_SRC) $(MODELS) tar cvvzf whizard-`date +%y%m%d`-`date +%H%M`.tar.gz $(FC_SRC) $(FC0_SRC) \ $(F77_SRC) $(CC_SRC) main_ut.f90 Makefile $(MODELS) clean: rm -f *.mod *.o whizard_test