Index: trunk/src/recola/recola.nw =================================================================== --- trunk/src/recola/recola.nw (revision 8251) +++ trunk/src/recola/recola.nw (revision 8252) @@ -1,3280 +1,3281 @@ % -*- ess-noweb-default-code-mode: f90-mode; noweb-default-code-mode: f90-mode; -*- % WHIZARD code as NOWEB source: interface to Recola 1-loop library @ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{Recola Interface} \section{Recola wrapper} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% <<[[recola_wrapper.f90]]>>= <> module recola_wrapper use recola !NODEP! use kinds <> <> use constants, only: zero use diagnostics, only: msg_fatal, msg_message, msg_debug, msg_debug2, D_ME_METHODS use io_units, only: given_output_unit <> <> <> <> <> contains <> end module recola_wrapper @ %def recola_wrapper @ <>= public :: rclwrap_is_active <>= logical, parameter :: rclwrap_is_active = .true. @ %def rclwrap_is_active @ Returns the particle string corresponding to a pdg code used in the Recola process definition <>= public :: get_recola_particle_string <>= elemental function get_recola_particle_string (pdg) result (name) type(string_t) :: name integer, intent(in) :: pdg select case (pdg) case (1) name = var_str ("d") case (-1) name = var_str ("d~") case (2) name = var_str ("u") case (-2) name = var_str ("u~") case (3) name = var_str ("s") case (-3) name = var_str ("s~") case (4) name = var_str ("c") case (-4) name = var_str ("c~") case (5) name = var_str ("b") case (-5) name = var_str ("b~") case (6) name = var_str ("t") case (-6) name = var_str ("t~") case (11) name = var_str ("e-") case (-11) name = var_str ("e+") case (12) name = var_str ("nu_e") case (-12) name = var_str ("nu_e~") case (13) name = var_str ("mu-") case (-13) name = var_str ("mu+") case (14) name = var_str ("nu_mu") case (-14) name = var_str ("nu_mu~") case (15) name = var_str ("tau-") case (-15) name = var_str ("tau+") case (16) name = var_str ("nu_tau") case (-16) name = var_str ("nu_tau~") case (21) name = var_str ("g") case (22) name = var_str ("A") case (23) name = var_str ("Z") case (24) name = var_str ("W+") case (-24) name = var_str ("W-") case (25) name = var_str ("H") end select end function get_recola_particle_string @ %def get_recola_particle_string @ <>= subroutine rclwrap_define_process (id, process_string, order) integer, intent(in) :: id type(string_t), intent(in) :: process_string type(string_t), intent(in) :: order if (debug_on) call msg_debug2 (D_ME_METHODS, "define_process_rcl") call define_process_rcl (id, char (process_string), char (order)) end subroutine rclwrap_define_process @ %def rclwrap_define_process @ This defines a wrapper for the information required to define a RECOLA process. It is used to collect the process definitions in an array. <>= type :: rcl_process_t private integer :: id type(string_t) :: process_string type(string_t) :: order contains <> end type rcl_process_t @ %def rcl_process_t @ <>= interface rcl_process_t module procedure new_rcl_process_t end interface @ %def rcl_process_t @ <>= function new_rcl_process_t (id, process_string, order) integer, intent(in) :: id type(string_t), intent(in) :: process_string, order type(rcl_process_t) :: new_rcl_process_t new_rcl_process_t%id = id new_rcl_process_t%process_string = process_string new_rcl_process_t%order = order end function new_rcl_process_t @ %def new_rcl_process_t <>= procedure :: get_params => rcl_process_get_params <>= subroutine rcl_process_get_params (prc, id, process_string, order) class(rcl_process_t), intent(in) :: prc integer, intent(out) :: id type(string_t), intent(out) :: process_string type(string_t), intent(out) :: order id = prc%id process_string = prc%process_string order = prc%order end subroutine rcl_process_get_params @ %def rcl_process_get_params @ Output. <>= procedure :: write => rcl_process_write <>= subroutine rcl_process_write (object, unit) class(rcl_process_t), intent(in) :: object integer, intent(in), optional :: unit integer :: u u = given_output_unit (unit) write (u, "(1x,A,I0,2(1x,A,1x))") "RECOLA process:", & "id=", object%id, "process_string=", char(object%process_string), & "order=", char(object%order) end subroutine rcl_process_write @ %def rcl_process_write @ This defines a singleton object, located in this module only, that controls RECOLA initialization and process management. When WHIZARD compiles processes, it should also run the RECOLA "`controller"', which actually initializes RECOLA for integration and manages process information in an array. The main complication is that this has to be done after all processes have been registered, and cannot be redone. We could work with module variables directly, but the singleton pattern, e.g., allows us to work with multiple RECOLA instances, if this becomes possible in the future. Type and object can be private. <>= type :: rcl_controller_t private logical :: active = .false. logical :: defined = .false. logical :: done = .false. integer :: recola_id = 0 type(rcl_process_t), dimension (:), allocatable :: processes integer :: n_processes = 0 contains <> end type rcl_controller_t @ %def rcl_controller_t <>= type(rcl_controller_t), target, save :: rcl_controller @ %def rcl_controller @ Add a RECOLA process to the controller. This will make sure that processes can be redefined if additional definitions are to be made after process generation. <>= procedure :: add_process => rcl_controller_add_process <>= subroutine rcl_controller_add_process (rcl, process) class(rcl_controller_t), intent(inout) :: rcl type(rcl_process_t), intent(in) :: process type(rcl_process_t), dimension (:), allocatable :: temp if (rcl%n_processes == size(rcl%processes)) then allocate( temp(2 * rcl%n_processes) ) temp(:rcl%n_processes) = rcl%processes call move_alloc(temp, rcl%processes) end if rcl%processes(rcl%n_processes + 1) = process rcl%n_processes = rcl%n_processes + 1 end subroutine rcl_controller_add_process @ %def rcl_controller_add_process @ Define all processes added to the controller, and only them. If processes have been defined before, RECOLA is reset. <>= procedure :: define_processes => rcl_controller_define_processes <>= subroutine rcl_controller_define_processes (rcl) class(rcl_controller_t), intent(inout) :: rcl integer :: id, i type(string_t) :: process_string type(string_t) :: order if (rcl%defined) then if (.not. rcl%done) call rclwrap_generate_processes () if (debug_on) call msg_debug2 (D_ME_METHODS, "reset_recola_rcl") call reset_recola_rcl () end if do i = 1, rcl%n_processes call rcl%processes(i)%get_params(id, process_string, order) call rclwrap_define_process (id, process_string, order) end do rcl%defined = .true. rcl%done = .false. end subroutine rcl_controller_define_processes @ %def rcl_controller_define_processes @ Revert to initial state. Also, reset RECOLA (only if it has already done something). <>= procedure :: reset => rcl_controller_reset <>= subroutine rcl_controller_reset (rcl) class(rcl_controller_t), intent(inout) :: rcl if (rcl%active .or. rcl%done) then if (debug_on) call msg_debug2 (D_ME_METHODS, "reset_recola_rcl") if (allocated (rcl%processes)) deallocate (rcl%processes) call reset_recola_rcl () end if rcl%active = .false. rcl%defined = .false. rcl%done = .false. rcl%recola_id = 0 rcl%n_processes = 0 end subroutine rcl_controller_reset @ %def rcl_controller_reset @ Output. <>= procedure :: write => rcl_controller_write <>= subroutine rcl_controller_write (object, unit) class(rcl_controller_t), intent(in) :: object integer, intent(in), optional :: unit integer :: u u = given_output_unit (unit) write (u, "(1x,A,2(1x,A,L1),2(1x,A,I0))") "RECOLA controller:", & "active=", object%active, "done=", object%done, & "id=", object%recola_id, "n_processes=", object%n_processes end subroutine rcl_controller_write @ %def rcl_controller_write @ Return a new numeric process ID, incrementing the counter once. <>= procedure :: get_new_id => rcl_controller_get_new_id <>= subroutine rcl_controller_get_new_id (object, id) class(rcl_controller_t), intent(inout) :: object integer, intent(out) :: id object%recola_id = object%recola_id + 1 id = object%recola_id end subroutine rcl_controller_get_new_id @ %def rcl_controller_get_new_id @ Return the current numeric process ID without incrementing the counter. <>= procedure :: get_current_id => rcl_controller_get_current_id <>= subroutine rcl_controller_get_current_id (object, id) class(rcl_controller_t), intent(inout) :: object integer, intent(out) :: id id = object%recola_id end subroutine rcl_controller_get_current_id @ %def rcl_controller_get_current_id @ Do not allow activation if processes have been calculated previously. Otherwise set the flag. <>= procedure :: activate => rcl_controller_activate <>= subroutine rcl_controller_activate (rcl) class(rcl_controller_t), intent(inout) :: rcl if ( .not. allocated(rcl%processes) ) allocate ( rcl%processes(10) ) rcl_controller%active = .true. end subroutine rcl_controller_activate @ %def rcl_controller_activate @ Start process initialization by calling the RECOLA API. Do not allow this twice (skip silently), and skip anyway if there is no activation. <>= procedure :: generate_processes => rcl_controller_generate_processes <>= subroutine rcl_controller_generate_processes (rcl) class(rcl_controller_t), intent(inout) :: rcl if (rcl_controller%active) then if (.not. rcl_controller%done) then call msg_message ("Recola: preparing processes for integration") call generate_processes_rcl () rcl_controller%done = .true. end if end if end subroutine rcl_controller_generate_processes @ %def rcl_controller_generate_processes @ Return a new numeric RECOLA process ID. The singleton nature of the controller guarantees that the ID is unique. <>= public :: rclwrap_get_new_recola_id <>= subroutine rclwrap_get_new_recola_id (id) integer, intent(out) :: id call rcl_controller%get_new_id (id) end subroutine rclwrap_get_new_recola_id @ %def rclwrap_get_new_recola_id @ Return the current numeric RECOLA process ID. This coincides with the amount of IDs currently in use. <>= public :: rclwrap_get_current_recola_id <>= function rclwrap_get_current_recola_id () result (n) integer :: n call rcl_controller%get_current_id (n) end function rclwrap_get_current_recola_id @ %def rclwrap_get_current_recola_id @ This procedure records the fact that there is a recola process pending, so we will have to call [[generate_processes]] before we can calculate anything with Recola. <>= public :: rclwrap_request_generate_processes <>= subroutine rclwrap_request_generate_processes () if (debug_on) call msg_debug2 (D_ME_METHODS, "request_generate_processes_rcl") call rcl_controller%activate () end subroutine rclwrap_request_generate_processes @ %def rclwrap_request_generate_processes @ Add a process to be defined later <>= public :: rclwrap_add_process <>= subroutine rclwrap_add_process (id, process_string, order) integer, intent(in) :: id type(string_t), intent(in) :: process_string, order type(rcl_process_t) :: prc if (debug_on) call msg_debug2 (D_ME_METHODS, "add_process_rcl: id", id) prc = rcl_process_t (id, process_string, order) call rcl_controller%add_process (prc) end subroutine rclwrap_add_process @ %def rclwrap_add_process @ Define all added processes. Reset if processes were already defined. <>= public :: rclwrap_define_processes <>= subroutine rclwrap_define_processes () if (debug_on) call msg_debug2 (D_ME_METHODS, "define_processes_rcl") call rcl_controller%define_processes () end subroutine rclwrap_define_processes @ %def rclwrap_define_processes @ We call this after all processes have been added and defined, so RECOLA can initialize itself for integration. <>= public :: rclwrap_generate_processes <>= subroutine rclwrap_generate_processes () if (debug_on) call msg_debug2 (D_ME_METHODS, "generate_processes_rcl") call rcl_controller%generate_processes () end subroutine rclwrap_generate_processes @ %def rclwrap_generate_processes @ <>= public :: rclwrap_compute_process <>= subroutine rclwrap_compute_process (id, p, order, sqme) integer, intent(in) :: id real(double), intent(in), dimension(:,:) :: p character(len=*), intent(in) :: order real(double), intent(out), dimension(0:1), optional :: sqme if (debug_on) call msg_debug2 (D_ME_METHODS, "compute_process_rcl") call compute_process_rcl (id, p, order, sqme) end subroutine rclwrap_compute_process @ %def rclwrap_compute_process @ <>= public :: rclwrap_get_amplitude <>= subroutine rclwrap_get_amplitude (id, g_power, order, col, hel, amp) integer, intent(in) :: id, g_power character(len=*), intent(in) :: order integer, dimension(:), intent(in) :: col, hel complex(double), intent(out) :: amp if (debug_on) call msg_debug2 (D_ME_METHODS, "get_amplitude_rcl") call get_amplitude_rcl (id, g_power, order, col, hel, amp) end subroutine rclwrap_get_amplitude @ %def rclwrap_get_amplitude @ <>= public :: rclwrap_get_squared_amplitude <>= subroutine rclwrap_get_squared_amplitude (id, alphas_power, order, sqme) integer, intent(in) :: id, alphas_power character(len=*), intent(in) :: order real(double), intent(out) :: sqme if (debug_on) call msg_debug2 (D_ME_METHODS, "get_squared_amplitude_rcl") call get_squared_amplitude_rcl (id, alphas_power, order, sqme) end subroutine rclwrap_get_squared_amplitude @ %def rclwrap_get_squared_amplitude @ <>= public :: rclwrap_set_pole_mass <>= subroutine rclwrap_set_pole_mass (pdg_id, mass, width) integer, intent(in) :: pdg_id real(double), intent(in) :: mass, width if (debug_on) call msg_debug2 (D_ME_METHODS, "rclwrap_set_pole_mass of ", pdg_id) select case (abs(pdg_id)) case (11) if (width > zero) & call msg_fatal ("Recola pole mass: Attempting to set non-zero electron width!") call set_pole_mass_electron_rcl (mass) case (13) call set_pole_mass_muon_rcl (mass, width) case (15) call set_pole_mass_tau_rcl (mass, width) case (1) if (width > zero) & call msg_fatal ("Recola pole mass: Attempting to set non-zero down-quark width!") call set_pole_mass_down_rcl (mass) case (2) if (width > zero) & call msg_fatal ("Recola pole mass: Attempting to set non-zero up-quark width!") call set_pole_mass_up_rcl (mass) case (3) if (width > zero) & call msg_fatal ("Recola pole mass: Attempting to set non-zero strange-quark width!") call set_pole_mass_strange_rcl (mass) case (4) call set_pole_mass_charm_rcl (mass, width) case (5) call set_pole_mass_bottom_rcl (mass, width) case (6) call set_pole_mass_top_rcl (mass, width) case (23) call set_pole_mass_z_rcl (mass, width) case (24) call set_pole_mass_w_rcl (mass, width) case (25) call set_pole_mass_h_rcl (mass, width) case default call msg_fatal ("Recola pole mass: Unsupported particle") end select end subroutine rclwrap_set_pole_mass @ %def rclwrap_set_pole_mass @ <>= public :: rclwrap_set_onshell_mass <>= subroutine rclwrap_set_onshell_mass (pdg_id, mass, width) integer, intent(in) :: pdg_id real(double), intent(in) :: mass, width if (debug_on) call msg_debug2 (D_ME_METHODS, "rclwrap_set_onshell_mass of ", pdg_id) select case (abs(pdg_id)) case (23) call set_onshell_mass_z_rcl (mass, width) case (24) call set_onshell_mass_w_rcl (mass, width) case default call msg_fatal ("Recola onshell mass: Only for W and Z") end select end subroutine rclwrap_set_onshell_mass @ %def rclwrap_set_onshell_mass @ <>= public :: rclwrap_use_gfermi_scheme <>= subroutine rclwrap_use_gfermi_scheme (gf) real(double), intent(in), optional :: gf if (debug_on) call msg_debug2 (D_ME_METHODS, "use_gfermi_scheme_rcl", & real(gf, kind=default)) call use_gfermi_scheme_rcl (gf) end subroutine rclwrap_use_gfermi_scheme @ %def rclwrap_use_gfermi_scheme @ <>= public :: rclwrap_set_light_fermions <>= subroutine rclwrap_set_light_fermions (m) real(double), intent(in) :: m if (debug_on) call msg_debug2 (D_ME_METHODS, "set_light_fermions_rcl", & real(m, kind=default)) call set_light_fermions_rcl (m) end subroutine rclwrap_set_light_fermions @ %def rclwrap_set_light_fermions @ <>= public :: rclwrap_set_light_fermion <>= subroutine rclwrap_set_light_fermion (pdg_id) integer, intent(in) :: pdg_id if (debug_on) call msg_debug2 (D_ME_METHODS, "rclwrap_set_light_fermion", pdg_id) select case (abs(pdg_id)) case (1) call set_light_down_rcl () case (2) call set_light_up_rcl () case (3) call set_light_strange_rcl () case (4) call set_light_charm_rcl () case (5) call set_light_bottom_rcl () case (6) call set_light_top_rcl () case (11) call set_light_electron_rcl () case (13) call set_light_muon_rcl () case (15) call set_light_tau_rcl () end select end subroutine rclwrap_set_light_fermion @ %def rclwrap_set_light_fermion @ <>= public :: rclwrap_unset_light_fermion <>= subroutine rclwrap_unset_light_fermion (pdg_id) integer, intent(in) :: pdg_id if (debug_on) call msg_debug2 (D_ME_METHODS, "rclwrap_unset_light_fermion", pdg_id) select case (abs(pdg_id)) case (1) call unset_light_down_rcl () case (2) call unset_light_up_rcl () case (3) call unset_light_strange_rcl () case (4) call unset_light_charm_rcl () case (5) call unset_light_bottom_rcl () case (6) call unset_light_top_rcl () case (11) call unset_light_electron_rcl () case (13) call unset_light_muon_rcl () case (15) call unset_light_tau_rcl () end select end subroutine rclwrap_unset_light_fermion @ %def rclwrap_unset_light_fermion @ <>= public :: rclwrap_set_onshell_scheme <>= subroutine rclwrap_set_onshell_scheme if (debug_on) call msg_debug2 (D_ME_METHODS, "set_on_shell_scheme_rcl") call set_on_shell_scheme_rcl () end subroutine rclwrap_set_onshell_scheme @ %def rclwrap_set_onshell_scheme @ <>= public :: rclwrap_set_alpha_s <>= subroutine rclwrap_set_alpha_s (alpha_s, mu, nf) real(double), intent(in) :: alpha_s, mu integer, intent(in) :: nf if (debug_on) call msg_debug2 (D_ME_METHODS, "set_alphas_rcl") call set_alphas_rcl (alpha_s, mu, nf) end subroutine rclwrap_set_alpha_s @ %def rclwrap_set_alpha_s @ <>= public :: rclwrap_get_alpha_s <>= function rclwrap_get_alpha_s () result (alpha_s) real(double) :: alpha_s if (debug_on) call msg_debug2 (D_ME_METHODS, "get_alphas_rcl") call get_alphas_rcl (alpha_s) end function rclwrap_get_alpha_s @ %def rclwrap_get_alpha_s @ <>= public :: rclwrap_get_helicity_configurations <>= subroutine rclwrap_get_helicity_configurations (id, hel) integer, intent(in) :: id integer, intent(inout), dimension(:,:), allocatable :: hel call get_helicity_configurations_rcl (id, hel) end subroutine rclwrap_get_helicity_configurations @ %def rclwrap_get_helicity_configurations @ <>= public :: rclwrap_get_color_configurations <>= subroutine rclwrap_get_color_configurations (id, col) integer, intent(in) :: id integer, intent(out), dimension(:,:), allocatable :: col call get_colour_configurations_rcl (id, col) end subroutine rclwrap_get_color_configurations @ %def rclwrap_get_color_configurations @ Selects dimensional regularization for soft singularities. <>= public :: rclwrap_use_dim_reg_soft <>= subroutine rclwrap_use_dim_reg_soft () if (debug_on) call msg_debug2 (D_ME_METHODS, "use_dim_reg_soft_rcl") call use_dim_reg_soft_rcl () end subroutine rclwrap_use_dim_reg_soft @ %def rclwrap_use_dim_reg_soft @ Selects mass regularization for soft singularities and sets the mass regulator in GeV to [[m]]. <>= public :: rclwrap_use_mass_reg_soft <>= subroutine rclwrap_use_mass_reg_soft (m) real(double), intent(in) :: m if (debug_on) call msg_debug2 (D_ME_METHODS, "use_mass_reg_soft_rcl") call use_mass_reg_soft_rcl (m) end subroutine rclwrap_use_mass_reg_soft @ %def rclwrap_use_mass_reg_soft @ Sets the UV pole parameterization $\Delta_{UV}$. <>= public :: rclwrap_set_delta_uv <>= subroutine rclwrap_set_delta_uv (d) real(double), intent(in) :: d if (debug_on) call msg_debug2 (D_ME_METHODS, "set_delta_uv_rcl") call set_delta_uv_rcl (d) end subroutine rclwrap_set_delta_uv @ %def rclwrap_set_delta_uv @ <>= public :: rclwrap_set_mu_uv <>= subroutine rclwrap_set_mu_uv (mu) real(double), intent(in) :: mu if (debug_on) call msg_debug2 (D_ME_METHODS, "set_mu_uv_rcl") call set_mu_uv_rcl (mu) end subroutine rclwrap_set_mu_uv @ %def rclwrap_set_mu_uv @ Sets the IR pole parameterizations $\Delta_{IR}$ and $\Delta_2$. <>= public :: rclwrap_set_delta_ir <>= subroutine rclwrap_set_delta_ir (d, d2) real(double), intent(in) :: d, d2 if (debug_on) call msg_debug2 (D_ME_METHODS, "set_delta_ir_rcl", & real(d, kind=default)) if (debug_on) call msg_debug2 (D_ME_METHODS, "set_delta_ir_rcl", & real(d2, kind=default)) call set_delta_ir_rcl (d, d2) end subroutine rclwrap_set_delta_ir @ %def rclwrap_set_delta_ir @ <>= public :: rclwrap_set_mu_ir <>= subroutine rclwrap_set_mu_ir (mu) real(double), intent(in) :: mu if (debug_on) call msg_debug2 (D_ME_METHODS, "set_mu_ir_rcl") call set_mu_ir_rcl (mu) end subroutine rclwrap_set_mu_ir @ %def rclwrap_set_mu_ir @ <>= public :: rclwrap_get_renormalization_scale <>= subroutine rclwrap_get_renormalization_scale (mu) real(double), intent(out) :: mu if (debug_on) call msg_debug2 (D_ME_METHODS, "get_renormalization_scale_rcl") call get_renormalization_scale_rcl (mu) end subroutine rclwrap_get_renormalization_scale @ %def rclwrap_get_renormalization_scale @ <>= public :: rclwrap_get_flavor_scheme <>= subroutine rclwrap_get_flavor_scheme (nf) integer, intent(out) :: nf if (debug_on) call msg_debug2 (D_ME_METHODS, "get_flavour_scheme_rcl") call get_flavour_scheme_rcl (nf) end subroutine rclwrap_get_flavor_scheme @ %def rclwrap_get_flavor_scheme @ <>= public :: rclwrap_use_alpha0_scheme <>= subroutine rclwrap_use_alpha0_scheme (al0) real(double), intent(in), optional :: al0 if (debug_on) call msg_debug2 (D_ME_METHODS, "use_alpha0_scheme_rcl") call use_alpha0_scheme_rcl (al0) end subroutine rclwrap_use_alpha0_scheme @ %def rclwrap_use_alpha0_scheme @ <>= public :: rclwrap_use_alphaz_scheme <>= subroutine rclwrap_use_alphaz_scheme (alz) real(double), intent(in), optional :: alz if (debug_on) call msg_debug2 (D_ME_METHODS, "use_alphaz_scheme_rcl") call use_alphaz_scheme_rcl (alz) end subroutine rclwrap_use_alphaz_scheme @ %def rclwrap_use_alphaz_scheme @ <>= public :: rclwrap_set_complex_mass_scheme <>= subroutine rclwrap_set_complex_mass_scheme () if (debug_on) call msg_debug2 (D_ME_METHODS, "set_complex_mass_scheme_rcl") call set_complex_mass_scheme_rcl () end subroutine rclwrap_set_complex_mass_scheme @ %def rclwrap_set_complex_mass_scheme @ <>= public :: rclwrap_set_resonant_particle <>= subroutine rclwrap_set_resonant_particle (pdg_id) integer, intent(in) :: pdg_id if (debug_on) call msg_debug2 (D_ME_METHODS, "set_resonant_particle_rcl") call set_resonant_particle_rcl (char(get_recola_particle_string (pdg_id))) end subroutine rclwrap_set_resonant_particle @ %def rclwrap_set_resonant_particle @ <>= public :: rclwrap_switch_on_resonant_self_energies <>= subroutine rclwrap_switch_on_resonant_self_energies () if (debug_on) call msg_debug2 (D_ME_METHODS, "switchon_resonant_selfenergies_rcl") call switchon_resonant_selfenergies_rcl () end subroutine rclwrap_switch_on_resonant_self_energies @ %def rclwrap_switch_on_resonant_self_energies @ <>= public :: rclwrap_switch_off_resonant_self_energies <>= subroutine rclwrap_switch_off_resonant_self_energies () if (debug_on) call msg_debug2 (D_ME_METHODS, "switchoff_resonant_selfenergies_rcl") call switchoff_resonant_selfenergies_rcl () end subroutine rclwrap_switch_off_resonant_self_energies @ %def rclwrap_switch_off_resonant_self_energies @ <>= public :: rclwrap_set_draw_level_branches <>= subroutine rclwrap_set_draw_level_branches (n) integer, intent(in) :: n if (debug_on) call msg_debug2 (D_ME_METHODS, "set_draw_level_branches_rcl") call set_draw_level_branches_rcl (n) end subroutine rclwrap_set_draw_level_branches @ %def rclwrap_set_draw_level_branches @ <>= public :: rclwrap_set_print_level_amplitude <>= subroutine rclwrap_set_print_level_amplitude (n) integer, intent(in) :: n if (debug_on) call msg_debug2 (D_ME_METHODS, "set_print_level_amplitude_rcl") call set_print_level_amplitude_rcl (n) end subroutine rclwrap_set_print_level_amplitude @ %def rclwrap_set_print_level_amplitude @ <>= public :: rclwrap_set_print_level_squared_amplitude <>= subroutine rclwrap_set_print_level_squared_amplitude (n) integer, intent(in) :: n if (debug_on) call msg_debug2 (D_ME_METHODS, "set_print_level_squared_amplitude_rcl") call set_print_level_squared_amplitude_rcl (n) end subroutine rclwrap_set_print_level_squared_amplitude @ %def rclwrap_set_print_level_squared_amplitude @ <>= public :: rclwrap_set_print_level_correlations <>= subroutine rclwrap_set_print_level_correlations (n) integer, intent(in) :: n if (debug_on) call msg_debug2 (D_ME_METHODS, "set_print_level_correlations_rcl") call set_print_level_correlations_rcl (n) end subroutine rclwrap_set_print_level_correlations @ %def rclwrap_set_print_level_correlations @ <>= public :: rclwrap_set_print_level_RAM <>= subroutine rclwrap_set_print_level_RAM (n) integer, intent(in) :: n if (debug_on) call msg_debug2 (D_ME_METHODS, "set_print_level_RAM_rcl") call set_print_level_RAM_rcl (n) end subroutine rclwrap_set_print_level_RAM @ %def rclwrap_set_print_level_RAM @ <>= public :: rclwrap_scale_coupling3 <>= subroutine rclwrap_scale_coupling3 (pdg_id1, pdg_id2, pdg_id3, factor) integer, intent(in) :: pdg_id1, pdg_id2, pdg_id3 complex(double), intent(in) :: factor if (debug_on) call msg_debug2 (D_ME_METHODS, "scale_coupling3_rcl") call scale_coupling3_rcl (factor, char(get_recola_particle_string (pdg_id1)), & char(get_recola_particle_string (pdg_id2)), char(get_recola_particle_string (pdg_id3))) end subroutine rclwrap_scale_coupling3 @ %def rclwrap_scale_coupling3 @ <>= public :: rclwrap_scale_coupling4 <>= subroutine rclwrap_scale_coupling4 (pdg_id1, pdg_id2, pdg_id3, pdg_id4, factor) integer, intent(in) :: pdg_id1, pdg_id2, pdg_id3, pdg_id4 complex(double), intent(in) :: factor if (debug_on) call msg_debug2 (D_ME_METHODS, "scale_coupling4_rcl") call scale_coupling4_rcl (factor, char(get_recola_particle_string (pdg_id1)), & char(get_recola_particle_string (pdg_id2)), char(get_recola_particle_string (pdg_id3)), & char(get_recola_particle_string (pdg_id4))) end subroutine rclwrap_scale_coupling4 @ %def rclwrap_scale_coupling4 @ <>= public :: rclwrap_switch_off_coupling3 <>= subroutine rclwrap_switch_off_coupling3 (pdg_id1, pdg_id2, pdg_id3) integer, intent(in) :: pdg_id1, pdg_id2, pdg_id3 if (debug_on) call msg_debug2 (D_ME_METHODS, "switchoff_coupling3_rcl") call switchoff_coupling3_rcl (char(get_recola_particle_string (pdg_id1)), & char(get_recola_particle_string (pdg_id2)), char(get_recola_particle_string (pdg_id3))) end subroutine rclwrap_switch_off_coupling3 @ %def rclwrap_switch_off_coupling3 @ <>= public :: rclwrap_switch_off_coupling4 <>= subroutine rclwrap_switch_off_coupling4 (pdg_id1, pdg_id2, pdg_id3, pdg_id4) integer, intent(in) :: pdg_id1, pdg_id2, pdg_id3, pdg_id4 if (debug_on) call msg_debug2 (D_ME_METHODS, "switchoff_coupling4_rcl") call switchoff_coupling4_rcl & (char(get_recola_particle_string (pdg_id1)), & char(get_recola_particle_string (pdg_id2)), & char(get_recola_particle_string (pdg_id3)), & char(get_recola_particle_string (pdg_id4))) end subroutine rclwrap_switch_off_coupling4 @ %def rclwrap_switch_off_coupling4 @ <>= public :: rclwrap_set_ifail <>= subroutine rclwrap_set_ifail (i) integer, intent(in) :: i if (debug_on) call msg_debug2 (D_ME_METHODS, "set_ifail_rcl") call set_ifail_rcl (i) end subroutine rclwrap_set_ifail @ %def rclwrap_set_ifail @ <>= public :: rclwrap_get_ifail <>= subroutine rclwrap_get_ifail (i) integer, intent(out) :: i if (debug_on) call msg_debug2 (D_ME_METHODS, "get_ifail_rcl") call get_ifail_rcl (i) end subroutine rclwrap_get_ifail @ %def rclwrap_get_ifail @ <>= public :: rclwrap_set_output_file <>= subroutine rclwrap_set_output_file (filename) character(len=*), intent(in) :: filename if (debug_on) call msg_debug2 (D_ME_METHODS, "set_output_file_rcl") call set_output_file_rcl (filename) end subroutine rclwrap_set_output_file @ %def rclwrap_set_output_file @ <>= public :: rclwrap_set_gs_power <>= subroutine rclwrap_set_gs_power (id, gs_array) integer, intent(in) :: id integer, dimension(:,:), intent(in) :: gs_array if (debug_on) call msg_debug2 (D_ME_METHODS, "set_gs_power_rcl") call set_gs_power_rcl (id, gs_array) end subroutine rclwrap_set_gs_power @ %def rclwrap_set_gs_power @ <>= public :: rclwrap_select_gs_power_born_amp <>= subroutine rclwrap_select_gs_power_born_amp (id, gs_power) integer, intent(in) :: id, gs_power if (debug_on) call msg_debug2 (D_ME_METHODS, "select_gs_power_BornAmpl_rcl") call select_gs_power_BornAmpl_rcl (id, gs_power) end subroutine rclwrap_select_gs_power_born_amp @ %def rclwrap_select_gs_power_born_amp @ <>= public :: rclwrap_unselect_gs_power_born_amp <>= subroutine rclwrap_unselect_gs_power_born_amp (id, gs_power) integer, intent(in) :: id, gs_power if (debug_on) call msg_debug2 (D_ME_METHODS, "unselect_gs_power_BornAmpl_rcl") call unselect_gs_power_BornAmpl_rcl (id, gs_power) end subroutine rclwrap_unselect_gs_power_born_amp @ %def rclwrap_unselect_gs_power_born_amp @ <>= public :: rclwrap_select_gs_power_loop_amp <>= subroutine rclwrap_select_gs_power_loop_amp (id, gs_power) integer, intent(in) :: id, gs_power if (debug_on) call msg_debug2 (D_ME_METHODS, "select_gs_power_LoopAmpl_rcl") call select_gs_power_LoopAmpl_rcl (id, gs_power) end subroutine rclwrap_select_gs_power_loop_amp @ %def rclwrap_select_gs_power_loop_amp @ <>= public :: rclwrap_unselect_gs_power_loop_amp <>= subroutine rclwrap_unselect_gs_power_loop_amp (id, gs_power) integer, intent(in) :: id, gs_power if (debug_on) call msg_debug2 (D_ME_METHODS, "unselect_gs_power_LoopAmpl_rcl") call unselect_gs_power_LoopAmpl_rcl (id, gs_power) end subroutine rclwrap_unselect_gs_power_loop_amp @ %def rclwrap_unselect_gs_power_loop_amp @ <>= public :: rclwrap_select_all_gs_powers_born_amp <>= subroutine rclwrap_select_all_gs_powers_born_amp (id) integer, intent(in) :: id if (debug_on) call msg_debug2 (D_ME_METHODS, "select_all_gs_powers_BornAmpl_rcl") call select_all_gs_powers_BornAmpl_rcl (id) end subroutine rclwrap_select_all_gs_powers_born_amp @ %def rclwrap_select_all_gs_powers_born_amp @ <>= public :: rclwrap_unselect_all_gs_powers_loop_amp <>= subroutine rclwrap_unselect_all_gs_powers_loop_amp (id) integer, intent(in) :: id if (debug_on) call msg_debug2 (D_ME_METHODS, "unselect_all_gs_powers_BornAmpl_rcl") call unselect_all_gs_powers_BornAmpl_rcl (id) end subroutine rclwrap_unselect_all_gs_powers_loop_amp @ %def rclwrap_unselect_all_gs_powers_loop_amp @ <>= public :: rclwrap_select_all_gs_powers_loop_amp <>= subroutine rclwrap_select_all_gs_powers_loop_amp (id) integer, intent(in) :: id if (debug_on) call msg_debug2 (D_ME_METHODS, "select_all_gs_powers_LoopAmpl_rcl") call select_all_gs_powers_LoopAmpl_rcl (id) end subroutine rclwrap_select_all_gs_powers_loop_amp @ %def rclwrap_select_all_gs_powers_loop_amp @ <>= public :: rclwrap_unselect_all_gs_powers_born_amp <>= subroutine rclwrap_unselect_all_gs_powers_born_amp (id) integer, intent(in) :: id if (debug_on) call msg_debug2 (D_ME_METHODS, "unselect_all_gs_powers_LoopAmpl_rcl") call unselect_all_gs_powers_LoopAmpl_rcl (id) end subroutine rclwrap_unselect_all_gs_powers_born_amp @ %def rclwrap_unselect_all_gs_powers_born_amp @ <>= public :: rclwrap_set_resonant_squared_momentum <>= subroutine rclwrap_set_resonant_squared_momentum (id, i_res, p2) integer, intent(in) :: id, i_res real(double), intent(in) :: p2 if (debug_on) call msg_debug2 (D_ME_METHODS, "set_resonant_squared_momentum_rcl") call set_resonant_squared_momentum_rcl (id, i_res, p2) end subroutine rclwrap_set_resonant_squared_momentum @ %def rclwrap_set_resonant_squared_momentum @ <>= public :: rclwrap_compute_running_alpha_s <>= subroutine rclwrap_compute_running_alpha_s (Q, nf, n_loops) real(double), intent(in) :: Q integer, intent(in) :: nf, n_loops if (debug_on) call msg_debug2 (D_ME_METHODS, "compute_running_alphas_rcl") call compute_running_alphas_rcl (Q, nf, n_loops) end subroutine rclwrap_compute_running_alpha_s @ %def rclwrap_compute_running_alpha_s @ <>= public :: rclwrap_set_dynamic_settings <>= subroutine rclwrap_set_dynamic_settings () if (debug_on) call msg_debug2 (D_ME_METHODS, "set_dynamic_settings_rcl") call set_dynamic_settings_rcl (1) end subroutine rclwrap_set_dynamic_settings @ %def rclwrap_set_dynamic_settings @ <>= public :: rclwrap_rescale_process <>= subroutine rclwrap_rescale_process (id, order, sqme) integer, intent(in) :: id character(len=*), intent(in) :: order real(double), dimension(0:1), intent(out), optional :: sqme if (debug_on) call msg_debug2 (D_ME_METHODS, "rescale_process_rcl") call rescale_process_rcl (id, order, sqme) end subroutine rclwrap_rescale_process @ %def rclwrap_rescale_process @ <>= public :: rclwrap_get_polarized_squared_amplitude <>= subroutine rclwrap_get_polarized_squared_amplitude (id, & alphas_power, order, hel, sqme) integer, intent(in) :: id, alphas_power character(len=*), intent(in) :: order integer, dimension(:), intent(in) :: hel real(double), intent(out) :: sqme if (debug_on) call msg_debug2 (D_ME_METHODS, "get_polarized_squared_amplitude_rcl") call get_polarized_squared_amplitude_rcl (id, alphas_power, & order, hel, sqme) end subroutine rclwrap_get_polarized_squared_amplitude @ %def rclwrap_get_polarized_squared_amplitude @ <>= public :: rclwrap_compute_color_correlation <>= subroutine rclwrap_compute_color_correlation (id, p, & i1, i2, sqme) integer, intent(in) :: id real(double), dimension(:,:), intent(in) :: p integer, intent(in) :: i1, i2 real(double), intent(out), optional :: sqme if (debug_on) call msg_debug2 (D_ME_METHODS, "compute_colour_correlation_rcl") call compute_colour_correlation_rcl (id, p, i1, i2, sqme) end subroutine rclwrap_compute_color_correlation @ %def rclwrap_compute_color_correlation @ <>= public :: rclwrap_compute_all_color_correlations <>= subroutine rclwrap_compute_all_color_correlations (id, p) integer, intent(in) :: id real(double), dimension(:,:), intent(in) :: p if (debug_on) call msg_debug2 (D_ME_METHODS, "compute_all_colour_correlations_rcl") call compute_all_colour_correlations_rcl (id, p) end subroutine rclwrap_compute_all_color_correlations @ %def rclwrap_compute_all_color_correlations @ <>= public :: rclwrap_rescale_color_correlation <>= subroutine rclwrap_rescale_color_correlation (id, i1, i2, sqme) integer, intent(in) :: id, i1, i2 real(double), intent(out), optional :: sqme if (debug_on) call msg_debug2 (D_ME_METHODS, "rescale_colour_correlation_rcl") call rescale_colour_correlation_rcl (id, i1, i2, sqme) end subroutine rclwrap_rescale_color_correlation @ %def rclwrap_rescale_color_correlation @ <>= public :: rclwrap_rescale_all_color_correlations <>= subroutine rclwrap_rescale_all_color_correlations (id) integer, intent(in) :: id if (debug_on) call msg_debug2 (D_ME_METHODS, "rescale_all_colour_correlations_rcl") call rescale_all_colour_correlations_rcl (id) end subroutine rclwrap_rescale_all_color_correlations @ %def rclwrap_rescale_all_color_correlations @ <>= public :: rclwrap_get_color_correlation <>= subroutine rclwrap_get_color_correlation (id, alphas_power, i1, i2, sqme) integer, intent(in) :: id, alphas_power, i1, i2 real(double), intent(out) :: sqme if (debug_on) call msg_debug2 (D_ME_METHODS, "get_colour_correlation_rcl") call get_colour_correlation_rcl (id, alphas_power, i1, i2, sqme) end subroutine rclwrap_get_color_correlation @ %def rclwrap_get_color_correlation @ <>= public :: rclwrap_compute_spin_correlation <>= subroutine rclwrap_compute_spin_correlation (id, p, i_photon, pol, sqme) integer, intent(in) :: id real(double), dimension(:,:), intent(in) :: p integer, intent(in) :: i_photon complex(double), dimension(:), intent(in) :: pol real(double), intent(out), optional :: sqme if (debug_on) call msg_debug2 (D_ME_METHODS, "compute_spin_correlation_rcl") call compute_spin_correlation_rcl (id, p, i_photon, pol, sqme) end subroutine rclwrap_compute_spin_correlation @ %def rclwrap_compute_spin_correlation @ <>= public :: rclwrap_rescale_spin_correlation <>= subroutine rclwrap_rescale_spin_correlation (id, i_photon, pol, sqme) integer, intent(in) :: id, i_photon complex(double), dimension(:), intent(in) :: pol real(double), intent(out), optional :: sqme if (debug_on) call msg_debug2 (D_ME_METHODS, "rescale_spin_correlation_rcl") call rescale_spin_correlation_rcl (id, i_photon, pol, sqme) end subroutine rclwrap_rescale_spin_correlation @ %def rclwrap_rescale_spin_correlation @ <>= public :: rclwrap_get_spin_correlation <>= subroutine rclwrap_get_spin_correlation (id, alphas_power, sqme) integer, intent(in) :: id, alphas_power real(double), intent(out) :: sqme if (debug_on) call msg_debug2 (D_ME_METHODS, "get_spin_correlation_rcl") call get_spin_correlation_rcl (id, alphas_power, sqme) end subroutine rclwrap_get_spin_correlation @ %def rclwrap_get_spin_correlation @ <>= public :: rclwrap_compute_spin_color_correlation <>= subroutine rclwrap_compute_spin_color_correlation (id, p, & i_gluon, i_spectator, pol, sqme) integer, intent(in) :: id real(double), dimension(:,:), intent(in) :: p integer, intent(in) :: i_gluon, i_spectator complex(double), dimension(:), intent(in) :: pol real(double), intent(out), optional :: sqme if (debug_on) call msg_debug2 (D_ME_METHODS, "compute_spin_colour_correlation_rcl") call compute_spin_colour_correlation_rcl (id, p, & i_gluon, i_spectator, pol, sqme) end subroutine rclwrap_compute_spin_color_correlation @ %def rclwrap_compute_spin_color_correlation @ <>= public :: rclwrap_rescale_spin_color_correlation <>= subroutine rclwrap_rescale_spin_color_correlation (id, i_gluon, & i_spectator, pol, sqme) integer, intent(in) :: id, i_gluon, i_spectator complex(double), dimension(:), intent(in) :: pol real(double), intent(out), optional :: sqme if (debug_on) call msg_debug2 (D_ME_METHODS, "rescale_spin_colour_correlation_rcl") call rescale_spin_colour_correlation_rcl (id, i_gluon, & i_spectator, pol, sqme) end subroutine rclwrap_rescale_spin_color_correlation @ %def rclwrap_rescale_spin_color_correlation @ <>= public :: rclwrap_get_spin_color_correlation <>= subroutine rclwrap_get_spin_color_correlation (id, alphas_power, & i_gluon, i_spectator, sqme) integer, intent(in) :: id, alphas_power, i_gluon, i_spectator real(double), intent(out) :: sqme if (debug_on) call msg_debug2 (D_ME_METHODS, "get_spin_colour_correlation_rcl") call get_spin_colour_correlation_rcl (id, alphas_power, & i_gluon, i_spectator, sqme) end subroutine rclwrap_get_spin_color_correlation @ %def rclwrap_get_spin_color_correlation @ <>= public :: rclwrap_get_momenta <>= subroutine rclwrap_get_momenta (id, p) integer, intent(in) :: id real(double), dimension(:,:), intent(out) :: p if (debug_on) call msg_debug2 (D_ME_METHODS, "get_momenta_rcl") call get_momenta_rcl (id, p) end subroutine rclwrap_get_momenta @ %def rclwrap_get_momenta @ The reset routine is essential. But note that it doesn't reset the Recola parameters, just the processes. For LOL, Recola's reset routine crashes the program if there was no process before. So, rather reset indirectly via the controller. <>= public :: rclwrap_reset_recola <>= subroutine rclwrap_reset_recola if (debug_on) call msg_debug (D_ME_METHODS, "rclwrap_reset_recola") call rcl_controller%reset () end subroutine rclwrap_reset_recola @ %def rclwrap_reset_recola @ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Recola dummy replacement module} <<[[recola_wrapper_dummy.f90]]>>= <> module recola_wrapper use kinds <> <> <> <> contains <> end module recola_wrapper @ %def recola_wrapper_dummy @ <>= public :: rclwrap_is_active <>= logical, parameter :: rclwrap_is_active = .false. @ %def rclwrap_is_active @ <>= public :: get_recola_particle_string <>= elemental function get_recola_particle_string (pdg) result (name) type(string_t) :: name integer, intent(in) :: pdg name = var_str ("?") end function get_recola_particle_string @ %def get_recola_paritcle_string @ <>= public :: rclwrap_get_new_recola_id <>= subroutine rclwrap_get_new_recola_id (id) integer, intent(out) :: id id = 0 end subroutine rclwrap_get_new_recola_id @ %def rclwrap_get_new_recola_id @ <>= public :: rclwrap_get_current_recola_id <>= function rclwrap_get_current_recola_id () result (n) integer :: n n = 0 end function rclwrap_get_current_recola_id @ %def rclwrap_get_current_recola_id @ <>= public :: rclwrap_request_generate_processes <>= subroutine rclwrap_request_generate_processes () end subroutine rclwrap_request_generate_processes @ %def rclwrap_request_generate_processes @ <>= public :: rclwrap_add_process <>= subroutine rclwrap_add_process (id, process_string, order) integer, intent(in) :: id type(string_t), intent(in) :: process_string, order end subroutine rclwrap_add_process @ %def rclwrap_add_process @ <>= public :: rclwrap_define_processes <>= subroutine rclwrap_define_processes () end subroutine rclwrap_define_processes @ %def rclwrap_define_processes @ <>= public :: rclwrap_generate_processes <>= subroutine rclwrap_generate_processes () end subroutine rclwrap_generate_processes @ %def rclwrap_generate_processes @ <>= public :: rclwrap_compute_process <>= subroutine rclwrap_compute_process (id, p, order, sqme) integer, intent(in) :: id real(double), intent(in), dimension(:,:) :: p character(len=*), intent(in) :: order real(double), intent(out), dimension(0:1), optional :: sqme end subroutine rclwrap_compute_process @ %def rclwrap_compute_process @ <>= public :: rclwrap_get_amplitude <>= subroutine rclwrap_get_amplitude (id, g_power, order, col, hel, amp) integer, intent(in) :: id, g_power character(len=*), intent(in) :: order integer, dimension(:), intent(in) :: col, hel complex(double), intent(out) :: amp end subroutine rclwrap_get_amplitude @ %def rclwrap_get_amplitude @ <>= public :: rclwrap_get_squared_amplitude <>= subroutine rclwrap_get_squared_amplitude (id, alphas_power, order, sqme) integer, intent(in) :: id, alphas_power character(len=*), intent(in) :: order real(double), intent(out) :: sqme end subroutine rclwrap_get_squared_amplitude @ %def rclwrap_get_squared_amplitude @ <>= public :: rclwrap_set_pole_mass <>= subroutine rclwrap_set_pole_mass (pdg_id, mass, width) integer, intent(in) :: pdg_id real(double), intent(in) :: mass, width end subroutine rclwrap_set_pole_mass @ %def rclwrap_set_pole_mass @ <>= public :: rclwrap_set_onshell_mass <>= subroutine rclwrap_set_onshell_mass (pdg_id, mass, width) integer, intent(in) :: pdg_id real(double), intent(in) :: mass, width end subroutine rclwrap_set_onshell_mass @ %def rclwrap_set_onshell_mass @ <>= public :: rclwrap_use_gfermi_scheme <>= subroutine rclwrap_use_gfermi_scheme (gf) real(double), intent(in), optional :: gf end subroutine rclwrap_use_gfermi_scheme @ %def rclwrap_use_gfermi_scheme @ <>= public :: rclwrap_set_light_fermions <>= subroutine rclwrap_set_light_fermions (m) real(double), intent(in) :: m end subroutine rclwrap_set_light_fermions @ %def rclwrap_set_light_fermions @ <>= public :: rclwrap_set_light_fermion <>= subroutine rclwrap_set_light_fermion (pdg_id) integer, intent(in) :: pdg_id end subroutine rclwrap_set_light_fermion @ %def rclwrap_set_light_fermion @ <>= public :: rclwrap_unset_light_fermion <>= subroutine rclwrap_unset_light_fermion (pdg_id) integer, intent(in) :: pdg_id end subroutine rclwrap_unset_light_fermion @ %def rclwrap_unset_light_fermion @ <>= public :: rclwrap_set_onshell_scheme <>= subroutine rclwrap_set_onshell_scheme end subroutine rclwrap_set_onshell_scheme @ %def rclwrap_set_onshell_scheme @ <>= public :: rclwrap_set_alpha_s <>= subroutine rclwrap_set_alpha_s (alpha_s, mu, nf) real(double), intent(in) :: alpha_s, mu integer, intent(in) :: nf end subroutine rclwrap_set_alpha_s @ %def rclwrap_set_alpha_s @ <>= public :: rclwrap_get_alpha_s <>= function rclwrap_get_alpha_s () result (alpha_s) real(double) :: alpha_s end function rclwrap_get_alpha_s @ %def rclwrap_get_alpha_s @ <>= public :: rclwrap_get_helicity_configurations <>= subroutine rclwrap_get_helicity_configurations (id, hel) integer, intent(in) :: id integer, intent(inout), dimension(:,:), allocatable :: hel end subroutine rclwrap_get_helicity_configurations @ %def rclwrap_get_helicity_configurations @ <>= public :: rclwrap_get_color_configurations <>= subroutine rclwrap_get_color_configurations (id, col) integer, intent(in) :: id integer, intent(out), dimension(:,:), allocatable :: col end subroutine rclwrap_get_color_configurations @ %def rclwrap_get_color_configurations @ <>= public :: rclwrap_use_dim_reg_soft <>= subroutine rclwrap_use_dim_reg_soft () end subroutine rclwrap_use_dim_reg_soft @ %def rclwrap_use_dim_reg_soft @ <>= public :: rclwrap_use_mass_reg_soft <>= subroutine rclwrap_use_mass_reg_soft (m) real(double), intent(in) :: m end subroutine rclwrap_use_mass_reg_soft @ %def rclwrap_use_mass_reg_soft @ <>= public :: rclwrap_set_delta_uv <>= subroutine rclwrap_set_delta_uv (d) real(double), intent(in) :: d end subroutine rclwrap_set_delta_uv @ %def rclwrap_set_delta_uv @ <>= public :: rclwrap_set_mu_uv <>= subroutine rclwrap_set_mu_uv (mu) real(double), intent(in) :: mu end subroutine rclwrap_set_mu_uv @ %def rclwrap_set_mu_uv @ <>= public :: rclwrap_set_delta_ir <>= subroutine rclwrap_set_delta_ir (d, d2) real(double), intent(in) :: d, d2 end subroutine rclwrap_set_delta_ir @ %def rclwrap_set_delta_ir @ <>= public :: rclwrap_set_mu_ir <>= subroutine rclwrap_set_mu_ir (mu) real(double), intent(in) :: mu end subroutine rclwrap_set_mu_ir @ %def rclwrap_set_mu_ir @ <>= public :: rclwrap_get_renormalization_scale <>= subroutine rclwrap_get_renormalization_scale (mu) real(double), intent(out) :: mu end subroutine rclwrap_get_renormalization_scale @ %def rclwrap_get_renormalization_scale @ <>= public :: rclwrap_get_flavor_scheme <>= subroutine rclwrap_get_flavor_scheme (nf) integer, intent(out) :: nf end subroutine rclwrap_get_flavor_scheme @ %def rclwrap_get_flavor_scheme @ <>= public :: rclwrap_use_alpha0_scheme <>= subroutine rclwrap_use_alpha0_scheme (al0) real(double), intent(in), optional :: al0 end subroutine rclwrap_use_alpha0_scheme @ %def rclwrap_use_alpha0_scheme @ <>= public :: rclwrap_use_alphaz_scheme <>= subroutine rclwrap_use_alphaz_scheme (alz) real(double), intent(in), optional :: alz end subroutine rclwrap_use_alphaz_scheme @ %def rclwrap_use_alphaz_scheme @ <>= public :: rclwrap_set_complex_mass_scheme <>= subroutine rclwrap_set_complex_mass_scheme () end subroutine rclwrap_set_complex_mass_scheme @ %def rclwrap_set_complex_mass_scheme @ <>= public :: rclwrap_set_resonant_particle <>= subroutine rclwrap_set_resonant_particle (pdg_id) integer, intent(in) :: pdg_id end subroutine rclwrap_set_resonant_particle @ %def rclwrap_set_resonant_particle @ <>= public :: rclwrap_switch_on_resonant_self_energies <>= subroutine rclwrap_switch_on_resonant_self_energies () end subroutine rclwrap_switch_on_resonant_self_energies @ %def rclwrap_switch_on_resonant_self_energies @ <>= public :: rclwrap_switch_off_resonant_self_energies <>= subroutine rclwrap_switch_off_resonant_self_energies () end subroutine rclwrap_switch_off_resonant_self_energies @ %def rclwrap_switch_off_resonant_self_energies @ <>= public :: rclwrap_set_draw_level_branches <>= subroutine rclwrap_set_draw_level_branches (n) integer, intent(in) :: n end subroutine rclwrap_set_draw_level_branches @ %def rclwrap_set_draw_level_branches @ <>= public :: rclwrap_set_print_level_amplitude <>= subroutine rclwrap_set_print_level_amplitude (n) integer, intent(in) :: n end subroutine rclwrap_set_print_level_amplitude @ %def rclwrap_set_print_level_amplitude @ <>= public :: rclwrap_set_print_level_squared_amplitude <>= subroutine rclwrap_set_print_level_squared_amplitude (n) integer, intent(in) :: n end subroutine rclwrap_set_print_level_squared_amplitude @ %def rclwrap_set_print_level_squared_amplitude @ <>= public :: rclwrap_set_print_level_correlations <>= subroutine rclwrap_set_print_level_correlations (n) integer, intent(in) :: n end subroutine rclwrap_set_print_level_correlations @ %def rclwrap_set_print_level_correlations @ <>= public :: rclwrap_set_print_level_RAM <>= subroutine rclwrap_set_print_level_RAM (n) integer, intent(in) :: n end subroutine rclwrap_set_print_level_RAM @ %def rclwrap_set_print_level_RAM @ <>= public :: rclwrap_scale_coupling3 <>= subroutine rclwrap_scale_coupling3 (pdg_id1, pdg_id2, pdg_id3, factor) integer, intent(in) :: pdg_id1, pdg_id2, pdg_id3 complex(double), intent(in) :: factor end subroutine rclwrap_scale_coupling3 @ %def rclwrap_scale_coupling3 @ <>= public :: rclwrap_scale_coupling4 <>= subroutine rclwrap_scale_coupling4 (pdg_id1, pdg_id2, pdg_id3, pdg_id4, factor) integer, intent(in) :: pdg_id1, pdg_id2, pdg_id3, pdg_id4 complex(double), intent(in) :: factor end subroutine rclwrap_scale_coupling4 @ %def rclwrap_scale_coupling4 @ <>= public :: rclwrap_switch_off_coupling3 <>= subroutine rclwrap_switch_off_coupling3 (pdg_id1, pdg_id2, pdg_id3) integer, intent(in) :: pdg_id1, pdg_id2, pdg_id3 end subroutine rclwrap_switch_off_coupling3 @ %def rclwrap_switch_off_coupling3 @ <>= public :: rclwrap_switch_off_coupling4 <>= subroutine rclwrap_switch_off_coupling4 (pdg_id1, pdg_id2, pdg_id3, pdg_id4) integer, intent(in) :: pdg_id1, pdg_id2, pdg_id3, pdg_id4 end subroutine rclwrap_switch_off_coupling4 @ %def rclwrap_switch_off_coupling4 @ <>= public :: rclwrap_set_ifail <>= subroutine rclwrap_set_ifail (i) integer, intent(in) :: i end subroutine rclwrap_set_ifail @ %def rclwrap_set_ifail @ <>= public :: rclwrap_get_ifail <>= subroutine rclwrap_get_ifail (i) integer, intent(out) :: i end subroutine rclwrap_get_ifail @ %def rclwrap_get_ifail @ <>= public :: rclwrap_set_output_file <>= subroutine rclwrap_set_output_file (filename) character(len=*), intent(in) :: filename end subroutine rclwrap_set_output_file @ %def rclwrap_set_output_file @ <>= public :: rclwrap_set_gs_power <>= subroutine rclwrap_set_gs_power (id, gs_array) integer, intent(in) :: id integer, dimension(:,:), intent(in) :: gs_array end subroutine rclwrap_set_gs_power @ %def rclwrap_set_gs_power @ <>= public :: rclwrap_select_gs_power_born_amp <>= subroutine rclwrap_select_gs_power_born_amp (id, gs_power) integer, intent(in) :: id, gs_power end subroutine rclwrap_select_gs_power_born_amp @ %def rclwrap_select_gs_power_born_amp @ <>= public :: rclwrap_unselect_gs_power_born_amp <>= subroutine rclwrap_unselect_gs_power_born_amp (id, gs_power) integer, intent(in) :: id, gs_power end subroutine rclwrap_unselect_gs_power_born_amp @ %def rclwrap_unselect_gs_power_born_amp @ <>= public :: rclwrap_select_gs_power_loop_amp <>= subroutine rclwrap_select_gs_power_loop_amp (id, gs_power) integer, intent(in) :: id, gs_power end subroutine rclwrap_select_gs_power_loop_amp @ %def rclwrap_select_gs_power_loop_amp @ <>= public :: rclwrap_unselect_gs_power_loop_amp <>= subroutine rclwrap_unselect_gs_power_loop_amp (id, gs_power) integer, intent(in) :: id, gs_power end subroutine rclwrap_unselect_gs_power_loop_amp @ %def rclwrap_unselect_gs_power_loop_amp @ <>= public :: rclwrap_select_all_gs_powers_born_amp <>= subroutine rclwrap_select_all_gs_powers_born_amp (id) integer, intent(in) :: id end subroutine rclwrap_select_all_gs_powers_born_amp @ %def rclwrap_select_all_gs_powers_born_amp @ <>= public :: rclwrap_unselect_all_gs_powers_loop_amp <>= subroutine rclwrap_unselect_all_gs_powers_loop_amp (id) integer, intent(in) :: id end subroutine rclwrap_unselect_all_gs_powers_loop_amp @ %def rclwrap_unselect_all_gs_powers_loop_amp @ <>= public :: rclwrap_select_all_gs_powers_loop_amp <>= subroutine rclwrap_select_all_gs_powers_loop_amp (id) integer, intent(in) :: id end subroutine rclwrap_select_all_gs_powers_loop_amp @ %def rclwrap_select_all_gs_powers_loop_amp @ <>= public :: rclwrap_unselect_all_gs_powers_born_amp <>= subroutine rclwrap_unselect_all_gs_powers_born_amp (id) integer, intent(in) :: id end subroutine rclwrap_unselect_all_gs_powers_born_amp @ %def rclwrap_unselect_all_gs_powers_born_amp @ <>= public :: rclwrap_set_resonant_squared_momentum <>= subroutine rclwrap_set_resonant_squared_momentum (id, i_res, p2) integer, intent(in) :: id, i_res real(double), intent(in) :: p2 end subroutine rclwrap_set_resonant_squared_momentum @ %def rclwrap_set_resonant_squared_momentum @ <>= public :: rclwrap_compute_running_alpha_s <>= subroutine rclwrap_compute_running_alpha_s (Q, nf, n_loops) real(double), intent(in) :: Q integer, intent(in) :: nf, n_loops end subroutine rclwrap_compute_running_alpha_s @ %def rclwrap_compute_running_alpha_s @ <>= public :: rclwrap_set_dynamic_settings <>= subroutine rclwrap_set_dynamic_settings () end subroutine rclwrap_set_dynamic_settings @ %def rclwrap_set_dynamic_settings @ <>= public :: rclwrap_rescale_process <>= subroutine rclwrap_rescale_process (id, order, sqme) integer, intent(in) :: id character(len=*), intent(in) :: order real(double), dimension(0:1), intent(out), optional :: sqme end subroutine rclwrap_rescale_process @ %def rclwrap_rescale_process @ <>= public :: rclwrap_get_polarized_squared_amplitude <>= subroutine rclwrap_get_polarized_squared_amplitude (id, & alphas_power, order, hel, sqme) integer, intent(in) :: id, alphas_power character(len=*), intent(in) :: order integer, dimension(:), intent(in) :: hel real(double), intent(out) :: sqme end subroutine rclwrap_get_polarized_squared_amplitude @ %def rclwrap_get_polarized_squared_amplitude @ <>= public :: rclwrap_compute_color_correlation <>= subroutine rclwrap_compute_color_correlation (id, p, & i1, i2, sqme) integer, intent(in) :: id real(double), dimension(:,:), intent(in) :: p integer, intent(in) :: i1, i2 real(double), intent(out), optional :: sqme end subroutine rclwrap_compute_color_correlation @ %def rclwrap_compute_color_correlation @ <>= public :: rclwrap_compute_all_color_correlations <>= subroutine rclwrap_compute_all_color_correlations (id, p) integer, intent(in) :: id real(double), dimension(:,:), intent(in) :: p end subroutine rclwrap_compute_all_color_correlations @ %def rclwrap_compute_all_color_correlations @ <>= public :: rclwrap_rescale_color_correlation <>= subroutine rclwrap_rescale_color_correlation (id, i1, i2, sqme) integer, intent(in) :: id, i1, i2 real(double), intent(out), optional :: sqme end subroutine rclwrap_rescale_color_correlation @ %def rclwrap_rescale_color_correlation @ <>= public :: rclwrap_rescale_all_color_correlations <>= subroutine rclwrap_rescale_all_color_correlations (id) integer, intent(in) :: id end subroutine rclwrap_rescale_all_color_correlations @ %def rclwrap_rescale_all_color_correlations @ <>= public :: rclwrap_get_color_correlation <>= subroutine rclwrap_get_color_correlation (id, alphas_power, i1, i2, sqme) integer, intent(in) :: id, alphas_power, i1, i2 real(double), intent(out) :: sqme end subroutine rclwrap_get_color_correlation @ %def rclwrap_get_color_correlation @ <>= public :: rclwrap_compute_spin_correlation <>= subroutine rclwrap_compute_spin_correlation (id, p, i_photon, pol, sqme) integer, intent(in) :: id real(double), dimension(:,:), intent(in) :: p integer, intent(in) :: i_photon complex(double), dimension(:), intent(in) :: pol real(double), intent(out), optional :: sqme end subroutine rclwrap_compute_spin_correlation @ %def rclwrap_compute_spin_correlation @ <>= public :: rclwrap_rescale_spin_correlation <>= subroutine rclwrap_rescale_spin_correlation (id, i_photon, pol, sqme) integer, intent(in) :: id, i_photon complex(double), dimension(:), intent(in) :: pol real(double), intent(out), optional :: sqme end subroutine rclwrap_rescale_spin_correlation @ %def rclwrap_rescale_spin_correlation @ <>= public :: rclwrap_get_spin_correlation <>= subroutine rclwrap_get_spin_correlation (id, alphas_power, sqme) integer, intent(in) :: id, alphas_power real(double), intent(out) :: sqme end subroutine rclwrap_get_spin_correlation @ %def rclwrap_get_spin_correlation @ <>= public :: rclwrap_compute_spin_color_correlation <>= subroutine rclwrap_compute_spin_color_correlation (id, p, & i_gluon, i_spectator, pol, sqme) integer, intent(in) :: id real(double), dimension(:,:), intent(in) :: p integer, intent(in) :: i_gluon, i_spectator complex(double), dimension(:), intent(in) :: pol real(double), intent(out), optional :: sqme end subroutine rclwrap_compute_spin_color_correlation @ %def rclwrap_compute_spin_color_correlation @ <>= public :: rclwrap_rescale_spin_color_correlation <>= subroutine rclwrap_rescale_spin_color_correlation (id, i_gluon, & i_spectator, pol, sqme) integer, intent(in) :: id, i_gluon, i_spectator complex(double), dimension(:), intent(in) :: pol real(double), intent(out), optional :: sqme end subroutine rclwrap_rescale_spin_color_correlation @ %def rclwrap_rescale_spin_color_correlation @ <>= public :: rclwrap_get_spin_color_correlation <>= subroutine rclwrap_get_spin_color_correlation (id, alphas_power, & i_gluon, i_spectator, sqme) integer, intent(in) :: id, alphas_power, i_gluon, i_spectator real(double), intent(out) :: sqme end subroutine rclwrap_get_spin_color_correlation @ %def rclwrap_get_spin_color_correlation @ <>= public :: rclwrap_get_momenta <>= subroutine rclwrap_get_momenta (id, p) integer, intent(in) :: id real(double), dimension(:,:), intent(out) :: p end subroutine rclwrap_get_momenta @ %def rclwrap_get_momenta @ <>= public :: rclwrap_reset_recola <>= subroutine rclwrap_reset_recola end subroutine rclwrap_reset_recola @ %def rclwrap_reset_recola @ \clearpage %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Recola Core} The recola core object and auxiliary types and objects. <<[[prc_recola.f90]]>>= <> module prc_recola use recola_wrapper !NODEP! use kinds use constants, only: pi, zero <> <> use string_utils, only: str use system_defs, only: TAB use diagnostics use io_units use lorentz use physics_defs use variables, only: var_list_t use os_interface, only: os_data_t use sm_qcd, only: qcd_t use model_data, only: model_data_t use prc_core, only: prc_core_state_t use prc_core_def, only: prc_core_driver_t, prc_core_def_t use prc_external use process_libraries, only: process_library_t <> <> <> <> <> contains <> end module prc_recola @ %def prc_recola @ \subsection{Sanity check} Checks the [[rclwrap_is_active]] flag and aborts the program if the dummy is used. <>= public :: abort_if_recola_not_active <>= subroutine abort_if_recola_not_active () if (.not. rclwrap_is_active) call msg_fatal ("You want to use Recola, ", & [var_str("but either the compiler with which Whizard has been build "), & var_str("is not supported by it, or you have not linked Recola "), & var_str("correctly to Whizard. Either reconfigure Whizard with a path to "), & var_str("a valid Recola installation (for details consult the manual), "), & var_str("or choose a different matrix-element method.")]) end subroutine abort_if_recola_not_active @ %def abort_if_recola_not_active @ \subsection{Process definition} When defining a RECOLA process, we store the process-specific flags and parameters. Correction types are either QCD, QED, or full SM. <>= integer, parameter :: RECOLA_UNDEFINED = 0, RECOLA_QCD = 1, & RECOLA_QED = 2, RECOLA_FULL = 3 @ %def RECOLA_QCD RECOLA_QED RECOLA_FULL @ <>= public :: recola_def_t <>= type, extends (prc_external_def_t) :: recola_def_t type(string_t) :: suffix type(string_t) :: order integer :: alpha_power = 0 integer :: alphas_power = 0 integer :: corr = RECOLA_UNDEFINED contains <> end type recola_def_t @ %def recola_def_t @ <>= procedure, nopass :: type_string => recola_def_type_string <>= function recola_def_type_string () result (string) type(string_t) :: string string = "recola" end function recola_def_type_string @ %def recola_def_type_string @ Not implemented yet. <>= procedure :: write => recola_def_write <>= subroutine recola_def_write (object, unit) class(recola_def_t), intent(in) :: object integer, intent(in) :: unit end subroutine recola_def_write @ %def recola_def_write @ <>= procedure :: read => recola_def_read <>= subroutine recola_def_read (object, unit) class(recola_def_t), intent(out) :: object integer, intent(in) :: unit end subroutine recola_def_read @ %def recola_def_read @ The initializer has the responsibility to store all process- and method-specific parameters, such that they can be used later by the writer and by the driver for this process. Also, it allocates the writer. For RECOLA, the writer (i) creates full-fledged \oMega\ matrix element code which we need for the interface. (ii) registers the process definition with the RECOLA library which has been linked. The latter task does not involve external code. Note that all management stuff is taken care of by the base type(s) methods. Here, we introduce only RECOLA-specific procedures, in addition. The NLO flag is true only for virtual matrix elements. <>= procedure :: init => recola_def_init <>= subroutine recola_def_init (object, basename, model_name, & prt_in, prt_out, nlo_type, alpha_power, alphas_power, & - correction_type) + correction_type, restrictions) class(recola_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 integer, intent(in) :: alpha_power integer, intent(in) :: alphas_power type(string_t), intent(in) :: correction_type + type(string_t), intent(in), optional :: restrictions if (debug_on) call msg_debug (D_ME_METHODS, "recola_def_init: " & // char (basename) // ", nlo_type", nlo_type) object%basename = basename object%alpha_power = alpha_power object%alphas_power = alphas_power select case (char (correction_type)) case ("QCD") object%corr = RECOLA_QCD case ("QED") object%corr = RECOLA_QED case ("Full") object%corr = RECOLA_FULL end select allocate (recola_writer_t :: object%writer) select case (nlo_type) case (BORN) object%suffix = '_BORN' object%order = "LO" case (NLO_REAL) object%suffix = '_REAL' object%order = "LO" if (object%corr == RECOLA_QCD) object%alphas_power = alphas_power + 1 if (object%corr == RECOLA_QED) object%alpha_power = alpha_power + 1 case (NLO_VIRTUAL) object%suffix = '_LOOP' object%order = "NLO" case (NLO_SUBTRACTION) object%suffix = '_SUB' object%order = "LO" case (NLO_MISMATCH) object%suffix = '_MISMATCH' object%order = "LO" case (NLO_DGLAP) object%suffix = '_DGLAP' object%order = "LO" end select select type (writer => object%writer) class is (recola_writer_t) - call writer%init (model_name, prt_in, prt_out) + call writer%init (model_name, prt_in, prt_out, restrictions) call writer%set_id (basename // object%suffix) call writer%set_order (object%order) call writer%set_coupling_powers (object%alpha_power, object%alphas_power) end select end subroutine recola_def_init @ %def recola_def_init @ \subsection{Writer object} The RECOLA writer takes the additional resposibility of transferring process information to RECOLA. <>= type, extends (prc_external_writer_t) :: recola_writer_t private type(string_t) :: id type(string_t) :: order integer :: alpha_power = 0 integer :: alphas_power = 0 contains <> end type recola_writer_t @ %def recola_writer_t @ <>= procedure, nopass :: type_name => recola_writer_type_name <>= function recola_writer_type_name () result (string) type(string_t) :: string string = "recola" end function recola_writer_type_name @ %def recola_writer_type_name @ Set the process ID string as used by WHIZARD. <>= procedure :: set_id => recola_writer_set_id <>= subroutine recola_writer_set_id (writer, id) class(recola_writer_t), intent(inout) :: writer type(string_t), intent(in) :: id if (debug_on) call msg_debug2 (D_ME_METHODS, "Recola writer: id = " // char (id)) writer%id = id end subroutine recola_writer_set_id @ %def recola_writer_set_id @ Set the NLO flag. <>= procedure :: set_order => recola_writer_set_order <>= subroutine recola_writer_set_order (writer, order) class(recola_writer_t), intent(inout) :: writer type(string_t), intent(in) :: order if (debug_on) call msg_debug2 (D_ME_METHODS, "Recola writer: order = " // char (order)) writer%order = order end subroutine recola_writer_set_order @ %def recola_writer_set_order @ Set coupling powers. <>= procedure :: set_coupling_powers => recola_writer_set_coupling_powers <>= subroutine recola_writer_set_coupling_powers (writer, alpha_power, alphas_power) class(recola_writer_t), intent(inout) :: writer integer, intent(in) :: alpha_power integer, intent(in) :: alphas_power if (debug_on) call msg_debug2 (D_ME_METHODS, "Recola writer: alphas_power", alphas_power) if (debug_on) call msg_debug2 (D_ME_METHODS, "Recola writer: alpha_power", alpha_power) writer%alpha_power = alpha_power writer%alphas_power = alphas_power end subroutine recola_writer_set_coupling_powers @ %def recola_writer_set_coupling_powers @ The Makefile code contains all of the code that the [[prc_external]] base method generates, plus an extra clause that extracts a shorthand listing of all flavor combinations for the current process. This list is required by [[make source]], so it can be read and used for declaring the RECOLA processes. There is one glitch here: we use the component-specific source file but write a flavor list for the process, without component extension. That is, we must not have more than one component at this stage. NB: We might actually extend \oMega\ to produce this shorthand listing. <>= procedure :: write_makefile_code => recola_writer_write_makefile_code <>= function flv_file_name (id) type(string_t), intent(in) :: id type(string_t) :: flv_file_name flv_file_name = id // ".flv.dat" end function flv_file_name subroutine recola_writer_write_makefile_code & (writer, unit, id, os_data, verbose, testflag) class(recola_writer_t), intent(in) :: writer integer, intent(in) :: unit type(string_t), intent(in) :: id type(os_data_t), intent(in) :: os_data logical, intent(in) :: verbose logical, intent(in), optional :: testflag type(string_t) :: src_file type(string_t) :: flv_file call writer%base_write_makefile_code (unit, id, os_data, verbose, testflag) src_file = trim (char(id)) // ".f90" flv_file = flv_file_name (writer%id) write (unit, *) write (unit, "(5A)") "# Flavor state listing for RECOLA process generation" write (unit, "(5A)") char (flv_file), ": ", char (src_file) if (verbose) then write (unit, "(5A)", advance="no") TAB else write (unit, "(5A)") TAB, '@echo " MAKE ', char (flv_file), '"' write (unit, "(5A)", advance="no") TAB, "@" end if write (unit, "(5A)") & "grep 'data table_flavor_states' $< ", & "| sed -e 's/.*\/\(.*\)\/.*/\1/' -e 's/,//g' > $@" write (unit, "(5A)") "SOURCES += ", char (flv_file) write (unit, "(5A)") "CLEAN_SOURCES += ", char (flv_file) end subroutine recola_writer_write_makefile_code @ %def recola_writer_write_makefile_code @ To communicate the process definition to RECOLA, we must know the following: the process definition, expanded in terms of flavor states, and the process order (LO/NLO). We will ask for a new numeric ID, create a process string using RECOLA conventions, and define the process. The [[request_generate_processes]] enables the RECOLA internal process compiler, which can be called only after all processes have been defined. <>= procedure :: register_processes => prc_recola_register_processes <>= subroutine prc_recola_register_processes (writer, recola_ids) class(recola_writer_t), intent(in) :: writer integer :: recola_id integer :: i_flv integer :: n_tot integer :: unit, iostat integer, dimension (:), intent(inout) :: recola_ids integer, dimension(:), allocatable :: pdg type(string_t), dimension(:), allocatable :: particle_names type(string_t) :: process_string integer :: i_part !!! TODO (cw-2016-08-08): Include helicities call msg_message ("Recola: registering processes for '" // char (writer%id) // "'") i_flv = 0 n_tot = writer%n_in + writer%n_out allocate (pdg (n_tot)) allocate (particle_names (n_tot)) call open_flv_list (writer%id, unit) call rclwrap_request_generate_processes () SCAN_FLV_LIST: do read (unit, *, iostat = iostat) pdg if (iostat < 0) then exit SCAN_FLV_LIST else if (iostat > 0) then call err_flv_list (writer%id) end if i_flv = i_flv + 1 call rclwrap_get_new_recola_id (recola_id) recola_ids(i_flv) = recola_id particle_names(:) = get_recola_particle_string (pdg) process_string = var_str ("") do i_part = 1, n_tot process_string = process_string // & particle_names (i_part) // var_str (" ") if (i_part == writer%n_in) then process_string = process_string // var_str ("-> ") end if end do call msg_message ("Recola: " & // "process #" // char (str (recola_id)) & // ": " // char (process_string) & // "(" // char (writer%order) // ")") call rclwrap_add_process (recola_id, process_string, writer%order) call rclwrap_define_processes () end do SCAN_FLV_LIST call close_flv_list (unit) if (debug_on) call msg_debug (D_ME_METHODS, "RECOLA: processes for '" & // char (writer%id) // "' registered") end subroutine prc_recola_register_processes @ %def prc_recola_register_processes @ Manage the list of flavor combinations for the current process. We rely on this being created along with the \oMega\ call. <>= subroutine open_flv_list (id, unit) type(string_t), intent(in) :: id integer, intent(out) :: unit type(string_t) :: flv_file integer :: iostat flv_file = flv_file_name (id) open (file = char (flv_file), newunit = unit, & status = "old", action = "read", & iostat = iostat) if (iostat /= 0) then call msg_fatal ("Recola: attempt to open flavor-list file '" & // char (flv_file) // "' failed") end if end subroutine open_flv_list subroutine err_flv_list (id) type(string_t), intent(in) :: id type(string_t) :: flv_file flv_file = flv_file_name (id) call msg_fatal ("Recola: error while reading from flavor-list file '" & // char (flv_file) // "'") end subroutine err_flv_list subroutine close_flv_list (unit) integer, intent(in) :: unit close (unit) end subroutine close_flv_list @ %def open_flv_list @ %def err_flv_list @ %def close_flv_list @ \subsection{Driver object} A core driver is required by design. However, we are not going to load any external dynamical libraries, so this is a dummy. <>= type, extends (prc_external_driver_t) :: recola_driver_t contains <> end type recola_driver_t @ %def recola_driver_t @ <>= procedure :: allocate_driver => recola_def_allocate_driver <>= subroutine recola_def_allocate_driver (object, driver, basename) class(recola_def_t), intent(in) :: object class(prc_core_driver_t), intent(out), allocatable :: driver type(string_t), intent(in) :: basename if (debug_on) call msg_debug2 (D_ME_METHODS, "recola_def_allocate_driver") allocate (recola_driver_t :: driver) end subroutine recola_def_allocate_driver @ %def recola_def_allocate_driver @ <>= procedure, nopass :: type_name => recola_driver_type_name <>= function recola_driver_type_name () result (type) type(string_t) :: type type = "Recola" end function recola_driver_type_name @ %def recola_driver_type_name @ \subsection{Process object} We create [[prc_recola_t]] as an extension of the [[prc_external_t]], which in turn inherits from [[prc_core_t]]. This way, we can use a lot of the existing interfaces in the actual code. However, we have to stick to the rules and implement the deferred type-bound procedures of [[prc_core_t]]. <>= public :: prc_recola_t <>= type, extends (prc_external_t) :: prc_recola_t integer, dimension(:), allocatable :: recola_ids integer, dimension(:,:), allocatable :: color_state integer :: n_f = 0 logical :: helicity_and_color_arrays_are_replaced = .false. contains <> end type prc_recola_t @ %def prc_recola_t @ <>= procedure :: write_name => prc_recola_write_name <>= subroutine prc_recola_write_name (object, unit) class(prc_recola_t), intent(in) :: object integer, intent(in), optional :: unit integer :: u u = given_output_unit (unit) write (u,"(1x,A)") "Core: Recola" end subroutine prc_recola_write_name @ %def prc_recola_write_name @ <>= procedure :: has_matrix_element => prc_recola_has_matrix_element <>= function prc_recola_has_matrix_element (object) result (flag) logical :: flag class(prc_recola_t), intent(in) :: object flag = .true. end function prc_recola_has_matrix_element @ %def prc_recola_has_matrix_element @ Not implemented yet. <>= procedure :: write => prc_recola_write <>= subroutine prc_recola_write (object, unit) class(prc_recola_t), intent(in) :: object integer, intent(in), optional :: unit end subroutine prc_recola_write @ %def prc_recola_write @ \subsection{Accompanying state object} This must be implemented, but is unused. <>= type, extends (prc_external_state_t) :: recola_state_t contains <> end type recola_state_t @ %def recola_state_t @ <>= procedure :: write => recola_state_write <>= subroutine recola_state_write (object, unit) class(recola_state_t), intent(in) :: object integer, intent(in), optional :: unit end subroutine recola_state_write @ %def recola_state_write @ <>= procedure :: allocate_workspace => prc_recola_allocate_workspace <>= subroutine prc_recola_allocate_workspace (object, core_state) class(prc_recola_t), intent(in) :: object class(prc_core_state_t), intent(inout), allocatable :: core_state allocate (recola_state_t :: core_state) end subroutine prc_recola_allocate_workspace @ %def prc_recola_allocate_workspace @ \subsection{Recola process data} This information is stored in the associated [[def]] object. To obtain it, we need a type cast. <>= procedure :: get_alpha_power => prc_recola_get_alpha_power procedure :: get_alphas_power => prc_recola_get_alphas_power <>= function prc_recola_get_alpha_power (object) result (p) class(prc_recola_t), intent(in) :: object integer :: p p = 0 if (associated (object%def)) then select type (def => object%def) type is (recola_def_t) p = def%alpha_power end select end if end function prc_recola_get_alpha_power function prc_recola_get_alphas_power (object) result (p) class(prc_recola_t), intent(in) :: object integer :: p p = 0 if (associated (object%def)) then select type (def => object%def) type is (recola_def_t) p = def%alphas_power end select end if end function prc_recola_get_alphas_power @ %def prc_recola_get_alpha_power @ %def prc_recola_get_alphas_power @ <>= procedure :: compute_alpha_s => prc_recola_compute_alpha_s <>= subroutine prc_recola_compute_alpha_s (object, core_state, ren_scale) class(prc_recola_t), intent(in) :: object class(prc_external_state_t), intent(inout) :: core_state real(default), intent(in) :: ren_scale core_state%alpha_qcd = object%qcd%alpha%get (ren_scale) end subroutine prc_recola_compute_alpha_s @ %def prc_recola_compute_alpha_s @ <>= procedure :: includes_polarization => prc_recola_includes_polarization <>= function prc_recola_includes_polarization (object) result (polarized) logical :: polarized class(prc_recola_t), intent(in) :: object polarized = .false. end function prc_recola_includes_polarization @ %def prc_recola_includes_polarization @ \subsection{Prepare for process evaluation} This has become obsolete and is empty. <>= procedure :: prepare_external_code => & prc_recola_prepare_external_code <>= subroutine prc_recola_prepare_external_code & (core, flv_states, var_list, os_data, libname, model, i_core, is_nlo) class(prc_recola_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 if (debug_on) call msg_debug (D_ME_METHODS, "prc_recola_prepare_external_code (no-op)") end subroutine prc_recola_prepare_external_code @ %def prc_recola_prepare_external_code @ Set all Recola parameters to their correct values. We use the model object for masses and such. Note that the QCD object provides the [[n_f]] parameter which affects $\alpha_s$ evaluation. Note that this is executed before the [[init]] method below, which defines and prepares the Recola process objects. This is in line with the Recola workflow, however. <>= procedure :: set_parameters => prc_recola_set_parameters <>= subroutine prc_recola_set_parameters (object, qcd, model) class(prc_recola_t), intent(inout) :: object type(qcd_t), intent(in) :: qcd class(model_data_t), intent(in), target, optional :: model if (debug_on) call msg_debug (D_ME_METHODS, "RECOLA: set_parameters") object%qcd = qcd call rclwrap_set_dynamic_settings () call rclwrap_set_pole_mass & (11, dble(model%get_real (var_str ('me'))), 0._double) call rclwrap_set_pole_mass & (13, dble(model%get_real (var_str ('mmu'))), 0._double) call rclwrap_set_pole_mass & (15, dble(model%get_real (var_str ('mtau'))), 0._double) call rclwrap_set_pole_mass (1, 0._double, 0._double) call rclwrap_set_pole_mass (2, 0._double, 0._double) call rclwrap_set_pole_mass (3, dble(model%get_real (var_str ('ms'))), 0._double) call rclwrap_set_pole_mass (4, dble(model%get_real (var_str ('mc'))), 0._double) call rclwrap_set_pole_mass (5, dble(model%get_real (var_str ('mb'))), 0._double) call rclwrap_set_pole_mass (6, dble(model%get_real (var_str ('mtop'))), & dble(model%get_real (var_str ('wtop')))) call rclwrap_set_pole_mass (23, dble(model%get_real (var_str ('mZ'))), & dble(model%get_real (var_str ('wZ')))) call rclwrap_set_pole_mass (24, dble(model%get_real (var_str ('mW'))), & dble(model%get_real (var_str ('wW')))) call rclwrap_set_pole_mass (25, dble(model%get_real (var_str ('mH'))), & dble(model%get_real (var_str ('wH')))) call rclwrap_use_gfermi_scheme (dble(model%get_real (var_str ('GF')))) call rclwrap_set_light_fermions (0._double) call rclwrap_set_delta_ir (0._double, dble(pi**2 / 6)) end subroutine prc_recola_set_parameters @ %def prc_recola_set_parameters @ <>= procedure :: set_mu_ir => prc_recola_set_mu_ir <>= subroutine prc_recola_set_mu_ir (object, mu) class(prc_recola_t), intent(inout) :: object real(default), intent(in) :: mu call rclwrap_set_mu_ir (dble(mu)) end subroutine prc_recola_set_mu_ir @ %def prc_recola_set_mu_ir @ Extend the base-type initialization method by Recola-specific initialization. We take the process definitions from the [[def]] object, which has been filled before. The [[writer]] component of the process-definition object can now complete its task and prepare the Recola processes. Sadly, we have to completely reset Recola first, since Recola does not allow to modify \emph{anything} after process definition. Also, we cannot really make use of Recola's multi-process capability without violating the Whizard convention that the parameter settings at process integration time apply, not at process definition time. Each new process (i.e., process-integration) object will thus trigger a complete new Recola instance. <>= procedure :: init => prc_recola_init <>= subroutine prc_recola_init (object, def, lib, id, i_component) class(prc_recola_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 integer :: n_flv if (debug_on) call msg_debug (D_ME_METHODS, "RECOLA: init process object") call object%base_init (def, lib, id, i_component) n_flv = object%get_n_flvs (1) allocate (object%recola_ids(n_flv)) select type (writer => object%def%writer) type is (recola_writer_t) call writer%register_processes (object%recola_ids) end select call rclwrap_generate_processes () call object%replace_helicity_and_color_arrays () end subroutine prc_recola_init @ %def prc_recola_init @ Recola can compute dressed amplitudes, but it needs helicity and color to be in its own format to do so. <>= procedure :: replace_helicity_and_color_arrays => & prc_recola_replace_helicity_and_color_arrays <>= subroutine prc_recola_replace_helicity_and_color_arrays (object) class(prc_recola_t), intent(inout) :: object integer, dimension(:,:), allocatable :: col_recola integer :: i if (debug_on) call msg_debug (D_ME_METHODS, "RECOLA: replace_helicity_and_color_arrays") deallocate (object%data%hel_state) call rclwrap_get_helicity_configurations & (object%recola_ids(1), object%data%hel_state) call rclwrap_get_color_configurations (object%recola_ids(1), col_recola) allocate (object%color_state (object%data%n_in + object%data%n_out, & size (col_recola, dim = 2))) do i = 1, size (col_recola, dim = 2) object%color_state (:, i) = col_recola (:, i) end do end subroutine prc_recola_replace_helicity_and_color_arrays @ %def prc_recola_replace_helicity_and_color_arrays @ \subsection{Compute matrix element} Computes the amplitude as a function of the phase space point, the flavor, helicity and color index. It is currently only used in the form by [[prc_omega_t]], all the other ones use different interfaces. H With RECOLA, we might be able to use this, too. The current implementation can fail due to missing helicity initialization. <>= procedure :: compute_amplitude => prc_recola_compute_amplitude <>= function prc_recola_compute_amplitude & (object, j, p, f, h, c, fac_scale, ren_scale, alpha_qcd_forced, & core_state) result (amp) complex(default) :: amp class(prc_recola_t), intent(in) :: object integer, intent(in) :: j type(vector4_t), intent(in), dimension(:) :: 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 real(double), dimension(0:3, object%data%n_in + object%data%n_out) :: & p_recola integer :: i logical :: new_event complex(double) :: amp_dble if (debug_on) call msg_debug2 (D_ME_METHODS, "prc_recola_compute_amplitude") if (present (core_state)) then if (allocated (core_state)) then select type (core_state) type is (recola_state_t) new_event = core_state%new_kinematics core_state%new_kinematics = .false. end select end if end if if (new_event) then do i = 1, object%data%n_in + object%data%n_out p_recola(:, i) = dble(p(i)%p) end do call rclwrap_compute_process (object%recola_ids(f), p_recola, 'LO') end if call rclwrap_get_amplitude (object%recola_ids(f), 0, 'LO', & object%color_state (:, c), object%data%hel_state (h, :), amp_dble) amp = amp_dble end function prc_recola_compute_amplitude @ %def prc_recola_compute_amplitude @ <>= procedure :: compute_sqme => prc_recola_compute_sqme <>= subroutine prc_recola_compute_sqme (object, i_flv, i_hel, p, & ren_scale, sqme, bad_point) class(prc_recola_t), intent(in) :: object integer, intent(in) :: i_flv, i_hel type(vector4_t), dimension(:), intent(in) :: p real(default), intent(in) :: ren_scale real(default), intent(out) :: sqme logical, intent(out) :: bad_point real(double) :: sqme_dble real(double), dimension(0:3, object%data%n_in + object%data%n_out) :: & p_recola real(default) :: alpha_s integer :: i integer :: alphas_power ! TODO sbrass: Helicity for RECOLA if (debug_on) call msg_debug2 (D_ME_METHODS, "prc_recola_compute_sqme") do i = 1, object%data%n_in + object%data%n_out p_recola(:, i) = dble(p(i)%p) end do alpha_s = object%qcd%alpha%get (ren_scale) if (debug_on) call msg_debug2 (D_ME_METHODS, "alpha_s", alpha_s) if (debug_on) call msg_debug2 (D_ME_METHODS, "ren_scale", ren_scale) call rclwrap_set_alpha_s (dble (alpha_s), dble (ren_scale), object%qcd%n_f) call rclwrap_set_mu_ir (dble (ren_scale)) call rclwrap_compute_process (object%recola_ids(i_flv), p_recola, 'LO') call rclwrap_get_squared_amplitude & (object%recola_ids(i_flv), object%get_alphas_power (), 'LO', sqme_dble) sqme = real(sqme_dble, kind=default) bad_point = .false. end subroutine prc_recola_compute_sqme @ %def prc_recola_compute_sqme @ <>= procedure :: compute_sqme_virt => prc_recola_compute_sqme_virt <>= subroutine prc_recola_compute_sqme_virt (object, i_flv, i_hel, & p, ren_scale, sqme, bad_point) class(prc_recola_t), intent(in) :: object integer, intent(in) :: i_flv, i_hel type(vector4_t), dimension(:), intent(in) :: p real(default), intent(in) :: ren_scale real(default), dimension(4), intent(out) :: sqme real(default) :: amp logical, intent(out) :: bad_point real(double), dimension(0:3, object%data%n_in + object%data%n_out) :: & p_recola real(double) :: sqme_dble real(default) :: alpha_s integer :: i ! TODO sbrass Helicity for RECOLA if (debug_on) call msg_debug2 (D_ME_METHODS, "prc_recola_compute_sqme_virt") sqme = zero do i = 1, object%data%n_in + object%data%n_out p_recola(:, i) = dble(p(i)%p) end do call rclwrap_set_mu_ir (dble (ren_scale)) alpha_s = object%qcd%alpha%get (ren_scale) call rclwrap_set_alpha_s (dble (alpha_s), dble (ren_scale), object%qcd%n_f) call rclwrap_compute_process (object%recola_ids(i_flv), p_recola, 'NLO') !!! JRR, TODO: generalize for QED corrections call rclwrap_get_squared_amplitude & (object%recola_ids(i_flv), object%get_alphas_power () + 1, 'NLO', sqme_dble) sqme(3) = sqme_dble call rclwrap_get_squared_amplitude & (object%recola_ids(i_flv), object%get_alphas_power (), 'LO', sqme_dble) sqme(4) = sqme_dble bad_point = .false. end subroutine prc_recola_compute_sqme_virt @ %def prc_recola_compute_sqme_virt @ For RECOLA, explicit color factors need to multiplied to the off-diagonal elements of the color correlation matrix. The factor 1/2 from the normalization accoring to the RECOLA manual is covered by the fact that we are taking only one half of the symmetric matrix. <>= procedure :: compute_sqme_color_c_raw => prc_recola_compute_sqme_color_c_raw <>= subroutine prc_recola_compute_sqme_color_c_raw (object, i_flv, i_hel, & p, ren_scale, sqme_color_c, bad_point) class(prc_recola_t), intent(in) :: object integer, intent(in) :: i_hel, i_flv type(vector4_t), dimension(:), intent(in) :: p real(double), dimension(0:3, object%data%n_in + object%data%n_out) :: & p_recola real(default), intent(in) :: ren_scale real(default), dimension(:), intent(out) :: sqme_color_c logical, intent(out) :: bad_point integer :: i1, i2, i, n_tot real(double) :: sqme_dble do i = 1, object%data%n_in + object%data%n_out p_recola(:, i) = dble(p(i)%p) end do n_tot = object%data%n_in + object%data%n_out i = 0 do i1 = 1, n_tot do i2 = 1, i1-1 i = i + 1 call rclwrap_compute_color_correlation & (object%recola_ids(i_flv), p_recola, i1, i2, sqme_dble) sqme_color_c(i) = real (sqme_dble, kind=default) select case (abs (object%data%flv_state (i1, i_flv))) case (1:6) sqme_color_c(i) = CF * sqme_color_c(i) case (9,21) sqme_color_c(i) = CA * sqme_color_c(i) end select end do end do end subroutine prc_recola_compute_sqme_color_c_raw @ %def prc_recola_compute_sqme_color_c_raw @ \subsection{Unit tests} <<[[prc_recola_ut.f90]]>>= <> module prc_recola_ut use unit_tests use prc_recola_uti <> <> contains <> end module prc_recola_ut @ %def prc_recola_ut @ <<[[prc_recola_uti.f90]]>>= <> module prc_recola_uti use recola_wrapper !NODEP! use, intrinsic :: iso_c_binding !NODEP! use kinds <> use constants use format_utils, only: write_separator use numeric_utils, only: assert_equal use os_interface use particle_specifiers, only: new_prt_spec use prc_core_def use process_constants use process_libraries use prc_core use prc_omega <> <> contains <> <> end module prc_recola_uti @ %def prc_recola_uti @ <>= public :: prc_recola_test <>= subroutine prc_recola_test (u, results) integer, intent(in) :: u type(test_results_t), intent(inout) :: results <> end subroutine prc_recola_test @ %def prc_recola_test @ \subsubsection{Testing a fixed flavor matrix element computation} <>= function get_omega_parameter_array () result (par) real(default), dimension(25) :: par par = zero par(1) = 1.16637d-5 ! gf par(2) = 91.153480619182744_default ! mZ par(3) = 80.357973609877547_default ! mW par(4) = 125._default ! mH par(5) = rclwrap_get_alpha_s () ! alpha_s par(12) = 173.2_default ! mt par(14) = 2.4942663787728243_default ! wZ par(15) = 2.0842989982782196_default ! wW par(22) = one / sqrt (sqrt (two) * par(1)) ! par%v - Higgs expectation value par(23) = par(3) / par(2) ! par%cw par(24) = sqrt (one - par(23)**2) ! par%sw par(25) = two * par(24) * par(3) / par(22) end function get_omega_parameter_array @ %def get_omega_parameter_array @ <>= call test (prc_recola_1, "prc_recola_1", & "Registering a RECOLA process and computing the amplitude", & u, results) <>= public :: prc_recola_1 <>= subroutine prc_recola_1 (u) integer, intent(in) :: u real(double) :: p(0:3,1:4) real(double) :: sqrts = 500._double real(double) :: m_e = 0._double real(double) :: m_mu = 0._double real(double) :: p_x_out, p_y_out, p_z_out, p_z_in integer :: h_e_p, h_e_m, h_mu_p, h_mu_m, counter real(double) :: sqme integer :: i integer, dimension(:), allocatable :: col_recola, hel_recola complex(double) :: amp_recola complex(default) :: amp_recola_default real(default), parameter :: ee = 0.3 !!! Electromagnetic coupling type(process_library_t) :: lib class(prc_core_def_t), allocatable :: def type(process_def_entry_t), pointer :: entry type(string_t), dimension(:), allocatable :: prt_in, prt_out type(os_data_t) :: os_data type(process_constants_t) :: data class(prc_core_driver_t), allocatable :: driver complex(default) :: amp integer, dimension(:,:), allocatable :: helicities write (u, "(A)") "* Test output: prc_recola_1" write (u, "(A)") "* Purpose: Test interface to RECOLA and compare matrix elements with O'Mega" write (u, "(A)") p_z_in = sqrt ((sqrts / 2)**2 - m_e**2) p_z_out = 0._double p_y_out = sqrts / 10._default p_x_out = sqrt ((sqrts / 2)**2 - p_y_out**2 - p_z_out**2 - m_mu**2) p(:,1) = [sqrts / 2, 0._double, 0._double, p_z_in] p(:,2) = [sqrts / 2, 0._double, 0._double, -p_z_in] p(:,3) = [sqrts / 2, p_x_out, p_y_out, p_z_out] p(:,4) = [sqrts / 2, -p_x_out, -p_y_out, -p_z_out] write (u, "(A)") "Use phase-space point: " do i = 1, 4 write (u, "(4(F12.3,1x))") p(:,1) end do write (u, "(A)") call write_separator (u) write (u, "(A)") write (u, "(A)") "* RECOLA: Evaluate process" counter = 1 call rclwrap_request_generate_processes () write (u, "(A)") "* RECOLA: Define process e+ e- -> mu+ mu- at leading order" call rclwrap_add_process (counter, var_str ('e+ e- -> mu+ mu-'), var_str ('LO')) call rclwrap_define_processes () write (u, "(A)") "* RECOLA: generate process" call rclwrap_generate_processes () call rclwrap_compute_process (1, p, 'LO') call rclwrap_get_helicity_configurations (1, helicities) allocate (hel_recola (4), col_recola (4)) col_recola = [0,0,0,0] write (u, "(A)") "* Setting up Omega to compute the same amplitude" call lib%init (var_str ("omega1")) allocate (prt_in (2), prt_out (2)) prt_in = [var_str ("e+"), var_str ("e-")] prt_out = [var_str ("mu+"), var_str ("mu-")] allocate (omega_def_t :: def) select type (def) type is (omega_def_t) call def%init (var_str ("SM"), prt_in, prt_out, & ufo = .false., ovm = .false., cms_scheme = .true.) end select allocate (entry) call entry%init (var_str ("omega1_a"), model_name = var_str ("SM"), & n_in = 2, n_components = 1) call entry%import_component (1, n_out = 2, & prt_in = new_prt_spec (prt_in), & prt_out = new_prt_spec (prt_out), & method = var_str ("omega"), & variant = def) call lib%append (entry) call os_data%init () call lib%configure (os_data) call lib%write_makefile (os_data, force = .true., verbose = .false.) call lib%clean (os_data, distclean = .false.) call lib%write_driver (force = .true.) call lib%load (os_data) call lib%connect_process (var_str ("omega1_a"), 1, data, driver) select type (driver) type is (omega_driver_t) call driver%init (get_omega_parameter_array (), 3) call driver%new_event (real(p, kind = default)) do i = 1, 6 call rclwrap_get_amplitude (1, 0, 'LO', col_recola, helicities (:,i), amp_recola) end do do i = 1, 16 call rclwrap_get_amplitude (1, 0, 'LO', col_recola, data%hel_state (:,i), amp_recola) amp_recola = amp_recola * cmplx (0, -1, double) amp_recola_default = amp_recola call driver%get_amplitude (1, i, 1, amp) write(u,"(A,4(I2),A)") "Helicity: [",data%hel_state (:,i),"]" call assert_equal (u, amp, amp_recola_default, rel_smallness = 1.E-7_default) end do end select call rclwrap_reset_recola () write (u, "(A)") write (u, "(A)") "* End of test output: prc_recola_1" end subroutine prc_recola_1 @ %def prc_recola_1 @ \subsubsection{Testing a fixed flavor matrix element computation for 2->3} <>= call test (prc_recola_2, "prc_recola_2", & "Registering a RECOLA process and computing the amplitude for 2->3 process", & u, results) <>= public :: prc_recola_2 <>= subroutine prc_recola_2 (u) integer, intent(in) :: u real(double) :: p(0:3,1:5) real(double) :: sqrts = 700._double real(double) :: m_e = 0._double real(double) :: m_mu = 0._double real(double) :: p_x_out, p_y_out, p_z_out, p_z_in real(double) :: sqme integer :: i integer, dimension(:), allocatable :: col_recola, hel_recola integer, dimension(:,:), allocatable :: helicities complex(double) :: amp_recola complex(default) :: amp_recola_default real(default), parameter :: ee = 0.3 !!! Electromagnetic coupling type(process_library_t) :: lib class(prc_core_def_t), allocatable :: def type(process_def_entry_t), pointer :: entry type(string_t), dimension(:), allocatable :: prt_in, prt_out type(os_data_t) :: os_data type(process_constants_t) :: data class(prc_core_driver_t), allocatable :: driver complex(default) :: amp integer :: n_allowed write (u, "(A)") "* Test output: prc_recola_2" write (u, "(A)") "* Purpose: Test interface to RECOLA and compare matrix elements with O'Mega for 2->3 process" write (u, "(A)") p_z_in = sqrt ((sqrts / 2)**2 - m_e**2) p(:,1) = [sqrts / 2, 0._double, 0._double, p_z_in] p(:,2) = [sqrts / 2, 0._double, 0._double, -p_z_in] p(:,3) = [243.49323116_double, -141.69619338_double, -108.30640321_double, 165.77353656_double] p(:,4) = [337.53250628_double, 143.95931207_double, 110.19717026_double, -284.71124482_double] p(:,5) = [118.97426257_double, -2.2631186860_double, -1.8907670459_double, 118.93770827_double] write (u, "(A)") "Use phase-space point: " do i = 1, 5 write (u, "(4(F12.3,1x))") p(:,1) end do write (u, "(A)") call write_separator (u) write (u, "(A)") write (u, "(A)") "* RECOLA: Evaluate process" call rclwrap_request_generate_processes () write (u, "(A)") "* RECOLA: Define process e+ e- -> mu+ mu- A at leading order" call rclwrap_add_process (2, var_str ('e+ e- -> mu+ mu- A'), var_str ('LO')) call rclwrap_define_processes () write (u, "(A)") "* RECOLA: generate process" call rclwrap_generate_processes () call rclwrap_compute_process (2, p, 'LO') call rclwrap_get_helicity_configurations (2, helicities) allocate (hel_recola (5), col_recola (5)) col_recola = [0,0,0,0,0] write (u, "(A)") "* Setting up Omega to compute the same amplitude" call lib%init (var_str ("omega2")) allocate (prt_in (2), prt_out (3)) prt_in = [var_str ("e+"), var_str ("e-")] prt_out = [var_str ("mu+"), var_str ("mu-"), var_str("A")] allocate (omega_def_t :: def) select type (def) type is (omega_def_t) call def%init (var_str ("SM"), prt_in, prt_out, & ufo = .false., ovm = .false.) end select allocate (entry) call entry%init (var_str ("omega2_a"), model_name = var_str ("SM"), & n_in = 2, n_components = 1) call entry%import_component (1, n_out = 3, & prt_in = new_prt_spec (prt_in), & prt_out = new_prt_spec (prt_out), & method = var_str ("omega"), & variant = def) call lib%append (entry) call os_data%init () call lib%configure (os_data) call lib%write_makefile (os_data, force = .true., verbose = .false.) call lib%clean (os_data, distclean = .false.) call lib%write_driver (force = .true.) call lib%load (os_data) call lib%connect_process (var_str ("omega2_a"), 1, data, driver) select type (driver) type is (omega_driver_t) call driver%init (get_omega_parameter_array (), 3) call driver%new_event (real(p, kind = default)) do i = 1, 32 call rclwrap_get_amplitude & (2, 0, 'LO', col_recola, data%hel_state (:,i), amp_recola) if (data%hel_state(3,i) * data%hel_state(4,i) * & data%hel_state(5,i) == -1) then amp_recola = amp_recola * cmplx (0, -1, double) else amp_recola = amp_recola * cmplx (0, 1, double) end if amp_recola_default = amp_recola call driver%get_amplitude (1, i, 1, amp) write(u,"(A,5(I2),A)") "Helicity: [", data%hel_state (:,i),"]" write(u,"(A,2(F12.7,1x),A,2(F12.7,1x))") "RECOLA:", & amp_recola,", O'MEGA:", amp call assert_equal & (u, amp, amp_recola_default, rel_smallness = 1.E-6_default) end do end select call rclwrap_reset_recola () write (u, "(A)") write (u, "(A)") "* End of test output: prc_recola_2" end subroutine prc_recola_2 @ %def prc_recola_2 @ Index: trunk/src/blha/blha.nw =================================================================== --- trunk/src/blha/blha.nw (revision 8251) +++ trunk/src/blha/blha.nw (revision 8252) @@ -1,3354 +1,3354 @@ % -*- ess-noweb-default-code-mode: f90-mode; noweb-default-code-mode: f90-mode; -*- % WHIZARD code as NOWEB source: matrix elements and process libraries %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{BLHA Interface} \includemodulegraph{blha} The code in this chapter implements support for the BLHA record that communicates data for NLO processes. These are the modules: \begin{description} \item[blha\_config] \item[blha\_interface] \item[blha\_driver] \end{description} @ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Module definition} These modules implement the communication with one loop matrix element providers according to the Binoth LesHouches Accord Interface. The actual matrix element(s) are loaded as a dynamic library. This module defines the common OLP-interfaces defined through the Binoth Les-Houches accord. <<[[blha_olp_interfaces.f90]]>>= <> module blha_olp_interfaces use, intrinsic :: iso_c_binding !NODEP! use, intrinsic :: iso_fortran_env use kinds <> <> use constants use numeric_utils, only: vanishes use numeric_utils, only: extend_integer_array, crop_integer_array use io_units use string_utils use physics_defs use diagnostics use os_interface use lorentz use sm_qcd use interactions use flavors use model_data use pdg_arrays, only: is_gluon, is_quark use prclib_interfaces use process_libraries use prc_core_def use prc_core use prc_external use blha_config <> <> <> <> <> <> contains <> end module blha_olp_interfaces @ %def module blha_olp_interfaces @ <>= public :: blha_template_t <>= type :: blha_template_t integer :: I_BORN = 0 integer :: I_REAL = 1 integer :: I_LOOP = 2 integer :: I_SUB = 3 integer :: I_DGLAP = 4 logical, dimension(0:4) :: compute_component logical :: include_polarizations = .false. logical :: switch_off_muon_yukawas = .false. logical :: use_internal_color_correlations = .true. real(default) :: external_top_yukawa = -1._default integer :: ew_scheme contains <> end type blha_template_t @ %def blha_template_t @ <>= procedure :: write => blha_template_write <>= subroutine blha_template_write (blha_template, unit) class(blha_template_t), intent(in) :: blha_template integer, intent(in), optional :: unit integer :: u u = given_output_unit (unit) write (u,"(A,4(L1))") "Compute components: ", & blha_template%compute_component write (u,"(A,L1)") "Include polarizations: ", & blha_template%include_polarizations write (u,"(A,L1)") "Switch off muon yukawas: ", & blha_template%switch_off_muon_yukawas write (u,"(A,L1)") "Use internal color correlations: ", & blha_template%use_internal_color_correlations end subroutine blha_template_write @ %def blha_template_write @ Compute the total number of used helicity states for the given particle PDG codes, given a model. Applies only if polarization is supported. This yields the [[n_hel]] value as required below. <>= procedure :: get_n_hel => blha_template_get_n_hel <>= function blha_template_get_n_hel (blha_template, pdg, model) result (n_hel) class(blha_template_t), intent(in) :: blha_template integer, dimension(:), intent(in) :: pdg class(model_data_t), intent(in), target :: model integer :: n_hel type(flavor_t) :: flv integer :: f n_hel = 1 if (blha_template%include_polarizations) then do f = 1, size (pdg) call flv%init (pdg(f), model) n_hel = n_hel * flv%get_multiplicity () end do end if end function blha_template_get_n_hel @ %def blha_template_get_n_hel @ <>= integer, parameter :: I_ALPHA = 1 integer, parameter :: I_GF = 2 integer, parameter :: I_SW2 = 3 <>= public :: prc_blha_t <>= type, abstract, extends (prc_external_t) :: prc_blha_t integer :: n_particles integer :: n_hel integer :: n_proc integer, dimension(:, :), allocatable :: i_tree, i_spin_c, i_color_c integer, dimension(:, :), allocatable :: i_virt integer, dimension(:, :), allocatable :: i_hel logical, dimension(3) :: ew_parameter_mask integer :: sqme_tree_pos contains <> end type prc_blha_t @ %def prc_blha_t @ Obviously, this process-core type uses the BLHA interface. <>= procedure, nopass :: uses_blha => prc_blha_uses_blha <>= function prc_blha_uses_blha () result (flag) logical :: flag flag = .true. end function prc_blha_uses_blha @ %def prc_blha_uses_blha @ <>= public :: blha_driver_t <>= type, abstract, extends (prc_external_driver_t) :: blha_driver_t type(string_t) :: contract_file type(string_t) :: nlo_suffix logical :: include_polarizations = .false. logical :: switch_off_muon_yukawas = .false. real(default) :: external_top_yukawa = -1.0 procedure(olp_start),nopass, pointer :: & blha_olp_start => null () procedure(olp_eval), nopass, pointer :: & blha_olp_eval => null() procedure(olp_info), nopass, pointer :: & blha_olp_info => null () procedure(olp_set_parameter), nopass, pointer :: & blha_olp_set_parameter => null () procedure(olp_eval2), nopass, pointer :: & blha_olp_eval2 => null () procedure(olp_option), nopass, pointer :: & blha_olp_option => null () procedure(olp_polvec), nopass, pointer :: & blha_olp_polvec => null () procedure(olp_finalize), nopass, pointer :: & blha_olp_finalize => null () procedure(olp_print_parameter), nopass, pointer :: & blha_olp_print_parameter => null () contains <> end type blha_driver_t @ @ %def blha_driver_t <>= public :: prc_blha_writer_t <>= type, abstract, extends (prc_external_writer_t) :: prc_blha_writer_t type(blha_configuration_t) :: blha_cfg contains <> end type prc_blha_writer_t @ @ %def prc_blha_writer_t <>= public :: blha_def_t <>= type, abstract, extends (prc_external_def_t) :: blha_def_t type(string_t) :: suffix contains <> end type blha_def_t @ %def blha_def_t @ <>= public :: blha_state_t <>= type, abstract, extends (prc_external_state_t) :: blha_state_t contains <> end type blha_state_t @ %def blha_state_t @ <>= procedure :: reset_new_kinematics => blha_state_reset_new_kinematics <>= subroutine blha_state_reset_new_kinematics (object) class(blha_state_t), intent(inout) :: object object%new_kinematics = .true. end subroutine blha_state_reset_new_kinematics @ %def blha_state_reset_new_kinematics @ <>= integer, parameter, public :: OLP_PARAMETER_LIMIT = 10 integer, parameter, public :: OLP_MOMENTUM_LIMIT = 50 integer, parameter, public :: OLP_RESULTS_LIMIT = 60 <>= public :: olp_start <>= interface subroutine olp_start (contract_file_name, ierr) bind (C,name = "OLP_Start") import character(kind = c_char, len = 1), intent(in) :: contract_file_name integer(kind = c_int), intent(out) :: ierr end subroutine olp_start end interface @ %def olp_start_interface @ <>= public :: olp_eval <>= interface subroutine olp_eval (label, momenta, mu, parameters, res) & bind (C, name = "OLP_EvalSubProcess") import integer(kind = c_int), value, intent(in) :: label real(kind = c_double), value, intent(in) :: mu real(kind = c_double), dimension(OLP_MOMENTUM_LIMIT), intent(in) :: & momenta real(kind = c_double), dimension(OLP_PARAMETER_LIMIT), intent(in) :: & parameters real(kind = c_double), dimension(OLP_RESULTS_LIMIT), intent(out) :: res end subroutine olp_eval end interface @ %def olp_eval interface @ <>= public :: olp_info <>= interface subroutine olp_info (olp_file, olp_version, message) bind(C) import character(kind = c_char), intent(inout), dimension(15) :: olp_file character(kind = c_char), intent(inout), dimension(15) :: olp_version character(kind = c_char), intent(inout), dimension(255) :: message end subroutine olp_info end interface @ %def olp_info interface @ <>= public :: olp_set_parameter <>= interface subroutine olp_set_parameter & (variable_name, real_part, complex_part, success) bind(C) import character(kind = c_char,len = 1), intent(in) :: variable_name real(kind = c_double), intent(in) :: real_part, complex_part integer(kind = c_int), intent(out) :: success end subroutine olp_set_parameter end interface @ %def olp_set_parameter_interface @ <>= public :: olp_eval2 <>= interface subroutine olp_eval2 (label, momenta, mu, res, acc) bind(C) import integer(kind = c_int), intent(in) :: label real(kind = c_double), intent(in) :: mu real(kind = c_double), dimension(OLP_MOMENTUM_LIMIT), intent(in) :: momenta real(kind = c_double), dimension(OLP_RESULTS_LIMIT), intent(out) :: res real(kind = c_double), intent(out) :: acc end subroutine olp_eval2 end interface @ %def olp_eval2 interface @ <>= public :: olp_option <>= interface subroutine olp_option (line, stat) bind(C) import character(kind = c_char, len=1), intent(in) :: line integer(kind = c_int), intent(out) :: stat end subroutine end interface @ %def olp_option_interface @ <>= public :: olp_polvec <>= interface subroutine olp_polvec (p, q, eps) bind(C) import real(kind = c_double), dimension(0:3), intent(in) :: p, q real(kind = c_double), dimension(0:7), intent(out) :: eps end subroutine end interface @ %def olp_polvec_interface @ <>= public :: olp_finalize <>= interface subroutine olp_finalize () bind(C) import end subroutine olp_finalize end interface @ %def olp_finalize_interface @ <>= public :: olp_print_parameter <>= interface subroutine olp_print_parameter (filename) bind(C) import character(kind = c_char, len = 1), intent(in) :: filename end subroutine olp_print_parameter end interface @ %def olp_print_parameter_interface @ <>= public :: blha_result_array_size <>= pure function blha_result_array_size (n_part, amp_type) result (rsize) integer, intent(in) :: n_part, amp_type integer :: rsize select case (amp_type) case (BLHA_AMP_TREE) rsize = 1 case (BLHA_AMP_LOOP) rsize = 4 case (BLHA_AMP_COLOR_C) rsize = n_part * (n_part - 1) / 2 case (BLHA_AMP_SPIN_C) rsize = 2 * n_part**2 case default rsize = 0 end select end function blha_result_array_size @ %def blha_result_array_size @ <>= procedure :: create_momentum_array => prc_blha_create_momentum_array <>= function prc_blha_create_momentum_array (object, p) result (mom) class(prc_blha_t), intent(in) :: object type(vector4_t), intent(in), dimension(:) :: p real(double), dimension(5*object%n_particles) :: mom integer :: n, i, k n = size (p) if (n > 10) call msg_fatal ("Number of external particles exceeds" & // "size of BLHA-internal momentum array") mom = zero k = 1 do i = 1, n mom(k : k + 3) = vector4_get_components (p(i)) mom(k + 4) = invariant_mass (p(i)) k = k + 5 end do end function prc_blha_create_momentum_array @ %def prc_blha_create_momentum_array @ <>= procedure :: init => blha_template_init <>= subroutine blha_template_init (template, requires_polarizations, & switch_off_muon_yukawas, external_top_yukawa, ew_scheme) class(blha_template_t), intent(inout) :: template logical, intent(in) :: requires_polarizations, switch_off_muon_yukawas real(default), intent(in) :: external_top_yukawa type(string_t), intent(in) :: ew_scheme template%compute_component = .false. template%include_polarizations = requires_polarizations template%switch_off_muon_yukawas = switch_off_muon_yukawas template%external_top_yukawa = external_top_yukawa template%ew_scheme = ew_scheme_string_to_int (ew_scheme) end subroutine blha_template_init @ %def blha_template_init @ <>= procedure :: set_born => blha_template_set_born procedure :: set_real_trees => blha_template_set_real_trees procedure :: set_loop => blha_template_set_loop procedure :: set_subtraction => blha_template_set_subtraction procedure :: set_dglap => blha_template_set_dglap <>= subroutine blha_template_set_born (template) class(blha_template_t), intent(inout) :: template template%compute_component (template%I_BORN) = .true. end subroutine blha_template_set_born subroutine blha_template_set_real_trees (template) class(blha_template_t), intent(inout) :: template template%compute_component (template%I_REAL) = .true. end subroutine blha_template_set_real_trees subroutine blha_template_set_loop (template) class(blha_template_t), intent(inout) :: template template%compute_component(template%I_LOOP) = .true. end subroutine blha_template_set_loop subroutine blha_template_set_subtraction (template) class(blha_template_t), intent(inout) :: template template%compute_component (template%I_SUB) = .true. end subroutine blha_template_set_subtraction subroutine blha_template_set_dglap (template) class(blha_template_t), intent(inout) :: template template%compute_component (template%I_DGLAP) = .true. end subroutine blha_template_set_dglap @ %def blha_template_set_components @ <>= procedure :: set_internal_color_correlations & => blha_template_set_internal_color_correlations <>= subroutine blha_template_set_internal_color_correlations (template) class(blha_template_t), intent(inout) :: template template%use_internal_color_correlations = .true. end subroutine blha_template_set_internal_color_correlations @ %def blha_template_set_internal_color_correlations @ <>= procedure :: get_internal_color_correlations & => blha_template_get_internal_color_correlations <>= pure function blha_template_get_internal_color_correlations (template) & result (val) logical :: val class(blha_template_t), intent(in) :: template val = template%use_internal_color_correlations end function blha_template_get_internal_color_correlations @ %def blha_template_use_internal_color_correlations @ <>= procedure :: compute_born => blha_template_compute_born procedure :: compute_real_trees => blha_template_compute_real_trees procedure :: compute_loop => blha_template_compute_loop procedure :: compute_subtraction => blha_template_compute_subtraction procedure :: compute_dglap => blha_template_compute_dglap <>= pure function blha_template_compute_born (template) result (val) class(blha_template_t), intent(in) :: template logical :: val val = template%compute_component (template%I_BORN) end function blha_template_compute_born pure function blha_template_compute_real_trees (template) result (val) class(blha_template_t), intent(in) :: template logical :: val val = template%compute_component (template%I_REAL) end function blha_template_compute_real_trees pure function blha_template_compute_loop (template) result (val) class(blha_template_t), intent(in) :: template logical :: val val = template%compute_component (template%I_LOOP) end function blha_template_compute_loop pure function blha_template_compute_subtraction (template) result (val) class(blha_template_t), intent(in) :: template logical :: val val = template%compute_component (template%I_SUB) end function blha_template_compute_subtraction pure function blha_template_compute_dglap (template) result (val) class(blha_template_t), intent(in) :: template logical :: val val = template%compute_component (template%I_DGLAP) end function blha_template_compute_dglap @ %def blha_template_compute @ <>= procedure :: check => blha_template_check <>= function blha_template_check (template) result (val) class(blha_template_t), intent(in) :: template logical :: val val = count (template%compute_component) == 1 end function blha_template_check @ %def blha_template_check @ <>= procedure :: reset => blha_template_reset <>= subroutine blha_template_reset (template) class(blha_template_t), intent(inout) :: template template%compute_component = .false. end subroutine blha_template_reset @ %def blha_template_reset @ <>= procedure :: write => prc_blha_writer_write <>= subroutine prc_blha_writer_write (writer, unit) class(prc_blha_writer_t), intent(in) :: writer integer, intent(in) :: unit write (unit, "(1x,A)") char (writer%get_process_string ()) end subroutine prc_blha_writer_write @ @ %def prc_blha_writer_write <>= procedure :: get_process_string => prc_blha_writer_get_process_string <>= function prc_blha_writer_get_process_string (writer) result (s_proc) class(prc_blha_writer_t), intent(in) :: writer type(string_t) :: s_proc s_proc = var_str ("") end function prc_blha_writer_get_process_string -@ %def gosam_writer_get_process_string +@ %def prc_blha_writer_get_process_string @ <>= procedure :: get_n_proc => prc_blha_writer_get_n_proc <>= function prc_blha_writer_get_n_proc (writer) result (n_proc) class(prc_blha_writer_t), intent(in) :: writer integer :: n_proc n_proc = blha_configuration_get_n_proc (writer%blha_cfg) end function prc_blha_writer_get_n_proc @ %def prc_blha_writer_get_n_proc @ <>= procedure(blha_driver_set_GF), deferred :: & set_GF <>= abstract interface subroutine blha_driver_set_GF (driver, GF) import class(blha_driver_t), intent(inout) :: driver real(default), intent(in) :: GF end subroutine blha_driver_set_GF end interface @ %def blha_driver_set_GF @ <>= procedure(blha_driver_set_alpha_s), deferred :: & set_alpha_s <>= abstract interface subroutine blha_driver_set_alpha_s (driver, alpha_s) import class(blha_driver_t), intent(in) :: driver real(default), intent(in) :: alpha_s end subroutine blha_driver_set_alpha_s end interface @ %def set_alpha_s interface @ <>= procedure(blha_driver_set_weinberg_angle), deferred :: & set_weinberg_angle <>= abstract interface subroutine blha_driver_set_weinberg_angle (driver, sw2) import class(blha_driver_t), intent(inout) :: driver real(default), intent(in) :: sw2 end subroutine blha_driver_set_weinberg_angle end interface @ %def blha_driver_set_weinberg_angle @ <>= procedure(blha_driver_set_alpha_qed), deferred :: set_alpha_qed <>= abstract interface subroutine blha_driver_set_alpha_qed (driver, alpha) import class(blha_driver_t), intent(inout) :: driver real(default), intent(in) :: alpha end subroutine blha_driver_set_alpha_qed end interface @ %def blha_driver_set_alpha_qed @ <>= procedure(blha_driver_print_alpha_s), deferred :: & print_alpha_s <>= abstract interface subroutine blha_driver_print_alpha_s (object) import class(blha_driver_t), intent(in) :: object end subroutine blha_driver_print_alpha_s end interface @ %def print_alpha_s interface @ <>= public :: parameter_error_message <>= subroutine parameter_error_message (par) type(string_t), intent(in) :: par type(string_t) :: message message = "Setting of parameter " // par & // "failed. This happens because the chosen " & // "EWScheme in the BLHA file does not fit " & // "your parameter choice" call msg_fatal (char (message)) end subroutine parameter_error_message @ %def parameter_error_message @ <>= procedure :: set_mass_and_width => blha_driver_set_mass_and_width <>= subroutine blha_driver_set_mass_and_width & (driver, i_pdg, mass, width) class(blha_driver_t), intent(inout) :: driver integer, intent(in) :: i_pdg real(default), intent(in), optional :: mass real(default), intent(in), optional :: width type(string_t) :: buf character(kind=c_char,len=20) :: c_string integer :: ierr if (present (mass)) then buf = 'mass(' // str (abs(i_pdg)) // ')' c_string = char(buf) // c_null_char call driver%blha_olp_set_parameter & (c_string, dble(mass), 0._double, ierr) if (ierr == 0) then buf = "BLHA driver: Attempt to set mass of particle " // & str (abs(i_pdg)) // "failed" call msg_fatal (char(buf)) end if end if if (present (width)) then buf = 'width(' // str (abs(i_pdg)) // ')' c_string = char(buf)//c_null_char call driver%blha_olp_set_parameter & (c_string, dble(width), 0._double, ierr) if (ierr == 0) then buf = "BLHA driver: Attempt to set width of particle " // & str (abs(i_pdg)) // "failed" call msg_fatal (char(buf)) end if end if end subroutine blha_driver_set_mass_and_width @ %def blha_driver_set_mass_and_width @ <>= procedure(blha_driver_init_dlaccess_to_library), deferred :: & init_dlaccess_to_library <>= abstract interface subroutine blha_driver_init_dlaccess_to_library & (object, os_data, dlaccess, success) import class(blha_driver_t), intent(in) :: object type(os_data_t), intent(in) :: os_data type(dlaccess_t), intent(out) :: dlaccess logical, intent(out) :: success end subroutine blha_driver_init_dlaccess_to_library end interface @ %def interface blha_driver_init_dlaccess_to_library @ <>= procedure :: load => blha_driver_load <>= subroutine blha_driver_load (object, os_data, success) class(blha_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 ("OLP_Start")) call c_f_procpointer (c_fptr, object%blha_olp_start) call check_for_error (var_str ("OLP_Start")) c_fptr = dlaccess_get_c_funptr (dlaccess, var_str ("OLP_EvalSubProcess")) call c_f_procpointer (c_fptr, object%blha_olp_eval) call check_for_error (var_str ("OLP_EvalSubProcess")) c_fptr = dlaccess_get_c_funptr (dlaccess, var_str ("OLP_Info")) call c_f_procpointer (c_fptr, object%blha_olp_info) call check_for_error (var_str ("OLP_Info")) c_fptr = dlaccess_get_c_funptr (dlaccess, var_str ("OLP_SetParameter")) call c_f_procpointer (c_fptr, object%blha_olp_set_parameter) call check_for_error (var_str ("OLP_SetParameter")) c_fptr = dlaccess_get_c_funptr (dlaccess, var_str ("OLP_EvalSubProcess2")) call c_f_procpointer (c_fptr, object%blha_olp_eval2) call check_for_error (var_str ("OLP_EvalSubProcess2")) !!! The following three functions are not implemented in OpenLoops. !!! In another BLHA provider, they need to be implemented separately. !!! c_fptr = dlaccess_get_c_funptr (dlaccess, var_str ("OLP_Option")) !!! call c_f_procpointer (c_fptr, object%blha_olp_option) !!! call check_for_error (var_str ("OLP_Option")) !!! c_fptr = dlaccess_get_c_funptr (dlaccess, var_str ("OLP_Polvec")) !!! call c_f_procpointer (c_fptr, object%blha_olp_polvec) !!! call check_for_error (var_str ("OLP_Polvec")) !!! c_fptr = dlaccess_get_c_funptr (dlaccess, var_str ("OLP_Finalize")) !!! call c_f_procpointer (c_fptr, object%blha_olp_finalize) !!! call check_for_error (var_str ("OLP_Finalize")) c_fptr = dlaccess_get_c_funptr (dlaccess, var_str ("OLP_PrintParameter")) call c_f_procpointer (c_fptr, object%blha_olp_print_parameter) call check_for_error (var_str ("OLP_PrintParameter")) 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 blha_driver_load @ %def blha_driver_load @ <>= integer, parameter :: LEN_MAX_FLAVOR_STRING = 100 integer, parameter :: N_MAX_FLAVORS = 100 <>= procedure :: read_contract_file => blha_driver_read_contract_file <>= subroutine blha_driver_read_contract_file (driver, flavors, & amp_type, flv_index, hel_index, label, helicities) class(blha_driver_t), intent(inout) :: driver integer, intent(in), dimension(:,:) :: flavors integer, intent(out), dimension(:), allocatable :: amp_type, & flv_index, hel_index, label integer, intent(out), dimension(:,:) :: helicities integer :: unit, filestat character(len=LEN_MAX_FLAVOR_STRING) :: rd_line logical :: read_flavor, give_warning integer :: label_count, i_flv integer :: i_hel, n_in integer :: i_next, n_entries integer, dimension(size(flavors, 1) + 2) :: i_array integer, dimension(size(flavors, 1) + 2) :: hel_array integer, parameter :: NO_NUMBER = -1000 integer, parameter :: PROC_NOT_FOUND = -1001 integer, parameter :: list_incr = 50 integer :: n_found allocate (amp_type (N_MAX_FLAVORS), flv_index (N_MAX_FLAVORS), & hel_index (N_MAX_FLAVORS), label (N_MAX_FLAVORS)) amp_type = -1; flv_index = -1; hel_index = -1; label = -1 helicities = 0 n_in = size (helicities, dim = 2) n_entries = size (flavors, 1) + 2 unit = free_unit () open (unit, file = char (driver%contract_file), status="old") read_flavor = .false. label_count = 1 i_hel = 1 n_found = 0 give_warning = .false. do read (unit, "(A)", iostat = filestat) rd_line if (filestat == iostat_end) then exit else if (rd_line(1:13) == 'AmplitudeType') then if (i_hel > 2 * n_in) i_hel = 1 i_next = find_next_word_index (rd_line, 13) if (label_count > size (amp_type)) & call extend_integer_array (amp_type, list_incr) if (rd_line(i_next : i_next + 4) == 'Loop') then amp_type(label_count) = BLHA_AMP_LOOP else if (rd_line(i_next : i_next + 4) == 'Tree') then amp_type(label_count) = BLHA_AMP_TREE else if (rd_line(i_next : i_next + 6) == 'ccTree') then amp_type(label_count) = BLHA_AMP_COLOR_C else if (rd_line(i_next : i_next + 6) == 'scTree' .or. & rd_line(i_next : i_next + 14) == 'sctree_polvect') then amp_type(label_count) = BLHA_AMP_SPIN_C else call msg_fatal ("AmplitudeType present but AmpType not known!") end if read_flavor = .true. else if (read_flavor) then i_array = create_flavor_string (rd_line, n_entries) if (driver%include_polarizations) then hel_array = create_helicity_string (rd_line, n_entries) call check_helicity_array (hel_array, n_entries, n_in) else hel_array = 0 end if if (.not. all (i_array == PROC_NOT_FOUND)) then do i_flv = 1, size (flavors, 2) if (all (i_array (1 : n_entries - 2) == flavors (:,i_flv))) then if (label_count > size (label)) & call extend_integer_array (label, list_incr) label(label_count) = i_array (n_entries) if (label_count > size (flv_index)) & call extend_integer_array (flv_index, list_incr) flv_index (label_count) = i_flv if (label_count > size (hel_index)) & call extend_integer_array (hel_index, list_incr) hel_index (label_count) = i_hel if (driver%include_polarizations) then helicities (label(label_count), :) = hel_array (1:n_in) i_hel = i_hel + 1 end if n_found = n_found + 1 label_count = label_count + 1 exit end if end do give_warning = .false. else give_warning = .true. end if read_flavor = .false. end if end if end do call crop_integer_array (amp_type, label_count-1) if (n_found == 0) then call msg_fatal ("The desired process has not been found ", & [var_str ("by the OLP-Provider. Maybe the value of alpha_power "), & var_str ("or alphas_power does not correspond to the process. "), & var_str ("If you are using OpenLoops, you can set the option "), & var_str ("openloops_verbosity to a value larger than 1 to obtain "), & var_str ("more information")]) else if (give_warning) then call msg_warning ("Some processes have not been found in the OLC file.", & [var_str ("This is because these processes do not fit the required "), & var_str ("coupling alpha_power and alphas_power. Be aware that the "), & var_str ("results of this calculation are not necessarily an accurate "), & var_str ("description of the physics of interest.")]) end if close(unit) contains function create_flavor_string (s, n_entries) result (i_array) character(len=LEN_MAX_FLAVOR_STRING), intent(in) :: s integer, intent(in) :: n_entries integer, dimension(n_entries) :: i_array integer :: k, current_position integer :: i_entry k = 1; current_position = 1 do if (current_position > LEN_MAX_FLAVOR_STRING) & call msg_fatal ("Read OLC File: Current position exceeds maximum value") if (s(current_position:current_position) /= " ") then call create_flavor (s, i_entry, current_position) if (i_entry /= NO_NUMBER .and. i_entry /= PROC_NOT_FOUND) then i_array(k) = i_entry k = k + 1 if (k > n_entries) then return else call increment_current_position (s, current_position) end if else if (i_entry == PROC_NOT_FOUND) then i_array = PROC_NOT_FOUND return else call increment_current_position (s, current_position) end if else call increment_current_position (s, current_position) end if end do end function create_flavor_string function create_helicity_string (s, n_entries) result (hel_array) character(len = LEN_MAX_FLAVOR_STRING), intent(in) :: s integer, intent(in) :: n_entries integer, dimension(n_entries) :: hel_array integer :: k, current_position integer :: hel k = 1; current_position = 1 do if (current_position > LEN_MAX_FLAVOR_STRING) & call msg_fatal ("Read OLC File: Current position exceeds maximum value") if (s(current_position:current_position) /= " ") then call create_helicity (s, hel, current_position) if (hel >= -1 .and. hel <= 1) then hel_array(k) = hel k = k + 1 if (k > n_entries) then return else call increment_current_position (s, current_position) end if else call increment_current_position (s, current_position) end if else call increment_current_position (s, current_position) end if end do end function create_helicity_string subroutine increment_current_position (s, current_position) character(len = LEN_MAX_FLAVOR_STRING), intent(in) :: s integer, intent(inout) :: current_position current_position = find_next_word_index (s, current_position) end subroutine increment_current_position subroutine get_next_buffer (s, current_position, buf, last_buffer_index) character(len = LEN_MAX_FLAVOR_STRING), intent(in) :: s integer, intent(inout) :: current_position character(len = 10), intent(out) :: buf integer, intent(out) :: last_buffer_index integer :: i i = 1; buf = "" do if (s(current_position:current_position) /= " ") then buf(i:i) = s(current_position:current_position) i = i + 1; current_position = current_position + 1 else exit end if end do last_buffer_index = i end subroutine get_next_buffer function is_particle_buffer (buf, i) result (valid) logical :: valid character(len = 10), intent(in) :: buf integer, intent(in) :: i valid = (buf(1 : i - 1) /= "->" .and. buf(1 : i - 1) /= "|" & .and. buf(1 : i - 1) /= "Process") end function is_particle_buffer subroutine create_flavor (s, i_particle, current_position) character(len=LEN_MAX_FLAVOR_STRING), intent(in) :: s integer, intent(out) :: i_particle integer, intent(inout) :: current_position character(len=10) :: buf integer :: i, last_buffer_index call get_next_buffer (s, current_position, buf, last_buffer_index) i = last_buffer_index if (is_particle_buffer (buf, i)) then call strip_helicity (buf, i) i_particle = read_ival (var_str (buf(1 : i - 1))) else if (buf(1 : i - 1) == "Process") then i_particle = PROC_NOT_FOUND else i_particle = NO_NUMBER end if end subroutine create_flavor subroutine create_helicity (s, helicity, current_position) character(len = LEN_MAX_FLAVOR_STRING), intent(in) :: s integer, intent(out) :: helicity integer, intent(inout) :: current_position character(len = 10) :: buf integer :: i, last_buffer_index logical :: success call get_next_buffer (s, current_position, buf, last_buffer_index) i = last_buffer_index if (is_particle_buffer (buf, i)) then call strip_flavor (buf, i, helicity, success) else helicity = 0 end if end subroutine create_helicity subroutine strip_helicity (buf, i) character(len = 10), intent(in) :: buf integer, intent(inout) :: i integer :: i_last i_last = i - 1 if (i_last < 4) return if (buf(i_last - 2 : i_last) == "(1)") then i = i - 3 else if (buf(i_last - 3 : i_last) == "(-1)") then i = i - 4 end if end subroutine strip_helicity subroutine strip_flavor (buf, i, helicity, success) character(len = 10), intent(in) :: buf integer, intent(in) :: i integer, intent(out) :: helicity logical, intent(out) :: success integer :: i_last i_last = i - 1 helicity = 0 if (i_last < 4) return if (buf(i_last - 2 : i_last) == "(1)") then helicity = 1 success = .true. else if (buf(i_last - 3 : i_last) == "(-1)") then helicity = -1 success = .true. else success = .false. end if end subroutine strip_flavor function find_next_word_index (word, i_start) result (i_next) character(len = LEN_MAX_FLAVOR_STRING), intent(in) :: word integer, intent(in) :: i_start integer :: i_next i_next = i_start + 1 do if (word(i_next : i_next) /= " ") then exit else i_next = i_next + 1 end if if (i_next > LEN_MAX_FLAVOR_STRING) & call msg_fatal ("Find next word: line limit exceeded") end do end function find_next_word_index subroutine check_helicity_array (hel_array, n_entries, n_in) integer, intent(in), dimension(:) :: hel_array integer, intent(in) :: n_entries, n_in integer :: n_particles, i logical :: valid n_particles = n_entries - 2 !!! only allow polarisations for incoming fermions for now valid = all (hel_array (n_in + 1 : n_particles) == 0) do i = 1, n_in valid = valid .and. (hel_array(i) == 1 .or. hel_array(i) == -1) end do if (.not. valid) & call msg_fatal ("Invalid helicities encountered!") end subroutine check_helicity_array end subroutine blha_driver_read_contract_file @ %def blha_driver_read_contract_file @ <>= procedure :: set_alpha_qed => prc_blha_set_alpha_qed <>= subroutine prc_blha_set_alpha_qed (object, model) class(prc_blha_t), intent(inout) :: object type(model_data_t), intent(in), target :: model real(default) :: alpha alpha = one / model%get_real (var_str ('alpha_em_i')) select type (driver => object%driver) class is (blha_driver_t) call driver%set_alpha_qed (alpha) end select end subroutine prc_blha_set_alpha_qed @ %def prc_blha_set_alpha_qed @ <>= procedure :: set_GF => prc_blha_set_GF <>= subroutine prc_blha_set_GF (object, model) class(prc_blha_t), intent(inout) :: object type(model_data_t), intent(in), target :: model real(default) :: GF GF = model%get_real (var_str ('GF')) select type (driver => object%driver) class is (blha_driver_t) call driver%set_GF (GF) end select end subroutine prc_blha_set_GF @ %def prc_blha_set_GF @ <>= procedure :: set_weinberg_angle => prc_blha_set_weinberg_angle <>= subroutine prc_blha_set_weinberg_angle (object, model) class(prc_blha_t), intent(inout) :: object type(model_data_t), intent(in), target :: model real(default) :: sw2 sw2 = model%get_real (var_str ('sw2')) select type (driver => object%driver) class is (blha_driver_t) call driver%set_weinberg_angle (sw2) end select end subroutine prc_blha_set_weinberg_angle @ %def prc_blha_set_weinberg_angle @ <>= procedure :: set_electroweak_parameters => & prc_blha_set_electroweak_parameters <>= subroutine prc_blha_set_electroweak_parameters (object, model) class(prc_blha_t), intent(inout) :: object type(model_data_t), intent(in), target :: model if (count (object%ew_parameter_mask) == 0) then call msg_fatal ("Cannot decide EW parameter setting: No scheme set!") else if (count (object%ew_parameter_mask) > 1) then call msg_fatal ("Cannot decide EW parameter setting: More than one scheme set!") end if if (object%ew_parameter_mask (I_ALPHA)) call object%set_alpha_qed (model) if (object%ew_parameter_mask (I_GF)) call object%set_GF (model) if (object%ew_parameter_mask (I_SW2)) call object%set_weinberg_angle (model) end subroutine prc_blha_set_electroweak_parameters @ %def prc_blha_set_electrweak_parameters @ <>= procedure :: read_contract_file => prc_blha_read_contract_file <>= subroutine prc_blha_read_contract_file (object, flavors) class(prc_blha_t), intent(inout) :: object integer, intent(in), dimension(:,:) :: flavors integer, dimension(:), allocatable :: amp_type, flv_index, hel_index, label integer, dimension(:,:), allocatable :: helicities integer :: i_proc, i_hel allocate (helicities (N_MAX_FLAVORS, object%data%n_in)) select type (driver => object%driver) class is (blha_driver_t) call driver%read_contract_file (flavors, amp_type, flv_index, & hel_index, label, helicities) end select object%n_proc = count (amp_type >= 0) do i_proc = 1, object%n_proc if (amp_type (i_proc) < 0) exit if (hel_index(i_proc) < 0 .and. object%includes_polarization ()) & call msg_bug ("Object includes polarization, but helicity index is undefined.") i_hel = hel_index (i_proc) select case (amp_type (i_proc)) case (BLHA_AMP_TREE) if (allocated (object%i_tree)) then object%i_tree(flv_index(i_proc), i_hel) = label(i_proc) else call msg_fatal ("Tree matrix element present, & &but neither Born nor real indices are allocated!") end if case (BLHA_AMP_COLOR_C) if (allocated (object%i_color_c)) then object%i_color_c(flv_index(i_proc), i_hel) = label(i_proc) else call msg_fatal ("Color-correlated matrix element present, & &but cc-indices are not allocated!") end if case (BLHA_AMP_SPIN_C) if (allocated (object%i_spin_c)) then object%i_spin_c(flv_index(i_proc), i_hel) = label(i_proc) else call msg_fatal ("Spin-correlated matrix element present, & &but sc-indices are not allocated!") end if case (BLHA_AMP_LOOP) if (allocated (object%i_virt)) then object%i_virt(flv_index(i_proc), i_hel) = label(i_proc) else call msg_fatal ("Loop matrix element present, & &but virt-indices are not allocated!") end if case default call msg_fatal ("Undefined amplitude type") end select if (allocated (object%i_hel)) & object%i_hel (i_proc, :) = helicities (label(i_proc), :) end do end subroutine prc_blha_read_contract_file @ %def prc_blha_read_contract_file @ <>= procedure :: print_parameter_file => prc_blha_print_parameter_file <>= subroutine prc_blha_print_parameter_file (object, i_component) class(prc_blha_t), intent(in) :: object integer, intent(in) :: i_component type(string_t) :: filename select type (def => object%def) class is (blha_def_t) filename = def%basename // '_' // str (i_component) // '.olp_parameters' end select select type (driver => object%driver) class is (blha_driver_t) call driver%blha_olp_print_parameter (char(filename)//c_null_char) end select end subroutine prc_blha_print_parameter_file @ %def prc_blha_print_parameter_file @ <>= procedure :: compute_amplitude => prc_blha_compute_amplitude <>= function prc_blha_compute_amplitude & (object, j, p, f, h, c, fac_scale, ren_scale, alpha_qcd_forced, & core_state) result (amp) class(prc_blha_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 select type (core_state) class is (blha_state_t) core_state%alpha_qcd = object%qcd%alpha%get (fac_scale) end select amp = zero end function prc_blha_compute_amplitude @ @ %def prc_blha_compute_amplitude <>= procedure :: init_blha => prc_blha_init_blha <>= subroutine prc_blha_init_blha (object, blha_template, n_in, & n_particles, n_flv, n_hel) class(prc_blha_t), intent(inout) :: object type(blha_template_t), intent(in) :: blha_template integer, intent(in) :: n_in, n_particles, n_flv, n_hel object%n_particles = n_particles object%n_flv = n_flv object%n_hel = n_hel if (blha_template%compute_loop ()) then if (blha_template%include_polarizations) then allocate (object%i_virt (n_flv, n_hel), & object%i_color_c (n_flv, n_hel)) if (blha_template%use_internal_color_correlations) then allocate (object%i_hel (n_flv * n_in * n_hel * 2, n_in)) else allocate (object%i_hel (n_flv * n_in * n_hel, n_in)) end if else allocate (object%i_virt (n_flv, 1), object%i_color_c (n_flv, 1)) end if object%i_virt = -1 object%i_color_c = -1 else if (blha_template%compute_subtraction ()) then if (blha_template%include_polarizations) then allocate (object%i_tree (n_flv, n_hel), & object%i_color_c (n_flv, n_hel), & object%i_spin_c (n_flv, n_hel), & object%i_hel (3 * (n_flv * n_hel * n_in), n_in)) object%i_hel = 0 else allocate (object%i_tree (n_flv, 1), object%i_color_c (n_flv, 1) , & object%i_spin_c (n_flv, 1)) end if object%i_tree = -1 object%i_color_c = -1 object%i_spin_c = -1 else if (blha_template%compute_real_trees () .or. blha_template%compute_born () & .or. blha_template%compute_dglap ()) then if (blha_template%include_polarizations) then allocate (object%i_tree (n_flv, n_hel)) allocate (object%i_hel (n_flv * n_hel * n_in, n_in)) object%i_hel = 0 else allocate (object%i_tree (n_flv, 1)) end if object%i_tree = -1 end if call object%init_ew_parameters (blha_template%ew_scheme) select type (driver => object%driver) class is (blha_driver_t) driver%include_polarizations = blha_template%include_polarizations driver%switch_off_muon_yukawas = blha_template%switch_off_muon_yukawas driver%external_top_yukawa = blha_template%external_top_yukawa end select end subroutine prc_blha_init_blha @ %def prc_blha_init_blha @ <>= procedure :: set_mass_and_width => prc_blha_set_mass_and_width <>= subroutine prc_blha_set_mass_and_width (object, i_pdg, mass, width) class(prc_blha_t), intent(inout) :: object integer, intent(in) :: i_pdg real(default), intent(in) :: mass, width select type (driver => object%driver) class is (blha_driver_t) call driver%set_mass_and_width (i_pdg, mass, width) end select end subroutine prc_blha_set_mass_and_width @ %def prc_blha_set_mass_and_width @ <>= procedure :: set_particle_properties => prc_blha_set_particle_properties <>= subroutine prc_blha_set_particle_properties (object, model) class(prc_blha_t), intent(inout) :: object class(model_data_t), intent(in), target :: model integer :: i, i_pdg type(flavor_t) :: flv real(default) :: mass, width integer :: ierr real(default) :: top_yukawa do i = 1, OLP_N_MASSIVE_PARTICLES i_pdg = OLP_MASSIVE_PARTICLES(i) if (i_pdg < 0) cycle call flv%init (i_pdg, model) mass = flv%get_mass (); width = flv%get_width () select type (driver => object%driver) class is (blha_driver_t) call driver%set_mass_and_width (i_pdg, mass = mass, width = width) if (i_pdg == 5) call driver%blha_olp_set_parameter & ('yuk(5)'//c_null_char, dble(mass), 0._double, ierr) if (i_pdg == 6) then if (driver%external_top_yukawa > 0._default) then top_yukawa = driver%external_top_yukawa else top_yukawa = mass end if call driver%blha_olp_set_parameter & ('yuk(6)'//c_null_char, dble(top_yukawa), 0._double, ierr) end if if (driver%switch_off_muon_yukawas) then if (i_pdg == 13) call driver%blha_olp_set_parameter & ('yuk(13)' //c_null_char, 0._double, 0._double, ierr) end if end select end do end subroutine prc_blha_set_particle_properties @ %def prc_blha_set_particle_properties @ This mask adapts which electroweak parameters are supposed to set according to the chosen BLHA EWScheme. This is only implemented for the default OLP method so far. <>= procedure :: init_ew_parameters => prc_blha_init_ew_parameters <>= subroutine prc_blha_init_ew_parameters (object, ew_scheme) class(prc_blha_t), intent(inout) :: object integer, intent(in) :: ew_scheme object%ew_parameter_mask = .false. select case (ew_scheme) case (BLHA_EW_QED) object%ew_parameter_mask (I_ALPHA) = .true. case (BLHA_EW_GF) object%ew_parameter_mask (I_GF) = .true. end select end subroutine prc_blha_init_ew_parameters @ %def prc_blha_init_ew_parameters @ Computes a virtual matrix element from an interface to an external one-loop provider. The output of [[blha_olp_eval2]] is an array of [[dimension(4)]], corresponding to the $\epsilon^2$-, $\epsilon^1$- and $\epsilon^0$-poles of the virtual matrix element at position [[r(1:3)]] and the Born matrix element at position [[r(4)]]. The matrix element is rejected if its accuracy is larger than the maximal allowed accuracy. OpenLoops includes a factor of 1 / [[n_hel]] in the amplitudes, which we have to undo if polarized matrix elements are requested (GoSam does not support polarized matrix elements). <>= procedure :: compute_sqme_virt => prc_blha_compute_sqme_virt <>= subroutine prc_blha_compute_sqme_virt (object, & i_flv, i_hel, p, ren_scale, sqme, bad_point) class(prc_blha_t), intent(in) :: object integer, intent(in) :: i_flv, i_hel type(vector4_t), dimension(:), intent(in) :: p real(default), intent(in) :: ren_scale real(default), dimension(4), intent(out) :: sqme logical, intent(out) :: bad_point real(double), dimension(5 * object%n_particles) :: mom real(double), dimension(:), allocatable :: r real(double) :: mu_dble real(double) :: acc_dble real(default) :: acc real(default) :: alpha_s if (object%i_virt(i_flv, i_hel) >= 0) then allocate (r (blha_result_array_size (object%n_particles, BLHA_AMP_LOOP))) if (debug_on) call msg_debug2 (D_VIRTUAL, "prc_blha_compute_sqme_virt") if (debug_on) call msg_debug2 (D_VIRTUAL, "i_flv", i_flv) if (debug_on) call msg_debug2 (D_VIRTUAL, "object%i_virt(i_flv, i_hel)", object%i_virt(i_flv, i_hel)) if (debug2_active (D_VIRTUAL)) then call msg_debug2 (D_VIRTUAL, "use momenta: ") call vector4_write_set (p, show_mass = .true., & check_conservation = .true.) end if mom = object%create_momentum_array (p) if (vanishes (ren_scale)) & call msg_fatal ("prc_blha_compute_sqme_virt: ren_scale vanishes") mu_dble = dble(ren_scale) alpha_s = object%qcd%alpha%get (ren_scale) select type (driver => object%driver) class is (blha_driver_t) call driver%set_alpha_s (alpha_s) call driver%blha_olp_eval2 (object%i_virt(i_flv, i_hel), mom, mu_dble, r, acc_dble) end select acc = acc_dble sqme = r(1:4) bad_point = acc > object%maximum_accuracy if (object%includes_polarization ()) sqme = object%n_hel * sqme else sqme = zero end if end subroutine prc_blha_compute_sqme_virt @ %def prc_blha_compute_sqme_virt @ Computes a tree-level matrix element from an interface to an external one-loop provider. The matrix element is rejected if its accuracy is larger than the maximal allowed accuracy. OpenLoops includes a factor of 1 / [[n_hel]] in the amplitudes, which we have to undo if polarized matrix elements are requested (GoSam does not support polarized matrix elements). <>= procedure :: compute_sqme => prc_blha_compute_sqme <>= subroutine prc_blha_compute_sqme (object, i_flv, i_hel, p, & ren_scale, sqme, bad_point) class(prc_blha_t), intent(in) :: object integer, intent(in) :: i_flv, i_hel type(vector4_t), intent(in), dimension(:) :: p real(default), intent(in) :: ren_scale real(default), intent(out) :: sqme logical, intent(out) :: bad_point real(double), dimension(5*object%n_particles) :: mom real(double), dimension(OLP_RESULTS_LIMIT) :: r real(double) :: mu_dble, acc_dble real(default) :: acc, alpha_s if (object%i_tree(i_flv, i_hel) >= 0) then mom = object%create_momentum_array (p) if (vanishes (ren_scale)) & call msg_fatal ("prc_blha_compute_sqme: ren_scale vanishes") mu_dble = dble(ren_scale) alpha_s = object%qcd%alpha%get (ren_scale) select type (driver => object%driver) class is (blha_driver_t) call driver%set_alpha_s (alpha_s) call driver%blha_olp_eval2 (object%i_tree(i_flv, i_hel), mom, & mu_dble, r, acc_dble) sqme = r(object%sqme_tree_pos) end select acc = acc_dble bad_point = acc > object%maximum_accuracy if (object%includes_polarization ()) sqme = object%n_hel * sqme else sqme = zero end if end subroutine prc_blha_compute_sqme @ %def prc_blha_compute_sqme @ <>= public :: blha_color_c_fill_diag <>= subroutine blha_color_c_fill_diag (sqme_born, flavors, sqme_color_c) real(default), intent(in) :: sqme_born integer, intent(in), dimension(:) :: flavors real(default), intent(inout), dimension(:,:) :: sqme_color_c integer :: i do i = 1, size (flavors) if (is_quark (flavors(i))) then sqme_color_c (i, i) = -cf * sqme_born else if (is_gluon (flavors(i))) then sqme_color_c (i, i) = -ca * sqme_born else sqme_color_c (i, i) = zero end if end do end subroutine blha_color_c_fill_diag @ %def blha_color_c_fill_diag <>= public :: blha_color_c_fill_offdiag <>= subroutine blha_color_c_fill_offdiag (n, r, sqme_color_c, offset, n_flv) integer, intent(in) :: n real(default), intent(in), dimension(:) :: r real(default), intent(inout), dimension(:,:) :: sqme_color_c integer, intent(in), optional :: offset, n_flv integer :: i, j, pos, incr if (present (offset)) then incr = offset else incr = 0 end if pos = 0 do j = 1, n do i = 1, j if (i /= j) then pos = (j - 1) * (j - 2) / 2 + i if (present (n_flv)) incr = incr + n_flv - 1 if (present (offset)) pos = pos + incr sqme_color_c (i, j) = -r (pos) sqme_color_c (j, i) = sqme_color_c (i, j) end if end do end do end subroutine blha_color_c_fill_offdiag @ %def blha_color_c_fill_offdiag @ Computes a color-correlated matrix element from an interface to an external one-loop provider. The output of [[blha_olp_eval2]] is an array of [[dimension(n * (n - 1) / 2)]]. The matrix element is rejected if its accuracy is larger than the maximal allowed accuracy. OpenLoops includes a factor of 1 / [[n_hel]] in the amplitudes, which we have to undo if polarized matrix elements are requested (GoSam does not support polarized matrix elements). <>= procedure :: compute_sqme_color_c_raw => prc_blha_compute_sqme_color_c_raw <>= subroutine prc_blha_compute_sqme_color_c_raw & (object, i_flv, i_hel, p, ren_scale, rr, bad_point) class(prc_blha_t), intent(in) :: object integer, intent(in) :: i_flv, i_hel type(vector4_t), intent(in), dimension(:) :: p real(default), intent(in) :: ren_scale real(default), intent(out), dimension(:) :: rr logical, intent(out) :: bad_point real(double), dimension(5 * object%n_particles) :: mom real(double), dimension(size(rr)) :: r real(default) :: alpha_s, acc real(double) :: mu_dble, acc_dble if (object%i_color_c(i_flv, i_hel) >= 0) then mom = object%create_momentum_array (p) if (vanishes (ren_scale)) & call msg_fatal ("prc_blha_compute_sqme_color_c: ren_scale vanishes") mu_dble = dble(ren_scale) alpha_s = object%qcd%alpha%get (ren_scale) select type (driver => object%driver) class is (blha_driver_t) call driver%set_alpha_s (alpha_s) call driver%blha_olp_eval2 (object%i_color_c(i_flv, i_hel), & mom, mu_dble, r, acc_dble) end select rr = r acc = acc_dble bad_point = acc > object%maximum_accuracy if (object%includes_polarization ()) rr = object%n_hel * rr else rr = zero end if end subroutine prc_blha_compute_sqme_color_c_raw @ %def prc_blha_compute_sqme_color_c_raw @ <>= procedure :: compute_sqme_color_c => prc_blha_compute_sqme_color_c <>= subroutine prc_blha_compute_sqme_color_c & (object, i_flv, i_hel, p, ren_scale, born_color_c, bad_point, born_out) class(prc_blha_t), intent(inout) :: object integer, intent(in) :: i_flv, i_hel type(vector4_t), intent(in), dimension(:) :: p real(default), intent(in) :: ren_scale real(default), intent(inout), dimension(:,:) :: born_color_c real(default), intent(out), optional :: born_out logical, intent(out) :: bad_point real(default), dimension(:), allocatable :: r logical :: bad_point2 real(default) :: born integer, dimension(:), allocatable :: flavors allocate (r (blha_result_array_size & (size(born_color_c, dim=1), BLHA_AMP_COLOR_C))) call object%compute_sqme_color_c_raw (i_flv, i_hel, p, ren_scale, r, bad_point) select type (driver => object%driver) class is (blha_driver_t) if (allocated (object%i_tree)) then call object%compute_sqme (i_flv, i_hel, p, ren_scale, born, bad_point2) else born = zero end if if (present (born_out)) born_out = born end select call blha_color_c_fill_offdiag (object%n_particles, r, born_color_c) flavors = object%get_flv_state (i_flv) call blha_color_c_fill_diag (born, flavors, born_color_c) bad_point = bad_point .or. bad_point2 end subroutine prc_blha_compute_sqme_color_c @ %def prc_blha_compute_sqme_color_c @ <>= generic :: get_beam_helicities => get_beam_helicities_single generic :: get_beam_helicities => get_beam_helicities_array procedure :: get_beam_helicities_single => prc_blha_get_beam_helicities_single procedure :: get_beam_helicities_array => prc_blha_get_beam_helicities_array <>= function prc_blha_get_beam_helicities_single (object, i, invert_second) result (hel) integer, dimension(:), allocatable :: hel class(prc_blha_t), intent(in) :: object logical, intent(in), optional :: invert_second integer, intent(in) :: i logical :: inv inv = .false.; if (present (invert_second)) inv = invert_second allocate (hel (object%data%n_in)) hel = object%i_hel (i, :) if (inv .and. object%data%n_in == 2) hel(2) = -hel(2) end function prc_blha_get_beam_helicities_single @ %def prc_blha_get_beam_helicities_single @ <>= procedure :: includes_polarization => prc_blha_includes_polarization <>= function prc_blha_includes_polarization (object) result (polarized) logical :: polarized class(prc_blha_t), intent(in) :: object select type (driver => object%driver) class is (blha_driver_t) polarized = driver%include_polarizations end select end function prc_blha_includes_polarization @ %def prc_blha_includes_polarization @ <>= function prc_blha_get_beam_helicities_array (object, invert_second) result (hel) integer, dimension(:,:), allocatable :: hel class(prc_blha_t), intent(in) :: object logical, intent(in), optional :: invert_second integer :: i allocate (hel (object%n_proc, object%data%n_in)) do i = 1, object%n_proc hel(i,:) = object%get_beam_helicities (i, invert_second) end do end function prc_blha_get_beam_helicities_array @ %def prc_blha_get_beam_helicities_array @ <>= procedure(prc_blha_init_driver), deferred :: & init_driver <>= abstract interface subroutine prc_blha_init_driver (object, os_data) import class(prc_blha_t), intent(inout) :: object type(os_data_t), intent(in) :: os_data end subroutine prc_blha_init_driver end interface @ %def prc_blha_init_driver interface @ In general, the BLHA consits of a virtual matrix element and $n_{\rm{sub}}$ subtraction terms. The subtractions terms can be pure Born matrix elements (to be used in collinear subtraction or in internal color-correlation), color-correlated matrix elements or spin-correlated matrix elements. The numbers should be ordered in such a way that $\mathcal{V}_{\rm{fin}}$ is first, followed by the pure Born, the color-correlated and the spin-correlated matrix elements. This repeats $n_{\rm{flv}}$ times. Let $\nu_i$ be the position of the $ith$ virtual matrix element. The next $\mathcal{V}_{\rm{fin}}$ is at position $\nu_i = \nu_{i - 1} + n_{\rm{sub}} + 1$. Obviously, $\nu_1 = 1$. This allows us to determine the virtual matrix element positions using the recursive function implemented below. <>= public :: blha_loop_positions <>= recursive function blha_loop_positions (i_flv, n_sub) result (index) integer :: index integer, intent(in) :: i_flv, n_sub index = 0 if (i_flv == 1) then index = 1 else index = blha_loop_positions (i_flv - 1, n_sub) + n_sub + 1 end if end function blha_loop_positions @ %def blha_loop_positions @ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% The module is split into a configuration interface which manages configuration and handles the request and contract files, a module which interfaces the OLP matrix elements and a driver. <<[[blha_config.f90]]>>= <> module blha_config use kinds <> use io_units use constants use string_utils use variables, only: var_list_t use diagnostics use md5 use model_data use flavors use quantum_numbers use pdg_arrays use sorting use lexers use parser use syntax_rules use ifiles use beam_structures, only: beam_structure_t <> <> <> <> <> <> <> contains <> end module blha_config @ %def blha_config @ \section{Configuration} Parameters to enumerate the different options in the order. <>= integer, public, parameter :: & BLHA_CT_QCD = 1, BLHA_CT_EW = 2, BLHA_CT_QED = 3, BLHA_CT_OTHER = 4 integer, public, parameter :: & BLHA_IRREG_CDR = 1, BLHA_IRREG_DRED = 2, BLHA_IRREG_THV = 3, & BLHA_IRREG_MREG = 4, BLHA_IRREG_OTHER = 5 integer, public, parameter :: & BLHA_MPS_ONSHELL = 1, BLHA_MPS_OTHER = 2 integer, public, parameter :: & BLHA_MODE_GOSAM = 1, BLHA_MODE_FEYNARTS = 2, BLHA_MODE_GENERIC = 3, & BLHA_MODE_OPENLOOPS = 4 integer, public, parameter :: & BLHA_VERSION_1 = 1, BLHA_VERSION_2 = 2 integer, public, parameter :: & BLHA_AMP_LOOP = 1, BLHA_AMP_COLOR_C = 2, BLHA_AMP_SPIN_C = 3, & BLHA_AMP_TREE = 4, BLHA_AMP_LOOPINDUCED = 5 integer, public, parameter :: & BLHA_EW_QED = 0, & BLHA_EW_GF = 1, BLHA_EW_MZ = 2, BLHA_EW_MSBAR = 3, & BLHA_EW_0 = 4, BLHA_EW_RUN = 5 integer, public, parameter :: & BLHA_WIDTH_COMPLEX = 1, BLHA_WIDTH_FIXED = 2, & BLHA_WIDTH_RUNNING = 3, BLHA_WIDTH_POLE = 4, & BLHA_WIDTH_DEFAULT = 5 @ %def blha_ct_qcd blha_ct_ew blha_ct_qed blha_ct_other @ %def blha_irreg_cdr blha_irreg_dred blha_irreg_thv blha_irreg_mreg blha_irreg_other @ %def blha_mps_onshell blha_mps_other @ %def blha_mode_gosam blha_mode_feynarts blha_mode_generic @ %def blha version blha_amp blha_ew blha_width @ Those are the default pdg codes for massive particles in BLHA programs <>= integer, parameter, public :: OLP_N_MASSIVE_PARTICLES = 12 integer, dimension(OLP_N_MASSIVE_PARTICLES), public :: & OLP_MASSIVE_PARTICLES = [5, -5, 6, -6, 13, -13, 15, -15, 23, 24, -24, 25] integer, parameter :: OLP_HEL_UNPOLARIZED = 0 @ %def OLP_MASSIVE_PARTICLES @ The user might provide an extra command string for OpenLoops to apply special libraries instead of the default ones, such as signal-only amplitudes for off-shell top production. We check in this subroutine that the provided string is valid and print out the possible options to ease the user's memory. <>= integer, parameter :: N_KNOWN_SPECIAL_OL_METHODS = 3 <>= subroutine check_extra_cmd (extra_cmd) type(string_t), intent(in) :: extra_cmd type(string_t), dimension(N_KNOWN_SPECIAL_OL_METHODS) :: known_methods integer :: i logical :: found known_methods(1) = 'top' known_methods(2) = 'not' known_methods(3) = 'stop' if (extra_cmd == var_str ("")) return found = .false. do i = 1, N_KNOWN_SPECIAL_OL_METHODS found = found .or. & (extra_cmd == var_str ('extra approx ') // known_methods(i)) end do if (.not. found) & call msg_fatal ("The given extra OpenLoops method is not kown ", & [var_str ("Available commands are: "), & var_str ("extra approx top (only WbWb signal),"), & var_str ("extra approx stop (only WbWb singletop),"), & var_str ("extra approx not (no top in WbWb).")]) end subroutine check_extra_cmd @ %def check_extra_cmd @ This type contains the pdg code of the particle to be written in the process specification string and an optional additional information about the polarization of the particles. Note that the output can only be processed by OpenLoops. <>= type :: blha_particle_string_element_t integer :: pdg = 0 integer :: hel = OLP_HEL_UNPOLARIZED logical :: polarized = .false. contains <> end type blha_particle_string_element_t @ %def blha_particle_string_element_t @ <>= generic :: init => init_default generic :: init => init_polarized procedure :: init_default => blha_particle_string_element_init_default procedure :: init_polarized => blha_particle_string_element_init_polarized <>= subroutine blha_particle_string_element_init_default (blha_p, id) class(blha_particle_string_element_t), intent(out) :: blha_p integer, intent(in) :: id blha_p%pdg = id end subroutine blha_particle_string_element_init_default @ %def blha_particle_string_element_init_default @ <>= subroutine blha_particle_string_element_init_polarized (blha_p, id, hel) class(blha_particle_string_element_t), intent(out) :: blha_p integer, intent(in) :: id, hel blha_p%polarized = .true. blha_p%pdg = id blha_p%hel = hel end subroutine blha_particle_string_element_init_polarized @ %def blha_particle_string_element_init_polarized @ <>= generic :: write_pdg => write_pdg_unit generic :: write_pdg => write_pdg_character procedure :: write_pdg_unit => blha_particle_string_element_write_pdg_unit procedure :: write_pdg_character & => blha_particle_string_element_write_pdg_character <>= subroutine blha_particle_string_element_write_pdg_unit (blha_p, unit) class(blha_particle_string_element_t), intent(in) :: blha_p integer, intent(in), optional :: unit integer :: u u = given_output_unit (unit) write (u, '(I3)') blha_p%pdg end subroutine blha_particle_string_element_write_pdg_unit @ %def blha_particle_string_element_write_pdg_unit @ <>= subroutine blha_particle_string_element_write_pdg_character (blha_p, c) class(blha_particle_string_element_t), intent(in) :: blha_p character(3), intent(inout) :: c write (c, '(I3)') blha_p%pdg end subroutine blha_particle_string_element_write_pdg_character @ %def blha_particle_string_element_write_pdg_character @ <>= generic :: write_helicity => write_helicity_unit generic :: write_helicity => write_helicity_character procedure :: write_helicity_unit & => blha_particle_string_element_write_helicity_unit procedure :: write_helicity_character & => blha_particle_string_element_write_helicity_character <>= subroutine blha_particle_string_element_write_helicity_unit (blha_p, unit) class(blha_particle_string_element_t), intent(in) :: blha_p integer, intent(in), optional :: unit integer :: u u = given_output_unit (unit) write (u, '(A1,I0,A1)') '(', blha_p%hel, ')' end subroutine blha_particle_string_element_write_helicity_unit @ %def blha_particle_string_element_write_helicity_unit @ <>= subroutine blha_particle_string_element_write_helicity_character (blha_p, c) class(blha_particle_string_element_t), intent(in) :: blha_p character(4), intent(inout) :: c write (c, '(A1,I0,A1)') '(', blha_p%hel, ')' end subroutine blha_particle_string_element_write_helicity_character @ %def blha_particle_string_element_write_helicity_character @ This type encapsulates a BLHA request. <>= public :: blha_configuration_t public :: blha_cfg_process_node_t <>= type :: blha_cfg_process_node_t type(blha_particle_string_element_t), dimension(:), allocatable :: pdg_in, pdg_out integer, dimension(:), allocatable :: fingerprint integer :: nsub integer, dimension(:), allocatable :: ids integer :: amplitude_type type(blha_cfg_process_node_t), pointer :: next => null () end type blha_cfg_process_node_t type :: blha_configuration_t type(string_t) :: name class(model_data_t), pointer :: model => null () type(string_t) :: md5 integer :: version = 2 logical :: dirty = .false. integer :: n_proc = 0 real(default) :: accuracy_target logical :: debug_unstable = .false. integer :: mode = BLHA_MODE_GENERIC logical :: polarized = .false. type(blha_cfg_process_node_t), pointer :: processes => null () !integer, dimension(2) :: matrix_element_square_type = BLHA_MEST_SUM integer :: correction_type type(string_t) :: correction_type_other integer :: irreg = BLHA_IRREG_THV type(string_t) :: irreg_other integer :: massive_particle_scheme = BLHA_MPS_ONSHELL type(string_t) :: massive_particle_scheme_other type(string_t) :: model_file logical :: subdivide_subprocesses = .false. integer :: alphas_power = -1, alpha_power = -1 integer :: ew_scheme = BLHA_EW_GF integer :: width_scheme = BLHA_WIDTH_DEFAULT logical :: openloops_use_cms = .false. integer :: openloops_phs_tolerance = 0 type(string_t) :: openloops_extra_cmd integer :: openloops_stability_log = 0 logical :: openloops_use_collier = .false. end type blha_configuration_t @ %def blha_cffg_process_node_t blha_configuration_t @ Translate the SINDARIN input string to the corresponding named integer. <>= public :: ew_scheme_string_to_int <>= function ew_scheme_string_to_int (ew_scheme_str) result (ew_scheme_int) integer :: ew_scheme_int type(string_t), intent(in) :: ew_scheme_str select case (char (ew_scheme_str)) case ('GF', 'Gmu') ew_scheme_int = BLHA_EW_GF case ('alpha_qed') ew_scheme_int = BLHA_EW_QED case ('alpha_mz') ew_scheme_int = BLHA_EW_MZ case ('alpha_0', 'alpha_thompson') ew_scheme_int = BLHA_EW_0 case default call msg_fatal ("ew_scheme: " // char (ew_scheme_str) // & " not supported. Try 'Gmu', 'alpha_qed', 'alpha_mz' or 'alpha_0'.") end select end function ew_scheme_string_to_int @ %def ew_scheme_string_to_int @ @ Translate the SINDARIN input string to the corresponding named integer. <>= public :: correction_type_string_to_int <>= function correction_type_string_to_int (correction_type_str) result (correction_type_int) integer :: correction_type_int type(string_t), intent(in) :: correction_type_str select case (char (correction_type_str)) case ('QCD') correction_type_int = BLHA_CT_QCD case ('QED', 'EW') correction_type_int = BLHA_CT_EW case default call msg_warning ("nlo_correction_type: " // char (correction_type_str) // & " not supported. Try setting it to 'QCD', 'QED' or 'EW'.") end select end function correction_type_string_to_int @ %def correction_type_string_to_int @ This types control the creation of BLHA-interface files <>= public :: blha_flv_state_t public :: blha_master_t <>= type:: blha_flv_state_t integer, dimension(:), allocatable :: flavors integer :: flv_mult logical :: flv_real = .false. end type blha_flv_state_t type :: blha_master_t integer, dimension(5) :: blha_mode = BLHA_MODE_GENERIC logical :: compute_borns = .false. logical :: compute_real_trees = .false. logical :: compute_loops = .true. logical :: compute_correlations = .false. logical :: compute_dglap = .false. integer :: ew_scheme type(string_t), dimension(:), allocatable :: suffix type(blha_configuration_t), dimension(:), allocatable :: blha_cfg integer :: n_files = 0 integer, dimension(:), allocatable :: i_file_to_nlo_index contains <> end type blha_master_t @ %def blha_flv_state_t, blha_master_t @ Master-Routines <>= procedure :: set_methods => blha_master_set_methods <>= subroutine blha_master_set_methods (master, is_nlo, var_list) class(blha_master_t), intent(inout) :: master logical, intent(in) :: is_nlo type(var_list_t), intent(in) :: var_list type(string_t) :: method, born_me_method, real_tree_me_method type(string_t) :: loop_me_method, correlation_me_method type(string_t) :: dglap_me_method type(string_t) :: default_method logical :: cmp_born, cmp_real logical :: cmp_loop, cmp_corr logical :: cmp_dglap if (is_nlo) then method = var_list%get_sval (var_str ("$method")) born_me_method = var_list%get_sval (var_str ("$born_me_method")) if (born_me_method == "") born_me_method = method real_tree_me_method = var_list%get_sval (var_str ("$real_tree_me_method")) if (real_tree_me_method == "") real_tree_me_method = method loop_me_method = var_list%get_sval (var_str ("$loop_me_method")) if (loop_me_method == "") loop_me_method = method correlation_me_method = var_list%get_sval (var_str ("$correlation_me_method")) if (correlation_me_method == "") correlation_me_method = method dglap_me_method = var_list%get_sval (var_str ("$dglap_me_method")) if (dglap_me_method == "") dglap_me_method = method cmp_born = born_me_method /= 'omega' cmp_real = is_nlo .and. (real_tree_me_method /= 'omega') cmp_loop = is_nlo .and. (loop_me_method /= 'omega') cmp_corr = is_nlo .and. (correlation_me_method /= 'omega') cmp_dglap = is_nlo .and. (dglap_me_method /= 'omega') call set_me_method (1, loop_me_method) call set_me_method (2, correlation_me_method) call set_me_method (3, real_tree_me_method) call set_me_method (4, born_me_method) call set_me_method (5, dglap_me_method) else default_method = var_list%get_sval (var_str ("$method")) cmp_born = default_method /= 'omega' cmp_real = .false.; cmp_loop = .false.; cmp_corr = .false. call set_me_method (4, default_method) end if master%n_files = count ([cmp_born, cmp_real, cmp_loop, cmp_corr, cmp_dglap]) call set_nlo_indices () master%compute_borns = cmp_born master%compute_real_trees = cmp_real master%compute_loops = cmp_loop master%compute_correlations = cmp_corr master%compute_dglap = cmp_dglap contains subroutine set_nlo_indices () integer :: i_file allocate (master%i_file_to_nlo_index (master%n_files)) master%i_file_to_nlo_index = 0 i_file = 0 if (cmp_loop) then i_file = i_file + 1 master%i_file_to_nlo_index(i_file) = 1 end if if (cmp_corr) then i_file = i_file + 1 master%i_file_to_nlo_index(i_file) = 2 end if if (cmp_real) then i_file = i_file + 1 master%i_file_to_nlo_index(i_file) = 3 end if if (cmp_born) then i_file = i_file + 1 master%i_file_to_nlo_index(i_file) = 4 end if if (cmp_dglap) then i_file = i_file + 1 master%i_file_to_nlo_index(i_file) = 5 end if end subroutine set_nlo_indices subroutine set_me_method (i, me_method) integer, intent(in) :: i type(string_t) :: me_method select case (char (me_method)) case ('gosam') call master%set_gosam (i) case ('openloops') call master%set_openloops (i) end select end subroutine set_me_method end subroutine blha_master_set_methods @ %def blha_master_set_methods @ <>= procedure :: allocate_config_files => blha_master_allocate_config_files <>= subroutine blha_master_allocate_config_files (master) class(blha_master_t), intent(inout) :: master allocate (master%blha_cfg (master%n_files)) allocate (master%suffix (master%n_files)) end subroutine blha_master_allocate_config_files @ %def blha_master_allocate_config_files @ <>= procedure :: set_ew_scheme => blha_master_set_ew_scheme <>= subroutine blha_master_set_ew_scheme (master, ew_scheme) class(blha_master_t), intent(inout) :: master type(string_t), intent(in) :: ew_scheme master%ew_scheme = ew_scheme_string_to_int (ew_scheme) end subroutine blha_master_set_ew_scheme @ %def blha_master_set_ew_scheme @ <>= procedure :: set_correction_type => blha_master_set_correction_type <>= subroutine blha_master_set_correction_type (master, correction_type_str) class(blha_master_t), intent(inout) :: master type(string_t), intent(in) :: correction_type_str master%blha_cfg(:)%correction_type = correction_type_string_to_int (correction_type_str) end subroutine blha_master_set_correction_type @ %def blha_master_set_correction_type @ <>= procedure :: generate => blha_master_generate <>= subroutine blha_master_generate (master, basename, model, & n_in, alpha_power, alphas_power, flv_born, flv_real) class(blha_master_t), intent(inout) :: master type(string_t), intent(in) :: basename class(model_data_t), intent(in), target :: model integer, intent(in) :: n_in integer, intent(in) :: alpha_power, alphas_power integer, intent(in), dimension(:,:), allocatable :: flv_born, flv_real integer :: i_file if (master%n_files < 1) & call msg_fatal ("Attempting to generate OLP-files, but none are specified!") i_file = 1 call master%generate_loop (basename, model, n_in, alpha_power, & alphas_power, flv_born, i_file) call master%generate_correlation (basename, model, n_in, alpha_power, & alphas_power, flv_born, i_file) call master%generate_real_tree (basename, model, n_in, alpha_power, & alphas_power, flv_real, i_file) call master%generate_born (basename, model, n_in, alpha_power, & alphas_power, flv_born, i_file) call master%generate_dglap (basename, model, n_in, alpha_power, & alphas_power, flv_born, i_file) end subroutine blha_master_generate @ %def blha_master_generate @ <>= procedure :: generate_loop => blha_master_generate_loop <>= subroutine blha_master_generate_loop (master, basename, model, n_in, & alpha_power, alphas_power, flv_born, i_file) class(blha_master_t), intent(inout) :: master type(string_t), intent(in) :: basename class(model_data_t), intent(in), target :: model integer, intent(in) :: n_in integer, intent(in) :: alpha_power, alphas_power integer, dimension(:,:), allocatable, intent(in) :: flv_born integer, intent(inout) :: i_file type(blha_flv_state_t), dimension(:), allocatable :: blha_flavor integer :: i_flv if (master%compute_loops) then if (allocated (flv_born)) then allocate (blha_flavor (size (flv_born, 2))) do i_flv = 1, size (flv_born, 2) allocate (blha_flavor(i_flv)%flavors (size (flv_born(:,i_flv)))) blha_flavor(i_flv)%flavors = flv_born(:,i_flv) blha_flavor(i_flv)%flv_mult = 2 end do master%suffix(i_file) = blha_get_additional_suffix (var_str ("_LOOP")) call blha_init_virtual (master%blha_cfg(i_file), blha_flavor, & n_in, alpha_power, alphas_power, master%ew_scheme, & basename, model, master%blha_mode(1), master%suffix(i_file)) i_file = i_file + 1 else call msg_fatal ("BLHA Loops requested but " & // "Born flavor not existing") end if end if end subroutine blha_master_generate_loop @ %def blha_master_generate_loop @ <>= procedure :: generate_correlation => blha_master_generate_correlation <>= subroutine blha_master_generate_correlation (master, basename, model, n_in, & alpha_power, alphas_power, flv_born, i_file) class(blha_master_t), intent(inout) :: master type(string_t), intent(in) :: basename class(model_data_t), intent(in), target :: model integer, intent(in) :: n_in integer, intent(in) :: alpha_power, alphas_power integer, dimension(:,:), allocatable, intent(in) :: flv_born integer, intent(inout) :: i_file type(blha_flv_state_t), dimension(:), allocatable :: blha_flavor integer :: i_flv if (master%compute_correlations) then if (allocated (flv_born)) then allocate (blha_flavor (size (flv_born, 2))) do i_flv = 1, size (flv_born, 2) allocate (blha_flavor(i_flv)%flavors (size (flv_born(:,i_flv)))) blha_flavor(i_flv)%flavors = flv_born(:,i_flv) blha_flavor(i_flv)%flv_mult = 3 end do master%suffix(i_file) = blha_get_additional_suffix (var_str ("_SUB")) call blha_init_subtraction (master%blha_cfg(i_file), blha_flavor, & n_in, alpha_power, alphas_power, master%ew_scheme, & basename, model, master%blha_mode(2), master%suffix(i_file)) i_file = i_file + 1 else call msg_fatal ("BLHA Correlations requested but "& // "Born flavor not existing") end if end if end subroutine blha_master_generate_correlation @ %def blha_master_generate_correlation @ <>= procedure :: generate_real_tree => blha_master_generate_real_tree <>= subroutine blha_master_generate_real_tree (master, basename, model, n_in, & alpha_power, alphas_power, flv_real, i_file) class(blha_master_t), intent(inout) :: master type(string_t), intent(in) :: basename class(model_data_t), intent(in), target :: model integer, intent(in) :: n_in integer, intent(in) :: alpha_power, alphas_power integer, dimension(:,:), allocatable, intent(in) :: flv_real integer, intent(inout) :: i_file type(blha_flv_state_t), dimension(:), allocatable :: blha_flavor integer :: i_flv if (master%compute_real_trees) then if (allocated (flv_real)) then allocate (blha_flavor (size (flv_real, 2))) do i_flv = 1, size (flv_real, 2) allocate (blha_flavor(i_flv)%flavors (size (flv_real(:,i_flv)))) blha_flavor(i_flv)%flavors = flv_real(:,i_flv) blha_flavor(i_flv)%flv_mult = 1 end do master%suffix(i_file) = blha_get_additional_suffix (var_str ("_REAL")) call blha_init_real (master%blha_cfg(i_file), blha_flavor, & n_in, alpha_power, alphas_power, master%ew_scheme, & basename, model, master%blha_mode(3), master%suffix(i_file)) i_file = i_file + 1 else call msg_fatal ("BLHA Trees requested but "& // "Real flavor not existing") end if end if end subroutine blha_master_generate_real_tree @ %def blha_master_generate_real_tree @ <>= procedure :: generate_born => blha_master_generate_born <>= subroutine blha_master_generate_born (master, basename, model, n_in, & alpha_power, alphas_power, flv_born, i_file) class(blha_master_t), intent(inout) :: master type(string_t), intent(in) :: basename class(model_data_t), intent(in), target :: model integer, intent(in) :: n_in integer, intent(in) :: alpha_power, alphas_power integer, dimension(:,:), allocatable, intent(in) :: flv_born integer, intent(inout) :: i_file type(blha_flv_state_t), dimension(:), allocatable :: blha_flavor integer :: i_flv if (master%compute_borns) then if (allocated (flv_born)) then allocate (blha_flavor (size (flv_born, 2))) do i_flv = 1, size (flv_born, 2) allocate (blha_flavor(i_flv)%flavors (size (flv_born(:,i_flv)))) blha_flavor(i_flv)%flavors = flv_born(:,i_flv) blha_flavor(i_flv)%flv_mult = 1 end do master%suffix(i_file) = blha_get_additional_suffix (var_str ("_BORN")) call blha_init_born (master%blha_cfg(i_file), blha_flavor, & n_in, alpha_power, alphas_power, master%ew_scheme, & basename, model, master%blha_mode(4), master%suffix(i_file)) i_file = i_file + 1 end if end if end subroutine blha_master_generate_born @ %def blha_master_generate_born @ <>= procedure :: generate_dglap => blha_master_generate_dglap <>= subroutine blha_master_generate_dglap (master, basename, model, n_in, & alpha_power, alphas_power, flv_born, i_file) class(blha_master_t), intent(inout) :: master type(string_t), intent(in) :: basename class(model_data_t), intent(in), target :: model integer, intent(in) :: n_in integer, intent(in) :: alpha_power, alphas_power integer, dimension(:,:), allocatable, intent(in) :: flv_born integer, intent(inout) :: i_file type(blha_flv_state_t), dimension(:), allocatable :: blha_flavor integer :: i_flv if (master%compute_dglap) then if (allocated (flv_born)) then allocate (blha_flavor (size (flv_born, 2))) do i_flv = 1, size (flv_born, 2) allocate (blha_flavor(i_flv)%flavors (size (flv_born(:,i_flv)))) blha_flavor(i_flv)%flavors = flv_born(:,i_flv) blha_flavor(i_flv)%flv_mult = 1 end do master%suffix(i_file) = blha_get_additional_suffix (var_str ("_DGLAP")) call blha_init_born (master%blha_cfg(i_file), blha_flavor, & n_in, alpha_power, alphas_power, master%ew_scheme, & basename, model, master%blha_mode(5), master%suffix(i_file)) i_file = i_file + 1 end if end if end subroutine blha_master_generate_dglap @ %def blha_master_generate_dglap @ <>= procedure :: setup_additional_features => blha_master_setup_additional_features <>= subroutine blha_master_setup_additional_features (master, & phs_tolerance, use_cms, stability_log, use_collier, extra_cmd, beam_structure) class(blha_master_t), intent(inout) :: master integer, intent(in) :: phs_tolerance logical, intent(in) :: use_cms type(string_t), intent(in), optional :: extra_cmd integer, intent(in) :: stability_log logical, intent(in) :: use_collier type(beam_structure_t), intent(in), optional :: beam_structure integer :: i_file logical :: polarized, throw_warning polarized = .false. if (present (beam_structure)) polarized = beam_structure%has_polarized_beams () throw_warning = .false. if (use_cms) then throw_warning = throw_warning .or. (master%compute_loops & .and. master%blha_mode(1) /= BLHA_MODE_OPENLOOPS) throw_warning = throw_warning .or. (master%compute_correlations & .and. master%blha_mode(2) /= BLHA_MODE_OPENLOOPS) throw_warning = throw_warning .or. (master%compute_real_trees & .and. master%blha_mode(3) /= BLHA_MODE_OPENLOOPS) throw_warning = throw_warning .or. (master%compute_borns & .and. master%blha_mode(4) /= BLHA_MODE_OPENLOOPS) throw_warning = throw_warning .or. (master%compute_dglap & .and. master%blha_mode(5) /= BLHA_MODE_OPENLOOPS) if (throw_warning) call cms_warning () end if do i_file = 1, master%n_files if (phs_tolerance > 0) then select case (master%blha_mode (master%i_file_to_nlo_index(i_file))) case (BLHA_MODE_GOSAM) if (polarized) call gosam_error_message () case (BLHA_MODE_OPENLOOPS) master%blha_cfg(i_file)%openloops_use_cms = use_cms master%blha_cfg(i_file)%openloops_phs_tolerance = phs_tolerance master%blha_cfg(i_file)%polarized = polarized if (present (extra_cmd)) then master%blha_cfg(i_file)%openloops_extra_cmd = extra_cmd else master%blha_cfg(i_file)%openloops_extra_cmd = var_str ('') end if master%blha_cfg(i_file)%openloops_stability_log = stability_log master%blha_cfg(i_file)%openloops_use_collier = use_collier end select end if end do contains subroutine cms_warning () call msg_warning ("You have set ?openloops_use_cms = true, but not all active matrix ", & [var_str ("element methods are set to OpenLoops. Note that other "), & var_str ("methods might not necessarily support the complex mass "), & var_str ("scheme. This can yield inconsistencies in your NLO results!")]) end subroutine cms_warning subroutine gosam_error_message () call msg_fatal ("You are trying to evaluate a process at NLO ", & [var_str ("which involves polarized beams using GoSam. "), & var_str ("This feature is not supported yet. "), & var_str ("Please use OpenLoops instead")]) end subroutine gosam_error_message end subroutine blha_master_setup_additional_features @ %def blha_master_setup_additional_features @ <>= procedure :: set_gosam => blha_master_set_gosam <>= subroutine blha_master_set_gosam (master, i) class(blha_master_t), intent(inout) :: master integer, intent(in) :: i master%blha_mode(i) = BLHA_MODE_GOSAM end subroutine blha_master_set_gosam @ %def blha_master_set_gosam @ <>= procedure :: set_openloops => blha_master_set_openloops <>= subroutine blha_master_set_openloops (master, i) class(blha_master_t), intent(inout) :: master integer, intent(in) :: i master%blha_mode(i) = BLHA_MODE_OPENLOOPS end subroutine blha_master_set_openloops @ %def blha_master_set_openloops @ <>= procedure :: set_polarization => blha_master_set_polarization <>= subroutine blha_master_set_polarization (master, i) class(blha_master_t), intent(inout) :: master integer, intent(in) :: i master%blha_cfg(i)%polarized = .true. end subroutine blha_master_set_polarization @ %def blha_master_set_polarization @ <>= subroutine blha_init_born (blha_cfg, blha_flavor, n_in, & ap, asp, ew_scheme, basename, model, blha_mode, suffix) type(blha_configuration_t), intent(inout) :: blha_cfg type(blha_flv_state_t), intent(in), dimension(:) :: blha_flavor integer, intent(in) :: n_in integer, intent(in) :: ap, asp integer, intent(in) :: ew_scheme type(string_t), intent(in) :: basename type(model_data_t), intent(in), target :: model integer, intent(in) :: blha_mode type(string_t), intent(in) :: suffix integer, dimension(:), allocatable :: amp_type integer :: i allocate (amp_type (size (blha_flavor))) do i = 1, size (blha_flavor) amp_type(i) = BLHA_AMP_TREE end do call blha_configuration_init (blha_cfg, basename // suffix , & model, blha_mode) call blha_configuration_append_processes (blha_cfg, n_in, & blha_flavor, amp_type) call blha_configuration_set (blha_cfg, BLHA_VERSION_2, & irreg = BLHA_IRREG_CDR, alphas_power = asp, & alpha_power = ap, ew_scheme = ew_scheme, & debug = blha_mode == BLHA_MODE_GOSAM) end subroutine blha_init_born subroutine blha_init_virtual (blha_cfg, blha_flavor, n_in, & ap, asp, ew_scheme, basename, model, blha_mode, suffix) type(blha_configuration_t), intent(inout) :: blha_cfg type(blha_flv_state_t), intent(in), dimension(:) :: blha_flavor integer, intent(in) :: n_in integer, intent(in) :: ap, asp integer, intent(in) :: ew_scheme type(string_t), intent(in) :: basename type(model_data_t), intent(in), target :: model integer, intent(in) :: blha_mode type(string_t), intent(in) :: suffix integer, dimension(:), allocatable :: amp_type integer :: i allocate (amp_type (size (blha_flavor) * 2)) do i = 1, size (blha_flavor) amp_type(2 * i - 1) = BLHA_AMP_LOOP amp_type(2 * i) = BLHA_AMP_COLOR_C end do call blha_configuration_init (blha_cfg, basename // suffix , & model, blha_mode) call blha_configuration_append_processes (blha_cfg, n_in, & blha_flavor, amp_type) call blha_configuration_set (blha_cfg, BLHA_VERSION_2, & irreg = BLHA_IRREG_CDR, & alphas_power = asp, & alpha_power = ap, & ew_scheme = ew_scheme, & debug = blha_mode == BLHA_MODE_GOSAM) end subroutine blha_init_virtual subroutine blha_init_subtraction (blha_cfg, blha_flavor, n_in, & ap, asp, ew_scheme, basename, model, blha_mode, suffix) type(blha_configuration_t), intent(inout) :: blha_cfg type(blha_flv_state_t), intent(in), dimension(:) :: blha_flavor integer, intent(in) :: n_in integer, intent(in) :: ap, asp integer, intent(in) :: ew_scheme type(string_t), intent(in) :: basename type(model_data_t), intent(in), target :: model integer, intent(in) :: blha_mode type(string_t), intent(in) :: suffix integer, dimension(:), allocatable :: amp_type integer :: i allocate (amp_type (size (blha_flavor) * 3)) do i = 1, size (blha_flavor) amp_type(3 * i - 2) = BLHA_AMP_TREE amp_type(3 * i - 1) = BLHA_AMP_COLOR_C amp_type(3 * i) = BLHA_AMP_SPIN_C end do call blha_configuration_init (blha_cfg, basename // suffix , & model, blha_mode) call blha_configuration_append_processes (blha_cfg, n_in, & blha_flavor, amp_type) call blha_configuration_set (blha_cfg, BLHA_VERSION_2, & irreg = BLHA_IRREG_CDR, & alphas_power = asp, & alpha_power = ap, & ew_scheme = ew_scheme, & debug = blha_mode == BLHA_MODE_GOSAM) end subroutine blha_init_subtraction subroutine blha_init_real (blha_cfg, blha_flavor, n_in, & ap, asp, ew_scheme, basename, model, blha_mode, suffix) type(blha_configuration_t), intent(inout) :: blha_cfg type(blha_flv_state_t), intent(in), dimension(:) :: blha_flavor integer, intent(in) :: n_in integer, intent(in) :: ap, asp integer :: ap_ew, ap_qcd integer, intent(in) :: ew_scheme type(string_t), intent(in) :: basename type(model_data_t), intent(in), target :: model integer, intent(in) :: blha_mode type(string_t), intent(in) :: suffix integer, dimension(:), allocatable :: amp_type integer :: i allocate (amp_type (size (blha_flavor))) do i = 1, size (blha_flavor) amp_type(i) = BLHA_AMP_TREE end do select case (blha_cfg%correction_type) case (BLHA_CT_QCD) ap_ew = ap ap_qcd = asp + 1 case (BLHA_CT_EW) ap_ew = ap + 1 ap_qcd = asp end select call blha_configuration_init (blha_cfg, basename // suffix , & model, blha_mode) call blha_configuration_append_processes (blha_cfg, n_in, & blha_flavor, amp_type) call blha_configuration_set (blha_cfg, BLHA_VERSION_2, & irreg = BLHA_IRREG_CDR, & alphas_power = ap_qcd, & alpha_power = ap_ew, & ew_scheme = ew_scheme, & debug = blha_mode == BLHA_MODE_GOSAM) end subroutine blha_init_real @ %def blha_init_virtual blha_init_real @ %def blha_init_subtraction @ <>= public :: blha_get_additional_suffix <>= function blha_get_additional_suffix (base_suffix) result (suffix) type(string_t) :: suffix type(string_t), intent(in) :: base_suffix <> suffix = base_suffix <> end function blha_get_additional_suffix @ %def blha_master_extend_suffixes @ <>= 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 suffix = suffix // var_str ("_") // str (rank) end if @ <>= procedure :: write_olp => blha_master_write_olp <>= subroutine blha_master_write_olp (master, basename) class(blha_master_t), intent(in) :: master type(string_t), intent(in) :: basename integer :: unit type(string_t) :: filename integer :: i_file do i_file = 1, master%n_files filename = basename // master%suffix(i_file) // ".olp" unit = free_unit () open (unit, file = char (filename), status = 'replace', action = 'write') call blha_configuration_write (master%blha_cfg(i_file), unit) close (unit) end do end subroutine blha_master_write_olp @ %def blha_master_write_olp @ <>= procedure :: final => blha_master_final <>= subroutine blha_master_final (master) class(blha_master_t), intent(inout) :: master master%n_files = 0 deallocate (master%suffix) deallocate (master%blha_cfg) deallocate (master%i_file_to_nlo_index) end subroutine blha_master_final @ %def blha_master_final @ <>= public :: blha_configuration_init <>= subroutine blha_configuration_init (cfg, name, model, mode) type(blha_configuration_t), intent(inout) :: cfg type(string_t), intent(in) :: name class(model_data_t), target, intent(in) :: model integer, intent(in), optional :: mode if (.not. associated (cfg%model)) then cfg%name = name cfg%model => model end if if (present (mode)) cfg%mode = mode end subroutine blha_configuration_init @ %def blha_configuration_init @ Create an array of massive particle indices, to be used by the "MassiveParticle"-statement of the order file. <>= subroutine blha_configuration_get_massive_particles & (cfg, massive, i_massive) type(blha_configuration_t), intent(in) :: cfg logical, intent(out) :: massive integer, intent(out), dimension(:), allocatable :: i_massive integer, parameter :: max_particles = 10 integer, dimension(max_particles) :: i_massive_tmp integer, dimension(max_particles) :: checked type(blha_cfg_process_node_t), pointer :: current_process integer :: k integer :: n_massive n_massive = 0; k = 1 checked = 0 if (associated (cfg%processes)) then current_process => cfg%processes else call msg_fatal ("BLHA, massive particles: " // & "No processes allocated!") end if do call check_pdg_list (current_process%pdg_in%pdg) call check_pdg_list (current_process%pdg_out%pdg) if (k > max_particles) & call msg_fatal ("BLHA, massive particles: " // & "Max. number of particles exceeded!") if (associated (current_process%next)) then current_process => current_process%next else exit end if end do if (n_massive > 0) then allocate (i_massive (n_massive)) i_massive = i_massive_tmp (1:n_massive) massive = .true. else massive = .false. end if contains subroutine check_pdg_list (pdg_list) integer, dimension(:), intent(in) :: pdg_list integer :: i, i_pdg type(flavor_t) :: flv do i = 1, size (pdg_list) i_pdg = abs (pdg_list(i)) call flv%init (i_pdg, cfg%model) if (flv%get_mass () > 0._default) then !!! Avoid duplicates in output if (.not. any (checked == i_pdg)) then i_massive_tmp(k) = i_pdg checked(k) = i_pdg k = k + 1 n_massive = n_massive + 1 end if end if end do end subroutine check_pdg_list end subroutine blha_configuration_get_massive_particles @ %def blha_configuration_get_massive_particles @ <>= public :: blha_configuration_append_processes <>= subroutine blha_configuration_append_processes (cfg, n_in, flavor, amp_type) type(blha_configuration_t), intent(inout) :: cfg integer, intent(in) :: n_in type(blha_flv_state_t), dimension(:), intent(in) :: flavor integer, dimension(:), intent(in), optional :: amp_type integer :: n_tot type(blha_cfg_process_node_t), pointer :: current_node integer :: i_process, i_flv integer, dimension(:), allocatable :: pdg_in, pdg_out integer, dimension(:), allocatable :: flavor_state integer :: proc_offset, n_proc_tot proc_offset = 0; n_proc_tot = 0 do i_flv = 1, size (flavor) n_proc_tot = n_proc_tot + flavor(i_flv)%flv_mult end do if (.not. associated (cfg%processes)) & allocate (cfg%processes) current_node => cfg%processes do i_flv = 1, size (flavor) n_tot = size (flavor(i_flv)%flavors) allocate (pdg_in (n_in), pdg_out (n_tot - n_in)) allocate (flavor_state (n_tot)) flavor_state = flavor(i_flv)%flavors do i_process = 1, flavor(i_flv)%flv_mult pdg_in = flavor_state (1 : n_in) pdg_out = flavor_state (n_in + 1 : ) if (cfg%polarized) then select case (cfg%mode) case (BLHA_MODE_OPENLOOPS) call allocate_and_init_pdg_and_helicities (current_node, & pdg_in, pdg_out, amp_type (proc_offset + i_process)) case (BLHA_MODE_GOSAM) !!! Nothing special for GoSam yet. This exception is already caught !!! in blha_master_setup_additional_features end select else call allocate_and_init_pdg (current_node, pdg_in, pdg_out, & amp_type (proc_offset + i_process)) end if if (proc_offset + i_process /= n_proc_tot) then allocate (current_node%next) current_node => current_node%next end if if (i_process == flavor(i_flv)%flv_mult) & proc_offset = proc_offset + flavor(i_flv)%flv_mult end do deallocate (pdg_in, pdg_out) deallocate (flavor_state) end do contains subroutine allocate_and_init_pdg (node, pdg_in, pdg_out, amp_type) type(blha_cfg_process_node_t), intent(inout), pointer :: node integer, intent(in), dimension(:), allocatable :: pdg_in, pdg_out integer, intent(in) :: amp_type allocate (node%pdg_in (size (pdg_in))) allocate (node%pdg_out (size (pdg_out))) node%pdg_in%pdg = pdg_in node%pdg_out%pdg = pdg_out node%amplitude_type = amp_type end subroutine allocate_and_init_pdg subroutine allocate_and_init_pdg_and_helicities (node, pdg_in, pdg_out, amp_type) type(blha_cfg_process_node_t), intent(inout), pointer :: node integer, intent(in), dimension(:), allocatable :: pdg_in, pdg_out integer, intent(in) :: amp_type integer :: h1, h2 if (size (pdg_in) == 2) then do h1 = -1, 1, 2 do h2 = -1, 1, 2 call allocate_and_init_pdg (current_node, pdg_in, pdg_out, amp_type) current_node%pdg_in(1)%polarized = .true. current_node%pdg_in(2)%polarized = .true. current_node%pdg_in(1)%hel = h1 current_node%pdg_in(2)%hel = h2 if (h1 + h2 /= 2) then !!! not end of loop allocate (current_node%next) current_node => current_node%next end if end do end do else do h1 = -1, 1, 2 call allocate_and_init_pdg (current_node, pdg_in, pdg_out, amp_type) current_node%pdg_in(1)%polarized = .true. current_node%pdg_in(1)%hel = h1 if (h1 /= 1) then !!! not end of loop allocate (current_node%next) current_node => current_node%next end if end do end if end subroutine allocate_and_init_pdg_and_helicities end subroutine blha_configuration_append_processes @ %def blha_configuration_append_processes @ Change parameter(s). <>= public :: blha_configuration_set <>= subroutine blha_configuration_set (cfg, & version, irreg, massive_particle_scheme, & model_file, alphas_power, alpha_power, ew_scheme, width_scheme, & accuracy, debug) type(blha_configuration_t), intent(inout) :: cfg integer, optional, intent(in) :: version integer, optional, intent(in) :: irreg integer, optional, intent(in) :: massive_particle_scheme type(string_t), optional, intent(in) :: model_file integer, optional, intent(in) :: alphas_power, alpha_power integer, optional, intent(in) :: ew_scheme integer, optional, intent(in) :: width_scheme real(default), optional, intent(in) :: accuracy logical, optional, intent(in) :: debug if (present (version)) & cfg%version = version if (present (irreg)) & cfg%irreg = irreg if (present (massive_particle_scheme)) & cfg%massive_particle_scheme = massive_particle_scheme if (present (model_file)) & cfg%model_file = model_file if (present (alphas_power)) & cfg%alphas_power = alphas_power if (present (alpha_power)) & cfg%alpha_power = alpha_power if (present (ew_scheme)) & cfg%ew_scheme = ew_scheme if (present (width_scheme)) & cfg%width_scheme = width_scheme if (present (accuracy)) & cfg%accuracy_target = accuracy if (present (debug)) & cfg%debug_unstable = debug cfg%dirty = .false. end subroutine blha_configuration_set @ %def blha_configuration_set @ <>= public :: blha_configuration_get_n_proc <>= function blha_configuration_get_n_proc (cfg) result (n_proc) type(blha_configuration_t), intent(in) :: cfg integer :: n_proc n_proc = cfg%n_proc end function blha_configuration_get_n_proc @ %def blha_configuration_get_n_proc @ Write the BLHA file. Internal mode is intented for md5summing only. <>= public :: blha_configuration_write <>= subroutine blha_configuration_write (cfg, unit, internal, no_version) type(blha_configuration_t), intent(in) :: cfg integer, intent(in), optional :: unit logical, intent(in), optional :: internal, no_version integer :: u logical :: full type(string_t) :: buf type(blha_cfg_process_node_t), pointer :: node integer :: i character(3) :: pdg_char character(4) :: hel_char character(len=25), parameter :: pad = "" logical :: write_process, no_v no_v = .false. ; if (present (no_version)) no_v = no_version u = given_output_unit (unit); if (u < 0) return full = .true.; if (present (internal)) full = .not. internal if (full .and. cfg%dirty) call msg_bug ( & "BUG: attempted to write out a dirty BLHA configuration") if (full) then if (no_v) then write (u, "(A)") "# BLHA order written by WHIZARD [version]" else write (u, "(A)") "# BLHA order written by WHIZARD <>" end if write (u, "(A)") end if select case (cfg%mode) case (BLHA_MODE_GOSAM); buf = "GoSam" case (BLHA_MODE_OPENLOOPS); buf = "OpenLoops" case default; buf = "vanilla" end select write (u, "(A)") "# BLHA interface mode: " // char (buf) write (u, "(A)") "# process: " // char (cfg%name) write (u, "(A)") "# model: " // char (cfg%model%get_name ()) select case (cfg%version) case (1); buf = "BLHA1" case (2); buf = "BLHA2" end select write (u, '(A25,A)') "InterfaceVersion " // pad, char (buf) select case (cfg%correction_type) case (BLHA_CT_QCD); buf = "QCD" case (BLHA_CT_EW); buf = "EW" case (BLHA_CT_QED); buf = "QED" case default; buf = cfg%correction_type_other end select write (u,'(A25,A)') "CorrectionType" // pad, char (buf) select case (cfg%mode) case (BLHA_MODE_OPENLOOPS) buf = cfg%name // '.olc' write (u, '(A25,A)') "Extra AnswerFile" // pad, char (buf) end select select case (cfg%irreg) case (BLHA_IRREG_CDR); buf = "CDR" case (BLHA_IRREG_DRED); buf = "DRED" case (BLHA_IRREG_THV); buf = "tHV" case (BLHA_IRREG_MREG); buf = "MassReg" case default; buf = cfg%irreg_other end select write (u,'(A25,A)') "IRregularisation" // pad, char (buf) select case (cfg%massive_particle_scheme) case (BLHA_MPS_ONSHELL); buf = "OnShell" case default; buf = cfg%massive_particle_scheme_other end select if (cfg%mode == BLHA_MODE_GOSAM) & write (u,'(A25,A)') "MassiveParticleScheme" // pad, char (buf) select case (cfg%version) case (1) if (cfg%alphas_power >= 0) write (u,'(A25,A)') & "AlphasPower" // pad, int2char (cfg%alphas_power) if (cfg%alpha_power >= 0) write (u,'(A25,A)') & "AlphaPower " // pad, int2char (cfg%alpha_power) case (2) if (cfg%alphas_power >= 0) write (u,'(A25,A)') & "CouplingPower QCD " // pad, int2char (cfg%alphas_power) if (cfg%alpha_power >= 0) write (u, '(A25,A)') & "CouplingPower QED " // pad, int2char (cfg%alpha_power) end select if (cfg%ew_scheme > BLHA_EW_QED) then select case (cfg%mode) case (BLHA_MODE_GOSAM) select case (cfg%ew_scheme) case (BLHA_EW_GF); buf = "alphaGF" case (BLHA_EW_MZ); buf = "alphaMZ" case (BLHA_EW_MSBAR); buf = "alphaMSbar" case (BLHA_EW_0); buf = "alpha0" case (BLHA_EW_RUN); buf = "alphaRUN" end select write (u, '(A25, A)') "EWScheme " // pad, char (buf) case (BLHA_MODE_OPENLOOPS) select case (cfg%ew_scheme) case (BLHA_EW_0); buf = "alpha0" case (BLHA_EW_GF); buf = "Gmu" case default call msg_fatal ("OpenLoops input: Only supported EW schemes are 'alpha0' and 'Gmu'") end select write (u, '(A25, A)') "ewscheme " // pad, char (buf) end select end if select case (cfg%mode) case (BLHA_MODE_GOSAM) write (u, '(A25)', advance='no') "MassiveParticles " // pad do i = 1, size (OLP_MASSIVE_PARTICLES) if (OLP_MASSIVE_PARTICLES(i) > 0) & write (u, '(I2,1X)', advance='no') OLP_MASSIVE_PARTICLES(i) end do write (u,*) case (BLHA_MODE_OPENLOOPS) if (cfg%openloops_use_cms) then write (u, '(A25,I1)') "extra use_cms " // pad, 1 else write (u, '(A25,I1)') "extra use_cms " // pad, 0 end if write (u, '(A25,I1)') "extra me_cache " // pad, 0 if (cfg%openloops_phs_tolerance > 0) then write (u, '(A25,A4,I0)') "extra psp_tolerance " // pad, "10e-", & cfg%openloops_phs_tolerance end if call check_extra_cmd (cfg%openloops_extra_cmd) write (u, '(A)') char (cfg%openloops_extra_cmd) if (cfg%openloops_stability_log > 0) & write (u, '(A25,I1)') "extra stability_log " // pad, & cfg%openloops_stability_log if (cfg%openloops_use_collier) then write (u, '(A25,I1)') "extra preset " // pad, 2 else write (u, '(A25,I1)') "extra preset " // pad, 1 end if end select if (full) then write (u, "(A)") write (u, "(A)") "# Process definitions" write (u, "(A)") end if if (cfg%debug_unstable) & write (u, '(A25,A)') "DebugUnstable " // pad, "True" write (u, *) node => cfg%processes do while (associated (node)) write_process = .true. select case (node%amplitude_type) case (BLHA_AMP_LOOP); buf = "Loop" case (BLHA_AMP_COLOR_C); buf = "ccTree" case (BLHA_AMP_SPIN_C) if (cfg%mode == BLHA_MODE_OPENLOOPS) then buf = "sctree_polvect" else buf = "scTree" end if case (BLHA_AMP_TREE); buf = "Tree" case (BLHA_AMP_LOOPINDUCED); buf = "LoopInduced" end select if (write_process) then write (u, '(A25, A)') "AmplitudeType " // pad, char (buf) buf = "" do i = 1, size (node%pdg_in) call node%pdg_in(i)%write_pdg (pdg_char) if (node%pdg_in(i)%polarized) then call node%pdg_in(i)%write_helicity (hel_char) buf = (buf // pdg_char // hel_char) // " " else buf = (buf // pdg_char) // " " end if end do buf = buf // "-> " do i = 1, size (node%pdg_out) call node%pdg_out(i)%write_pdg (pdg_char) buf = (buf // pdg_char) // " " end do write (u, "(A)") char (trim (buf)) write (u, *) end if node => node%next end do end subroutine blha_configuration_write @ %def blha_configuration_write @ \subsection{Unit tests} Test module, followed by the corresponding implementation module. <<[[blha_ut.f90]]>>= <> module blha_ut use unit_tests use blha_uti <> <> contains <> end module blha_ut @ %def blha_ut @ <<[[blha_uti.f90]]>>= <> module blha_uti <> use format_utils, only: write_separator use variables, only: var_list_t use os_interface use models use blha_config <> <> contains <> <> end module blha_uti @ %def blha_uti @ API: driver for the unit tests below. <>= public :: blha_test <>= subroutine blha_test (u, results) integer, intent(in) :: u type(test_results_t), intent(inout) :: results call test(blha_1, "blha_1", "Test the creation of BLHA-OLP files", u, results) call test(blha_2, "blha_2", "Test the creation of BLHA-OLP files for "& &"multiple flavor structures", u, results) call test(blha_3, "blha_3", "Test helicity-information in OpenLoops OLP files", & u, results) end subroutine blha_test @ %def blha_test @ <>= subroutine setup_and_write_blha_configuration (u, single, polarized) integer, intent(in) :: u logical, intent(in), optional :: single logical, intent(in), optional :: polarized logical :: polrzd, singl type(blha_master_t) :: blha_master integer :: i integer :: n_in, n_out integer :: alpha_power, alphas_power integer, dimension(:,:), allocatable :: flv_born, flv_real type(string_t) :: proc_id, method, correction_type type(os_data_t) :: os_data type(model_list_t) :: model_list type(var_list_t) :: var_list type(model_t), pointer :: model => null () integer :: openloops_phs_tolerance polrzd = .false.; if (present (polarized)) polrzd = polarized singl = .true.; if (present (single)) singl = single if (singl) then write (u, "(A)") "* Process: e+ e- -> W+ W- b b~" n_in = 2; n_out = 4 alpha_power = 4; alphas_power = 0 allocate (flv_born (n_in + n_out, 1)) allocate (flv_real (n_in + n_out + 1, 1)) flv_born(1,1) = 11; flv_born(2,1) = -11 flv_born(3,1) = 24; flv_born(4,1) = -24 flv_born(5,1) = 5; flv_born(6,1) = -5 flv_real(1:6,1) = flv_born(:,1) flv_real(7,1) = 21 else write (u, "(A)") "* Process: e+ e- -> u:d:s U:D:S" n_in = 2; n_out = 2 alpha_power = 2; alphas_power = 0 allocate (flv_born (n_in + n_out, 3)) allocate (flv_real (n_in + n_out + 1, 3)) flv_born(1,:) = 11; flv_born(2,:) = -11 flv_born(3,1) = 1; flv_born(4,1) = -1 flv_born(3,2) = 2; flv_born(4,2) = -2 flv_born(3,3) = 3; flv_born(4,3) = -3 flv_real(1:4,:) = flv_born flv_real(5,:) = 21 end if proc_id = var_str ("BLHA_Test") call syntax_model_file_init () call os_data%init () call model_list%read_model & (var_str ("SM"), var_str ("SM.mdl"), os_data, model) write (u, "(A)") "* BLHA matrix elements assumed for all process components" write (u, "(A)") "* Mode: GoSam" method = var_str ("gosam") correction_type = var_str ("QCD") call var_list%append_string (var_str ("$born_me_method"), method) call var_list%append_string (var_str ("$real_tree_me_method"), method) call var_list%append_string (var_str ("$loop_me_method"), method) call var_list%append_string (var_str ("$correlation_me_method"), method) call blha_master%set_ew_scheme (var_str ("GF")) call blha_master%set_methods (.true., var_list) call blha_master%allocate_config_files () call blha_master%set_correction_type (correction_type) call blha_master%generate (proc_id, model, n_in, & alpha_power, alphas_power, flv_born, flv_real) call test_output (u) call blha_master%final () call var_list%final () write (u, "(A)") "* Switch to OpenLoops" openloops_phs_tolerance = 7 method = var_str ("openloops") correction_type = var_str ("QCD") call var_list%append_string (var_str ("$born_me_method"), method) call var_list%append_string (var_str ("$real_tree_me_method"), method) call var_list%append_string (var_str ("$loop_me_method"), method) call var_list%append_string (var_str ("$correlation_me_method"), method) call blha_master%set_methods (.true., var_list) call blha_master%allocate_config_files () call blha_master%set_correction_type (correction_type) call blha_master%generate (proc_id, model, n_in, & alpha_power, alphas_power, flv_born, flv_real) if (polrzd) then do i = 1, 4 call blha_master%set_polarization (i) end do end if call blha_master%setup_additional_features & (openloops_phs_tolerance, .false., 0, .false.) call test_output (u) contains subroutine test_output (u) integer, intent(in) :: u do i = 1, 4 call write_separator (u) call write_component_type (i, u) call write_separator (u) call blha_configuration_write & (blha_master%blha_cfg(i), u, no_version = .true.) end do end subroutine test_output subroutine write_component_type (i, u) integer, intent(in) :: i, u type(string_t) :: message, component_type message = var_str ("OLP-File content for ") select case (i) case (1) component_type = var_str ("loop") case (2) component_type = var_str ("subtraction") case (3) component_type = var_str ("real") case (4) component_type = var_str ("born") end select message = message // component_type // " matrix elements" write (u, "(A)") char (message) end subroutine write_component_type end subroutine setup_and_write_blha_configuration @ %def setup_and_write_blha_configuration @ <>= public :: blha_1 <>= subroutine blha_1 (u) integer, intent(in) :: u write (u, "(A)") "* Test output: blha_1" write (u, "(A)") "* Purpose: Test the creation of olp-files for single "& &"and unpolarized flavor structures" write (u, "(A)") call setup_and_write_blha_configuration (u, single = .true., polarized = .false.) end subroutine blha_1 @ %def blha_1 @ <>= public :: blha_2 <>= subroutine blha_2 (u) integer, intent(in) :: u write (u, "(A)") "* Test output: blha_2" write (u, "(A)") "* Purpose: Test the creation of olp-files for multiple "& &"and unpolarized flavor structures" write (u, "(A)") call setup_and_write_blha_configuration (u, single = .false., polarized = .false.) end subroutine blha_2 @ %def blha_2 @ <>= public :: blha_3 <>= subroutine blha_3 (u) integer, intent(in) :: u write (u, "(A)") "* Test output: blha_3" write (u, "(A)") "* Purpose: Test the creation of olp-files for single "& &"and polarized flavor structures" write (u, "(A)") call setup_and_write_blha_configuration (u, single = .true., polarized = .true.) end subroutine blha_3 @ %def blha_3 @ Index: trunk/src/gosam/gosam.nw =================================================================== --- trunk/src/gosam/gosam.nw (revision 8251) +++ trunk/src/gosam/gosam.nw (revision 8252) @@ -1,806 +1,807 @@ % -*- ess-noweb-default-code-mode: f90-mode; noweb-default-code-mode: f90-mode; -*- % WHIZARD code as NOWEB source: GoSam interface %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{GoSam Interface} \includemodulegraph{gosam} The code in this chapter makes amplitudes accessible to \whizard\ that are generated and computed by the GoSam package. These are the modules: \begin{description} \item[loop\_archive] Provide some useful extra functionality. \item[prc\_gosam] The actual interface, following the \whizard\ conventions for matrix-element generator methods. \end{description} \clearpage %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Gosam Interface} <<[[prc_gosam.f90]]>>= <> module prc_gosam use, intrinsic :: iso_c_binding !NODEP! use, intrinsic :: iso_fortran_env use kinds <> use io_units use constants use numeric_utils use system_defs, only: TAB use system_dependencies use file_utils use string_utils use physics_defs use diagnostics use os_interface use lorentz use interactions use pdg_arrays use sm_qcd use flavors use model_data use variables use process_constants use prclib_interfaces use process_libraries use prc_core_def use prc_core use blha_config use blha_olp_interfaces <> <> <> <> <> contains <> end module prc_gosam @ @ %def module prc_gosam @ <>= type, extends (prc_blha_writer_t) :: gosam_writer_t type(string_t) :: gosam_dir type(string_t) :: golem_dir type(string_t) :: samurai_dir type(string_t) :: ninja_dir type(string_t) :: form_dir type(string_t) :: qgraf_dir type(string_t) :: filter_lo, filter_nlo type(string_t) :: symmetries integer :: form_threads integer :: form_workspace type(string_t) :: fc contains <> end type gosam_writer_t @ @ %def gosam_writer_t <>= public :: gosam_def_t <>= type, extends (blha_def_t) :: gosam_def_t logical :: execute_olp = .true. contains <> end type gosam_def_t @ @ %def gosam_def_t <>= type, extends (blha_driver_t) :: gosam_driver_t type(string_t) :: gosam_dir type(string_t) :: olp_file type(string_t) :: olc_file type(string_t) :: olp_dir type(string_t) :: olp_lib contains <> end type gosam_driver_t @ @ %def gosam_driver_t <>= public :: prc_gosam_t <>= type, extends (prc_blha_t) :: prc_gosam_t logical :: initialized = .false. contains <> end type prc_gosam_t @ @ %def prc_gosam_t <>= type, extends (blha_state_t) :: gosam_state_t contains <> end type gosam_state_t @ %def gosam_state_t @ <>= procedure :: init => gosam_def_init <>= subroutine gosam_def_init (object, basename, model_name, & - prt_in, prt_out, nlo_type, var_list) + prt_in, prt_out, nlo_type, restrictions, var_list) class(gosam_def_t), intent(inout) :: object type(string_t), intent(in) :: basename type(string_t), intent(in) :: 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 (gosam_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) object%suffix = '_SUB' end select select type (writer => object%writer) type is (gosam_writer_t) - call writer%init (model_name, prt_in, prt_out) + call writer%init (model_name, prt_in, prt_out, restrictions) writer%filter_lo = var_list%get_sval (var_str ("$gosam_filter_lo")) writer%filter_nlo = var_list%get_sval (var_str ("$gosam_filter_nlo")) writer%symmetries = & var_list%get_sval (var_str ("$gosam_symmetries")) writer%form_threads = & var_list%get_ival (var_str ("form_threads")) writer%form_workspace = & var_list%get_ival (var_str ("form_workspace")) writer%fc = & var_list%get_sval (var_str ("$gosam_fc")) end select end subroutine gosam_def_init @ %def gosam_def_init @ <>= procedure :: write_config => gosam_writer_write_config <>= subroutine gosam_writer_write_config (gosam_writer) class(gosam_writer_t), intent(in) :: gosam_writer integer :: unit unit = free_unit () open (unit, file = "golem.in", status = "replace", action = "write") call gosam_writer%generate_configuration_file (unit) close(unit) end subroutine gosam_writer_write_config @ %def gosam_writer_write_config @ <>= procedure, nopass :: type_string => gosam_def_type_string <>= function gosam_def_type_string () result (string) type(string_t) :: string string = "gosam" end function gosam_def_type_string @ @ %def gosam_def_type_string <>= procedure :: write => gosam_def_write <>= subroutine gosam_def_write (object, unit) class(gosam_def_t), intent(in) :: object integer, intent(in) :: unit select type (writer => object%writer) type is (gosam_writer_t) call writer%write (unit) end select end subroutine gosam_def_write @ @ %def gosam_def_write <>= procedure :: read => gosam_def_read <>= subroutine gosam_def_read (object, unit) class(gosam_def_t), intent(out) :: object integer, intent(in) :: unit end subroutine gosam_def_read @ %def gosam_def_read @ <>= procedure :: allocate_driver => gosam_def_allocate_driver <>= subroutine gosam_def_allocate_driver (object, driver, basename) class(gosam_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 (gosam_driver_t :: driver) end subroutine gosam_def_allocate_driver @ @ %def gosam_def_allocate_driver <>= procedure, nopass :: type_name => gosam_writer_type_name <>= function gosam_writer_type_name () result (string) type(string_t) :: string string = "gosam" end function gosam_writer_type_name @ @ %def gosam_writer_type_name <>= procedure :: init => gosam_writer_init <>= pure subroutine gosam_writer_init (writer, model_name, prt_in, prt_out, restrictions) class(gosam_writer_t), intent(inout) :: writer type(string_t), intent(in) :: model_name type(string_t), dimension(:), intent(in) :: prt_in, prt_out type(string_t), intent(in), optional :: restrictions writer%gosam_dir = GOSAM_DIR writer%golem_dir = GOLEM_DIR writer%samurai_dir = SAMURAI_DIR writer%ninja_dir = NINJA_DIR writer%form_dir = FORM_DIR writer%qgraf_dir = QGRAF_DIR call writer%base_init (model_name, prt_in, prt_out) end subroutine gosam_writer_init @ %def gosam_writer_init @ <>= procedure, nopass :: type_name => gosam_driver_type_name <>= function gosam_driver_type_name () result (string) type(string_t) :: string string = "gosam" end function gosam_driver_type_name @ @ %def gosam_driver_type_name <>= procedure :: init_gosam => gosam_driver_init_gosam <>= subroutine gosam_driver_init_gosam (object, os_data, olp_file, & olc_file, olp_dir, olp_lib) class(gosam_driver_t), intent(inout) :: object type(os_data_t), intent(in) :: os_data type(string_t), intent(in) :: olp_file, olc_file, olp_dir, olp_lib object%gosam_dir = GOSAM_DIR object%olp_file = olp_file object%contract_file = olc_file object%olp_dir = olp_dir object%olp_lib = olp_lib end subroutine gosam_driver_init_gosam @ %def gosam_driver_init @ <>= procedure :: init_dlaccess_to_library => gosam_driver_init_dlaccess_to_library <>= subroutine gosam_driver_init_dlaccess_to_library & (object, os_data, dlaccess, success) class(gosam_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) :: libname, msg_buffer libname = object%olp_dir // '/.libs/libgolem_olp.' // & os_data%shrlib_ext msg_buffer = "One-Loop-Provider: Using Gosam" call msg_message (char(msg_buffer)) msg_buffer = "Loading library: " // libname call msg_message (char(msg_buffer)) call dlaccess_init (dlaccess, var_str ("."), libname, os_data) success = .not. dlaccess_has_error (dlaccess) end subroutine gosam_driver_init_dlaccess_to_library @ %def gosam_driver_init_dlaccess_to_library @ <>= procedure :: generate_configuration_file => & gosam_writer_generate_configuration_file <>= subroutine gosam_writer_generate_configuration_file & (object, unit) class(gosam_writer_t), intent(in) :: object integer, intent(in) :: unit type(string_t) :: fc_bin type(string_t) :: form_bin, qgraf_bin, haggies_bin type(string_t) :: fcflags_golem, ldflags_golem type(string_t) :: fcflags_samurai, ldflags_samurai type(string_t) :: fcflags_ninja, ldflags_ninja type(string_t) :: ldflags_avh_olo, ldflags_qcdloop fc_bin = DEFAULT_FC form_bin = object%form_dir // '/bin/tform' qgraf_bin = object%qgraf_dir // '/bin/qgraf' if (object%gosam_dir /= "") then haggies_bin = '/usr/bin/java -jar ' // object%gosam_dir // & '/share/golem/haggies/haggies.jar' else call msg_fatal ("generate_configuration_file: At least " // & "the GoSam Directory has to be specified!") end if if (object%golem_dir /= "") then fcflags_golem = "-I" // object%golem_dir // "/include/golem95" ldflags_golem = "-L" // object%golem_dir // "/lib -lgolem" end if if (object%samurai_dir /= "") then fcflags_samurai = "-I" // object%samurai_dir // "/include/samurai" ldflags_samurai = "-L" // object%samurai_dir // "/lib -lsamurai" ldflags_avh_olo = "-L" // object%samurai_dir // "/lib -lavh_olo" ldflags_qcdloop = "-L" // object%samurai_dir // "/lib -lqcdloop" end if if (object%ninja_dir /= "") then fcflags_ninja = "-I" // object%ninja_dir // "/include/ninja " & // "-I" // object%ninja_dir // "/include" ldflags_ninja = "-L" // object%ninja_dir // "/lib -lninja" end if write (unit, "(A)") "#+avh_olo.ldflags=" & // char (ldflags_avh_olo) write (unit, "(A)") "reduction_programs=golem95, samurai, ninja" write (unit, "(A)") "extensions=autotools" write (unit, "(A)") "#+qcdloop.ldflags=" & // char (ldflags_qcdloop) write (unit, "(A)") "#+zzz.extensions=qcdloop, avh_olo" write (unit, "(A)") "#fc.bin=" // char (fc_bin) write (unit, "(A)") "form.bin=" // char (form_bin) write (unit, "(A)") "qgraf.bin=" // char (qgraf_bin) write (unit, "(A)") "#golem95.fcflags=" // char (fcflags_golem) write (unit, "(A)") "#golem95.ldflags=" // char (ldflags_golem) write (unit, "(A)") "haggies.bin=" // char (haggies_bin) write (unit, "(A)") "#samurai.fcflags=" // char (fcflags_samurai) write (unit, "(A)") "#samurai.ldflags=" // char (ldflags_samurai) write (unit, "(A)") "#ninja.fcflags=" // char (fcflags_ninja) write (unit, "(A)") "#ninja.ldflags=" // char (ldflags_ninja) !!! This might collide with the mass-setup in the order-file !!! write (unit, "(A)") "zero=mU,mD,mC,mS,mB" !!! This is covered by the BLHA2 interface write (unit, "(A)") "PSP_check=False" if (char (object%filter_lo) /= "") & write (unit, "(A)") "filter.lo=" // char (object%filter_lo) if (char (object%filter_nlo) /= "") & write (unit, "(A)") "filter.nlo=" // char (object%filter_nlo) if (char (object%symmetries) /= "") & write (unit, "(A)") "symmetries=" // char(object%symmetries) write (unit, "(A,I0)") "form.threads=", object%form_threads write (unit, "(A,I0)") "form.workspace=", object%form_workspace if (char (object%fc) /= "") & write (unit, "(A)") "fc.bin=" // char(object%fc) end subroutine gosam_writer_generate_configuration_file @ %def gosam_writer_generate_configuration_file @ We have to assure that all files necessary for the configure process in the GoSam code are ready. This is done with a stamp mechanism. <>= procedure :: write_makefile => gosam_driver_write_makefile <>= subroutine gosam_driver_write_makefile (object, unit, libname) class(gosam_driver_t), intent(in) :: object integer, intent(in) :: unit type(string_t), intent(in) :: libname write (unit, "(2A)") "OLP_FILE = ", char (object%olp_file) write (unit, "(2A)") "OLP_DIR = ", char (object%olp_dir) write (unit, "(A)") write (unit, "(A)") "all: $(OLP_DIR)/config.log" write (unit, "(2A)") TAB, "make -C $(OLP_DIR) install" write (unit, "(A)") write (unit, "(3A)") "$(OLP_DIR)/config.log: " write (unit, "(4A)") TAB, char (object%gosam_dir // "/bin/gosam.py "), & "--olp $(OLP_FILE) --destination=$(OLP_DIR)", & " -f -z" write (unit, "(3A)") TAB, "cd $(OLP_DIR); ./autogen.sh --prefix=", & "$(dir $(abspath $(lastword $(MAKEFILE_LIST))))" end subroutine gosam_driver_write_makefile @ %def gosam_driver_write_makefile @ <>= procedure :: set_alpha_s => gosam_driver_set_alpha_s <>= subroutine gosam_driver_set_alpha_s (driver, alpha_s) class(gosam_driver_t), intent(in) :: driver real(default), intent(in) :: alpha_s integer :: ierr call driver%blha_olp_set_parameter & (c_char_'alphaS'//c_null_char, & dble (alpha_s), 0._double, ierr) end subroutine gosam_driver_set_alpha_s @ %def gosam_driver_set_alpha_s @ <>= procedure :: set_alpha_qed => gosam_driver_set_alpha_qed <>= subroutine gosam_driver_set_alpha_qed (driver, alpha) class(gosam_driver_t), intent(inout) :: driver real(default), intent(in) :: alpha integer :: ierr call driver%blha_olp_set_parameter & (c_char_'alpha'//c_null_char, & dble (alpha), 0._double, ierr) if (ierr == 0) call parameter_error_message (var_str ('alpha')) end subroutine gosam_driver_set_alpha_qed @ %def gosam_driver_set_alpha_qed @ <>= procedure :: set_GF => gosam_driver_set_GF <>= subroutine gosam_driver_set_GF (driver, GF) class(gosam_driver_t), intent(inout) :: driver real(default), intent(in) :: GF integer :: ierr call driver%blha_olp_set_parameter & (c_char_'GF'//c_null_char, & dble(GF), 0._double, ierr) if (ierr == 0) call parameter_error_message (var_str ('GF')) end subroutine gosam_driver_set_GF @ %def gosam_driver_set_GF @ <>= procedure :: set_weinberg_angle => gosam_driver_set_weinberg_angle <>= subroutine gosam_driver_set_weinberg_angle (driver, sw2) class(gosam_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')) end subroutine gosam_driver_set_weinberg_angle @ %def gosam_driver_set_weinberg_angle @ <>= procedure :: print_alpha_s => gosam_driver_print_alpha_s <>= subroutine gosam_driver_print_alpha_s (object) class(gosam_driver_t), intent(in) :: object call object%blha_olp_print_parameter (c_char_'alphaS'//c_null_char) end subroutine gosam_driver_print_alpha_s @ %def gosam_driver_print_alpha_s @ <>= procedure :: prepare_library => prc_gosam_prepare_library <>= subroutine prc_gosam_prepare_library (object, os_data, libname) class(prc_gosam_t), intent(inout) :: object type(os_data_t), intent(in) :: os_data type(string_t), intent(in) :: libname select type (writer => object%def%writer) type is (gosam_writer_t) call writer%write_config () end select call object%create_olp_library (libname) call object%load_driver (os_data) end subroutine prc_gosam_prepare_library @ %def prc_gosam_prepare_library @ <>= procedure :: prepare_external_code => & prc_gosam_prepare_external_code <>= subroutine prc_gosam_prepare_external_code & (core, flv_states, var_list, os_data, libname, model, i_core, is_nlo) class(prc_gosam_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 core%sqme_tree_pos = 4 call core%prepare_library (os_data, libname) call core%start () call core%read_contract_file (flv_states) call core%set_particle_properties (model) call core%set_electroweak_parameters (model) call core%print_parameter_file (i_core) end subroutine prc_gosam_prepare_external_code @ %def prc_gosam_prepare_external_code @ <>= procedure :: write_makefile => prc_gosam_write_makefile <>= subroutine prc_gosam_write_makefile (object, unit, libname) class(prc_gosam_t), intent(in) :: object integer, intent(in) :: unit type(string_t), intent(in) :: libname select type (driver => object%driver) type is (gosam_driver_t) call driver%write_makefile (unit, libname) end select end subroutine prc_gosam_write_makefile @ %def prc_gosam_write_makefile @ <>= procedure :: execute_makefile => prc_gosam_execute_makefile <>= subroutine prc_gosam_execute_makefile (object, libname) class(prc_gosam_t), intent(in) :: object type(string_t), intent(in) :: libname select type (driver => object%driver) type is (gosam_driver_t) call os_system_call ("make -f " // & libname // "_gosam.makefile") end select end subroutine prc_gosam_execute_makefile @ %def prc_gosam_execute_makefile @ <>= procedure :: create_olp_library => prc_gosam_create_olp_library <>= subroutine prc_gosam_create_olp_library (object, libname) class(prc_gosam_t), intent(inout) :: object type(string_t), intent(in) :: libname integer :: unit select type (driver => object%driver) type is (gosam_driver_t) unit = free_unit () open (unit, file = char (libname // "_gosam.makefile"), & status = "replace", action= "write") call object%write_makefile (unit, libname) close (unit) call object%execute_makefile (libname) end select end subroutine prc_gosam_create_olp_library @ %def prc_gosam_create_olp_library @ <>= procedure :: load_driver => prc_gosam_load_driver <>= subroutine prc_gosam_load_driver (object, os_data) class(prc_gosam_t), intent(inout) :: object type(os_data_t), intent(in) :: os_data logical :: dl_success select type (driver => object%driver) type is (gosam_driver_t) call driver%load (os_data, dl_success) if (.not. dl_success) & call msg_fatal ("GoSam Libraries could not be loaded") end select end subroutine prc_gosam_load_driver @ %def prc_gosam_load_driver @ <>= procedure :: start => prc_gosam_start <>= subroutine prc_gosam_start (object) class(prc_gosam_t), intent(inout) :: object integer :: ierr if (object%includes_polarization()) & call msg_fatal ('GoSam does not support polarized beams!') select type (driver => object%driver) type is (gosam_driver_t) call driver%blha_olp_start (string_f2c (driver%contract_file), ierr) end select end subroutine prc_gosam_start @ %def prc_gosam_start @ <>= procedure :: write => prc_gosam_write <>= subroutine prc_gosam_write (object, unit) class(prc_gosam_t), intent(in) :: object integer, intent(in), optional :: unit call msg_message (unit = unit, string = "GOSAM") end subroutine prc_gosam_write @ @ %def prc_gosam_write <>= procedure :: write_name => prc_gosam_write_name <>= subroutine prc_gosam_write_name (object, unit) class(prc_gosam_t), intent(in) :: object integer, intent(in), optional :: unit integer :: u u = given_output_unit (unit) write (u,"(1x,A)") "Core: GoSam" end subroutine prc_gosam_write_name @ %def prc_gosam_write_name @ <>= procedure :: init_driver => prc_gosam_init_driver <>= subroutine prc_gosam_init_driver (object, os_data) class(prc_gosam_t), intent(inout) :: object type(os_data_t), intent(in) :: os_data type(string_t) :: olp_file, olc_file, olp_dir type(string_t) :: suffix select type (def => object%def) type is (gosam_def_t) suffix = def%suffix olp_file = def%basename // suffix // '.olp' olc_file = def%basename // suffix // '.olc' olp_dir = def%basename // suffix // '_olp_modules' class default call msg_bug ("prc_gosam_init_driver: core_def should be of gosam-type") end select select type(driver => object%driver) type is (gosam_driver_t) driver%nlo_suffix = suffix call driver%init_gosam (os_data, olp_file, olc_file, olp_dir, & var_str ("libgolem_olp")) end select end subroutine prc_gosam_init_driver @ %def prc_gosam_init_driver @ <>= procedure :: set_initialized => prc_gosam_set_initialized <>= subroutine prc_gosam_set_initialized (prc_gosam) class(prc_gosam_t), intent(inout) :: prc_gosam prc_gosam%initialized = .true. end subroutine prc_gosam_set_initialized @ %def prc_gosam_set_initialized @ The BLHA-interface conventions require the quantity $S_{ij} = \langle M_{i,+}|T_iT_j|M_{i,-}\rangle$ to be produced, where $i$ is the position of the splitting gluon. However, $\tilde{M} = \langle M_{i,-}|M_{i,+}\rangle$ is needed. This can be obtained using color conservation, $\sum_{j} T_j|M\rangle = 0$, so that \begin{equation*} \sum_{j \neq i} S_{ij} = -\langle M_{i,+}|T_i^2|M_{i,-}\rangle = -C_A \langle M_{i,+}|M_{i,-}\rangle = -C_A \tilde{M}^* \end{equation*} According to BLHA conventions, the real part of $S_{ij}$ is located at positions $2i + 2nj$ in the output array, where $n$ denotes the number of external particles and the enumeration of particles starts at zero. The subsequent position, i.e. $2i + 2nj + 1$ is designated to the imaginary part of $S_{ij}$. Note that, since the first array position is 1, the implemented position association deviates from the above one in the addition of 1. <>= <>= procedure :: compute_sqme_spin_c => prc_gosam_compute_sqme_spin_c <>= subroutine prc_gosam_compute_sqme_spin_c (object, & i_flv, i_hel, em, p, ren_scale, me_sc, bad_point) class(prc_gosam_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 complex(default), intent(out) :: me_sc logical, intent(out) :: bad_point real(double), dimension(5*object%n_particles) :: mom real(double), dimension(OLP_RESULTS_LIMIT) :: r real(double) :: ren_scale_dble integer :: i, igm1, n integer :: pos_real, pos_imag real(double) :: acc_dble real(default) :: acc, alpha_s if (object%i_spin_c(i_flv, i_hel) >= 0) then me_sc = cmplx (zero ,zero, kind=default) mom = object%create_momentum_array (p) if (vanishes (ren_scale)) & call msg_fatal ("prc_gosam_compute_sqme_spin_c: ren_scale vanishes") alpha_s = object%qcd%alpha%get (ren_scale) ren_scale_dble = dble (ren_scale) select type (driver => object%driver) type is (gosam_driver_t) call driver%set_alpha_s (alpha_s) call driver%blha_olp_eval2 (object%i_spin_c(i_flv, i_hel), & mom, ren_scale_dble, r, acc_dble) end select igm1 = em - 1 n = size(p) do i = 0, n - 1 pos_real = 2 * igm1 + 2 * n * i + 1 pos_imag = pos_real + 1 me_sc = me_sc + cmplx (r(pos_real), r(pos_imag), default) end do me_sc = - conjg(me_sc) / CA acc = acc_dble if (acc > object%maximum_accuracy) bad_point = .true. else r = 0._double end if end subroutine prc_gosam_compute_sqme_spin_c @ %def prc_gosam_compute_sqme_spin_c @ <>= procedure :: allocate_workspace => prc_gosam_allocate_workspace <>= subroutine prc_gosam_allocate_workspace (object, core_state) class(prc_gosam_t), intent(in) :: object class(prc_core_state_t), intent(inout), allocatable :: core_state allocate (gosam_state_t :: core_state) end subroutine prc_gosam_allocate_workspace @ %def prc_gosam_allocate_workspace @ <>= procedure :: write => gosam_state_write <>= subroutine gosam_state_write (object, unit) class(gosam_state_t), intent(in) :: object integer, intent(in), optional :: unit call msg_warning (unit = unit, string = "gosam_state_write: What to write?") end subroutine gosam_state_write @ %def prc_gosam_state_write @ <>= procedure :: set_particle_properties => prc_gosam_set_particle_properties <>= subroutine prc_gosam_set_particle_properties (object, model) class(prc_gosam_t), intent(inout) :: object class(model_data_t), intent(in), target :: model integer :: i, i_pdg type(flavor_t) :: flv real(default) :: mass, width integer :: ierr real(default) :: top_yukawa do i = 1, OLP_N_MASSIVE_PARTICLES i_pdg = OLP_MASSIVE_PARTICLES(i) if (i_pdg < 0) cycle call flv%init (i_pdg, model) mass = flv%get_mass (); width = flv%get_width () select type (driver => object%driver) class is (blha_driver_t) if (i_pdg == 13) then call driver%set_mass_and_width (i_pdg, mass = mass) else call driver%set_mass_and_width (i_pdg, mass = mass, width = width) end if if (i_pdg == 5) call driver%blha_olp_set_parameter & ('yuk(5)'//c_null_char, dble(mass), 0._double, ierr) if (i_pdg == 6) then if (driver%external_top_yukawa > 0._default) then top_yukawa = driver%external_top_yukawa else top_yukawa = mass end if call driver%blha_olp_set_parameter & ('yuk(6)'//c_null_char, dble(top_yukawa), 0._double, ierr) end if if (driver%switch_off_muon_yukawas) then if (i_pdg == 13) call driver%blha_olp_set_parameter & ('yuk(13)' //c_null_char, 0._double, 0._double, ierr) end if end select end do end subroutine prc_gosam_set_particle_properties @ %def prc_gosam_set_particle_properties Index: trunk/src/openloops/openloops.nw =================================================================== --- trunk/src/openloops/openloops.nw (revision 8251) +++ trunk/src/openloops/openloops.nw (revision 8252) @@ -1,684 +1,685 @@ % -*- 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 @ <>= 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 () 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, var_list) + 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) + 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')) 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')) 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')) 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')) 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 <>= subroutine openloops_driver_load_sc_procedure (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 end subroutine openloops_driver_load_sc_procedure @ %def openloops_driver_load_sc_procedure @ <>= 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) 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 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%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 else sqme_spin_c = zero end if end subroutine prc_openloops_compute_sqme_spin_c @ %def prc_openloops_compute_sqme_spin_c @ Index: trunk/src/me_methods/me_methods.nw =================================================================== --- trunk/src/me_methods/me_methods.nw (revision 8251) +++ trunk/src/me_methods/me_methods.nw (revision 8252) @@ -1,6903 +1,6924 @@ % -*- 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 <> <