Index: trunk/src/openloops/openloops.nw =================================================================== --- trunk/src/openloops/openloops.nw (revision 8325) +++ trunk/src/openloops/openloops.nw (revision 8326) @@ -1,691 +1,738 @@ % -*- ess-noweb-default-code-mode: f90-mode; noweb-default-code-mode: f90-mode; -*- % WHIZARD code as NOWEB source: interface to OpenLoops 1-loop library %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{OpenLoops Interface} \includemodulegraph{openloops} The interface to OpenLoops. \clearpage %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% <<[[prc_openloops.f90]]>>= <> module prc_openloops use, intrinsic :: iso_c_binding !NODEP! use kinds use io_units <> use string_utils, only: str use constants use numeric_utils use diagnostics <> use system_dependencies use physics_defs use variables use os_interface use lorentz use interactions use sm_qcd use sm_physics, only: top_width_sm_lo, top_width_sm_qcd_nlo_jk use model_data use prclib_interfaces use prc_core_def use prc_core use blha_config use blha_olp_interfaces <> <> <> <> <> <> contains <> end module prc_openloops @ %def module prc_openloops @ <>= real(default), parameter :: openloops_default_bmass = 0._default real(default), parameter :: openloops_default_topmass = 172._default real(default), parameter :: openloops_default_topwidth = 0._default real(default), parameter :: openloops_default_wmass = 80.399_default real(default), parameter :: openloops_default_wwidth = 0._default real(default), parameter :: openloops_default_zmass = 91.1876_default real(default), parameter :: openloops_default_zwidth = 0._default real(default), parameter :: openloops_default_higgsmass = 125._default real(default), parameter :: openloops_default_higgswidth = 0._default integer :: N_EXTERNAL = 0 @ %def openloops default parameter @ <>= abstract interface subroutine ol_evaluate_scpowheg (id, pp, emitter, res, resmunu) bind(C) import integer(kind = c_int), value :: id, emitter real(kind = c_double), intent(in) :: pp(5 * N_EXTERNAL) real(kind = c_double), intent(out) :: res, resmunu(16) end subroutine ol_evaluate_scpowheg end interface @ %def ol_evaluate_scpowheg interface @ +<>= + abstract interface + subroutine ol_getparameter_double (variable_name, value) bind(C) + import + character(kind = c_char,len = 1), intent(in) :: variable_name + real(kind = c_double), intent(out) :: value + end subroutine ol_getparameter_double + end interface + +@ %def ol_getparameter_double interface +@ <>= type, extends (prc_blha_writer_t) :: openloops_writer_t contains <> end type openloops_writer_t @ %def openloops_writer_t @ <>= public :: openloops_def_t <>= type, extends (blha_def_t) :: openloops_def_t integer :: verbosity contains <> end type openloops_def_t @ %def openloops_def_t @ <>= type, extends (blha_driver_t) :: openloops_driver_t integer :: n_external = 0 type(string_t) :: olp_file procedure(ol_evaluate_scpowheg), nopass, pointer :: & evaluate_spin_correlations_powheg => null () + procedure(ol_getparameter_double), nopass, pointer :: & + get_parameter_double => null () contains <> end type openloops_driver_t @ %def openloops_driver_t @ <>= type :: openloops_threshold_data_t logical :: nlo = .true. real(default) :: alpha_ew real(default) :: sinthw real(default) :: m_b, m_W real(default) :: vtb contains <> end type openloops_threshold_data_t @ %def openloops_threshold_data_t @ <>= procedure :: compute_top_width => & openloops_threshold_data_compute_top_width <>= function openloops_threshold_data_compute_top_width & (data, mtop, alpha_s) result (wtop) real(default) :: wtop class(openloops_threshold_data_t), intent(in) :: data real(default), intent(in) :: mtop, alpha_s if (data%nlo) then wtop = top_width_sm_qcd_nlo_jk (data%alpha_ew, data%sinthw, & data%vtb, mtop, data%m_W, data%m_b, alpha_s) else wtop = top_width_sm_lo (data%alpha_ew, data%sinthw, data%vtb, & mtop, data%m_W, data%m_b) end if end function openloops_threshold_data_compute_top_width @ %def openloops_threshold_data_compute_top_width @ <>= public :: openloops_state_t <>= type, extends (blha_state_t) :: openloops_state_t type(openloops_threshold_data_t), allocatable :: threshold_data contains <> end type openloops_state_t @ %def openloops_state_t @ <>= procedure :: init_threshold => openloops_state_init_threshold <>= subroutine openloops_state_init_threshold (object, model) class(openloops_state_t), intent(inout) :: object type(model_data_t), intent(in) :: model if (model%get_name () == "SM_tt_threshold") then allocate (object%threshold_data) associate (data => object%threshold_data) data%nlo = btest (int (model%get_real (var_str ('offshell_strategy'))), 0) data%alpha_ew = one / model%get_real (var_str ('alpha_em_i')) data%sinthw = model%get_real (var_str ('sw')) data%m_b = model%get_real (var_str ('mb')) data%m_W = model%get_real (var_str ('mW')) data%vtb = model%get_real (var_str ('Vtb')) end associate end if end subroutine openloops_state_init_threshold @ %def openloops_state_init_threshold @ <>= public :: prc_openloops_t <>= type, extends (prc_blha_t) :: prc_openloops_t contains <> end type prc_openloops_t @ %def prc_openloops_t @ <>= procedure, nopass :: type_name => openloops_writer_type_name <>= function openloops_writer_type_name () result (string) type(string_t) :: string string = "openloops" end function openloops_writer_type_name @ @ %def openloops_writer_type_name <>= procedure :: init => openloops_def_init <>= subroutine openloops_def_init (object, basename, model_name, & prt_in, prt_out, nlo_type, restrictions, var_list) class(openloops_def_t), intent(inout) :: object type(string_t), intent(in) :: basename, model_name type(string_t), dimension(:), intent(in) :: prt_in, prt_out integer, intent(in) :: nlo_type type(string_t), intent(in), optional :: restrictions type(var_list_t), intent(in) :: var_list <> object%basename = basename allocate (openloops_writer_t :: object%writer) select case (nlo_type) case (BORN) object%suffix = '_BORN' case (NLO_REAL) object%suffix = '_REAL' case (NLO_VIRTUAL) object%suffix = '_LOOP' case (NLO_SUBTRACTION, NLO_MISMATCH) object%suffix = '_SUB' case (NLO_DGLAP) object%suffix = '_DGLAP' end select <> select type (writer => object%writer) class is (prc_blha_writer_t) call writer%init (model_name, prt_in, prt_out, restrictions) end select object%verbosity = var_list%get_ival (var_str ("openloops_verbosity")) end subroutine openloops_def_init @ %def openloops_def_init @ Add additional suffix for each rank of the communicator, such that the filenames do not clash. <>= integer :: n_size, rank <>= call MPI_comm_rank (MPI_COMM_WORLD, rank) call MPI_Comm_size (MPI_COMM_WORLD, n_size) if (n_size > 1) then object%suffix = object%suffix // var_str ("_") // str (rank) end if @ <>= procedure, nopass :: type_string => openloops_def_type_string <>= function openloops_def_type_string () result (string) type(string_t) :: string string = "openloops" end function openloops_def_type_string @ @ %def openloops_def_type_string <>= procedure :: write => openloops_def_write <>= subroutine openloops_def_write (object, unit) class(openloops_def_t), intent(in) :: object integer, intent(in) :: unit select type (writer => object%writer) type is (openloops_writer_t) call writer%write (unit) end select end subroutine openloops_def_write @ @ %def openloops_def_write <>= procedure :: init_dlaccess_to_library => openloops_driver_init_dlaccess_to_library <>= subroutine openloops_driver_init_dlaccess_to_library & (object, os_data, dlaccess, success) class(openloops_driver_t), intent(in) :: object type(os_data_t), intent(in) :: os_data type(dlaccess_t), intent(out) :: dlaccess logical, intent(out) :: success type(string_t) :: ol_library, msg_buffer ol_library = OPENLOOPS_DIR // '/lib/libopenloops.' // & os_data%shrlib_ext msg_buffer = "One-Loop-Provider: Using OpenLoops" call msg_message (char(msg_buffer)) msg_buffer = "Loading library: " // ol_library call msg_message (char(msg_buffer)) if (os_file_exist (ol_library)) then call dlaccess_init (dlaccess, var_str (""), ol_library, os_data) else call msg_fatal ("Link OpenLoops: library not found") end if success = .not. dlaccess_has_error (dlaccess) end subroutine openloops_driver_init_dlaccess_to_library @ %def openloops_driver_init_dlaccess_to_library @ <>= procedure :: set_alpha_s => openloops_driver_set_alpha_s <>= subroutine openloops_driver_set_alpha_s (driver, alpha_s) class(openloops_driver_t), intent(in) :: driver real(default), intent(in) :: alpha_s integer :: ierr if (associated (driver%blha_olp_set_parameter)) then call driver%blha_olp_set_parameter & (c_char_'alphas'//c_null_char, & dble (alpha_s), 0._double, ierr) else call msg_fatal ("blha_olp_set_parameter not associated!") end if - if (ierr == 0) call parameter_error_message (var_str ('alphas')) + if (ierr == 0) call parameter_error_message (var_str ('alphas'), & + var_str ('openloops_driver_set_alpha_s')) end subroutine openloops_driver_set_alpha_s @ %def openloops_driver_set_alpha_s @ <>= procedure :: set_alpha_qed => openloops_driver_set_alpha_qed <>= subroutine openloops_driver_set_alpha_qed (driver, alpha) class(openloops_driver_t), intent(inout) :: driver real(default), intent(in) :: alpha integer :: ierr call driver%blha_olp_set_parameter & (c_char_'alpha_qed'//c_null_char, & dble (alpha), 0._double, ierr) - if (ierr == 0) call parameter_error_message (var_str ('alpha_qed')) + if (ierr == 0) call ew_parameter_error_message (var_str ('alpha_qed')) end subroutine openloops_driver_set_alpha_qed @ %def openloops_driver_set_alpha_qed @ <>= procedure :: set_GF => openloops_driver_set_GF <>= subroutine openloops_driver_set_GF (driver, GF) class(openloops_driver_t), intent(inout) :: driver real(default), intent(in) :: GF integer :: ierr call driver%blha_olp_set_parameter & (c_char_'Gmu'//c_null_char, & dble(GF), 0._double, ierr) - if (ierr == 0) call parameter_error_message (var_str ('Gmu')) + if (ierr == 0) call ew_parameter_error_message (var_str ('Gmu')) end subroutine openloops_driver_set_GF @ %def openloops_driver_set_GF @ <>= procedure :: set_weinberg_angle => openloops_driver_set_weinberg_angle <>= subroutine openloops_driver_set_weinberg_angle (driver, sw2) class(openloops_driver_t), intent(inout) :: driver real(default), intent(in) :: sw2 integer :: ierr call driver%blha_olp_set_parameter & (c_char_'sw2'//c_null_char, & dble(sw2), 0._double, ierr) - if (ierr == 0) call parameter_error_message (var_str ('sw2')) + if (ierr == 0) call ew_parameter_error_message (var_str ('sw2')) end subroutine openloops_driver_set_weinberg_angle @ %def openloops_driver_set_weinberg_angle @ <>= procedure :: print_alpha_s => openloops_driver_print_alpha_s <>= subroutine openloops_driver_print_alpha_s (object) class(openloops_driver_t), intent(in) :: object call object%blha_olp_print_parameter (c_char_'alphas'//c_null_char) end subroutine openloops_driver_print_alpha_s @ %def openloops_driver_print_alpha_s @ <>= procedure, nopass :: type_name => openloops_driver_type_name <>= function openloops_driver_type_name () result (type) type(string_t) :: type type = "OpenLoops" end function openloops_driver_type_name @ %def openloops_driver_type_name @ <>= - procedure :: load_sc_procedure => openloops_driver_load_sc_procedure + procedure :: load_procedures => openloops_driver_load_procedures <>= - subroutine openloops_driver_load_sc_procedure (object, os_data, success) + subroutine openloops_driver_load_procedures (object, os_data, success) class(openloops_driver_t), intent(inout) :: object type(os_data_t), intent(in) :: os_data logical, intent(out) :: success type(dlaccess_t) :: dlaccess type(c_funptr) :: c_fptr logical :: init_success call object%init_dlaccess_to_library (os_data, dlaccess, init_success) c_fptr = dlaccess_get_c_funptr (dlaccess, var_str ("ol_evaluate_scpowheg")) call c_f_procpointer (c_fptr, object%evaluate_spin_correlations_powheg) - if (dlaccess_has_error (dlaccess)) then - call msg_fatal ("Could not load Openloops-powheg spin correlations!") - else - success = .true. - end if + call check_for_error (var_str ("ol_evaluate_scpowheg")) + c_fptr = dlaccess_get_c_funptr (dlaccess, var_str ("ol_getparameter_double")) + call c_f_procpointer (c_fptr, object%get_parameter_double) + call check_for_error (var_str ("ol_getparameter_double")) + success = .true. + + contains + subroutine check_for_error (function_name) + type(string_t), intent(in) :: function_name + if (dlaccess_has_error (dlaccess)) & + call msg_fatal (char ("Loading of " // function_name // " failed!")) + end subroutine check_for_error + end subroutine openloops_driver_load_procedures - end subroutine openloops_driver_load_sc_procedure - -@ %def openloops_driver_load_sc_procedure +@ %def openloops_driver_load_procedures @ <>= procedure :: read => openloops_def_read <>= subroutine openloops_def_read (object, unit) class(openloops_def_t), intent(out) :: object integer, intent(in) :: unit end subroutine openloops_def_read @ %def openloops_def_read @ <>= procedure :: allocate_driver => openloops_def_allocate_driver <>= subroutine openloops_def_allocate_driver (object, driver, basename) class(openloops_def_t), intent(in) :: object class(prc_core_driver_t), intent(out), allocatable :: driver type(string_t), intent(in) :: basename if (.not. allocated (driver)) allocate (openloops_driver_t :: driver) end subroutine openloops_def_allocate_driver @ @ %def openloops_def_allocate_driver <>= procedure :: write => openloops_state_write <>= subroutine openloops_state_write (object, unit) class(openloops_state_t), intent(in) :: object integer, intent(in), optional :: unit end subroutine openloops_state_write @ %def prc_openloops_state_write @ <>= procedure :: allocate_workspace => prc_openloops_allocate_workspace <>= subroutine prc_openloops_allocate_workspace (object, core_state) class(prc_openloops_t), intent(in) :: object class(prc_core_state_t), intent(inout), allocatable :: core_state allocate (openloops_state_t :: core_state) end subroutine prc_openloops_allocate_workspace @ %def prc_openloops_allocate_workspace @ <>= procedure :: init_driver => prc_openloops_init_driver <>= subroutine prc_openloops_init_driver (object, os_data) class(prc_openloops_t), intent(inout) :: object type(os_data_t), intent(in) :: os_data type(string_t) :: olp_file, olc_file type(string_t) :: suffix select type (def => object%def) type is (openloops_def_t) suffix = def%suffix olp_file = def%basename // suffix // '.olp' olc_file = def%basename // suffix // '.olc' class default call msg_bug ("prc_openloops_init_driver: core_def should be openloops-type") end select select type (driver => object%driver) type is (openloops_driver_t) driver%olp_file = olp_file driver%contract_file = olc_file driver%nlo_suffix = suffix end select end subroutine prc_openloops_init_driver @ %def prc_openloops_init_driver @ <>= procedure :: write => prc_openloops_write <>= subroutine prc_openloops_write (object, unit) class(prc_openloops_t), intent(in) :: object integer, intent(in), optional :: unit call msg_message (unit = unit, string = "OpenLoops") end subroutine prc_openloops_write @ @ %def prc_openloops_write <>= procedure :: write_name => prc_openloops_write_name <>= subroutine prc_openloops_write_name (object, unit) class(prc_openloops_t), intent(in) :: object integer, intent(in), optional :: unit integer :: u u = given_output_unit (unit) write (u,"(1x,A)") "Core: OpenLoops" end subroutine prc_openloops_write_name @ @ %def prc_openloops_write_name <>= procedure :: prepare_library => prc_openloops_prepare_library <>= subroutine prc_openloops_prepare_library (object, os_data, model) class(prc_openloops_t), intent(inout) :: object type(os_data_t), intent(in) :: os_data type(model_data_t), intent(in), target :: model call object%load_driver (os_data) call object%reset_parameters () call object%set_particle_properties (model) - call object%set_electroweak_parameters (model) select type(def => object%def) type is (openloops_def_t) call object%set_verbosity (def%verbosity) end select end subroutine prc_openloops_prepare_library @ %def prc_openloops_prepare_library @ <>= procedure :: load_driver => prc_openloops_load_driver <>= subroutine prc_openloops_load_driver (object, os_data) class(prc_openloops_t), intent(inout) :: object type(os_data_t), intent(in) :: os_data logical :: success select type (driver => object%driver) type is (openloops_driver_t) call driver%load (os_data, success) - call driver%load_sc_procedure (os_data, success) + call driver%load_procedures (os_data, success) end select end subroutine prc_openloops_load_driver @ %def prc_openloops_load_driver @ <>= procedure :: start => prc_openloops_start <>= subroutine prc_openloops_start (object) class(prc_openloops_t), intent(inout) :: object integer :: ierr select type (driver => object%driver) type is (openloops_driver_t) call driver%blha_olp_start (char (driver%olp_file)//c_null_char, ierr) end select end subroutine prc_openloops_start @ %def prc_openloops_start @ <>= procedure :: set_n_external => prc_openloops_set_n_external <>= subroutine prc_openloops_set_n_external (object, n) class(prc_openloops_t), intent(inout) :: object integer, intent(in) :: n N_EXTERNAL = n end subroutine prc_openloops_set_n_external @ %def prc_openloops_set_n_external @ <>= procedure :: reset_parameters => prc_openloops_reset_parameters <>= subroutine prc_openloops_reset_parameters (object) class(prc_openloops_t), intent(inout) :: object integer :: ierr select type (driver => object%driver) type is (openloops_driver_t) call driver%blha_olp_set_parameter ('mass(5)'//c_null_char, & dble(openloops_default_bmass), 0._double, ierr) call driver%blha_olp_set_parameter ('mass(6)'//c_null_char, & dble(openloops_default_topmass), 0._double, ierr) call driver%blha_olp_set_parameter ('width(6)'//c_null_char, & dble(openloops_default_topwidth), 0._double, ierr) call driver%blha_olp_set_parameter ('mass(23)'//c_null_char, & dble(openloops_default_zmass), 0._double, ierr) call driver%blha_olp_set_parameter ('width(23)'//c_null_char, & dble(openloops_default_zwidth), 0._double, ierr) call driver%blha_olp_set_parameter ('mass(24)'//c_null_char, & dble(openloops_default_wmass), 0._double, ierr) call driver%blha_olp_set_parameter ('width(24)'//c_null_char, & dble(openloops_default_wwidth), 0._double, ierr) call driver%blha_olp_set_parameter ('mass(25)'//c_null_char, & dble(openloops_default_higgsmass), 0._double, ierr) call driver%blha_olp_set_parameter ('width(25)'//c_null_char, & dble(openloops_default_higgswidth), 0._double, ierr) end select end subroutine prc_openloops_reset_parameters @ %def prc_openloops_reset_parameters @ Set the verbosity level for openloops. The different levels are as follows: \begin{itemize} \item[0] minimal output (startup message et.al.) \item[1] show which libraries are loaded \item[2] show debug information of the library loader, but not during run time \item[3] show debug information during run time \item[4] output for each call of [[set_parameters]]. \end{itemize} <>= procedure :: set_verbosity => prc_openloops_set_verbosity <>= subroutine prc_openloops_set_verbosity (object, verbose) class(prc_openloops_t), intent(inout) :: object integer, intent(in) :: verbose integer :: ierr select type (driver => object%driver) type is (openloops_driver_t) call driver%blha_olp_set_parameter ('verbose'//c_null_char, & dble(verbose), 0._double, ierr) end select end subroutine prc_openloops_set_verbosity @ %def prc_openloops_set_verbosity @ <>= procedure :: prepare_external_code => & prc_openloops_prepare_external_code <>= subroutine prc_openloops_prepare_external_code & (core, flv_states, var_list, os_data, libname, model, i_core, is_nlo) class(prc_openloops_t), intent(inout) :: core integer, intent(in), dimension(:,:), allocatable :: flv_states type(var_list_t), intent(in) :: var_list type(os_data_t), intent(in) :: os_data type(string_t), intent(in) :: libname type(model_data_t), intent(in), target :: model integer, intent(in) :: i_core logical, intent(in) :: is_nlo + integer :: ierr core%sqme_tree_pos = 1 call core%set_n_external (core%data%get_n_tot ()) call core%prepare_library (os_data, model) call core%start () + call core%set_electroweak_parameters (model) + select type (driver => core%driver) + type is (openloops_driver_t) + !!! We have to set the external vector boson wavefunction to the MadGraph convention + !!! in order to be consistent with our calculation of the spin correlated contributions. + call driver%blha_olp_set_parameter ('wf_v_select'//c_null_char, & + 3._double, 0._double, ierr) + if (ierr == 0) call parameter_error_message (var_str ('wf_v_select'), & + var_str ('prc_openloops_prepare_external_code')) + end select call core%read_contract_file (flv_states) call core%print_parameter_file (i_core) end subroutine prc_openloops_prepare_external_code @ %def prc_openloops_prepare_external_code @ Computes a spin-correlated matrix element from an interface to an external one-loop provider. The output of [[blha_olp_eval2]] is an array of [[dimension(16)]]. The current interface does not give out an accuracy, so that [[bad_point]] is always [[.false.]]. OpenLoops includes a factor of 1 / [[n_hel]] in the amplitudes, which we have to undo if polarized matrix elements are requested. <>= procedure :: compute_sqme_spin_c => prc_openloops_compute_sqme_spin_c <>= subroutine prc_openloops_compute_sqme_spin_c (object, & i_flv, i_hel, em, p, ren_scale, sqme_spin_c, bad_point) class(prc_openloops_t), intent(inout) :: object integer, intent(in) :: i_flv, i_hel integer, intent(in) :: em type(vector4_t), intent(in), dimension(:) :: p real(default), intent(in) :: ren_scale real(default), intent(out), dimension(16) :: sqme_spin_c logical, intent(out) :: bad_point real(double), dimension(5*N_EXTERNAL) :: mom real(double) :: res real(double), dimension(16) :: res_munu real(default) :: alpha_s if (object%i_spin_c(i_flv, i_hel) >= 0) then mom = object%create_momentum_array (p) if (vanishes (ren_scale)) call msg_fatal & ("prc_openloops_compute_sqme_spin_c: ren_scale vanishes") alpha_s = object%qcd%alpha%get (ren_scale) select type (driver => object%driver) type is (openloops_driver_t) call driver%set_alpha_s (alpha_s) call driver%evaluate_spin_correlations_powheg & (object%i_spin_c(i_flv, i_hel), mom, em, res, res_munu) end select sqme_spin_c = res_munu bad_point = .false. if (object%includes_polarization ()) & sqme_spin_c = object%n_hel * sqme_spin_c if (debug_on) then if (sum(sqme_spin_c) == 0) then call msg_debug(D_SUBTRACTION,'Spin-correlated matrix elements provided by OpenLoops are zero!') end if end if else sqme_spin_c = zero end if end subroutine prc_openloops_compute_sqme_spin_c @ %def prc_openloops_compute_sqme_spin_c @ +<>= + procedure :: get_alpha_qed => prc_openloops_get_alpha_qed +<>= + function prc_openloops_get_alpha_qed (object) result (alpha_qed) + class(prc_openloops_t), intent(in) :: object + real(double) :: value + real(default) :: alpha_qed + select type (driver => object%driver) + type is (openloops_driver_t) + call driver%get_parameter_double ('alpha_qed'//c_null_char, value) + alpha_qed = value + return + end select + call msg_fatal ("prc_openloops_get_alpha_qed: called by wrong driver, only supported for OpenLoops!") + end function prc_openloops_get_alpha_qed + +@ %def prc_openloops_get_alpha_qed Index: trunk/src/me_methods/me_methods.nw =================================================================== --- trunk/src/me_methods/me_methods.nw (revision 8325) +++ trunk/src/me_methods/me_methods.nw (revision 8326) @@ -1,6924 +1,6926 @@ % -*- ess-noweb-default-code-mode: f90-mode; noweb-default-code-mode: f90-mode; -*- % WHIZARD code as NOWEB source: matrix elements and process libraries %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{Interface for Matrix Element Objects} \includemodulegraph{me_methods} These modules manage internal and, in particular, external matrix-element code. \begin{description} \item[prc\_core] We define the abstract [[prc_core_t]] type which handles all specific features of kinematics matrix-element evaluation that depend on a particular class of processes. This abstract type supplements the [[prc_core_def_t]] type and related types in another module. Together, they provide a complete set of matrix-element handlers that are implemented in the concrete types below. \end{description} These are the implementations: \begin{description} \item[prc\_template\_me] Implements matrix-element code without actual content (trivial value), but full-fledged interface. This can be used for injecting user-defined matrix-element code. \item[prc\_omega] Matrix elements calculated by \oMega\ are the default for WHIZARD. Here, we provide all necessary support. \item[prc\_external] Matrix elements provided or using external (not \oMega) code or libraries. This is an abstract type, with concrete extensions below: \item[prc\_external\_test] Concrete implementation of the external-code type, actually using some pre-defined test matrix elements. \item[prc\_gosam] Interface for matrix elements computed using \gosam. \item[prc\_openloops] Interface for matrix elements computed using OpenLoops. \item[prc\_recola] Interface for matrix elements computed using Recola. \item[prc\_threshold] Interface for matrix elements for the top-pair threshold, that use external libraries. \end{description} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Abstract process core} In this module we provide abstract data types for process classes. Each process class represents a set of processes which are handled by a common ``method'', e.g., by the \oMega\ matrix-element generator. The process class is also able to select a particular implementation for the phase-space and integration modules. For a complete implementation of a process class, we have to provide extensions of the following abstract types: \begin{description} \item[prc\_core\_def\_t] process and matrix-element configuration \item[prc\_writer\_t] (optional) writing external matrix-element code \item[prc\_driver\_t] accessing the matrix element (internal or external) \item[prc\_core\_t] evaluating kinematics and matrix element. The process core also selects phase-space and integrator implementations as appropriate for the process class and configuration. \end{description} In the actual process-handling data structures, each process component contains an instance of such a process class as its core. This allows us to keep the [[processes]] module below, which supervises matrix-element evaluation, integration, and event generation, free of any reference to concrete implementations (for the process class, phase space, and integrator). There are no unit tests, these are deferred to the [[processes]] module. <<[[prc_core.f90]]>>= <> module prc_core <> <> use io_units use diagnostics use os_interface, only: os_data_t use lorentz use interactions use variables, only: var_list_t use model_data, only: model_data_t use process_constants use prc_core_def use process_libraries use sf_base <> <> <> <> contains <> end module prc_core @ %def prc_core @ \subsection{The process core} The process core is of abstract data type. Different types of matrix element will be represented by different implementations. <>= public :: prc_core_t <>= type, abstract :: prc_core_t class(prc_core_def_t), pointer :: def => null () logical :: data_known = .false. type(process_constants_t) :: data class(prc_core_driver_t), allocatable :: driver logical :: use_color_factors = .false. integer :: nc = 3 contains <> end type prc_core_t @ %def prc_core_t @ In any case there must be an output routine. <>= procedure(prc_core_write), deferred :: write <>= abstract interface subroutine prc_core_write (object, unit) import class(prc_core_t), intent(in) :: object integer, intent(in), optional :: unit end subroutine prc_core_write end interface @ %def prc_core_write @ Just type the name of the actual core method. <>= procedure(prc_core_write_name), deferred :: write_name <>= abstract interface subroutine prc_core_write_name (object, unit) import class(prc_core_t), intent(in) :: object integer, intent(in), optional :: unit end subroutine prc_core_write_name end interface @ %def prc_core_write_name @ For initialization, we assign a pointer to the process entry in the relevant library. This allows us to access all process functions via the implementation of [[prc_core_t]]. We declare the [[object]] as [[intent(inout)]], since just after allocation it may be useful to store some extra data in the object, which we can then use in the actual initialization. This applies to extensions of [[prc_core]] which override the [[init]] method. <>= procedure :: init => prc_core_init procedure :: base_init => prc_core_init <>= subroutine prc_core_init (object, def, lib, id, i_component) class(prc_core_t), intent(inout) :: object class(prc_core_def_t), intent(in), target :: def type(process_library_t), intent(in), target :: lib type(string_t), intent(in) :: id integer, intent(in) :: i_component object%def => def call lib%connect_process (id, i_component, object%data, object%driver) object%data_known = .true. end subroutine prc_core_init @ %def prc_core_init @ Return true if the matrix element generation was successful. This can be tested by looking at the number of generated flavor states, which should be nonzero. <>= procedure :: has_matrix_element => prc_core_has_matrix_element <>= function prc_core_has_matrix_element (object) result (flag) class(prc_core_t), intent(in) :: object logical :: flag flag = object%data%n_flv /= 0 end function prc_core_has_matrix_element @ %def prc_core_has_matrix_element @ Return true if this process-core type needs extra code that has to be compiled and/or linked, beyond the default \oMega\ framework. This depends only on the concrete type, and the default is no. <>= procedure, nopass :: needs_external_code => prc_core_needs_external_code <>= function prc_core_needs_external_code () result (flag) logical :: flag flag = .false. end function prc_core_needs_external_code @ %def prc_core_needs_external_code @ The corresponding procedure to create and load extra libraries. The base procedure must not be called but has to be overridden, if extra code is required. <>= procedure :: prepare_external_code => & prc_core_prepare_external_code <>= subroutine prc_core_prepare_external_code & (core, flv_states, var_list, os_data, libname, model, i_core, is_nlo) class(prc_core_t), intent(inout) :: core integer, intent(in), dimension(:,:), allocatable :: flv_states type(var_list_t), intent(in) :: var_list type(os_data_t), intent(in) :: os_data type(string_t), intent(in) :: libname type(model_data_t), intent(in), target :: model integer, intent(in) :: i_core logical, intent(in) :: is_nlo call core%write () call msg_bug ("prc_core_prepare_external_code called & &but not overridden") end subroutine prc_core_prepare_external_code @ %def prc_core_prepare_external_code @ Return true if this process-core type uses the BLHA interface. This depends only on the concrete type, and the default is no. <>= procedure, nopass :: uses_blha => prc_core_uses_blha <>= function prc_core_uses_blha () result (flag) logical :: flag flag = .false. end function prc_core_uses_blha @ %def prc_core_uses_blha @ Tell whether a particular combination of flavor/helicity/color state is allowed for the matrix element. <>= procedure(prc_core_is_allowed), deferred :: is_allowed <>= abstract interface function prc_core_is_allowed (object, i_term, f, h, c) result (flag) import class(prc_core_t), intent(in) :: object integer, intent(in) :: i_term, f, h, c logical :: flag end function prc_core_is_allowed end interface @ %def prc_core_is_allowed @ Set the constant process data for a specific term. By default, these are the constants stored inside the object, ignoring the term index. Type extensions may override this and provide term-specific data. <>= procedure :: get_constants => prc_core_get_constants <>= subroutine prc_core_get_constants (object, data, i_term) class(prc_core_t), intent(in) :: object type(process_constants_t), intent(out) :: data integer, intent(in) :: i_term data = object%data end subroutine prc_core_get_constants @ %def prc_core_get_constants @ The strong coupling is not among the process constants. The default implementation is to return a negative number, which indicates that $\alpha_s$ is not available. This may be overridden by an implementation that provides an (event-specific) value. The value can be stored in the process-specific workspace. <>= procedure :: get_alpha_s => prc_core_get_alpha_s <>= function prc_core_get_alpha_s (object, core_state) result (alpha) class(prc_core_t), intent(in) :: object class(prc_core_state_t), intent(in), allocatable :: core_state real(default) :: alpha alpha = -1 end function prc_core_get_alpha_s @ %def prc_core_get_alpha_s @ Allocate the workspace associated to a process component. The default is that there is no workspace, so we do nothing. A type extension may override this and allocate a workspace object of appropriate type, which can be used in further calculations. In any case, the [[intent(out)]] attribute deletes any previously allocated workspace. <>= procedure :: allocate_workspace => prc_core_ignore_workspace <>= subroutine prc_core_ignore_workspace (object, core_state) class(prc_core_t), intent(in) :: object class(prc_core_state_t), intent(inout), allocatable :: core_state end subroutine prc_core_ignore_workspace @ %def prc_core_ignore_workspace @ Compute the momenta in the hard interaction, taking the seed kinematics as input. The [[i_term]] index tells us which term we want to compute. (The standard method is to just transfer the momenta to the hard interaction.) <>= procedure(prc_core_compute_hard_kinematics), deferred :: & compute_hard_kinematics <>= abstract interface subroutine prc_core_compute_hard_kinematics & (object, p_seed, i_term, int_hard, core_state) import class(prc_core_t), intent(in) :: object type(vector4_t), dimension(:), intent(in) :: p_seed integer, intent(in) :: i_term type(interaction_t), intent(inout) :: int_hard class(prc_core_state_t), intent(inout), allocatable :: core_state end subroutine prc_core_compute_hard_kinematics end interface @ %def prc_core_compute_hard_kinematics @ Compute the momenta in the effective interaction, taking the hard kinematics as input. (This is called only if parton recombination is to be applied for the process variant.) <>= procedure(prc_core_compute_eff_kinematics), deferred :: & compute_eff_kinematics <>= abstract interface subroutine prc_core_compute_eff_kinematics & (object, i_term, int_hard, int_eff, core_state) import class(prc_core_t), intent(in) :: object integer, intent(in) :: i_term type(interaction_t), intent(in) :: int_hard type(interaction_t), intent(inout) :: int_eff class(prc_core_state_t), intent(inout), allocatable :: core_state end subroutine prc_core_compute_eff_kinematics end interface @ %def prc_core_compute_eff_kinematics @ The process core must implement this function. Here, [[j]] is the index of the particular term we want to compute. The amplitude may depend on the factorization and renormalization scales. The [[core_state]] (workspace) argument may be used if it is provided by the caller. Otherwise, the routine should compute the result directly. <>= procedure(prc_core_compute_amplitude), deferred :: compute_amplitude <>= abstract interface function prc_core_compute_amplitude & (object, j, p, f, h, c, fac_scale, ren_scale, alpha_qcd_forced, & core_state) result (amp) import complex(default) :: amp class(prc_core_t), intent(in) :: object integer, intent(in) :: j type(vector4_t), dimension(:), intent(in) :: p integer, intent(in) :: f, h, c real(default), intent(in) :: fac_scale, ren_scale real(default), intent(in), allocatable :: alpha_qcd_forced class(prc_core_state_t), intent(inout), allocatable, optional :: & core_state end function prc_core_compute_amplitude end interface @ %def prc_core_compute_amplitude @ \subsection{Storage for intermediate results} The abstract [[prc_core_state_t]] type allows process cores to set up temporary workspace. The object is an extra argument for each of the individual calculations between kinematics setup and matrix-element evaluation. <>= public :: prc_core_state_t <>= type, abstract :: prc_core_state_t contains procedure(workspace_write), deferred :: write procedure(workspace_reset_new_kinematics), deferred :: reset_new_kinematics end type prc_core_state_t @ %def prc_core_state_t @ For debugging, we should at least have an output routine. <>= abstract interface subroutine workspace_write (object, unit) import class(prc_core_state_t), intent(in) :: object integer, intent(in), optional :: unit end subroutine workspace_write end interface @ %def workspace_write @ This is used during the NLO calculation, see there for more information. <>= abstract interface subroutine workspace_reset_new_kinematics (object) import class(prc_core_state_t), intent(inout) :: object end subroutine workspace_reset_new_kinematics end interface @ %def workspace_reset_new_kinematics @ \subsection{Helicity selection data} This is intended for use with \oMega, but may also be made available to other process methods. We set thresholds for counting the times a specific helicity amplitude is zero. When the threshold is reached, we skip this amplitude in subsequent calls. For initializing the helicity counters, we need an object that holds the two parameters, the threshold (large real number) and the cutoff (integer). A helicity value suppressed by more than [[threshold]] (a value which multiplies [[epsilon]], to be compared with the average of the current amplitude, default is $10^{10}$) is treated as zero. A matrix element is assumed to be zero and not called again if it has been zero [[cutoff]] times. <>= public :: helicity_selection_t <>= type :: helicity_selection_t logical :: active = .false. real(default) :: threshold = 0 integer :: cutoff = 0 contains <> end type helicity_selection_t @ %def helicity_selection_t @ Output. If the selection is inactive, print nothing. <>= procedure :: write => helicity_selection_write <>= subroutine helicity_selection_write (object, unit) class(helicity_selection_t), intent(in) :: object integer, intent(in), optional :: unit integer :: u u = given_output_unit (unit) if (object%active) then write (u, "(3x,A)") "Helicity selection data:" write (u, "(5x,A,ES17.10)") & "threshold =", object%threshold write (u, "(5x,A,I0)") & "cutoff = ", object%cutoff end if end subroutine helicity_selection_write @ %def helicity_selection_write @ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Test process type} For the following tests, we define a simple implementation of the abstract [[prc_core_t]], designed such as to complement the [[prc_test_t]] process definition type. Note that it is not given that the actual process is defined as [[prc_test_t]] type. We enforce this by calling [[prc_test_create_library]]. The driver component in the process core will then become of type [[prc_test_t]]. @ <<[[prc_test_core.f90]]>>= <> module prc_test_core <> use io_units use lorentz use interactions use prc_test use prc_core <> <> <> contains <> end module prc_test_core @ %def prc_test_core <>= public :: test_t <>= type, extends (prc_core_t) :: test_t contains <> end type test_t @ %def test_t <>= procedure :: write => test_write <>= subroutine test_write (object, unit) class(test_t), intent(in) :: object integer, intent(in), optional :: unit integer :: u u = given_output_unit (unit) write (u, "(3x,A)") "test type implementing prc_test" end subroutine test_write @ %def test_write <>= procedure :: write_name => test_write_name <>= subroutine test_write_name (object, unit) class(test_t), intent(in) :: object integer, intent(in), optional :: unit integer :: u u = given_output_unit (unit) write (u,"(1x,A)") "Core: prc_test" end subroutine test_write_name @ %def test_write_name @ This process type always needs a MC parameter set and a single term. This only state is always allowed. <>= procedure :: needs_mcset => test_needs_mcset procedure :: get_n_terms => test_get_n_terms procedure :: is_allowed => test_is_allowed <>= function test_needs_mcset (object) result (flag) class(test_t), intent(in) :: object logical :: flag flag = .true. end function test_needs_mcset function test_get_n_terms (object) result (n) class(test_t), intent(in) :: object integer :: n n = 1 end function test_get_n_terms function test_is_allowed (object, i_term, f, h, c) result (flag) class(test_t), intent(in) :: object integer, intent(in) :: i_term, f, h, c logical :: flag flag = .true. end function test_is_allowed @ %def test_needs_mcset @ %def test_get_n_terms @ %def test_is_allowed @ Transfer the generated momenta directly to the hard interaction in the (only) term. We assume that everything has been set up correctly, so the array fits. <>= procedure :: compute_hard_kinematics => test_compute_hard_kinematics <>= subroutine test_compute_hard_kinematics & (object, p_seed, i_term, int_hard, core_state) class(test_t), intent(in) :: object type(vector4_t), dimension(:), intent(in) :: p_seed integer, intent(in) :: i_term type(interaction_t), intent(inout) :: int_hard class(prc_core_state_t), intent(inout), allocatable :: core_state call int_hard%set_momenta (p_seed) end subroutine test_compute_hard_kinematics @ %def test_compute_hard_kinematics @ This procedure is not called for [[test_t]], just a placeholder. <>= procedure :: compute_eff_kinematics => test_compute_eff_kinematics <>= subroutine test_compute_eff_kinematics & (object, i_term, int_hard, int_eff, core_state) class(test_t), intent(in) :: object integer, intent(in) :: i_term type(interaction_t), intent(in) :: int_hard type(interaction_t), intent(inout) :: int_eff class(prc_core_state_t), intent(inout), allocatable :: core_state end subroutine test_compute_eff_kinematics @ %def test_compute_eff_kinematics @ Transfer the incoming momenta of [[p_seed]] directly to the effective interaction, and vice versa for the outgoing momenta. [[int_hard]] is left untouched since [[int_eff]] is an alias (via pointer) to it. <>= procedure :: recover_kinematics => test_recover_kinematics <>= subroutine test_recover_kinematics & (object, p_seed, int_hard, int_eff, core_state) class(test_t), intent(in) :: object type(vector4_t), dimension(:), intent(inout) :: p_seed type(interaction_t), intent(inout) :: int_hard type(interaction_t), intent(inout) :: int_eff class(prc_core_state_t), intent(inout), allocatable :: core_state integer :: n_in n_in = int_eff%get_n_in () call int_eff%set_momenta (p_seed(1:n_in), outgoing = .false.) p_seed(n_in+1:) = int_eff%get_momenta (outgoing = .true.) end subroutine test_recover_kinematics @ %def test_recover_kinematics @ Compute the amplitude. The driver ignores all quantum numbers and, in fact, returns a constant. Nevertheless, we properly transfer the momentum vectors. <>= procedure :: compute_amplitude => test_compute_amplitude <>= function test_compute_amplitude & (object, j, p, f, h, c, fac_scale, ren_scale, alpha_qcd_forced, core_state) & result (amp) class(test_t), intent(in) :: object integer, intent(in) :: j type(vector4_t), dimension(:), intent(in) :: p integer, intent(in) :: f, h, c real(default), intent(in) :: fac_scale, ren_scale real(default), intent(in), allocatable :: alpha_qcd_forced class(prc_core_state_t), intent(inout), allocatable, optional :: core_state complex(default) :: amp real(default), dimension(:,:), allocatable :: parray integer :: i, n_tot select type (driver => object%driver) type is (prc_test_t) if (driver%scattering) then n_tot = 4 else n_tot = 3 end if allocate (parray (0:3,n_tot)) forall (i = 1:n_tot) parray(:,i) = vector4_get_components (p(i)) amp = driver%get_amplitude (parray) end select end function test_compute_amplitude @ %def test_compute_amplitude @ @ \section{Template matrix elements} Here, we provide template matrix elements that are in structure very similar to \oMega\ matrix elements, but do not need its infrastructure. Per default, the matrix elements are flat, i.e. they have the constant value one. Analogous to the \oMega\ implementation, this section implements the interface to the templates (via the makefile) and the driver. <<[[prc_template_me.f90]]>>= <> module prc_template_me use, intrinsic :: iso_c_binding !NODEP! use kinds <> use io_units use system_defs, only: TAB use diagnostics use os_interface use lorentz use flavors use interactions use model_data use particle_specifiers, only: new_prt_spec use process_constants use prclib_interfaces use prc_core_def use process_libraries use prc_core <> <