diff --git a/CMakeLists.txt b/CMakeLists.txt
index 4a7bde1..ae387d6 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,435 +1,443 @@
cmake_minimum_required(VERSION 3.1 FATAL_ERROR)
set(CMAKE_LEGACY_CYGWIN_WIN32 0)
set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
project("HEJ" VERSION 2.0.5 LANGUAGES C CXX)
# Set a default build type if none was specified
set(default_build_type "RelWithDebInfo")
if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES)
message(STATUS "Setting build type to '${default_build_type}' as none was specified.")
set(CMAKE_BUILD_TYPE "${default_build_type}" CACHE
STRING "Choose the type of build." FORCE)
# Set the possible values of build type for cmake-gui
set_property(CACHE CMAKE_BUILD_TYPE PROPERTY STRINGS
"Debug" "Release" "MinSizeRel" "RelWithDebInfo")
endif()
## Flags for the compiler. No warning allowed.
if (CMAKE_COMPILER_IS_GNUCC OR CMAKE_COMPILER_IS_GNUCXX)
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wextra")
elseif (MSVC)
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /W4 /WX /EHsc")
endif()
set(CMAKE_CXX_STANDARD_REQUIRED ON)
set(CMAKE_CXX_STANDARD 14)
## Create Version
set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}/cmake/Modules/")
# Get the latest abbreviated commit hash of the working branch
execute_process(
COMMAND git rev-parse HEAD
WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
OUTPUT_VARIABLE PROJECT_GIT_REVISION
OUTPUT_STRIP_TRAILING_WHITESPACE
)
# Get the current working branch
execute_process(
COMMAND git rev-parse --abbrev-ref HEAD
WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
OUTPUT_VARIABLE PROJECT_GIT_BRANCH
OUTPUT_STRIP_TRAILING_WHITESPACE
)
## target directories for install
set(INSTALL_INCLUDE_DIR_BASE include)
set(INSTALL_INCLUDE_DIR ${INSTALL_INCLUDE_DIR_BASE}/HEJ)
set(INSTALL_BIN_DIR bin)
set(INSTALL_LIB_DIR lib)
set(INSTALL_CONFIG_DIR lib/cmake/HEJ)
## Template files
configure_file( ${CMAKE_CURRENT_SOURCE_DIR}/cmake/Templates/Version.hh.in
${PROJECT_BINARY_DIR}/include/HEJ/Version.hh @ONLY )
configure_file( ${CMAKE_CURRENT_SOURCE_DIR}/cmake/Templates/HEJ-config.cc.in
${PROJECT_BINARY_DIR}/src/bin/HEJ-config.cc @ONLY )
# Generate CMake config file
include(CMakePackageConfigHelpers)
configure_package_config_file(
cmake/Templates/hej-config.cmake.in
${CMAKE_CURRENT_BINARY_DIR}/${INSTALL_CONFIG_DIR}/hej-config.cmake
INSTALL_DESTINATION ${INSTALL_CONFIG_DIR}
PATH_VARS INSTALL_INCLUDE_DIR_BASE INSTALL_LIB_DIR
)
write_basic_package_version_file(
${CMAKE_CURRENT_BINARY_DIR}/${INSTALL_CONFIG_DIR}/hej-config-version.cmake
COMPATIBILITY SameMajorVersion
)
install(
FILES ${CMAKE_CURRENT_BINARY_DIR}/${INSTALL_CONFIG_DIR}/hej-config.cmake
${CMAKE_CURRENT_BINARY_DIR}/${INSTALL_CONFIG_DIR}/hej-config-version.cmake
DESTINATION ${INSTALL_CONFIG_DIR})
## Add directories and find dependences
include_directories(${CMAKE_CURRENT_SOURCE_DIR}/include ${PROJECT_BINARY_DIR}/include)
find_package(fastjet REQUIRED)
include_directories(${fastjet_INCLUDE_DIRS})
find_package(CLHEP 2.3 REQUIRED)
include_directories(${CLHEP_INCLUDE_DIRS})
find_package(LHAPDF REQUIRED)
include_directories(${LHAPDF_INCLUDE_DIRS})
## Amend unintuitive behaviour of FindBoost.cmake
## Priority of BOOST_ROOT over BOOSTROOT matches FindBoost.cmake
## at least for cmake 3.12
if(DEFINED BOOST_ROOT)
message("BOOST_ROOT set - only looking for Boost in ${BOOST_ROOT}")
set(Boost_NO_BOOST_CMAKE ON)
elseif(DEFINED BOOSTROOT)
message("BOOSTROOT set - only looking for Boost in ${BOOSTROOT}")
set(Boost_NO_BOOST_CMAKE ON)
endif()
find_package(Boost REQUIRED COMPONENTS iostreams)
include_directories(${Boost_INCLUDE_DIRS})
find_package(yaml-cpp) # requiring yaml does not work with fedora
include_directories(${YAML_CPP_INCLUDE_DIR})
if(${EXCLUDE_HepMC})
message(STATUS "Skipping HepMC")
# avoid "unused variable" warning if EXCLUDE_rivet is set by user
set(EXCLUDE_rivet TRUE)
else()
find_package(HepMC 2 EXACT)
endif()
if(${HepMC_FOUND})
set(
CMAKE_CXX_FLAGS
"${CMAKE_CXX_FLAGS} -DHEJ_BUILD_WITH_HepMC2"
)
include_directories(${HepMC_INCLUDE_DIRS})
if(${EXCLUDE_rivet})
message(STATUS "Skipping rivet")
else()
find_package(rivet)
endif()
if(${rivet_FOUND})
include_directories(${rivet_INCLUDE_DIRS})
set(
CMAKE_CXX_FLAGS
"${CMAKE_CXX_FLAGS} -DHEJ_BUILD_WITH_RIVET"
)
endif()
endif()
if(${EXCLUDE_HepMC3})
message(STATUS "Skipping HepMC3")
else()
## version 3.1: Finding version not possible in 3.1.1, see HepMC3 git 29e7a6c3
find_package(HepMC3)
## print standard message
find_package_handle_standard_args( HepMC3
FOUND_VAR
HepMC3_FOUND
REQUIRED_VARS
HEPMC3_INCLUDE_DIR
HEPMC3_LIBRARIES
# VERSION_VAR
# HEPMC3_VERSION
)
if(${HepMC3_FOUND})
set(
CMAKE_CXX_FLAGS
"${CMAKE_CXX_FLAGS} -DHEJ_BUILD_WITH_HepMC3"
)
include_directories(${HEPMC3_INCLUDE_DIR})
endif()
endif()
if(${EXCLUDE_QCDloop})
message(STATUS "Skipping QCDloop")
else()
find_package(QCDloop 2)
endif()
if(${QCDloop_FOUND})
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DHEJ_BUILD_WITH_QCDLOOP")
include_directories(SYSTEM ${QCDloop_INCLUDE_DIRS})
endif()
if(${EXCLUDE_HighFive})
message(STATUS "Skipping HighFive")
else()
find_package(HighFive QUIET)
find_package_handle_standard_args( HighFive CONFIG_MODE )
endif()
if(${HighFive_FOUND})
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DHEJ_BUILD_WITH_HDF5")
include_directories($<TARGET_PROPERTY:HighFive,INTERFACE_INCLUDE_DIRECTORIES>)
endif()
add_subdirectory(src)
## define executable
add_executable(HEJ src/bin/HEJ.cc)
## link libraries
target_link_libraries(HEJ hejlib)
add_executable(HEJ-config src/bin/HEJ-config.cc)
file(GLOB hej_headers
${CMAKE_CURRENT_SOURCE_DIR}/include/HEJ/*.hh
${PROJECT_BINARY_DIR}/include/HEJ/*.hh
)
file(GLOB hej_detail_headers ${CMAKE_CURRENT_SOURCE_DIR}/include/HEJ/detail/*.hh)
file(GLOB lhef_headers ${CMAKE_CURRENT_SOURCE_DIR}/include/LHEF/*.h)
install(FILES ${hej_headers} DESTINATION ${INSTALL_INCLUDE_DIR})
install(FILES ${hej_detail_headers} DESTINATION ${INSTALL_INCLUDE_DIR}/detail/)
install(FILES ${lhef_headers} DESTINATION ${INSTALL_INCLUDE_DIR_BASE}/LHEF/)
install(TARGETS HEJ HEJ-config DESTINATION ${INSTALL_BIN_DIR})
## tests
enable_testing()
set(tst_dir "${CMAKE_CURRENT_SOURCE_DIR}/t")
set(tst_ME_data_dir "${tst_dir}/ME_data")
# test event classification
add_executable(test_classify ${tst_dir}/test_classify.cc)
target_link_libraries(test_classify hejlib)
add_test(
NAME t_classify
COMMAND test_classify
)
add_executable(test_classify_ref ${tst_dir}/test_classify_ref.cc)
target_link_libraries(test_classify_ref hejlib)
add_test(
NAME t_classify_ref
COMMAND test_classify_ref ${tst_dir}/classify_ref ${tst_dir}/classify.lhe.gz
)
add_test(
NAME t_classify_ref_4j
COMMAND test_classify_ref ${tst_dir}/classify_ref_4j ${tst_dir}/4j.lhe.gz
)
# test phase space point
add_executable(test_psp ${tst_dir}/test_psp.cc)
target_link_libraries(test_psp hejlib)
add_test(
NAME t_psp
COMMAND test_psp ${tst_dir}/psp_gen.lhe.gz
)
# test importing scales
add_library(scales SHARED ${tst_dir}/scales.cc)
add_executable(test_scale_import ${tst_dir}/test_scale_import)
target_link_libraries(test_scale_import hejlib)
add_test(
NAME t_scale_import
COMMAND test_scale_import ${tst_dir}/jet_config_with_import.yml
)
# test scale arithmetic (e.g. 2*H_T/4)
add_executable(test_scale_arithmetics ${tst_dir}/test_scale_arithmetics)
target_link_libraries(test_scale_arithmetics hejlib)
add_test(
NAME t_scale_arithmetics
COMMAND test_scale_arithmetics ${tst_dir}/jet_config.yml ${tst_dir}/2j.lhe.gz
)
# test "ParameterDescription"
add_executable(test_descriptions ${tst_dir}/test_descriptions)
target_link_libraries(test_descriptions hejlib)
add_test(
NAME t_descriptions
COMMAND test_descriptions
)
# test "EventParameters*Weight"
add_executable(test_parameters ${tst_dir}/test_parameters)
target_link_libraries(test_parameters hejlib)
add_test(
NAME test_parameters
COMMAND test_parameters
)
# test colour generation
add_executable(test_colours ${tst_dir}/test_colours)
target_link_libraries(test_colours hejlib)
add_test(
NAME t_colour_flow
COMMAND test_colours
)
# test matrix elements
add_executable(test_ME_generic ${tst_dir}/test_ME_generic.cc)
target_link_libraries(test_ME_generic hejlib)
add_test(
NAME t_ME_j
COMMAND test_ME_generic ${tst_ME_data_dir}/config_mtinf.yml ${tst_ME_data_dir}/ME_jets_tree.dat ${tst_ME_data_dir}/PSP_jets.lhe.gz
)
add_test(
NAME t_ME_h
COMMAND test_ME_generic ${tst_ME_data_dir}/config_mtinf.yml ${tst_ME_data_dir}/ME_h_mtinf_tree.dat ${tst_ME_data_dir}/PSP_h.lhe.gz
)
if(${QCDloop_FOUND})
add_test(
NAME t_ME_h_mt
COMMAND test_ME_generic ${tst_ME_data_dir}/config_mt.yml ${tst_ME_data_dir}/ME_h_mt_tree.dat ${tst_ME_data_dir}/PSP_h.lhe.gz
)
add_test(
NAME t_ME_h_mtmb
COMMAND test_ME_generic ${tst_ME_data_dir}/config_mtmb.yml ${tst_ME_data_dir}/ME_h_mtmb_tree.dat ${tst_ME_data_dir}/PSP_h.lhe.gz
)
endif()
add_test(
NAME t_ME_w_FKL
COMMAND test_ME_generic ${tst_ME_data_dir}/config_w_ME.yml ${tst_ME_data_dir}/ME_w_FKL_tree.dat ${tst_ME_data_dir}/PSP_w_FKL.lhe.gz
)
add_test(
NAME t_ME_w_FKL_virt
COMMAND test_ME_generic ${tst_ME_data_dir}/config_w_ME.yml ${tst_ME_data_dir}/ME_w_FKL_virt.dat ${tst_ME_data_dir}/PSP_w_FKL.lhe.gz
)
add_test(
NAME t_ME_Wp
COMMAND test_ME_generic ${tst_ME_data_dir}/config_w_ME.yml ${tst_ME_data_dir}/ME_Wp.dat ${tst_ME_data_dir}/PSP_Wp.lhe.gz
)
add_test(
NAME t_ME_Wm
COMMAND test_ME_generic ${tst_ME_data_dir}/config_w_ME.yml ${tst_ME_data_dir}/ME_Wm.dat ${tst_ME_data_dir}/PSP_Wm.lhe.gz
)
# test main executable
file(COPY "${tst_dir}/jet_config.yml" DESTINATION "${CMAKE_BINARY_DIR}")
set(test_config "${CMAKE_BINARY_DIR}/jet_config.yml")
if(${HepMC3_FOUND})
file(READ ${test_config} config)
file(WRITE ${test_config} "${config} - tst.hepmc\n")
endif()
if(${HepMC_FOUND})
file(READ ${test_config} config)
file(WRITE ${test_config} "${config} - tst.hepmc2\n")
if(${rivet_FOUND})
file(READ ${test_config} config)
file(WRITE ${test_config} "${config}\nanalysis:\n rivet: MC_XS\n output: tst\n")
endif()
endif()
set(test_cmd_main "$<TARGET_FILE:HEJ>\\\;${test_config}\\\;${tst_dir}/2j.lhe.gz")
# check that HepMC3 output is correct
if(${HepMC3_FOUND})
add_executable(check_hepmc ${tst_dir}/check_hepmc.cc)
target_link_libraries(check_hepmc hejlib)
set(test_cmd_hepmc "$<TARGET_FILE:check_hepmc>\\\;tst.hepmc")
else()
set(test_cmd_hepmc "")
endif()
# check that LHEF output is correct
add_executable(check_lhe ${tst_dir}/check_lhe.cc)
target_link_libraries(check_lhe hejlib)
set(tes_cmd_lhe "$<TARGET_FILE:check_lhe>\\\;tst.lhe")
# Run dependent tests in one command to ensure correct execution order
# Note: The commands are concatenated with "\;" to escape CMake lists.
# Thus arguments have to be escaped twice "\\\;".
# e.g. "cmd1\;cmd2\\\;arg1\\\;arg2" is executed like "cmd1 && cmd2 arg1 arg2"
add_test(
NAME t_main
COMMAND ${CMAKE_COMMAND}
-DCMDS=${test_cmd_main}\;${test_cmd_hepmc}\;${test_cmd_lhe}
-P ${CMAKE_CURRENT_SOURCE_DIR}/cmake/run_multiple_tests.cmake
)
# check HDF5 reader
if(${HighFive_FOUND})
add_executable(test_hdf5 ${tst_dir}/test_hdf5.cc)
target_link_libraries(test_hdf5 hejlib)
add_test(
NAME t_hdf5
COMMAND test_hdf5 ${tst_dir}/Wm9-g4-repack.hdf5
)
endif()
# check rivet interface
if(${RIVET_FOUND})
add_executable(check_rivet ${tst_dir}/check_rivet.cc)
target_link_libraries(check_rivet hejlib)
add_test(
NAME t_rivet
COMMAND check_rivet
)
endif()
# test boson reconstruction
add_executable(cmp_events ${tst_dir}/cmp_events.cc)
target_link_libraries(cmp_events hejlib)
add_test(
NAME t_epnu_2j_noW
COMMAND cmp_events ${tst_dir}/epnu2jLOFKL_unweight.lhe.tar.gz ${tst_dir}/epnu2jLOFKL_unweight_noW.lhe.gz
)
# test resummed result
add_executable(check_res ${tst_dir}/check_res.cc)
target_link_libraries(check_res hejlib)
if(${TEST_ALL}) # deactivate long tests by default
add_test(
NAME t_2j
COMMAND check_res ${tst_dir}/2j.lhe.gz 3.49391e+07 419684
)
add_test(
NAME t_3j
COMMAND check_res ${tst_dir}/3j.lhe.gz 2.37902e+06 25746.6
)
add_test(
+ NAME t_3j_unof
+ COMMAND check_res ${tst_dir}/3j.lhe.gz 133399 4688.83 unof
+ )
+ add_test(
+ NAME t_3j_unob
+ COMMAND check_res ${tst_dir}/3j.lhe.gz 105247 3449.45 unob
+ )
+ add_test(
NAME t_4j
COMMAND check_res ${tst_dir}/4j.lhe.gz 603713 72822.6
)
add_test(
NAME t_h_3j
COMMAND check_res ${tst_dir}/h_3j.lhe.gz 0.821622 0.0220334
)
add_test(
NAME t_h_3j_unof
COMMAND check_res ${tst_dir}/h_3j_uno.lhe.gz 0.0127362 0.000271555 unof
)
add_test(
NAME t_h_3j_unob
COMMAND check_res ${tst_dir}/h_3j_uno.lhe.gz 0.0130615 0.000224793 unob
)
add_test(
NAME t_epnu_2j
COMMAND check_res ${tst_dir}/epnu2jLOFKL_unweight.lhe.tar.gz 262.7 3
)
add_test(
NAME t_MGepnu_3j
COMMAND check_res ${tst_dir}/MGepnu3j_unweighted.tar.gz 38.9512 1
)
add_test(
NAME t_MGemnubar_3j
COMMAND check_res ${tst_dir}/MGemnubar3j_unweighted.tar.gz 24.1575 1
)
add_test(
NAME t_MGepnu_3j_unof
COMMAND check_res ${tst_dir}/MGepnu3j_unweighted.tar.gz 9.63702 0.128355 unof
)
add_test(
NAME t_MGepnu_3j_unob
COMMAND check_res ${tst_dir}/MGepnu3j_unweighted.tar.gz 9.70119 0.108436 unob
)
add_test(
NAME t_MGepnu_3j_splitf
COMMAND check_res ${tst_dir}/MGepnu3j_unweighted.tar.gz 2.91995 0.0463182 splitf
)
add_test(
NAME t_MGepnu_3j_splitb
COMMAND check_res ${tst_dir}/MGepnu3j_unweighted.tar.gz 3.40708 0.0550975 splitb
)
add_test(
NAME t_MGepnu_4j
COMMAND check_res ${tst_dir}/MGepnu4j_unweighted.tar.gz 10.2542 0.135106
)
add_test(
NAME t_MGemnubar_4j
COMMAND check_res ${tst_dir}/MGemnubar4j_unweighted.tar.gz 5.57909 0.0300496
)
add_test(
NAME t_MGepnu_4j_qqxmid
COMMAND check_res ${tst_dir}/MGepnu4j_unweighted.tar.gz 0.732084 0.005 qqxmid
)
endif()
diff --git a/Changes.md b/Changes.md
index 712d775..ed7de34 100644
--- a/Changes.md
+++ b/Changes.md
@@ -1,58 +1,60 @@
# Changelog
This is the log for changes to the HEJ program. Further changes to the HEJ API
are documented in `Changes-API.md`. If you are using HEJ as a library, please
also read the changes there.
## Version 2.X
### 2.X.0
* Resummation for W bosons with jets
- New subleading processes `extremal qqx` & `central qqx` for a quark and
anti-quark in the final state, e.g. `g g => u d_bar Wm g` (the other
subleading processes also work with W's)
- `HEJFOG` can generate mutliple jets together with a (off-shell) W bosons
decaying into lepton & neutrino
+* Resummation can now be performed on `unordered` subleading processes
+ in pure jets.
* Allow multiplication and division of multiple scale functions e.g.
`H_T/2*m_j1j2`
* Print cross sections at end of run
* Follow HepMC convention for particle Status codes: incoming = 11,
decaying = 2, outgoing = 1 (unchanged)
* Partons now have a Colour charge
- Colours are read from and written to LHE files
- For reweighted events the colours are created according to leading colour in
the FKL limit
* Allow changing the regulator lambda in input (`regulator parameter`, only for
advanced users)
* Use `git-lfs` for raw data in test (`make test` now requires `git-lfs`)
* Added support to read `hdf5` event files suggested in
[arXiv:1905.05120](https://arxiv.org/abs/1905.05120) (needs
[HighFive](https://github.com/BlueBrain/HighFive))
* Support input with avarage weight equal to the cross section (`IDWTUP=1 or 4`)
* Rename `non-HEJ` Processes to `non-resummable`
* Dropped support for HepMC 3.0.0, either HepMC version 2 or >3.1 is required
- It is now possible to write out both HepMC 2 and HepMC 3 events at the same
time
## 2.0.5
* Fixed event classification for input not ordered in rapidity
### 2.0.4
* Fixed wrong path of `HEJ_INCLUDE_DIR` in `hej-config.cmake`
### 2.0.3
* Fixed parsing of (numerical factor) * (base scale) in configuration
* Don't change scale names, but sanitise Rivet output file names instead
### 2.0.2
* Changed scale names to `"_over_"` and `"_times_"` for proper file names (was
`"/"` and `"*"` before)
### 2.0.1
* Fixed name of fixed-order generator in error message.
diff --git a/FixedOrderGen/include/Subleading.hh b/FixedOrderGen/include/Subleading.hh
index c209b2d..36f2976 100644
--- a/FixedOrderGen/include/Subleading.hh
+++ b/FixedOrderGen/include/Subleading.hh
@@ -1,20 +1,24 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#pragma once
namespace HEJFOG{
- /**
- * Bit position of different subleading channels
- * e.g. (unsigned int) 1 => only unordered
- */
- enum Subleading: unsigned {
- none = 0u,
- all = ~0u,
- uno = 1u,
- unordered = uno,
- qqx = 2u
- };
+ namespace channels{
+ /**
+ * Bit position of different subleading channels
+ * e.g. (unsigned int) 1 => only unordered
+ */
+ enum Subleading: unsigned {
+ none = 0u,
+ all = ~0u,
+ uno = 1u,
+ unordered = uno,
+ qqx = 2u
+ };
+ }
+
+ using Subleading = channels::Subleading;
}
diff --git a/FixedOrderGen/src/config.cc b/FixedOrderGen/src/config.cc
index 2a87eb6..48aeaff 100644
--- a/FixedOrderGen/src/config.cc
+++ b/FixedOrderGen/src/config.cc
@@ -1,411 +1,409 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "config.hh"
#include <cctype>
#include "Subleading.hh"
#include "HEJ/config.hh"
#include "HEJ/YAMLreader.hh"
namespace HEJFOG{
using HEJ::set_from_yaml;
using HEJ::set_from_yaml_if_defined;
using HEJ::pid::ParticleID;
namespace{
//! Get YAML tree of supported options
/**
* The configuration file is checked against this tree of options
* in assert_all_options_known.
*/
YAML::Node const & get_supported_options(){
const static YAML::Node supported = [](){
YAML::Node supported;
static const auto opts = {
"process", "events", "subleading fraction","subleading channels",
"scales", "scale factors", "max scale ratio", "pdf",
"event output", "analysis", "import scales"
};
// add subnodes to "supported" - the assigned value is irrelevant
for(auto && opt: opts) supported[opt] = "";
for(auto && jet_opt: {"min pt", "peak pt", "algorithm", "R", "max rapidity"}){
supported["jets"][jet_opt] = "";
}
for(auto && particle_type: {"Higgs", "Wp", "W+", "Wm", "W-", "Z"}){
for(auto && particle_opt: {"mass", "width"}){
supported["particle properties"][particle_type][particle_opt] = "";
}
supported["particle properties"][particle_type]["decays"]["into"] = "";
supported["particle properties"][particle_type]["decays"]["branching ratio"] = "";
}
for(auto && opt: {"mt", "use impact factors", "include bottom", "mb"}){
supported["Higgs coupling"][opt] = "";
}
for(auto && beam_opt: {"energy", "particles"}){
supported["beam"][beam_opt] = "";
}
for(auto && unweight_opt: {"sample size", "max deviation"}){
supported["unweight"][unweight_opt] = "";
}
for(auto && opt: {"name", "seed"}){
supported["random generator"][opt] = "";
}
return supported;
}();
return supported;
}
JetParameters get_jet_parameters(
YAML::Node const & node, std::string const & entry
){
const auto p = HEJ::get_jet_parameters(node, entry);
JetParameters result;
result.def = p.def;
result.min_pt = p.min_pt;
set_from_yaml(result.max_y, node, entry, "max rapidity");
set_from_yaml_if_defined(result.peak_pt, node, entry, "peak pt");
if(result.peak_pt && *result.peak_pt <= result.min_pt)
throw std::invalid_argument{
"Value of option 'peak pt' has to be larger than 'min pt'."
};
return result;
}
Beam get_Beam(
YAML::Node const & node, std::string const & entry
){
Beam beam;
std::vector<HEJ::ParticleID> particles;
set_from_yaml(beam.energy, node, entry, "energy");
set_from_yaml_if_defined(particles, node, entry, "particles");
if(! particles.empty()){
for(HEJ::ParticleID particle: particles){
if(particle != HEJ::pid::p && particle != HEJ::pid::p_bar){
throw std::invalid_argument{
"Unsupported value in option " + entry + ": particles:"
" only proton ('p') and antiproton ('p_bar') beams are supported"
};
}
}
if(particles.size() != 2){
throw std::invalid_argument{"Not exactly two beam particles"};
}
beam.particles.front() = particles.front();
beam.particles.back() = particles.back();
}
return beam;
}
std::vector<std::string> split(
std::string const & str, std::string const & delims
){
std::vector<std::string> result;
for(size_t begin, end = 0; end != str.npos;){
begin = str.find_first_not_of(delims, end);
if(begin == str.npos) break;
end = str.find_first_of(delims, begin + 1);
result.emplace_back(str.substr(begin, end - begin));
}
return result;
}
std::invalid_argument invalid_incoming(std::string const & what){
return std::invalid_argument{
"Incoming particle type " + what + " not supported,"
" incoming particles have to be 'p', 'p_bar' or partons"
};
}
std::invalid_argument invalid_outgoing(std::string const & what){
return std::invalid_argument{
"Outgoing particle type " + what + " not supported,"
" outgoing particles have to be 'j', 'photon', 'H', 'Wm', 'Wp', 'e-', 'e+', 'nu_e', 'nu_e_bar'"
};
}
HEJ::ParticleID reconstruct_boson_id(
std::vector<HEJ::ParticleID> const & ids
){
assert(ids.size()==2);
const int pidsum = ids[0] + ids[1];
if(pidsum == +1) {
assert(HEJ::is_antilepton(ids[0]));
if(HEJ::is_antineutrino(ids[0])) {
throw HEJ::not_implemented{"lepton-flavour violating final state"};
}
assert(HEJ::is_neutrino(ids[1]));
// charged antilepton + neutrino means we had a W+
return HEJ::pid::Wp;
}
if(pidsum == -1) {
assert(HEJ::is_antilepton(ids[0]));
if(HEJ::is_neutrino(ids[1])) {
throw HEJ::not_implemented{"lepton-flavour violating final state"};
}
assert(HEJ::is_antineutrino(ids[0]));
// charged lepton + antineutrino means we had a W-
return HEJ::pid::Wm;
}
throw HEJ::not_implemented{
"final state with leptons "+HEJ::name(ids[0])+" and "+HEJ::name(ids[1])
+" not supported"
};
}
Process get_process(
YAML::Node const & node, std::string const & entry
){
Process result;
std::string process_string;
set_from_yaml(process_string, node, entry);
assert(! process_string.empty());
const auto particles = split(process_string, " \n\t\v=>");
if(particles.size() < 3){
throw std::invalid_argument{
"Bad format in option process: '" + process_string
+ "', expected format is 'in1 in2 => out1 ...'"
};
}
result.incoming.front() = HEJ::to_ParticleID(particles[0]);
result.incoming.back() = HEJ::to_ParticleID(particles[1]);
for(size_t i = 0; i < result.incoming.size(); ++i){
const HEJ::ParticleID in = result.incoming[i];
if(
in != HEJ::pid::proton && in != HEJ::pid::p_bar
&& !HEJ::is_parton(in)
){
throw invalid_incoming(particles[i]);
}
}
result.njets = 0;
for(size_t i = result.incoming.size(); i < particles.size(); ++i){
assert(! particles[i].empty());
if(particles[i] == "j") ++result.njets;
else if(std::isdigit(particles[i].front())
&& particles[i].back() == 'j')
result.njets += std::stoi(particles[i]);
else{
const auto pid = HEJ::to_ParticleID(particles[i]);
if(pid==HEJ::pid::Higgs || pid==HEJ::pid::Wp || pid==HEJ::pid::Wm){
if(result.boson)
throw std::invalid_argument{
"More than one outgoing boson is not supported"
};
if(!result.boson_decay.empty())
throw std::invalid_argument{
"Production of a boson together with a lepton is not supported"
};
result.boson = pid;
} else if (HEJ::is_anylepton(pid)){
// Do not accept more leptons, if two leptons are already mentioned
if( result.boson_decay.size()>=2 )
throw std::invalid_argument{"Too many leptons provided"};
if(result.boson)
throw std::invalid_argument{
"Production of a lepton together with a boson is not supported"
};
result.boson_decay.emplace_back(pid);
} else {
throw invalid_outgoing(particles[i]);
}
}
}
if(result.njets < 2){
throw std::invalid_argument{
"Process has to include at least two jets ('j')"
};
}
if(!result.boson_decay.empty()){
std::sort(begin(result.boson_decay),end(result.boson_decay));
assert(!result.boson);
result.boson = reconstruct_boson_id(result.boson_decay);
}
return result;
}
HEJFOG::Subleading to_subleading_channel(YAML::Node const & yaml){
std::string name;
- using HEJFOG::Subleading;
+ using namespace HEJFOG::channels;
set_from_yaml(name, yaml);
if(name == "none")
return none;
if(name == "all")
return all;
if(name == "unordered" || name == "uno")
return uno;
if(name == "qqx")
return qqx;
throw HEJ::unknown_option("Unknown subleading channel '"+name+"'");
}
unsigned int get_subleading_channels(YAML::Node const & node){
using YAML::NodeType;
- using HEJFOG::Subleading;
+ using namespace HEJFOG::channels;
// all channels allowed by default
if(!node) return all;
switch(node.Type()){
case NodeType::Undefined:
return all;
case NodeType::Null:
return none;
case NodeType::Scalar:
return to_subleading_channel(node);
case NodeType::Map:
throw HEJ::invalid_type{"map is not a valid option for subleading channels"};
case NodeType::Sequence:
unsigned int channels = HEJFOG::Subleading::none;
for(auto && channel_node: node){
channels |= get_subleading_channels(channel_node);
}
return channels;
}
throw std::logic_error{"unreachable"};
}
Decay get_decay(YAML::Node const & node){
Decay decay;
set_from_yaml(decay.products, node, "into");
decay.branching_ratio=1;
set_from_yaml_if_defined(decay.branching_ratio, node, "branching ratio");
return decay;
}
std::vector<Decay> get_decays(YAML::Node const & node){
using YAML::NodeType;
if(!node) return {};
switch(node.Type()){
case NodeType::Null:
case NodeType::Undefined:
return {};
case NodeType::Scalar:
throw HEJ::invalid_type{"value is not a list of decays"};
case NodeType::Map:
return {get_decay(node)};
case NodeType::Sequence:
std::vector<Decay> result;
for(auto && decay_str: node){
result.emplace_back(get_decay(decay_str));
}
return result;
}
throw std::logic_error{"unreachable"};
}
ParticleProperties get_particle_properties(
YAML::Node const & node, std::string const & entry
){
ParticleProperties result;
set_from_yaml(result.mass, node, entry, "mass");
set_from_yaml(result.width, node, entry, "width");
try{
result.decays = get_decays(node[entry]["decays"]);
}
catch(HEJ::missing_option const & ex){
throw HEJ::missing_option{entry + ": decays: " + ex.what()};
}
catch(HEJ::invalid_type const & ex){
throw HEJ::invalid_type{entry + ": decays: " + ex.what()};
}
return result;
}
ParticlesPropMap get_all_particles_properties(YAML::Node const & node){
ParticlesPropMap result;
for(auto const & entry: node) {
const auto name = entry.first.as<std::string>();
const auto id = HEJ::to_ParticleID(name);
result.emplace(id, get_particle_properties(node,name));
}
return result;
}
UnweightSettings get_unweight(
YAML::Node const & node, std::string const & entry
){
UnweightSettings result;
set_from_yaml(result.sample_size, node, entry, "sample size");
if(result.sample_size <= 0){
throw std::invalid_argument{
"negative sample size " + std::to_string(result.sample_size)
};
}
set_from_yaml(result.max_dev, node, entry, "max deviation");
return result;
}
Config to_Config(YAML::Node const & yaml){
try{
HEJ::assert_all_options_known(yaml, get_supported_options());
}
catch(HEJ::unknown_option const & ex){
throw HEJ::unknown_option{std::string{"Unknown option '"} + ex.what() + "'"};
}
Config config;
config.process = get_process(yaml, "process");
set_from_yaml(config.events, yaml, "events");
config.jets = get_jet_parameters(yaml, "jets");
config.beam = get_Beam(yaml, "beam");
for(size_t i = 0; i < config.process.incoming.size(); ++i){
const auto & in = config.process.incoming[i];
using namespace HEJ::pid;
if( (in == p || in == p_bar) && in != config.beam.particles[i]){
throw std::invalid_argument{
"Particle type of beam " + std::to_string(i+1) + " incompatible"
+ " with type of incoming particle " + std::to_string(i+1)
};
}
}
set_from_yaml(config.pdf_id, yaml, "pdf");
set_from_yaml(config.subleading_fraction, yaml, "subleading fraction");
if(config.subleading_fraction < 0 || config.subleading_fraction > 1){
throw std::invalid_argument{
"subleading fraction has to be between 0 and 1"
};
}
if(config.subleading_fraction == 0)
config.subleading_channels = Subleading::none;
else
config.subleading_channels = get_subleading_channels(yaml["subleading channels"]);
- if(!config.process.boson && config.subleading_channels != Subleading::none)
- throw HEJ::not_implemented("Subleading processes for pure Jet production not implemented yet");
if(yaml["particle properties"]){
config.particles_properties = get_all_particles_properties(
yaml["particle properties"]);
}
if(config.process.boson
&& config.particles_properties.find(*(config.process.boson))
== config.particles_properties.end())
throw HEJ::missing_option("Process wants to generate boson "
+std::to_string(*(config.process.boson))+", but particle properties are missing");
set_from_yaml_if_defined(config.analysis_parameters, yaml, "analysis");
config.scales = HEJ::to_ScaleConfig(yaml);
set_from_yaml_if_defined(config.output, yaml, "event output");
config.rng = HEJ::to_RNGConfig(yaml, "random generator");
config.Higgs_coupling = HEJ::get_Higgs_coupling(yaml, "Higgs coupling");
if(yaml["unweight"]) config.unweight = get_unweight(yaml, "unweight");
return config;
}
} // namespace anonymous
Config load_config(std::string const & config_file){
try{
return to_Config(YAML::LoadFile(config_file));
}
catch(...){
std::cerr << "Error reading " << config_file << ":\n ";
throw;
}
}
}
diff --git a/doc/developer_manual/currents.tex b/doc/developer_manual/currents.tex
index 8892772..16e795c 100644
--- a/doc/developer_manual/currents.tex
+++ b/doc/developer_manual/currents.tex
@@ -1,592 +1,592 @@
\section{Currents}
\label{sec:currents_impl}
The following section contains a list of all the currents implemented
in \HEJ. Clean up of the code structure is ongoing. All $W$+Jet currents
are located in \texttt{src/Wjets.cc}, all Higgs+Jets currents are
defined in \texttt{src/Hjets.cc}, and pure jet currents are defined in
in \texttt{src/jets.cc}. All of these have their own separate header
files: \texttt{include/HEJ/Wjets.hh}, \texttt{include/HEJ/Hjets.hh} and
\texttt{include/HEJ/jets.hh} respectively.
The naming convention for the current contraction $\left\|S_{f_1 f_2\to f_1
f_2}\right\|^2$ is \lstinline!ME_[Boson]_[subleading-type]_[incoming]!. For
example \lstinline!ME_W_unob_qq! corresponds to the contraction $j_W^\mu
j_{\text{uno}, \mu}$ ($qQ\to \bar{q}WQg$). For bosons on the same side as the
subleading we drop the connecting underscore, e.g. \lstinline!ME_Wuno_qq!
gives $j_{W,\text{uno}}^\mu j_\mu$ ($qQ\to g\bar{q}WQ$).
\subsection{Pure Jets}
\subsubsection{Quark}
\label{sec:current_quark}
\begin{align}
j_\mu(p_i,p_j)=\bar{u}(p_i)\gamma_\mu u(p_j)
\end{align}
The exact for depends on the helicity and direction (forward/backwards) for the quarks. Currently all different contractions of incoming and outgoing states are defined in \lstinline!joi!, \lstinline!jio! and \lstinline!joo!.
\subsubsection{Gluon}
In \HEJ the currents for gluons and quarks are the same, up to a colour factor $K_g/C_F$, where
\begin{align}
K_g(p_1^-, p_a^-) = \frac{1}{2}\left(\frac{p_1^-}{p_a^-} + \frac{p_a^-}{p_1^-}\right)\left(C_A -
\frac{1}{C_A}\right)+\frac{1}{C_A}.
\end{align}
Thus we can just reuse the results from sec.~\ref{sec:current_quark}.
\subsubsection{Single unordered gluon}
Configuration $q(p_a) \to g(p_1) q(p_2) g^*(\tilde{q}_2)$~\cite{Andersen:2017kfc}
\begin{align}
\label{eq:juno}
\begin{split}
&j^{{\rm uno}\; \mu\ cd}(p_2,p_1,p_a) = i \varepsilon_{1\nu} \left( T_{2i}^{c}T_{ia}^d\
\left(U_1^{\mu\nu}-L^{\mu\nu} \right) + T_{2i}^{d}T_{ia}^c\ \left(U_2^{\mu\nu} +
L^{\mu\nu} \right) \right). \\
U_1^{\mu\nu} &= \frac1{s_{21}} \left( j_{21}^\nu j_{1a}^\mu + 2 p_2^\nu
j_{2a}^\mu \right) \qquad \qquad U_2^{\mu\nu} = \frac1{t_{a1}} \left( 2
j_{2a}^\mu p_a^\nu - j_{21}^\mu j_{1a}^\nu \right) \\
L^{\mu\nu} &= \frac1{t_{a2}} \left(-2p_1^\mu j_{2a}^\nu+2p_1.j_{2a}
g^{\mu\nu} + (\tilde{q}_1+\tilde{q}_2)^\nu j_{2a}^\mu + \frac{t_{b2}}{2} j_{2a}^\mu \left(
\frac{p_2^\nu}{p_1.p_2} + \frac{p_b^\nu}{p_1.p_b} \right) \right) ,
\end{split}
\end{align}
-$j^{{\rm uno}\; \mu}$ is currently not calculated as a separate current, but always as needed for the ME (i.e. in \lstinline!jM2unoXXX!).
+$j^{{\rm uno}\; \mu}$ is currently not calculated as a separate current, but always as needed for the ME (i.e. in \lstinline!ME_unob_XX!).
\subsection{Higgs}
Different rapidity orderings \todo{give name of functions}
\begin{enumerate}
\item $qQ\to HqQ/qHQ/qQH$ (any rapidity order, full LO ME) $\Rightarrow$ see~\ref{sec:V_H}
\item $qg\to Hqg$ (Higgs outside quark) $\Rightarrow$ see~\ref{sec:V_H}
\item $qg\to qHg$ (central Higgs) $\Rightarrow$ see~\ref{sec:V_H}
\item $qg\to qgH$ (Higgs outside gluon) $\Rightarrow$ see~\ref{sec:jH_mt}
\item $gg\to gHg$ (central Higgs) $\Rightarrow$ see~\ref{sec:V_H}
\item $gg\to ggH$ (Higgs outside gluon) $\Rightarrow$ see~\ref{sec:jH_mt}
\end{enumerate}
\subsubsection{Higgs gluon vertex}
\label{sec:V_H}
The coupling of the Higgs boson to gluons via a virtual quark loop can be written as
\begin{align}
\label{eq:VH}
V^{\mu\nu}_H(q_1, q_2) = \mathgraphics{build/figures/V_H.pdf} &=
\frac{\alpha_s m^2}{\pi v}\big[
g^{\mu\nu} T_1(q_1, q_2) - q_2^\mu q_1^\nu T_2(q_1, q_2)
\big]\, \\
&\xrightarrow{m \to \infty} \frac{\alpha_s}{3\pi
v} \left(g^{\mu\nu} q_1\cdot q_2 - q_2^\mu q_1^\nu\right).
\end{align}
The outgoing momentum of the Higgs boson is $p_H = q_1 - q_2$.
As a contraction with two currents this by implemented in \lstinline!cHdot! inside \texttt{src/Hjets.cc}.
The form factors $T_1$ and $T_2$ are then given by~\cite{DelDuca:2003ba}
\begin{align}
\label{eq:T_1}
T_1(q_1, q_2) ={}& -C_0(q_1, q_2)\*\left[2\*m^2+\frac{1}{2}\*\left(q_1^2+q_2^2-p_H^2\right)+\frac{2\*q_1^2\*q_2^2\*p_H^2}{\lambda}\right]\notag\\
&-\left[B_0(q_2)-B_0(p_H)\right]\*\frac{q_2^2}{\lambda}\*\left(q_2^2-q_1^2-p_H^2\right)\notag\\
&-\left[B_0(q_1)-B_0(p_H)\right]\*\frac{q_1^2}{\lambda}\*\left(q_1^2-q_2^2-p_H^2\right)-1\,,\displaybreak[0]\\
\label{eq:T_2}
T_2(q_1, q_2) ={}& C_0(q_1,
q_2)\*\left[\frac{4\*m^2}{\lambda}\*\left(p_H^2-q_1^2-q_2^2\right)-1-\frac{4\*q_1^2\*q_2^2}{\lambda}
-
\frac{12\*q_1^2\*q_2^2\*p_H^2}{\lambda^2}\*\left(q_1^2+q_2^2-p_H^2\right)\right]\notag\\
&-\left[B_0(q_2)-B_0(p_H)\right]\*\left[\frac{2\*q_2^2}{\lambda}+\frac{12\*q_1^2\*q_2^2}{\lambda^2}\*\left(q_2^2-q_1^2+p_H^2\right)\right]\notag\\
&-\left[B_0(q_1)-B_0(p_H)\right]\*\left[\frac{2\*q_1^2}{\lambda}+\frac{12\*q_1^2\*q_2^2}{\lambda^2}\*\left(q_1^2-q_2^2+p_H^2\right)\right]\notag\\
&-\frac{2}{\lambda}\*\left(q_1^2+q_2^2-p_H^2\right)\,,
\end{align}
where we have used the scalar bubble and triangle integrals
\begin{align}
\label{eq:B0}
B_0\left(p\right) ={}& \int \frac{d^dl}{i\pi^{\frac{d}{2}}}
\frac{1}{\left(l^2-m^2\right)\left((l+p)^2-m^2\right)}\,,\\
\label{eq:C0}
C_0\left(p,q\right) ={}& \int \frac{d^dl}{i\pi^{\frac{d}{2}}} \frac{1}{\left(l^2-m^2\right)\left((l+p)^2-m^2\right)\left((l+p-q)^2-m^2\right)}\,,
\end{align}
and the K\"{a}ll\'{e}n function
\begin{equation}
\label{eq:lambda}
\lambda = q_1^4 + q_2^4 + p_H^4 - 2\*q_1^2\*q_2^2 - 2\*q_1^2\*p_H^2- 2\*q_2^2\*p_H^2\,.
\end{equation}
The Integrals as such are provided by \QCDloop{} (see wrapper functions \lstinline!B0DD! and \lstinline!C0DD! in \texttt{src/Hjets.cc}).
In the code we are sticking to the convention of~\cite{DelDuca:2003ba}, thus instead of the $T_{1/2}$ we implement (in the functions \lstinline!A1! and \lstinline!A2!)
\begin{align}
\label{eq:A_1}
A_1(q_1, q_2) ={}& \frac{i}{16\pi^2}\*T_2(-q_1, q_2)\,,\\
\label{eq:A_2}
A_2(q_1, q_2) ={}& -\frac{i}{16\pi^2}\*T_1(-q_1, q_2)\,.
\end{align}
\subsubsection{Peripheral Higgs emission - Finite quark mass}
\label{sec:jH_mt}
We describe the emission of a peripheral Higgs boson close to a
scattering gluon with an effective current. In the following we consider
a lightcone decomposition of the gluon momenta, i.e. $p^\pm = E \pm p_z$
and $p_\perp = p_x + i p_y$. The incoming gluon momentum $p_a$ defines
the $-$ direction, so that $p_a^+ = p_{a\perp} = 0$. The outgoing
momenta are $p_1$ for the gluon and $p_H$ for the Higgs boson. We choose
the following polarisation vectors:
\begin{equation}
\label{eq:pol_vectors}
\epsilon_\mu^\pm(p_a) = \frac{j_\mu^\pm(p_1, p_a)}{\sqrt{2}
\bar{u}^\pm(p_a)u^\mp(p_1)}\,, \quad \epsilon_\mu^{\pm,*}(p_1) = -\frac{j_\mu^\pm(p_1, p_a)}{\sqrt{2}
\bar{u}^\mp(p_1)u^\pm(p_a)}\,.
\end{equation}
Following~\cite{DelDuca:2001fn}, we introduce effective polarisation
vectors to describe the contraction with the Higgs-boson production
vertex eq.~\eqref{eq:VH}:
\begin{align}
\label{eq:eps_H}
\epsilon_{H,\mu}(p_a) = \frac{T_2(p_a, p_a-p_H)}{(p_a-p_H)^2}\big[p_a\cdot
p_H\epsilon_\mu(p_a) - p_H\cdot\epsilon(p_a) p_{a,\mu}\big]\,,\\
\epsilon_{H,\mu}^*(p_1) = -\frac{T_2(p_1+p_H, p_1)}{(p_1+p_H)^2}\big[p_1\cdot
p_H\epsilon_\mu^*(p_1) - p_H\cdot\epsilon^*(p_1) p_{1,\mu}\big]\,,
\end{align}
We also employ the usual short-hand notation
\begin{equation}
\label{eq:spinor_helicity}
\spa i.j = \bar{u}^-(p_i)u^+(p_j)\,,\qquad \spb i.j =
\bar{u}^+(p_i)u^-(p_j)\,, \qquad[ i | H | j\rangle = j_\mu^+(p_i, p_j)p_H^\mu\,.
\end{equation}
Without loss of generality, we consider only the case where the incoming
gluon has positive helicity. The remaining helicity configurations can
be obtained through parity transformation.
Labelling the effective current by the helicities of the gluons we obtain
for the same-helicity case
\begin{equation}
\label{eq:jH_same_helicity}
\begin{split}
j_{H,\mu}^{++}{}&(p_1,p_a,p_H) =
\frac{m^2}{\pi v}\bigg[\\
&-\sqrt{\frac{2p_1^-}{p_a^-}}\frac{p_{1\perp}^*}{|p_{1\perp}|}\frac{t_2}{\spb a.1}\epsilon^{+,*}_{H,\mu}(p_1)
+\sqrt{\frac{2p_a^-}{p_1^-}}\frac{p_{1\perp}^*}{|p_{1\perp}|}\frac{t_2}{\spa 1.a}\epsilon^{+}_{H,\mu}(p_a)\\
&+ [1|H|a\rangle \bigg(
\frac{\sqrt{2}}{\spa 1.a}\epsilon^{+}_{H,\mu}(p_a)
+ \frac{\sqrt{2}}{\spb a.1}\epsilon^{+,*}_{H,\mu}(p_1)-\frac{\spa 1.a T_2(p_a, p_a-p_H)}{\sqrt{2}(p_a-p_H)^2}\epsilon^{+,*}_{\mu}(p_1)\\
&
\qquad
-\frac{\spb a.1 T_2(p_1+p_H,
p_1)}{\sqrt{2}(p_1+p_H)^2}\epsilon^{+}_{\mu}(p_a)-\frac{RH_4}{\sqrt{2}\spb a.1}\epsilon^{+,*}_{\mu}(p_1)+\frac{RH_5}{\sqrt{2}\spa 1.a}\epsilon^{+}_{\mu}(p_a)
\bigg)\\
&
- \frac{[1|H|a\rangle^2}{2 t_1}(p_{a,\mu} RH_{10} - p_{1,\mu} RH_{12})\bigg]
\end{split}
\end{equation}
with $t_1 = (p_a-p_1)^2$, $t_2 = (p_a-p_1-p_H)^2$ and $R = 8 \pi^2$. Eq.~\eqref{eq:jH_same_helicity}
is implemented by \lstinline!g_gH_HC! in \texttt{src/Hjets.cc}
\footnote{\lstinline!g_gH_HC! and \lstinline!g_gH_HNC! includes an additional
$1/t_2$ factor, which should be in the Matrix element instead.}.
The currents with a helicity flip is given through
\begin{equation}
\label{eq:jH_helicity_flip}
\begin{split}
j_{H,\mu}^{+-}{}&(p_1,p_a,p_H) =
\frac{m^2}{\pi v}\bigg[\\
&-\sqrt{\frac{2p_1^-}{p_a^-}}\frac{p_{1\perp}^*}{|p_{1\perp}|}\frac{t_2}{\spb
a.1}\epsilon^{-,*}_{H,\mu}(p_1)
+\sqrt{\frac{2p_a^-}{p_1^-}}\frac{p_{1\perp}}{|p_{1\perp}|}\frac{t_2}{\spb a.1}\epsilon^{+}_{H,\mu}(p_a)\\
&+ [1|H|a\rangle \left(
\frac{\sqrt{2}}{\spb a.1} \epsilon^{-,*}_{H,\mu}(p_1)
-\frac{\spa 1.a T_2(p_a, p_a-p_H)}{\sqrt{2}(p_a-p_H)^2}\epsilon^{-,*}_{\mu}(p_1)
- \frac{RH_4}{\sqrt{2}\spb a.1}\epsilon^{-,*}_{\mu}(p_1)\right)
\\
&+ [a|H|1\rangle \left(
\frac{\sqrt{2}}{\spb a.1}\epsilon^{+}_{H,\mu}(p_a)
-\frac{\spa 1.a
T_2(p_1+p_H,p_1)}{\sqrt{2}(p_1+p_H)^2}\epsilon^{+}_{\mu}(p_a)
+\frac{RH_5}{\sqrt{2}\spb a.1}\epsilon^{+}_{\mu}(p_a)
\right)\\
& - \frac{[1|H|a\rangle [a|H|1\rangle}{2 \spb a.1 ^2}(p_{a,\mu} RH_{10} - p_{1,\mu}
RH_{12})\\
&+ \frac{\spa 1.a}{\spb a.1}\bigg(RH_1p_{1,\mu}-RH_2p_{a,\mu}+2
p_1\cdot p_H \frac{T_2(p_1+p_H, p_1)}{(p_1+p_H)^2} p_{a,\mu}
\\
&
\qquad- 2p_a \cdot p_H \frac{T_2(p_a, p_a-p_H)}{(p_a-p_H)^2} p_{1,\mu}+ T_1(p_a-p_1, p_a-p_1-p_H)\frac{(p_1 + p_a)_\mu}{t_1}\\
&\qquad-\frac{(p_1+p_a)\cdot p_H}{t_1} T_2(p_a-p_1, p_a-p_1-p_H)(p_1 - p_a)_\mu
\bigg)
\bigg]\,,
\end{split}
\end{equation}
and implemented by \lstinline!g_gH_HNC! again in \texttt{src/Hjets.cc}.
If we instead choose the gluon momentum in the $+$ direction, so that
$p_a^- = p_{a\perp} = 0$, the corresponding currents are obtained by
replacing $p_1^- \to p_1^+, p_a^- \to p_a^+,
\frac{p_{1\perp}}{|p_{1\perp}|} \to -1$ in the second line of
eq.~\eqref{eq:jH_same_helicity} and eq.~\eqref{eq:jH_helicity_flip} (see variables \lstinline!ang1a! and \lstinline!sqa1! in the implementation).
The form factors $H_1,H_2,H_4,H_5, H_{10}, H_{12}$ are given in~\cite{DelDuca:2003ba}, and are implemented under their name in \texttt{src/Hjets.cc}. They reduce down to fundamental QCD integrals, which are again provided by \QCDloop.
\subsubsection{Peripheral Higgs emission - Infinite top mass}
\label{sec:jH_eff}
To get the result with infinite top mass we could either take the limit $m_t\to \infty$ in~\eqref{eq:jH_helicity_flip} and~\eqref{eq:jH_same_helicity}, or use the \textit{impact factors} as given in~\cite{DelDuca:2003ba}. Both methods are equivalent, and lead to the same result. For the first one would find
\begin{align}
\lim_{m_t\to\infty} m_t^2 H_1 &= i \frac{1}{24 \pi^2}\\
\lim_{m_t\to\infty} m_t^2 H_2 &=-i \frac{1}{24 \pi^2}\\
\lim_{m_t\to\infty} m_t^2 H_4 &= i \frac{1}{24 \pi^2}\\
\lim_{m_t\to\infty} m_t^2 H_5 &=-i \frac{1}{24 \pi^2}\\
\lim_{m_t\to\infty} m_t^2 H_{10} &= 0 \\
\lim_{m_t\to\infty} m_t^2 H_{12} &= 0.
\end{align}
\todo{double check this, see James thesis eq. 4.33}
However only the second method is implemented in the code through \lstinline!C2gHgp!
and \lstinline!C2gHgm! inside \texttt{src/Hjets.cc}, each function
calculates the square of eq. (4.23) and (4.22) from~\cite{DelDuca:2003ba} respectively.
\subsection{$W$+Jets}
\label{sec:currents_W}
\subsubsection{Quark+$W$}
\begin{figure}
\centering
\begin{minipage}[b]{0.2\textwidth}
\includegraphics[width=\textwidth]{figures/Wbits.pdf}
\end{minipage}
\begin{minipage}[b]{0.1\textwidth}
\centering{=}
\vspace{0.7cm}
\end{minipage}
\begin{minipage}[b]{0.2\textwidth}
\includegraphics[width=\textwidth]{figures/Wbits2.pdf}
\end{minipage}
\begin{minipage}[b]{0.1\textwidth}
\centering{+}
\vspace{0.7cm}
\end{minipage}
\begin{minipage}[b]{0.2\textwidth}
\includegraphics[width=\textwidth]{figures/Wbits3.pdf}
\end{minipage}
\caption{The $j_W$ current is constructed from the two diagrams which
contribute to the emission of a $W$-boson from a given quark line.}
\label{fig:jW}
\end{figure}
For a $W$ emission we require a fermion. The $j_W$ current is actually a sum of
two separate contributions, see figure~\ref{fig:jW}, one with a $W$-emission
from the initial state, and one with the $W$-emission from the final state.
Mathematically this can be seen as the following two
terms~\cite{Andersen:2012gk}\todo{cite W subleading paper}:
\begin{align}
\label{eq:Weffcur1}
j_W^\mu(p_a,p_{\ell},p_{\bar{\ell}}, p_1) =&\ \frac{g_W^2}{2}\
\frac1{p_W^2-M_W^2+i\ \Gamma_W M_W}\ \bar{u}^-(p_\ell) \gamma_\alpha
v^-(p_{\bar\ell})\nonumber \\
& \cdot \left( \frac{ \bar{u}^-(p_1) \gamma^\alpha (\slashed{p}_W +
\slashed{p}_1)\gamma^\mu u^-(p_a)}{(p_W+p_1)^2} +
\frac{ \bar{u}^-(p_1)\gamma^\mu (\slashed{p}_a + \slashed{p}_W)\gamma^\alpha u^-(p_a)}{(p_a-p_W)^2} \right).
\end{align}
There are a couple of subtleties here. There is a minus sign
distinction between the quark-anti-quark cases due to the fermion flow
of the propagator in the current. Note that the type of $W$ emission
(+ or -) will depend on the quark flavour, and that the handedness of
the quark-line is given by whether its a quark or anti-quark.
\subsubsection{$W$+uno}
\begin{figure}
\centering
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics{figures/wuno1}
\caption{}
\label{fig:U1diags}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics{figures/wuno2}
\caption{}
\label{fig:U2diags}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics{figures/wuno3}
\caption{}
\label{fig:Cdiags}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics{figures/wuno4}
\caption{}
\label{fig:Ddiags}
\end{subfigure}
\vspace{0.4cm}
\caption{Examples of each of the four categories of Feynman diagram which
contribute to at leading-order; there are twelve in total. (a) is an example where the gluon and $W$
boson are emitted from the same quark line and the gluon comes after the
$t$-channel propagator. In (b), the gluon and $W$ boson are emitted from
the same quark line and the gluon comes before the $t$-channel proagator.
In (c) the gluon is emitted from the $t$-channel gluon and in (d) the gluon
is emitted from the $b$--$3$ quark line.}
\label{fig:Wunodiags}
\end{figure}
It is necessary to include subleading processes in $W$+Jets also. All of
these currents have been built for the \lstinline!Tensor! Class detailed in
section~\ref{sec:tensor}. Similarly to the pure jet case, the uno currents are
not calculated separately, and only in the ME functions when required
in the \texttt{src/Wjets.cc} file. For unordered emissions a new
current is required, $j_{W,{\rm uno}}$, it is only non-zero for
$h_a=h_1=-$ and hence we have suppressed its helicity indices. It is
derived from the 12 leading-order Feynman diagrams in the QMRK
limit (see figure~\ref{fig:Wunodiags}). Using $T^m_{ij}$ represent fundamental
colour matrices between quark state $i$ and $j$ with adjoint index $m$ we find
\begin{align}\label{eq:wunocurrent}
\begin{split}
j^{d\,\mu}_{\rm W,uno}(p_a,p_1,p_2,p_\ell,p_{\bar{\ell}}) =& \ i \varepsilon_{\nu}(p_1)\
\bar{u}^-(p_\ell) \gamma_\rho v^-(p_{\bar \ell}) \\ & \quad \times\
\left(T^1_{2i} T^d_{ia} (\tilde U_1^{\nu\mu\rho}-\tilde L^{\nu\mu\rho}) +
T^d_{2i} T^1_{ia} (\tilde U_2^{\nu\mu\rho}+\tilde L^{\nu\mu\rho}) \right),
\end{split}
\end{align}
where expressions for $\tilde U_{1,2}^{\nu\mu\rho}$ and $\tilde L^{\nu\mu\rho}$
are given as:
\begin{align}
\label{eq:U1tensor}
\begin{split}
\tilde U_1^{\nu\mu\rho} =&\frac{\langle 2|\nu (\slashed{p}_2+ \slashed{p}_1)\mu (\slashed{p}_a - \slashed{p}_W)\rho P_L |a\rangle }{s_{12}t_{aW}} + \frac{\langle 2|\nu (\slashed{p}_2+ \slashed{p}_1)\rho P_L (\slashed{p}_2+\slashed{p}_1 + \slashed{p}_W)\mu |a\rangle }{s_{12}s_{12W}} \\
&+ \frac{\langle 2|\rho P_L (\slashed{p}_2+ \slashed{p}_W) \nu
(\slashed{p}_1 + \slashed{p}_2+\slashed{p}_W)\mu |a\rangle}{s_{2W}s_{12W}}.
\end{split}
\end{align}
\begin{align}
\label{eq:U2tensor}
\begin{split}
\tilde U_2^{\nu\mu\rho} =&\frac{\langle 2|\mu (\slashed{p}_a-\slashed{p}_W-\slashed{p}_1)\nu (\slashed{p}_a - \slashed{p}_W)\rho P_L |a\rangle }{t_{aW1}t_{aW}} + \frac{\langle 2|\mu (\slashed{p}_a-\slashed{p}_W- \slashed{p}_1)\rho P_L (\slashed{p}_a-\slashed{p}_1) \nu |a\rangle }{t_{a1W}t_{a1}} \\
&+ \frac{\langle 2|\rho P_L (\slashed{p}_2+ \slashed{p}_W) \mu
(\slashed{p}_a-\slashed{p}_1)\nu |a\rangle}{s_{2W}t_{a1}}.
\end{split}
\end{align}
\begin{align}
\label{eq:Ltensor}
\tilde L^{\nu\mu\rho} &= \frac{q_2^2}{2t_{aW2}} \left[\frac{\langle 2 |\mu (\slashed{p}_a - \slashed{p}_W) \rho P_L | a\rangle}{t_{aW}} + \frac{\langle 2 |\rho P_L (\slashed{p}_2 + \slashed{p}_W) \mu | a\rangle }{s_{2W}} \right]
\cdot \left( \frac{p_b^{\nu}}{p_b\cdot p_1} + \frac{p_3^{\nu}}{p_3\cdot p_1} \right) \nonumber \\
&\quad+\frac{1}{t_{aW2}}\left[\frac{\langle 2 |\sigma (\slashed{p}_a - \slashed{p}_W) \rho P_L | a\rangle}{t_{aW}} + \frac{\langle 2 |\rho P_L (\slashed{p}_2 + \slashed{p}_W) \sigma | a\rangle }{s_{2W}} \right] \nonumber \\
&\vphantom{+\frac{1}{t_{aW2}}}\quad\cdot \left( g^{\sigma \mu} (q_1 +q_2)^\nu + g^{\mu \nu}(-q_2 +p_1)^\sigma+ g^{\nu \sigma}(-p_1 -q_1)^\mu \right).
\end{align}
\subsubsection{$W$+Extremal $\mathbf{q\bar{q}}$}
The $W$+Jet sub-leading processes containing an extremal $q\bar{q}$ are
related by crossing symmetry to the $W$+Jet unordered processes. This
means that one can simply perform a crossing symmetry argument on
eq.~\ref{eq:wunocurrent} to arrive at the extremal $q\bar{q}$ current
required.We show the basic structure of the extremal $q\bar{q}$
current in figure~\ref{fig:qgimp}, neglecting the $W$-emission for
simplicity.
\begin{figure}
\centering
\includegraphics[width=0.3\textwidth]{{qqbarex_schem}}
\caption{Schematic structure of the $gq \to \bar{Q}Qq$ amplitude in the limit
$y_1 \sim y_2 \ll y_3$}
\label{fig:qgimp}
\end{figure}
\begin{figure}
\centering
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics{qqbarex1}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics{qqbarex2}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics{qqbarex4}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics{qqbarex5}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics{qqbarex3}
\end{subfigure}
\caption{The five tree-level graphs which contribute to the process $gq \to \bar{Q}Qq$.}
\label{fig:qg_qQQ_graphs}
\end{figure}
We can obtain the current for $g\rightarrow W q \bar{q}$ by evaluating
the current for $W$ plus unordered emissions with the normal arguments
$p_a \leftrightarrow -p_1 $ interchanged. This is a non-trivial
statement: due to the minimality of the approximations made, the
crossing symmetry normally present in the full amplitude may be
extended to the factorised current.
We must again note that swapping $p_a \leftrightarrow -p_1$ will lead
to $u$-spinors with momenta with negative energy. These are identical
to $v$-spinors with momenta with positive energy, up to an overall
phase which is common to all terms, and can therefore be neglected.
Mathematically, this is given by:
\begin{align}\label{eq:crossedJ}
j^\mu_{\rm W,g\to Q\bar{Q}}(p_a,p_1,p_2,p_\ell,p_{\bar{\ell}}) =i \varepsilon_{g\nu}
\langle \ell | \rho | \bar \ell \rangle_L
\left(T^1_{2i} T^d_{ia} (\tilde U_{1,X}^{\nu\mu\rho}-\tilde L^{\nu\mu\rho}_X) + T^d_{2i} T^1_{ia} (\tilde U_{2,X}^{\nu\mu\rho}+\tilde L_X^{\nu\mu\rho}) \right),
\end{align}
where the components are now given by
\begin{align}
\label{eq:U1tensorX}
\begin{split}
\tilde U_{1,X}^{\nu\mu\rho} =&\frac{\langle 2|\nu (\slashed{p}_a- \slashed{p}_2)\mu (\slashed{p}_1 + \slashed{p}_W)\rho P_L |1\rangle }{t_{a2}s_{1W}} + \frac{\langle 2|\nu (\slashed{p}_a- \slashed{p}_2)\rho P_L (\slashed{p}_a-\slashed{p}_2 - \slashed{p}_W)\mu |1\rangle }{t_{a2}t_{a2W}} \\
&- \frac{\langle 2|\rho P_L (\slashed{p}_2+ \slashed{p}_W) \nu
(\slashed{p}_a - \slashed{p}_2-\slashed{p}_W)\mu
|1\rangle}{s_{2W}t_{a2W}}, \\
\tilde U_{2,X}^{\nu\mu\rho} =&-\frac{\langle 2|\mu (\slashed{p}_a-\slashed{p}_W-\slashed{p}_1)\nu (\slashed{p}_1 + \slashed{p}_W)\rho P_L |1\rangle }{t_{aW1}s_{1W}} + \frac{\langle 2|\mu (\slashed{p}_a-\slashed{p}_W- \slashed{p}_1)\rho P_L (\slashed{p}_a-\slashed{p}_1) \nu |1\rangle }{t_{a1W}t_{a1}} \\
&+ \frac{\langle 2|\rho P_L (\slashed{p}_2+ \slashed{p}_W) \mu
(\slashed{p}_a-\slashed{p}_1)\nu |1\rangle}{s_{2W}t_{a1}}, \\
\tilde L^{\nu\mu\rho}_X &= \frac{q_2^2}{2s_{1W2}} \left[\frac{\langle 2 |\mu (\slashed{p}_1 + \slashed{p}_W) \rho P_L | 1\rangle}{s_{1W}} + \frac{\langle 2 |\rho P_L (\slashed{p}_2 + \slashed{p}_W) \mu | 1\rangle }{s_{2W}} \right]
\cdot \left( \frac{p_b^{\nu}}{p_a\cdot p_b} + \frac{p_3^{\nu}}{p_a\cdot p_3} \right) \\
&\quad+\frac{1}{s_{W12}}\left[-\frac{\langle 2 |\sigma (\slashed{p}_1 + \slashed{p}_W) \rho P_L | a\rangle}{t_{aW}} + \frac{\langle 2 |\rho P_L (\slashed{p}_2 + \slashed{p}_W) \sigma | 1\rangle }{s_{2W}} \right] \\
&\vphantom{+\frac{1}{t_{aW2}}}\quad\cdot \left( g^{\sigma \mu} (q_1 +q_2)^\nu + g^{\mu \nu}(-q_2 +p_1)^\sigma+ g^{\nu \sigma}(-p_1 -q_1)^\mu \right).
\end{split}
\end{align}
Notice in particular the similarity to the $W$+uno scenario (from which
this has been derived).
\subsubsection{Central $\mathbf{q\bar{q}}$ Vertex}
The final subleading process in the $W$+Jet case is the Central
$q\bar{q}$ vertex. This subleading process does not require an altered
current, but an effective vertex which is contracted with two regular
\HEJ currents. This complexity is dealt with nicely by the \lstinline!Tensor!
class, which is detailed in section~\ref{sec:tensor}.
The $W$-emission can be from the central effective vertex (scenario
dealt with by the function \texttt{jM2WqqtoqQQq()} in the file
\texttt{src/Wjets.cc}); or from either of the external quark legs
(scenario dealt with by \texttt{jM2WqqtoqQQqW()} in same file). In
the pure jets case, there are 7 separate diagrams which contribute to
this, which can be seen in figure~\ref{fig:qq_qQQq_graphs}. In the $W$+Jets
case, there are then 45 separate contributions.
\begin{figure}
\centering
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics{figures/qqbarcen1}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics{figures/qqbarcen2}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics{figures/qqbarcen3}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics{figures/qqbarcen4}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics{figures/qqbarcen5}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics{figures/qqbarcen6}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics{figures/qqbarcen7}
\end{subfigure}
\caption{All Feynman diagrams which contribute to $qq' \to qQ\bar{Q}q'$ at
leading order.}
\label{fig:qq_qQQq_graphs}
\end{figure}
The end result is of the effective vertex, after derivation, is:
\begin{align}
\label{eq:EffectiveVertexqqbar}
\begin{split}
V^{\mu\nu}_{\text{Eff}}=&
\frac{C_1}{s_{23AB}}\left(X^{\mu\nu\sigma}_{1a}\hat{t_1} + X^{\mu\nu\sigma}_{4b}\hat{t_3} +V^{\mu\nu\sigma}_{3g}\right)J_{\text{V} \sigma}(p_2,p_A,p_B,p_3)
\\
&\quad +iC_2X^{\mu\nu}_{Unc}+iC_3X^{\mu\nu}_{Cro},
\end{split}
\end{align}
where:
\begin{eqnarray}
\begin{split}
C_1=&T^e_{1q}T^g_{qa}T^e_{23}T^g_{4b} -
T^g_{1q}T^e_{qa}T^e_{23}T^g_{4b} = f^{egc}T^c_{1a}T^e_{23}T^g_{4b},
\\
C_2=&T^g_{1a}T^g_{2q}T^{g'}_{q3}T^{g'}_{4b}
\\
C_3=&T^g_{1a}T^{g'}_{2q}T^g_{q3}T^{g'}_{4b}
\end{split}
\end{eqnarray}
are the colour factors of different contributions. The following
tensor structures correspond to groupings of diagrams in
figure~\ref{fig:qq_qQQq_graphs}.
\begin{eqnarray}
\label{eq:1aFixed}
X_{1a}^{\mu\nu\sigma} &=
\frac{-g^{\mu\nu}}{s_{23AB}\hat{t_3}}\left(\frac{p^\sigma_a}{(s_{a23AB})} +
\frac{p^\sigma_1}{(s_{123AB})}\right)
\\
\label{eq:4bFixed}
X_{4b}^{\mu\nu\sigma}
&=\frac{g^{\mu\nu}}{s_{23AB}\hat{t_1}}\left(\frac{p^\sigma_b}{(s_{23bAB})}+\frac{p^\sigma_4}{(s_{234AB}}\right)
\end{eqnarray}
correspond to the first and second row of diagrams in figure~\ref{fig:qq_qQQq_graphs}.
\begin{align}
\label{eq:3GluonWEmit}
\begin{split}
X^{\mu\nu}_{3g}=\frac{1}{
\hat{t}_1s_{23AB}\,\hat{t}_3}
\bigg[&\left(q_1+p_2+p_3+p_A+p_B\right)^\nu
g^{\mu\sigma}+
\\
&\quad\left(q_3-p_2-p_3-p_A-p_B\right)^\mu g^{\sigma\nu}-
\\
& \qquad\left(q_1+q_3\right)^\sigma g^{\mu\nu}\bigg]J_{\text{V} \sigma}(p_2,p_A,p_B,p_3)
\end{split}
\end{align}
corresponds to the left diagram on the third row in
figure~\ref{fig:qq_qQQq_graphs}. One notes that all of these contributions have
the same colour factor, and as such we can group them together nicely
before summing over helicities etc. As such, the function
\texttt{MSymW()} in \texttt{src/Wjets.cc} returns a \lstinline!Tensor!
containing the information from these 5 groupings of contributions (30 diagrams
in total).
\begin{align}
\begin{split}
X^{\mu\nu}_{Unc}=\frac{\langle A|\sigma P_L|B\rangle}{\hat{t_1}\hat{t_3}} \bar{u}_2&\left[-\frac{
\gamma^\sigma P_L(\slashed{p}_2+\slashed{p}_A+\slashed{p}_B)\gamma^\mu
(\slashed{q}_3+ \slashed{p}_3)\gamma^\nu}{(s_{2AB})(t_{unc_{2}})}\right.+
\\
&\qquad\left. \frac{\gamma^\mu(\slashed{q}_1-\slashed{p}_2)\gamma^\sigma
P_L(\slashed{q}_3+\slashed{p}_3)\gamma^\nu}{(t_{unc_{1}})(t_{unc_{2}})}\right. +
\\
&\qquad\qquad\left. \frac{\gamma^\mu(\slashed{q}_1-\slashed{p}_2)\gamma^\nu(\slashed{p}_3+\slashed{p}_A+\slashed{p}_B)\gamma^\sigma P_L
}{(t_{unc_1})(s_{3AB})}\right]v_3
\end{split}
\end{align}
corresponds to the diagram on the right of row three in
figure~\ref{fig:qq_qQQq_graphs}. This contribution to the effective vertex can
be obtained in the code with the function \texttt{MUncW()} in file
\texttt{src/Wjets.cc}.
\begin{align}
\begin{split}
X^{\mu\nu}_{Cro}=\frac{\langle
A|\sigma P_L|B\rangle}{\hat{t_1}\hat{t_3}} \bar{u}_2&\left[-\frac{
\gamma^\nu(\slashed{q}_3+\slashed{p}_2)\gamma^\mu
(\slashed{p}_3+\slashed{p}_A+\slashed{p}_B)\gamma^\sigma P_L}{(t_{cro_1})(s_{3AB})}\right.+
\\
&\qquad\left. \frac{\gamma^\nu(\slashed{q}_3+\slashed{p}_2)\gamma^\sigma
P_L(\slashed{q}_1-\slashed{p}_3)\gamma^\mu}{(t_{cro_{1}})(t_{cro_{2}})}\right.+
\\ &\qquad\qquad\left
. \frac{\gamma^\sigma P_L(\slashed{p}_2+\slashed{p}_A+\slashed{p}_B)\gamma^\nu(\slashed{q}_1-\slashed{p}_3)\gamma^\mu
}{(s_{2AB})(t_{cro_2})}\right]v_3
\end{split}
\end{align}
corresponds to the last diagram in figure~\ref{fig:qq_qQQq_graphs}. This
contribution to the effective vertex can be obtained in the code with
the function \texttt{MCroW()} in file \texttt{src/Wjets.cc}.
%%% Local Variables:
%%% mode: latex
%%% TeX-master: "developer_manual"
%%% End:
diff --git a/doc/sphinx/HEJ.rst b/doc/sphinx/HEJ.rst
index cc56fe2..de02a72 100644
--- a/doc/sphinx/HEJ.rst
+++ b/doc/sphinx/HEJ.rst
@@ -1,327 +1,325 @@
.. _`Running HEJ 2`:
Running HEJ 2
=============
Quick start
-----------
In order to run HEJ 2, you need a configuration file and a file
containing fixed-order events. A sample configuration is given by the
:file:`config.yml` file distributed together with HEJ 2. Events in the
Les Houches Event File format can be generated with standard Monte Carlo
generators like `MadGraph5_aMC@NLO <https://launchpad.net/mg5amcnlo>`_
or `Sherpa <https://sherpa.hepforge.org/trac/wiki>`_. If HEJ 2 was
compiled with `HDF5 <https://www.hdfgroup.org/>`_ support, it can also
read event files in the format suggested in
`arXiv:1905.05120 <https://arxiv.org/abs/1905.05120>`_.
HEJ 2 assumes that the cross section is given by the sum of the event
weights. Depending on the fixed-order generator it may be necessary to
adjust the weights in the Les Houches Event File accordingly.
The processes supported by HEJ 2 are
- Pure multijet production
- Production of a Higgs boson with jets
- Production of a W boson with jets
..
- *TODO* Production of a Z boson or photon with jets
where at least two jets are required in each case. For the time being,
only leading-order events are supported.
After generating an event file :file:`events.lhe` adjust the parameters
under the `fixed order jets`_ setting in :file:`config.yml` to the
settings in the fixed-order generation. Resummation can then be added by
running::
HEJ config.yml events.lhe
Using the default settings, this will produce an output event file
:file:`HEJ.lhe` with events including high-energy resummation.
When using the `Docker image <https://hub.docker.com/r/hejdock/hej>`_,
HEJ can be run with
.. code-block:: bash
docker run -v $PWD:$PWD -w $PWD hejdock/hej HEJ config.yml events.lhe
.. _`HEJ 2 settings`:
Settings
--------
HEJ 2 configuration files follow the `YAML <http://yaml.org/>`_
format. The following configuration parameters are supported:
.. _`trials`:
**trials**
High-energy resummation is performed by generating a number of
resummation phase space configurations corresponding to an input
fixed-order event. This parameter specifies how many such
configurations HEJ 2 should try to generate for each input
event. Typical values vary between 10 and 100.
.. _`min extparton pt`:
**min extparton pt**
Specifies the minimum transverse momentum in GeV of the most forward
and the most backward parton. This setting is needed to regulate an
otherwise uncancelled divergence. Its value should be slightly below
the minimum transverse momentum of jets specified by `resummation
jets: min pt`_. See also the `max ext soft pt fraction`_ setting.
.. _`max ext soft pt fraction`:
**max ext soft pt fraction**
Specifies the maximum fraction that soft radiation can contribute to
the transverse momentum of each the most forward and the most backward
jet. Values between around 0.05 and 0.1 are recommended. See also the
`min extparton pt`_ setting.
.. _`fixed order jets`:
**fixed order jets**
This tag collects a number of settings specifying the jet definition
in the event input. The settings should correspond to the ones used in
the fixed-order Monte Carlo that generated the input events.
.. _`fixed order jets: min pt`:
**min pt**
Minimum transverse momentum in GeV of fixed-order jets.
.. _`fixed order jets: algorithm`:
**algorithm**
The algorithm used to define jets. Allowed settings are
:code:`kt`, :code:`cambridge`, :code:`antikt`,
:code:`cambridge for passive`. See the `FastJet
<http://fastjet.fr/>`_ documentation for a description of these
algorithms.
.. _`fixed order jets: R`:
**R**
The R parameter used in the jet algorithm, roughly corresponding
to the jet radius in the plane spanned by the rapidity and the
azimuthal angle.
.. _`resummation jets`:
**resummation jets**
This tag collects a number of settings specifying the jet definition
in the observed, i.e. resummed events. These settings are optional, by
default the same values as for the `fixed order jets`_ are assumed.
.. _`resummation jets: min pt`:
**min pt**
Minimum transverse momentum in GeV of resummation jets. This
should be between 25% and 50% larger than the minimum transverse
momentum of fixed order jets set by `fixed order jets: min pt`_.
.. _`resummation jets: algorithm`:
**algorithm**
The algorithm used to define jets. The HEJ 2 approach to
resummation relies on properties of :code:`antikt` jets, so this
value is strongly recommended. For a list of possible other
values, see the `fixed order jets: algorithm`_ setting.
.. _`resummation jets: R`:
**R**
The R parameter used in the jet algorithm.
.. _`FKL`:
**FKL**
Specifies how to treat events respecting FKL rapidity ordering. These
configurations are dominant in the high-energy limit. The possible
values are :code:`reweight` to enable resummation, :code:`keep` to
keep the events as they are up to a possible change of
renormalisation and factorisation scale, and :code:`discard` to
discard these events.
.. _`unordered`:
**unordered**
Specifies how to treat events with one emission that does not respect FKL
ordering, e.g. :code:`u d => g u d`. In the high-energy limit, such
configurations are logarithmically suppressed compared to FKL configurations.
- The possible values are the same as for the `FKL`_ setting, but
- :code:`reweight` is currently only supported for Higgs or W bosons plus jets
- production.
+ The possible values are the same as for the `FKL`_ setting.
.. _`extremal qqx`:
**extremal qqx**
Specifies how to treat events with a quark-antiquark pair as extremal partons
in rapidity, e.g. :code:`g d => u u_bar d`. In the high-energy limit, such
configurations are logarithmically suppressed compared to FKL configurations.
The possible values are the same as for the `FKL`_ setting, but
:code:`reweight` is currently only supported for W boson plus jets
production.
.. _`central qqx`:
**central qqx**
Specifies how to treat events with a quark-antiquark pair central in
rapidity, e.g. :code:`g g => g u u_bar g`. In the high-energy limit, such
configurations are logarithmically suppressed compared to FKL configurations.
The possible values are the same as for the `FKL`_ setting, but
:code:`reweight` is currently only supported for W boson plus jets
production.
.. _`non-resummable`:
**non-resummable**
Specifies how to treat events where no resummation is possible. The
allowed values are :code:`keep` to keep the events as they are up to
a possible change of renormalisation and factorisation scale and
:code:`discard` to discard these events.
.. _`scales`:
**scales**
Specifies the renormalisation and factorisation scales for the output
events. This can either be a single entry or a list :code:`[scale1,
scale2, ...]`. For the case of a list the first entry defines the
central scale. Possible values are fixed numbers to set the scale in
GeV or the following:
- :code:`H_T`: The sum of the scalar transverse momenta of all
final-state particles
- :code:`max jet pperp`: The maximum transverse momentum of all jets
- :code:`jet invariant mass`: Sum of the invariant masses of all jets
- :code:`m_j1j2`: Invariant mass between the two hardest jets.
Scales can be multiplied or divided by overall factors, e.g. :code:`H_T/2`.
It is also possible to import scales from an external library, see
:ref:`Custom scales`
.. _`scale factors`:
**scale factors**
A list of numeric factors by which each of the `scales`_ should be
multiplied. Renormalisation and factorisation scales are varied
independently. For example, a list with entries :code:`[0.5, 2]`
would give the four scale choices (0.5μ\ :sub:`r`, 0.5μ\ :sub:`f`);
(0.5μ\ :sub:`r`, 2μ\ :sub:`f`); (2μ\ :sub:`r`, 0.5μ\ :sub:`f`); (2μ\
:sub:`r`, 2μ\ :sub:`f`) in this order. The ordering corresponds to
the order of the final event weights.
.. _`max scale ratio`:
**max scale ratio**
Specifies the maximum factor by which renormalisation and
factorisation scales may difer. For a value of :code:`2` and the
example given for the `scale factors`_ the scale choices
(0.5μ\ :sub:`r`, 2μ\ :sub:`f`) and (2μ\ :sub:`r`, 0.5μ\ :sub:`f`)
will be discarded.
.. _`log correction`:
**log correction**
Whether to include corrections due to the evolution of the strong
coupling constant in the virtual corrections. Allowed values are
:code:`true` and :code:`false`.
.. _`event output`:
**event output**
Specifies the name of a single event output file or a list of such
files. The file format is either specified explicitly or derived from
the suffix. For example, :code:`events.lhe` or, equivalently
:code:`Les Houches: events.lhe` generates an output event file
:code:`events.lhe` in the Les Houches format. The supported formats
are
- :code:`file.lhe` or :code:`Les Houches: file`: The Les Houches
event file format.
- :code:`file.hepmc2` or :code:`HepMC2: file`: HepMC format version 2.
- :code:`file.hepmc3` or :code:`HepMC3: file`: HepMC format version 3.
- :code:`file.hepmc` or :code:`HepMC: file`: The latest supported
version of the HepMC format, currently version 3.
.. _`random generator`:
**random generator**
Sets parameters for random number generation.
.. _`random generator: name`:
**name**
Which random number generator to use. Currently, :code:`mixmax`
and :code:`ranlux64` are supported. Mixmax is recommended. See
the `CLHEP documentation
<http://proj-clhep.web.cern.ch/proj-clhep/index.html#docu>`_ for
details on the generators.
.. _`random generator: seed`:
**seed**
The seed for random generation. This should be a single number for
:code:`mixmax` and the name of a state file for :code:`ranlux64`.
.. _`analysis`:
**analysis**
Name and Setting for the event analyses; either a custom
analysis plugin or Rivet. For the first the :code:`plugin` sub-entry
should be set to the analysis file path. All further entries are passed on
to the analysis. To use Rivet a list of Rivet analyses have to be
given in :code:`rivet` and prefix for the yoda file has to be set
through :code:`output`. See :ref:`Writing custom analyses` for details.
.. _`Higgs coupling`:
**Higgs coupling**
This collects a number of settings concerning the effective coupling
of the Higgs boson to gluons. This is only relevant for the
production process of a Higgs boson with jets and only supported if
HEJ 2 was compiled with `QCDLoop
<https://github.com/scarrazza/qcdloop>`_ support.
.. _`Higgs coupling: use impact factors`:
**use impact factors**
Whether to use impact factors for the coupling to the most forward
and most backward partons. Impact factors imply the infinite
top-quark mass limit.
.. _`Higgs coupling: mt`:
**mt**
The value of the top-quark mass in GeV. If this is not specified,
the limit of an infinite mass is taken.
.. _`Higgs coupling: include bottom`:
**include bottom**
Whether to include the Higgs coupling to bottom quarks.
.. _`Higgs coupling: mb`:
**mb**
The value of the bottom-quark mass in GeV. Only used for the Higgs
coupling, external bottom-quarks are always assumed to be massless.
Advanced Settings
~~~~~~~~~~~~~~~~~
All of the following settings are optional. Please **do not set** any of the
following options, unless you know exactly what you are doing. The default
behaviour gives the most reliable results for a wide range of observables.
.. _`regulator parameter`:
**regulator parameter**
Slicing parameter to regularise the subtraction term, called :math:`\lambda`
in `arxiv:1706.01002 <https://arxiv.org/abs/1706.01002>`_. Default is 0.2
diff --git a/doc/sphinx/HEJFOG.rst b/doc/sphinx/HEJFOG.rst
index 70046d7..4a22e7d 100644
--- a/doc/sphinx/HEJFOG.rst
+++ b/doc/sphinx/HEJFOG.rst
@@ -1,311 +1,310 @@
The HEJ Fixed Order Generator
=============================
For high jet multiplicities event generation with standard fixed-order
generators becomes increasingly cumbersome. For example, the
leading-order production of a Higgs Boson with five or more jets is
computationally prohibitively expensive.
To this end, HEJ 2 provides the ``HEJFOG`` fixed-order generator
that allows to generate events with high jet multiplicities. To
facilitate the computation the limit of Multi-Regge Kinematics with
large invariant masses between all outgoing particles is assumed in the
matrix elements. The typical use of the ``HEJFOG`` is to supplement
low-multiplicity events from standard generators with high-multiplicity
events before using the HEJ 2 program to add high-energy
resummation.
Installation
------------
The ``HEJFOG`` comes bundled together with HEJ 2 and the
installation is very similar. After downloading HEJ 2 and
installing the prerequisites as described in :ref:`Installation` the
``HEJFOG`` can be installed with::
cmake /path/to/FixedOrderGen -DCMAKE_INSTALL_PREFIX=target/directory
make install
where :file:`/path/to/FixedOrderGen` refers to the :file:`FixedOrderGen`
subdirectory in the HEJ 2 directory. If HEJ 2 was
installed to a non-standard location, it may be necessary to specify the
directory containing :file:`HEJ-config.cmake`. If the base installation
directory is :file:`/path/to/HEJ`, :file:`HEJ-config.cmake` should be
found in :file:`/path/to/HEJ/lib/cmake/HEJ` and the commands for
installing the ``HEJFOG`` would read::
cmake /path/to/FixedOrderGen -DHEJ_DIR=/path/to/HEJ/lib/cmake/HEJ -DCMAKE_INSTALL_PREFIX=target/directory
make install
The installation can be tested with::
make test
provided that the CT10nlo PDF set is installed.
Running the fixed-order generator
---------------------------------
After installing the ``HEJFOG`` you can modify the provided
configuration file :file:`configFO.yml` and run the generator with::
HEJFOG configFO.yml
The resulting event file, by default named :file:`HEJFO.lhe`, can then be
fed into HEJ 2 like any event file generated from a standard
fixed-order generator, see :ref:`Running HEJ 2`.
Settings
--------
Similar to HEJ 2, the ``HEJFOG`` uses a `YAML
<http://yaml.org/>`_ configuration file. The settings are
.. _`process`:
**process**
The scattering process for which events are being generated. The
format is
:code:`in1 in2 => out1 out2 ...`
The incoming particles, :code:`in1`, :code:`in2` can be
- quarks: :code:`u`, :code:`d`, :code:`u_bar`, and so on
- gluons: :code:`g`
- protons :code:`p` or antiprotons :code:`p_bar`
The outgoing particles, can be
- jets: :code:`j`, multiple jets can be grouped together e.g. :code:`2j`
- bosons: Any gauge boson, they require `particle properties`_
- leptons: if they can be reconstructed to a :code:`W+` or :code:`W-`,
e.g. :code:`e+ nu_e`
At most one of the outgoing particles can be a boson, the rest has to be
partonic. At the moment only Higgs :code:`h`, :code:`W+` and :code:`W-`
bosons are supported. All other outgoing particles are jets. There have to be
at least two jets. Instead of the leptons decays of the bosons can be added
through the :ref:`particle properties<particle properties: particle:
decays>`. The latter is required to decay a Higgs boson. So :code:`p p => Wp
j j` is the same as :code:`p p => e+ nu_e 2j`, if the decay :code:`Wp => e+
nu_e` is specified in `particle properties`_.
.. _`events`:
**events**
Specifies the number of events to generate.
.. _`jets`:
**jets**
Defines the properties of the generated jets.
.. _`jets: min pt`:
**min pt**
Minimum jet transverse momentum in GeV.
.. _`jets: peak pt`:
**peak pt**
Optional setting to specify the dominant jet transverse momentum
in GeV. If the generated events are used as input for HEJ
resummation, this should be set to the minimum transverse momentum
of the resummation jets. In all cases it has to be larger than
:code:`min pt`. The effect is that only a small fraction of jets
will be generated with a transverse momentum below the value of
this setting.
.. _`jets: algorithm`:
**algorithm**
The algorithm used to define jets. Allowed settings are
:code:`kt`, :code:`cambridge`, :code:`antikt`,
:code:`cambridge for passive`. See the `FastJet
<http://fastjet.fr/>`_ documentation for a description of these
algorithms.
.. _`jets: R`:
**R**
The R parameter used in the jet algorithm.
.. _`jets: max rapidity`:
**max rapidity**
Maximum absolute value of the jet rapidity.
.. _`beam`:
**beam**
Defines various properties of the collider beam.
.. _`beam: energy`:
**energy**
The beam energy in GeV. For example, the 13
TeV LHC corresponds to a value of 6500.
.. _`beam: particles`:
**particles**
A list :code:`[p1, p2]` of two beam particles. The only supported
entries are protons :code:`p` and antiprotons :code:`p_bar`.
.. _`pdf`:
**pdf**
The `LHAPDF number <https://lhapdf.hepforge.org/pdfsets>`_ of the PDF set.
For example, 230000 corresponds to an NNPDF 2.3 NLO PDF set.
.. _`subleading fraction`:
**subleading fraction**
This setting is related to the fraction of events that are not a FKL
configuration. All possible subleading process are listed in
:ref:`subleading channels<subleading channels>`. This value must be
positive and not larger than 1. It should typically be chosen between
0.01 and 0.1. Note that while this parameter influences the chance of
generating subleading configurations, it generally does not
correspond to the actual fraction of subleading events.
.. _`subleading channels`:
**subleading channels**
Optional parameters to select a specific subleading process, multiple
channels can be selected at once. If multiple subleading configurations
are possible one will be selected at random for each event, thus each
final state will include at most one subleading process, e.g. either
:code:`unordered` or :code:`qqx`. Only has an effect if
:code:`subleading fraction` is non-zero. The following values are allowed:
- :code:`all`: All channels allowed. This is the default.
- :code:`none`: No subleading contribution, only the leading (FKL)
configurations are allowed. This is equivalent to
:code:`subleading fraction: 0`.
- :code:`unordered`: Unordered emission allowed.
Unordered emission are any rapidity ordering where exactly one gluon is
emitted outside the rapidity ordering required in FKL events. More
precisely, if at least one of the incoming particles is a quark or
antiquark and there are more than two jets in the final state,
:code:`subleading fraction` states the probability that the flavours
of the outgoing particles are assigned in such a way that an unordered
- configuration arises. Unordered emissions are currently not implemented
- for pure jet production.
+ configuration arises.
- :code:`qqx`: Production of a quark-antiquark pair. In the leading
(FKL) configurations all partons except for the most backward and
forward ones are gluons. If the :code:`qqx` channel is allowed,
:code:`subleading fraction` gives the probability of emitting a
quark-antiquark pair in place of two gluons that are adjacent in
rapidity. Quark-antiquark pairs are currently only implemented for
W boson production.
.. _`unweight`:
**unweight**
This setting defines the parameters for the partial unweighting of
events. You can disable unweighting by removing this entry from the
configuration file.
.. _`unweight: sample size`:
**sample size**
The number of weighted events used to calibrate the unweighting.
A good default is to set this to the number of target
`events`_. If the number of `events`_ is large this can
lead to significant memory consumption and a lower value should be
chosen. Contrarily, for large multiplicities the unweighting
efficiency becomes worse and the sample size should be increased.
.. _`unweight: max deviation`:
**max deviation**
Controls the range of events to which unweighting is applied. A
larger value means that a larger fraction of events are unweighted.
Typical values are between -1 and 1.
.. _`particle properties`:
**particle properties**
Specifies various properties of the different particles (Higgs, W or Z).
This is only relevant if the chosen `process`_ is the production of the
corresponding particles with jets. E.g. for the `process`_
:code:`p p => h 2j` the :code:`mass`, :code:`width` and (optionally)
:code:`decays` of the :code:`Higgs` boson are required, while all other
particle properties will be ignored.
.. _`particle properties: particle`:
**Higgs, W+, W- or Z**
The particle (Higgs, |W+|, |W-|, Z) for which the following
properties are defined.
.. |W+| replace:: W\ :sup:`+`
.. |W-| replace:: W\ :sup:`-`
.. _`particle properties: particle: mass`:
**mass**
The mass of the particle in GeV.
.. _`particle properties: particle: width`:
**width**
The total decay width of the particle in GeV.
.. _`particle properties: particle: decays`:
**decays**
Optional setting specifying the decays of the particle. Only the decay
into two particles is implemented. Each decay has the form
:code:`{into: [p1,p2], branching ratio: r}`
where :code:`p1` and :code:`p2` are the particle names of the
decay product (e.g. :code:`photon`) and :code:`r` is the branching
ratio. The branching ratio is optional (:code:`1` by default).
Decays of a Higgs boson are treated as the production and subsequent
decay of an on-shell Higgs boson, so decays into e.g. Z bosons are not
supported. In contrast W bosons are decayed off-shell, thus the
branching ratio should not be changed or set to :code:`1`.
.. _`scales`:
**scales**
Specifies the renormalisation and factorisation scales for the output
events. For details, see the corresponding entry in the HEJ 2
:ref:`HEJ 2 settings`. Note that this should usually be a
single value, as the weights resulting from additional scale choices
will simply be ignored when adding high-energy resummation with HEJ 2.
.. _`event output`:
**event output**
Specifies the name of a single event output file or a list of such
files. See the corresponding entry in the HEJ 2
:ref:`HEJ 2 settings` for details.
.. _`RanLux init`:
.. _`random generator`:
**random generator**
Sets parameters for random number generation. See the HEJ 2
:ref:`HEJ 2 settings` for details.
.. _`analysis`:
**analysis**
Specifies the name and settings for a custom analysis library. This
can be useful to specify cuts at the fixed-order level. See the
corresponding entry in the HEJ 2 :ref:`HEJ 2 settings`
for details.
.. _`Higgs coupling`:
**Higgs coupling**
This collects a number of settings concerning the effective coupling
of the Higgs boson to gluons. See the corresponding entry in the
HEJ 2 :ref:`HEJ 2 settings` for details.
diff --git a/include/HEJ/jets.hh b/include/HEJ/jets.hh
index 94ba2c5..178a631 100644
--- a/include/HEJ/jets.hh
+++ b/include/HEJ/jets.hh
@@ -1,257 +1,339 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
/** \file
* \brief Functions computing the square of current contractions in pure jets.
*
* This file contains all the necessary functions to compute the
- * current contractions for all valid pure jet HEJ processes. It will
- * also contain some pure jet ME components used in other process ME
- * calculations
+ * current contractions for all valid pure jet HEJ processes, which
+ * so far is FKL and unordered processes. It will also contain some
+ * pure jet ME components used in other process ME calculations
*
* @TODO add a namespace
*/
#pragma once
#include <complex>
#include <ostream>
#include <CLHEP/Vector/LorentzVector.h>
typedef std::complex<double> COM;
typedef COM current[4];
typedef CLHEP::HepLorentzVector HLV;
//! Square of qQ->qQ Pure Jets Scattering Current
/**
* @param p1out Momentum of final state quark
* @param p1in Momentum of initial state quark
* @param p2out Momentum of final state quark
* @param p2in Momentum of intial state quark
* @returns Square of the current contractions for qQ->qQ Scattering
*/
double ME_qQ (HLV p1out, HLV p1in, HLV p2out, HLV p2in);
//! Square of qQbar->qQbar Pure Jets Scattering Current
/**
* @param p1out Momentum of final state quark
* @param p1in Momentum of initial state quark
* @param p2out Momentum of final state anti-quark
* @param p2in Momentum of intial state anti-quark
* @returns Square of the current contractions for qQbar->qQbar Scattering
*
* @note this can be used for qbarQ->qbarQ Scattering by inputting arguments
* appropriately.
*/
double ME_qQbar (HLV p1out, HLV p1in, HLV p2out, HLV p2in);
//! Square of qbarQbar->qbarQbar Pure Jets Scattering Current
/**
* @param p1out Momentum of final state anti-quark
* @param p1in Momentum of initial state anti-quark
* @param p2out Momentum of final state anti-quark
* @param p2in Momentum of intial state anti-quark
* @returns Square of the current contractions for qbarQbar->qbarQbar Scattering
*/
double ME_qbarQbar (HLV p1out, HLV p1in, HLV p2out, HLV p2in);
//! Square of qg->qg Pure Jets Scattering Current
/**
* @param p1out Momentum of final state quark
* @param p1in Momentum of initial state quark
* @param p2out Momentum of final state gluon
* @param p2in Momentum of intial state gluon
* @returns Square of the current contractions for qg->qg Scattering
*
* @note this can be used for gq->gq Scattering by inputting arguments
* appropriately.
*/
double ME_qg (HLV p1out, HLV p1in, HLV p2out, HLV p2in);
//! Square of qbarg->qbarg Pure Jets Scattering Current
/**
* @param p1out Momentum of final state anti-quark
* @param p1in Momentum of initial state anti-quark
* @param p2out Momentum of final state gluon
* @param p2in Momentum of intial state gluon
* @returns Square of the current contractions for qbarg->qbarg Scattering
*
* @note this can be used for gqbar->gqbar Scattering by inputting arguments
* appropriately.
*/
double ME_qbarg (HLV p1out, HLV p1in, HLV p2out, HLV p2in);
//! Square of gg->gg Pure Jets Scattering Current
/**
* @param p1out Momentum of final state gluon
* @param p1in Momentum of initial state gluon
* @param p2out Momentum of final state gluon
* @param p2in Momentum of intial state gluon
* @returns Square of the current contractions for gg->gg Scattering
*/
double ME_gg (HLV p1out, HLV p1in, HLV p2out, HLV p2in);
+//Unordered Backwards contributions:
+//! Square of qQ->qQ Pure Jets Scattering Current
+/**
+ * @param p1out Momentum of final state quark
+ * @param p1in Momentum of initial state quark
+ * @param pg Momentum of unordered gluon
+ * @param p2out Momentum of final state quark
+ * @param p2in Momentum of intial state quark
+ * @returns Square of the current contractions for qQ->qQ Scattering
+ */
+double ME_unob_qQ (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in);
+
+//! Square of qbarQ->qbarQ Pure Jets Unordered backwards Scattering Current
+/**
+ * @param p1out Momentum of final state anti-quark
+ * @param p1in Momentum of initial state anti-quark
+ * @param pg Momentum of unordered gluon
+ * @param p2out Momentum of final state quark
+ * @param p2in Momentum of intial state quark
+ * @returns Square of the current contractions for qbarQ->qbarQ Scattering
+ *
+ * @note this can be used for unof contributions by inputting
+ * arguments appropriately.
+ */
+double ME_unob_qbarQ (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in);
+
+//! Square of qQbar->qQbar Pure Jets Unordered backwards Scattering Current
+/**
+ * @param p1out Momentum of final state quark
+ * @param p1in Momentum of initial state quark
+ * @param pg Momentum of unordered gluon
+ * @param p2out Momentum of final state anti-quark
+ * @param p2in Momentum of intial state anti-quark
+ * @returns Square of the current contractions for qQbar->qQbar Scattering
+ *
+ * @note this can be used for unof contributions by inputting
+ * arguments appropriately.
+ */
+double ME_unob_qQbar (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in);
+
+//! Square of qbarQbar->qbarQbar Pure Jets Unordered backwards Scattering Current
+/**
+ * @param p1out Momentum of final state anti-quark
+ * @param p1in Momentum of initial state anti-quark
+ * @param pg Momentum of unordered gluon
+ * @param p2out Momentum of final state anti-quark
+ * @param p2in Momentum of intial state anti-quark
+ * @returns Square of the current contractions for qbarQbar->qbarQbar Scattering
+ *
+ * @note this can be used for unof contributions by inputting
+ * arguments appropriately.
+ */
+double ME_unob_qbarQbar (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in);
+
+//! Square of qg->qg Pure Jets Unordered backwards Scattering Current
+/**
+ * @param p1out Momentum of final state gluon
+ * @param p1in Momentum of initial state gluon
+ * @param pg Momentum of unordered gluon
+ * @param p2out Momentum of final state quark
+ * @param p2in Momentum of intial state quark
+ * @returns Square of the current contractions for qg->qg Scattering
+ *
+ * @note this can be used for unof contributions by inputting
+ * arguments appropriately.
+ */
+double ME_unob_qg (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in);
+
+//! Square of qbarg->qbarg Pure Jets Unordered backwards Scattering Current
+/**
+ * @param p1out Momentum of final state gluon
+ * @param p1in Momentum of initial state gluon
+ * @param pg Momentum of unordered gluon
+ * @param p2out Momentum of final state anti-quark
+ * @param p2in Momentum of intial state anti-quark
+ * @returns Square of the current contractions for qbarg->qbarg Scattering
+ *
+ * @note this can be used for unof contributions by inputting
+ * arguments appropriately.
+ */
+double ME_unob_qbarg (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in);
+
/** \class CCurrent jets.hh "include/HEJ/jets.hh"
* \brief This is the a new class structure for currents.
*/
class CCurrent
{
public:
CCurrent(COM sc0, COM sc1, COM sc2, COM sc3)
:c0(sc0),c1(sc1),c2(sc2),c3(sc3)
{};
CCurrent(const HLV p)
{
c0=p.e();
c1=p.px();
c2=p.py();
c3=p.pz();
};
CCurrent()
{};
CCurrent operator+(const CCurrent& other);
CCurrent operator-(const CCurrent& other);
CCurrent operator*(const double x);
CCurrent operator*(const COM x);
CCurrent operator/(const double x);
CCurrent operator/(const COM x);
friend std::ostream& operator<<(std::ostream& os, const CCurrent& cur);
COM dot(HLV p1);
COM dot(CCurrent p1);
COM c0,c1,c2,c3;
};
/* std::ostream& operator <<(std::ostream& os, const CCurrent& cur); */
CCurrent operator * ( double x, CCurrent& m);
CCurrent operator * ( COM x, CCurrent& m);
CCurrent operator / ( double x, CCurrent& m);
CCurrent operator / ( COM x, CCurrent& m);
//! Current <incoming state | mu | outgoing state>
/**
* This is a wrapper function around \see joi() note helicity flip to
* give same answer.
*/
void jio(HLV pin, bool helin, HLV pout, bool helout, current &cur);
//! Current <outgoing state | mu | outgoing state>
/**
* @param pi bra state momentum
* @param heli helicity of pi
* @param pj ket state momentum
* @param helj helicity of pj. (must be same as heli)
* @param cur reference to current which is saved.
*
* This function is for building <i (out)| mu |j (out)> currents. It
* must be called with pi as the bra, and pj as the ket.
*
* @TODO Remove heli/helj and just have helicity of current as argument.
*/
void joo(HLV pi, bool heli, HLV pj, bool helj, current &cur);
//! Current <outgoing state | mu | incoming state>
/**
* @param pout bra state momentum
* @param helout helicity of pout
* @param pin ket state momentum
* @param helin helicity of pin. (must be same as helout)
* @param cur reference to current which is saved.
*
* This function is for building <out| mu |in> currents. It must be
* called with pout as the bra, and pin as the ket. jio calls this
* with flipped helicity
*
* @TODO Remove helout/helin and just have helicity of current as argument.
*/
void joi(HLV pout, bool helout, HLV pin, bool helin, current &cur);
//! Current <outgoing state | mu | incoming state>
/**
* This is a wrapper function around the void function of the same name. \see joi
*/
CCurrent joi (HLV pout, bool helout, HLV pin, bool helin);
//! Current <incoming state | mu | outgoing state>
/**
* This is a wrapper function around the void function of the same name. \see jio
*/
CCurrent jio (HLV pout, bool helout, HLV pin, bool helin);
//! Current <outgoing state | mu | outgoing state>
/**
* This is a wrapper function around the void function of the same name. \see joo
*/
CCurrent joo (HLV pout, bool helout, HLV pin, bool helin);
inline COM cdot(const current & j1, const current & j2)
{
return j1[0]*j2[0]-j1[1]*j2[1]-j1[2]*j2[2]-j1[3]*j2[3];
}
inline COM cdot(const HLV & p, const current & j1) {
return j1[0]*p.e()-j1[1]*p.x()-j1[2]*p.y()-j1[3]*p.z();
}
inline void cmult(const COM & factor, const current & j1, current &cur)
{
cur[0]=factor*j1[0];
cur[1]=factor*j1[1];
cur[2]=factor*j1[2];
cur[3]=factor*j1[3];
}
// WHY!?!
inline void cadd(const current & j1, const current & j2, const current & j3,
const current & j4, const current & j5, current &sum)
{
sum[0]=j1[0]+j2[0]+j3[0]+j4[0]+j5[0];
sum[1]=j1[1]+j2[1]+j3[1]+j4[1]+j5[1];
sum[2]=j1[2]+j2[2]+j3[2]+j4[2]+j5[2];
sum[3]=j1[3]+j2[3]+j3[3]+j4[3]+j5[3];
}
inline void cadd(const current & j1, const current & j2, const current & j3,
const current & j4, current &sum) {
sum[0] = j1[0] + j2[0] + j3[0] + j4[0];
sum[1] = j1[1] + j2[1] + j3[1] + j4[1];
sum[2] = j1[2] + j2[2] + j3[2] + j4[2];
sum[3] = j1[3] + j2[3] + j3[3] + j4[3];
}
inline void cadd(const current & j1, const current & j2, const current & j3,
current &sum)
{
sum[0]=j1[0]+j2[0]+j3[0];
sum[1]=j1[1]+j2[1]+j3[1];
sum[2]=j1[2]+j2[2]+j3[2];
sum[3]=j1[3]+j2[3]+j3[3];
}
inline void cadd(const current & j1, const current & j2, current &sum)
{
sum[0]=j1[0]+j2[0];
sum[1]=j1[1]+j2[1];
sum[2]=j1[2]+j2[2];
sum[3]=j1[3]+j2[3];
}
inline double abs2(const COM & a)
{
return (a*conj(a)).real();
}
inline double vabs2(const CCurrent & cur)
{
return abs2(cur.c0)-abs2(cur.c1)-abs2(cur.c2)-abs2(cur.c3);
}
inline double vre(const CCurrent & a, const CCurrent & b)
{
return real(a.c0*conj(b.c0)-a.c1*conj(b.c1)-a.c2*conj(b.c2)-a.c3*conj(b.c3));
}
//! @TODO These are not currents and should be moved elsewhere.
double K_g(double p1minus, double paminus);
double K_g(HLV const & pout, HLV const & pin);
diff --git a/src/MatrixElement.cc b/src/MatrixElement.cc
index cbfba22..3fbdbe7 100644
--- a/src/MatrixElement.cc
+++ b/src/MatrixElement.cc
@@ -1,1704 +1,1773 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "HEJ/MatrixElement.hh"
#include <algorithm>
#include <assert.h>
#include <limits>
#include <math.h>
#include <stddef.h>
#include <unordered_map>
#include <utility>
#include "CLHEP/Vector/LorentzVector.h"
#include "fastjet/ClusterSequence.hh"
#include "HEJ/Constants.hh"
#include "HEJ/Wjets.hh"
#include "HEJ/Hjets.hh"
#include "HEJ/jets.hh"
#include "HEJ/PDG_codes.hh"
#include "HEJ/event_types.hh"
#include "HEJ/Event.hh"
#include "HEJ/exceptions.hh"
#include "HEJ/Particle.hh"
#include "HEJ/utility.hh"
namespace HEJ{
double MatrixElement::omega0(
double alpha_s, double mur,
fastjet::PseudoJet const & q_j
) const {
const double lambda = param_.regulator_lambda;
const double result = - alpha_s*N_C/M_PI*log(q_j.perp2()/(lambda*lambda));
if(! param_.log_correction) return result;
// use alpha_s(sqrt(q_j*lambda)), evolved to mur
return (
1. + alpha_s/(4.*M_PI)*beta0*log(mur*mur/(q_j.perp()*lambda))
)*result;
}
Weights MatrixElement::operator()(Event const & event) const {
return tree(event)*virtual_corrections(event);
}
Weights MatrixElement::tree(Event const & event) const {
return tree_param(event)*tree_kin(event);
}
Weights MatrixElement::tree_param(Event const & event) const {
if(! is_resummable(event.type())) {
return Weights{0., std::vector<double>(event.variations().size(), 0.)};
}
Weights result;
// only compute once for each renormalisation scale
std::unordered_map<double, double> known;
result.central = tree_param(event, event.central().mur);
known.emplace(event.central().mur, result.central);
for(auto const & var: event.variations()) {
const auto ME_it = known.find(var.mur);
if(ME_it == end(known)) {
const double wt = tree_param(event, var.mur);
result.variations.emplace_back(wt);
known.emplace(var.mur, wt);
}
else {
result.variations.emplace_back(ME_it->second);
}
}
return result;
}
Weights MatrixElement::virtual_corrections(Event const & event) const {
if(! is_resummable(event.type())) {
return Weights{0., std::vector<double>(event.variations().size(), 0.)};
}
Weights result;
// only compute once for each renormalisation scale
std::unordered_map<double, double> known;
result.central = virtual_corrections(event, event.central().mur);
known.emplace(event.central().mur, result.central);
for(auto const & var: event.variations()) {
const auto ME_it = known.find(var.mur);
if(ME_it == end(known)) {
const double wt = virtual_corrections(event, var.mur);
result.variations.emplace_back(wt);
known.emplace(var.mur, wt);
}
else {
result.variations.emplace_back(ME_it->second);
}
}
return result;
}
double MatrixElement::virtual_corrections_W(
Event const & event,
double mur,
Particle const & WBoson
) const{
auto const & in = event.incoming();
const auto partons = filter_partons(event.outgoing());
fastjet::PseudoJet const & pa = in.front().p;
#ifndef NDEBUG
fastjet::PseudoJet const & pb = in.back().p;
double const norm = (in.front().p + in.back().p).E();
#endif
assert(std::is_sorted(partons.begin(), partons.end(), rapidity_less{}));
assert(partons.size() >= 2);
assert(pa.pz() < pb.pz());
fastjet::PseudoJet q = pa - partons[0].p;
size_t first_idx = 0;
size_t last_idx = partons.size() - 1;
bool wc = true;
bool wqq = false;
// With extremal qqx or unordered gluon outside the extremal
// partons then it is not part of the FKL ladder and does not
// contribute to the virtual corrections. W emitted from the
// most backward leg must be taken into account in t-channel
if (event.type() == event_type::FKL) {
if (in[0].type != partons[0].type ){
q -= WBoson.p;
wc = false;
}
}
else if (event.type() == event_type::unob) {
q -= partons[1].p;
++first_idx;
if (in[0].type != partons[1].type ){
q -= WBoson.p;
wc = false;
}
}
else if (event.type() == event_type::qqxexb) {
q -= partons[1].p;
++first_idx;
if (abs(partons[0].type) != abs(partons[1].type)){
q -= WBoson.p;
wc = false;
}
}
if(event.type() == event_type::unof
|| event.type() == event_type::qqxexf){
--last_idx;
}
size_t first_idx_qqx = last_idx;
size_t last_idx_qqx = last_idx;
//if qqxMid event, virtual correction do not occur between
//qqx pair.
if(event.type() == event_type::qqxmid){
const auto backquark = std::find_if(
begin(partons) + 1, end(partons) - 1 ,
[](Particle const & s){ return (s.type != pid::gluon); }
);
if(backquark == end(partons) || (backquark+1)->type==pid::gluon) return 0;
if(abs(backquark->type) != abs((backquark+1)->type)) {
wqq=true;
wc=false;
}
last_idx = std::distance(begin(partons), backquark);
first_idx_qqx = last_idx+1;
}
double exponent = 0;
const double alpha_s = alpha_s_(mur);
for(size_t j = first_idx; j < last_idx; ++j){
exponent += omega0(alpha_s, mur, q)*(
partons[j+1].rapidity() - partons[j].rapidity()
);
q -=partons[j+1].p;
} // End Loop one
if (last_idx != first_idx_qqx) q -= partons[last_idx+1].p;
if (wqq) q -= WBoson.p;
for(size_t j = first_idx_qqx; j < last_idx_qqx; ++j){
exponent += omega0(alpha_s, mur, q)*(
partons[j+1].rapidity() - partons[j].rapidity()
);
q -= partons[j+1].p;
}
if (wc) q -= WBoson.p;
assert(
nearby(q, -1*pb, norm)
|| is_AWZH_boson(partons.back().type)
|| event.type() == event_type::unof
|| event.type() == event_type::qqxexf
);
return exp(exponent);
}
double MatrixElement::virtual_corrections(
Event const & event,
double mur
) const{
auto const & in = event.incoming();
auto const & out = event.outgoing();
fastjet::PseudoJet const & pa = in.front().p;
#ifndef NDEBUG
fastjet::PseudoJet const & pb = in.back().p;
double const norm = (in.front().p + in.back().p).E();
#endif
const auto AWZH_boson = std::find_if(
begin(out), end(out),
[](Particle const & p){ return is_AWZH_boson(p); }
);
if(AWZH_boson != end(out) && abs(AWZH_boson->type) == pid::Wp){
return virtual_corrections_W(event, mur, *AWZH_boson);
}
assert(std::is_sorted(out.begin(), out.end(), rapidity_less{}));
assert(out.size() >= 2);
assert(pa.pz() < pb.pz());
fastjet::PseudoJet q = pa - out[0].p;
size_t first_idx = 0;
size_t last_idx = out.size() - 1;
// if there is a Higgs boson, extremal qqx or unordered gluon
// outside the extremal partons then it is not part of the FKL
// ladder and does not contribute to the virtual corrections
if((out.front().type == pid::Higgs)
|| event.type() == event_type::unob
|| event.type() == event_type::qqxexb){
q -= out[1].p;
++first_idx;
}
if((out.back().type == pid::Higgs)
|| event.type() == event_type::unof
|| event.type() == event_type::qqxexf){
--last_idx;
}
size_t first_idx_qqx = last_idx;
size_t last_idx_qqx = last_idx;
//if qqxMid event, virtual correction do not occur between
//qqx pair.
if(event.type() == event_type::qqxmid){
const auto backquark = std::find_if(
begin(out) + 1, end(out) - 1 ,
[](Particle const & s){ return (s.type != pid::gluon && is_parton(s.type)); }
);
if(backquark == end(out) || (backquark+1)->type==pid::gluon) return 0;
last_idx = std::distance(begin(out), backquark);
first_idx_qqx = last_idx+1;
}
double exponent = 0;
const double alpha_s = alpha_s_(mur);
for(size_t j = first_idx; j < last_idx; ++j){
exponent += omega0(alpha_s, mur, q)*(
out[j+1].rapidity() - out[j].rapidity()
);
q -= out[j+1].p;
}
if (last_idx != first_idx_qqx) q -= out[last_idx+1].p;
for(size_t j = first_idx_qqx; j < last_idx_qqx; ++j){
exponent += omega0(alpha_s, mur, q)*(
out[j+1].rapidity() - out[j].rapidity()
);
q -= out[j+1].p;
}
assert(
nearby(q, -1*pb, norm)
|| out.back().type == pid::Higgs
|| event.type() == event_type::unof
|| event.type() == event_type::qqxexf
);
return exp(exponent);
}
} // namespace HEJ
namespace {
//! Lipatov vertex for partons emitted into extremal jets
double C2Lipatov(
CLHEP::HepLorentzVector const & qav,
CLHEP::HepLorentzVector const & qbv,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & p2
){
CLHEP::HepLorentzVector temptrans=-(qav+qbv);
CLHEP::HepLorentzVector p5=qav-qbv;
CLHEP::HepLorentzVector CL=temptrans
+ p1*(qav.m2()/p5.dot(p1) + 2.*p5.dot(p2)/p1.dot(p2))
- p2*(qbv.m2()/p5.dot(p2) + 2.*p5.dot(p1)/p1.dot(p2));
return -CL.dot(CL);
}
//! Lipatov vertex with soft subtraction for partons emitted into extremal jets
double C2Lipatovots(
CLHEP::HepLorentzVector const & qav,
CLHEP::HepLorentzVector const & qbv,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & p2,
double lambda
) {
double kperp=(qav-qbv).perp();
if (kperp>lambda)
return C2Lipatov(qav, qbv, p1, p2)/(qav.m2()*qbv.m2());
double Cls=(C2Lipatov(qav, qbv, p1, p2)/(qav.m2()*qbv.m2()));
return Cls-4./(kperp*kperp);
}
//! Lipatov vertex
double C2Lipatov( // B
CLHEP::HepLorentzVector const & qav,
CLHEP::HepLorentzVector const & qbv,
CLHEP::HepLorentzVector const & pim,
CLHEP::HepLorentzVector const & pip,
CLHEP::HepLorentzVector const & pom,
CLHEP::HepLorentzVector const & pop
){
CLHEP::HepLorentzVector temptrans=-(qav+qbv);
CLHEP::HepLorentzVector p5=qav-qbv;
CLHEP::HepLorentzVector CL=temptrans
+ qav.m2()*(1./p5.dot(pip)*pip + 1./p5.dot(pop)*pop)/2.
- qbv.m2()*(1./p5.dot(pim)*pim + 1./p5.dot(pom)*pom)/2.
+ ( pip*(p5.dot(pim)/pip.dot(pim) + p5.dot(pom)/pip.dot(pom))
+ pop*(p5.dot(pim)/pop.dot(pim) + p5.dot(pom)/pop.dot(pom))
- pim*(p5.dot(pip)/pip.dot(pim) + p5.dot(pop)/pop.dot(pim))
- pom*(p5.dot(pip)/pip.dot(pom) + p5.dot(pop)/pop.dot(pom)) )/2.;
return -CL.dot(CL);
}
//! Lipatov vertex with soft subtraction
double C2Lipatovots(
CLHEP::HepLorentzVector const & qav,
CLHEP::HepLorentzVector const & qbv,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & p2,
double lambda
) {
double kperp=(qav-qbv).perp();
if (kperp>lambda)
return C2Lipatov(qav, qbv, pa, pb, p1, p2)/(qav.m2()*qbv.m2());
double Cls=(C2Lipatov(qav, qbv, pa, pb, p1, p2)/(qav.m2()*qbv.m2()));
double temp=Cls-4./(kperp*kperp);
return temp;
}
/** Matrix element squared for tree-level current-current scattering
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
+ * @param pg Unordered gluon momentum
+ * @param pn Particle n Momentum
+ * @param pb Particle b Momentum
+ * @param p1 Particle 1 Momentum
+ * @param pa Particle a Momentum
+ * @returns ME Squared for Tree-Level Current-Current Scattering
+ */
+ double ME_uno_current(
+ int aptype, int bptype,
+ CLHEP::HepLorentzVector const & pg,
+ CLHEP::HepLorentzVector const & pn,
+ CLHEP::HepLorentzVector const & pb,
+ CLHEP::HepLorentzVector const & p1,
+ CLHEP::HepLorentzVector const & pa
+ ){
+ assert(aptype!=21); // aptype cannot be gluon
+ if (bptype==21) {
+ if (aptype > 0)
+ return ME_unob_qg(pg,p1,pa,pn,pb);
+ else
+ return ME_unob_qbarg(pg,p1,pa,pn,pb);
+ }
+ else if (bptype<0) { // ----- || -----
+ if (aptype > 0)
+ return ME_unob_qQbar(pg,p1,pa,pn,pb);
+ else
+ return ME_unob_qbarQbar(pg,p1,pa,pn,pb);
+ }
+ else { //bptype == quark
+ if (aptype > 0)
+ return ME_unob_qQ(pg,p1,pa,pn,pb);
+ else
+ return ME_unob_qbarQ(pg,p1,pa,pn,pb);
+ }
+ }
+
+ /** Matrix element squared for tree-level current-current scattering
+ * @param aptype Particle a PDG ID
+ * @param bptype Particle b PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @returns ME Squared for Tree-Level Current-Current Scattering
*/
double ME_current(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa
){
if (aptype==21&&bptype==21) {
return ME_gg(pn,pb,p1,pa);
} else if (aptype==21&&bptype!=21) {
if (bptype > 0)
return ME_qg(pn,pb,p1,pa);
else
return ME_qbarg(pn,pb,p1,pa);
}
else if (bptype==21&&aptype!=21) { // ----- || -----
if (aptype > 0)
return ME_qg(p1,pa,pn,pb);
else
return ME_qbarg(p1,pa,pn,pb);
}
else { // they are both quark
if (bptype>0) {
if (aptype>0)
return ME_qQ(pn,pb,p1,pa);
else
return ME_qQbar(pn,pb,p1,pa);
}
else {
if (aptype>0)
return ME_qQbar(p1,pa,pn,pb);
else
return ME_qbarQbar(pn,pb,p1,pa);
}
}
throw std::logic_error("unknown particle types");
}
/** Matrix element squared for tree-level current-current scattering With W+Jets
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param wc Boolean. True->W Emitted from b. Else; emitted from leg a
* @returns ME Squared for Tree-Level Current-Current Scattering
*/
double ME_W_current(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const wc
){
// We know it cannot be gg incoming.
assert(!(aptype==21 && bptype==21));
if (aptype==21&&bptype!=21) {
if (bptype > 0)
return ME_W_qg(pn,plbar,pl,pb,p1,pa);
else
return ME_W_qbarg(pn,plbar,pl,pb,p1,pa);
}
else if (bptype==21&&aptype!=21) { // ----- || -----
if (aptype > 0)
return ME_W_qg(p1,plbar,pl,pa,pn,pb);
else
return ME_W_qbarg(p1,plbar,pl,pa,pn,pb);
}
else { // they are both quark
if (wc==true){ // emission off b, (first argument pbout)
if (bptype>0) {
if (aptype>0)
return ME_W_qQ(pn,plbar,pl,pb,p1,pa);
else
return ME_W_qQbar(pn,plbar,pl,pb,p1,pa);
}
else {
if (aptype>0)
return ME_W_qbarQ(pn,plbar,pl,pb,p1,pa);
else
return ME_W_qbarQbar(pn,plbar,pl,pb,p1,pa);
}
}
else{ // emission off a, (first argument paout)
if (aptype > 0) {
if (bptype > 0)
return ME_W_qQ(p1,plbar,pl,pa,pn,pb);
else
return ME_W_qQbar(p1,plbar,pl,pa,pn,pb);
}
else { // a is anti-quark
if (bptype > 0)
return ME_W_qbarQ(p1,plbar,pl,pa,pn,pb);
else
return ME_W_qbarQbar(p1,plbar,pl,pa,pn,pb);
}
}
}
throw std::logic_error("unknown particle types");
}
/** Matrix element squared for backwards uno tree-level current-current
* scattering With W+Jets
*
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param pg Unordered gluon momentum
* @param wc Boolean. True->W Emitted from b. Else; emitted from leg a
* @returns ME Squared for unob Tree-Level Current-Current Scattering
*/
double ME_W_unob_current(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pg,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const wc
){
// we know they are not both gluons
if (bptype == 21 && aptype != 21) { // b gluon => W emission off a
if (aptype > 0)
return ME_Wuno_qg(p1,pa,pn,pb,pg,plbar,pl);
else
return ME_Wuno_qbarg(p1,pa,pn,pb,pg,plbar,pl);
}
else { // they are both quark
if (wc) {// emission off b, i.e. b is first current
if (bptype>0){
if (aptype>0)
return ME_W_unob_qQ(p1,pa,pn,pb,pg,plbar,pl);
else
return ME_W_unob_qQbar(p1,pa,pn,pb,pg,plbar,pl);
}
else{
if (aptype>0)
return ME_W_unob_qbarQ(p1,pa,pn,pb,pg,plbar,pl);
else
return ME_W_unob_qbarQbar(p1,pa,pn,pb,pg,plbar,pl);
}
}
else {// wc == false, emission off a, i.e. a is first current
if (aptype > 0) {
if (bptype > 0) //qq
return ME_Wuno_qQ(p1,pa,pn,pb,pg,plbar,pl);
else //qqbar
return ME_Wuno_qQbar(p1,pa,pn,pb,pg,plbar,pl);
}
else { // a is anti-quark
if (bptype > 0) //qbarq
return ME_Wuno_qbarQ(p1,pa,pn,pb,pg,plbar,pl);
else //qbarqbar
return ME_Wuno_qbarQbar(p1,pa,pn,pb,pg,plbar,pl);
}
}
}
throw std::logic_error("unknown particle types");
}
/** Matrix element squared for uno forward tree-level current-current
* scattering With W+Jets
*
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param pg Unordered gluon momentum
* @param wc Boolean. True->W Emitted from b. Else; emitted from leg a
* @returns ME Squared for unof Tree-Level Current-Current Scattering
*/
double ME_W_unof_current(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pg,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const wc
){
// we know they are not both gluons
if (aptype==21 && bptype!=21) {//a gluon => W emission off b
if (bptype > 0)
return ME_Wuno_qg(pn, pb, p1, pa, pg, plbar, pl);
else
return ME_Wuno_qbarg(pn, pb, p1, pa, pg, plbar, pl);
}
else { // they are both quark
if (wc) {// emission off b, i.e. b is first current
if (bptype>0){
if (aptype>0)
return ME_Wuno_qQ(pn, pb, p1, pa, pg, plbar, pl);
else
return ME_Wuno_qQbar(pn, pb, p1, pa, pg, plbar, pl);
}
else{
if (aptype>0)
return ME_Wuno_qbarQ(pn, pb, p1, pa, pg, plbar, pl);
else
return ME_Wuno_qbarQbar(pn, pb, p1, pa, pg, plbar, pl);
}
}
else {// wc == false, emission off a, i.e. a is first current
if (aptype > 0) {
if (bptype > 0) //qq
return ME_W_unof_qQ(p1, pa, pn, pb, pg, plbar, pl);
- // return ME_W_unof_qQ(pn,pb,p1,pa,pg,plbar,pl);
else //qqbar
return ME_W_unof_qQbar(p1, pa, pn, pb, pg, plbar, pl);
}
else { // a is anti-quark
if (bptype > 0) //qbarq
return ME_W_unof_qbarQ(p1, pa, pn, pb, pg, plbar, pl);
else //qbarqbar
return ME_W_unof_qbarQbar(p1, pa, pn, pb, pg, plbar, pl);
}
}
}
throw std::logic_error("unknown particle types");
}
/** \brief Matrix element squared for backward qqx tree-level current-current
* scattering With W+Jets
*
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pa Initial state a Momentum
* @param pb Initial state b Momentum
* @param pq Final state q Momentum
* @param pqbar Final state qbar Momentum
* @param pn Final state n Momentum
* @param plbar Final state anti-lepton momentum
* @param pl Final state lepton momentum
* @param wc Boolean. True->W Emitted from b. Else; emitted from leg a
* @returns ME Squared for qqxb Tree-Level Current-Current Scattering
*/
double ME_W_qqxb_current(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & pq,
CLHEP::HepLorentzVector const & pqbar,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const wc
){
// CAM factors for the qqx amps, and qqbar ordering (default, qbar extremal)
bool swapQuarkAntiquark=false;
double CFbackward;
if (pqbar.rapidity() > pq.rapidity()){
swapQuarkAntiquark=true;
CFbackward = (0.5*(3.-1./3.)*(pa.minus()/(pq.minus())+(pq.minus())/pa.minus())+1./3.)*3./4.;
}
else{
CFbackward = (0.5*(3.-1./3.)*(pa.minus()/(pqbar.minus())+(pqbar.minus())/pa.minus())+1./3.)*3./4.;
}
// With qqbar we could have 2 incoming gluons and W Emission
if (aptype==21&&bptype==21) {//a gluon, b gluon gg->qqbarWg
// This will be a wqqx emission as there is no other possible W Emission Site.
if (swapQuarkAntiquark)
return ME_WExqqx_qqbarg(pa, pqbar, plbar, pl, pq, pn, pb)*CFbackward;
else
return ME_WExqqx_qbarqg(pa, pq, plbar, pl, pqbar, pn, pb)*CFbackward;
}
else if (aptype==21&&bptype!=21 ) {//a gluon => W emission off b leg or qqx
if (!wc){ // W Emitted from backwards qqx
if (swapQuarkAntiquark)
return ME_WExqqx_qqbarQ(pa, pqbar, plbar, pl, pq, pn, pb)*CFbackward;
else
return ME_WExqqx_qbarqQ(pa, pq, plbar, pl, pqbar, pn, pb)*CFbackward;
}
else { // W Must be emitted from forwards leg.
if(bptype > 0){
if (swapQuarkAntiquark)
return ME_W_Exqqx_QQq(pb, pa, pn, pqbar, pq, plbar, pl, false)*CFbackward;
else
return ME_W_Exqqx_QQq(pb, pa, pn, pq, pqbar, plbar, pl, false)*CFbackward;
} else {
if (swapQuarkAntiquark)
return ME_W_Exqqx_QQq(pb, pa, pn, pqbar, pq, plbar, pl, true)*CFbackward;
else
return ME_W_Exqqx_QQq(pb, pa, pn, pq, pqbar, plbar, pl, true)*CFbackward;
}
}
}
throw std::logic_error("Incompatible incoming particle types with qqxb");
}
/* \brief Matrix element squared for forward qqx tree-level current-current
* scattering With W+Jets
*
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pa Initial state a Momentum
* @param pb Initial state b Momentum
* @param pq Final state q Momentum
* @param pqbar Final state qbar Momentum
* @param p1 Final state 1 Momentum
* @param plbar Final state anti-lepton momentum
* @param pl Final state lepton momentum
* @param wc Boolean. True->W Emitted from b. Else; emitted from leg a
* @returns ME Squared for qqxf Tree-Level Current-Current Scattering
*/
double ME_W_qqxf_current(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & pq,
CLHEP::HepLorentzVector const & pqbar,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const wc
){
// CAM factors for the qqx amps, and qqbar ordering (default, qbar extremal)
bool swapQuarkAntiquark=false;
double CFforward;
if (pqbar.rapidity() < pq.rapidity()){
swapQuarkAntiquark=true;
CFforward = (0.5*(3.-1./3.)*(pb.plus()/(pq.plus())+(pq.plus())/pb.plus())+1./3.)*3./4.;
}
else{
CFforward = (0.5*(3.-1./3.)*(pb.plus()/(pqbar.plus())+(pqbar.plus())/pb.plus())+1./3.)*3./4.;
}
// With qqbar we could have 2 incoming gluons and W Emission
if (aptype==21&&bptype==21) {//a gluon, b gluon gg->qqbarWg
// This will be a wqqx emission as there is no other possible W Emission Site.
if (swapQuarkAntiquark)
return ME_WExqqx_qqbarg(pb, pqbar, plbar, pl, pq, p1, pa)*CFforward;
else
return ME_WExqqx_qbarqg(pb, pq, plbar, pl, pqbar, p1, pa)*CFforward;
}
else if (bptype==21&&aptype!=21) {// b gluon => W emission off a or qqx
if (wc){ // W Emitted from forwards qqx
if (swapQuarkAntiquark)
return ME_WExqqx_qqbarQ(pb, pqbar, plbar, pl, pq, p1, pa)*CFforward;
else
return ME_WExqqx_qbarqQ(pb, pq, plbar, pl, pqbar, p1, pa)*CFforward;
}
// W Must be emitted from backwards leg.
if (aptype > 0){
if (swapQuarkAntiquark)
return ME_W_Exqqx_QQq(pa, pb, p1, pqbar, pq, plbar, pl, false)*CFforward;
else
return ME_W_Exqqx_QQq(pa, pb, p1, pq, pqbar, plbar, pl, false)*CFforward;
} else {
if (swapQuarkAntiquark)
return ME_W_Exqqx_QQq(pa, pb, p1, pqbar, pq, plbar, pl, true)*CFforward;
else
return ME_W_Exqqx_QQq(pa, pb, p1, pq, pqbar, plbar, pl, true)*CFforward;
}
}
throw std::logic_error("Incompatible incoming particle types with qqxf");
}
/* \brief Matrix element squared for central qqx tree-level current-current
* scattering With W+Jets
*
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param nabove Number of gluons emitted before central qqxpair
* @param nbelow Number of gluons emitted after central qqxpair
* @param pa Initial state a Momentum
* @param pb Initial state b Momentum\
* @param pq Final state qbar Momentum
* @param pqbar Final state q Momentum
* @param partons Vector of all outgoing partons
* @param plbar Final state anti-lepton momentum
* @param pl Final state lepton momentum
* @param wqq Boolean. True siginfies W boson is emitted from Central qqx
* @param wc Boolean. wc=true signifies w boson emitted from leg b; if wqq=false.
* @returns ME Squared for qqxmid Tree-Level Current-Current Scattering
*/
double ME_W_qqxmid_current(
int aptype, int bptype,
int nabove, int nbelow,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & pq,
CLHEP::HepLorentzVector const & pqbar,
std::vector<HLV> const & partons,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const wqq, bool const wc
){
// CAM factors for the qqx amps, and qqbar ordering (default, pq backwards)
bool swapQuarkAntiquark=false;
if (pqbar.rapidity() < pq.rapidity()){
swapQuarkAntiquark=true;
}
double wt=1.;
if (aptype==21) wt*=K_g(partons.front(),pa)/HEJ::C_F;
if (bptype==21) wt*=K_g(partons.back(),pb)/HEJ::C_F;
if (aptype <=0 && bptype <=0){ // Both External AntiQuark
if (wqq==1){//emission from central qqbar
return wt*ME_WCenqqx_qq(pa, pb, pl, plbar, partons,true,true,
swapQuarkAntiquark, nabove);
}
else if (wc==1){//emission from b leg
return wt*ME_W_Cenqqx_qq(pa, pb, pl, plbar, partons, true,true,
swapQuarkAntiquark, nabove, nbelow, true);
}
else { // emission from a leg
return wt*ME_W_Cenqqx_qq(pa, pb, pl,plbar, partons, true,true,
swapQuarkAntiquark, nabove, nbelow, false);
}
} // end both antiquark
else if (aptype<=0){ // a is antiquark
if (wqq==1){//emission from central qqbar
return wt*ME_WCenqqx_qq(pa, pb, pl, plbar, partons, false, true,
swapQuarkAntiquark, nabove);
}
else if (wc==1){//emission from b leg
return wt*ME_W_Cenqqx_qq(pa, pb, pl, plbar, partons,false,true,
swapQuarkAntiquark, nabove, nbelow, true);
}
else { // emission from a leg
return wt*ME_W_Cenqqx_qq(pa, pb, pl, plbar, partons, false, true,
swapQuarkAntiquark, nabove, nbelow, false);
}
} // end a is antiquark
else if (bptype<=0){ // b is antiquark
if (wqq==1){//emission from central qqbar
return wt*ME_WCenqqx_qq(pa, pb, pl, plbar, partons, true, false,
swapQuarkAntiquark, nabove);
}
else if (wc==1){//emission from b leg
return wt*ME_W_Cenqqx_qq(pa, pb, pl, plbar, partons, true, false,
swapQuarkAntiquark, nabove, nbelow, true);
}
else { // emission from a leg
return wt*ME_W_Cenqqx_qq(pa, pb, pl, plbar, partons, true, false,
swapQuarkAntiquark, nabove, nbelow, false);
}
} //end b is antiquark
else{ //Both Quark or gluon
if (wqq==1){//emission from central qqbar
return wt*ME_WCenqqx_qq(pa, pb, pl, plbar, partons, false, false,
swapQuarkAntiquark, nabove);}
else if (wc==1){//emission from b leg
return wt*ME_W_Cenqqx_qq(pa, pb, pl, plbar, partons, false, false,
swapQuarkAntiquark, nabove, nbelow, true);
}
else { // emission from a leg
return wt*ME_W_Cenqqx_qq(pa, pb, pl, plbar, partons, false, false,
swapQuarkAntiquark, nabove, nbelow, false);
}
}
throw std::logic_error("unknown particle types");
}
/** \brief Matrix element squared for tree-level current-current scattering with Higgs
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param qH t-channel momentum before Higgs
* @param qHp1 t-channel momentum after Higgs
* @returns ME Squared for Tree-Level Current-Current Scattering with Higgs
*/
double ME_Higgs_current(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & qH, // t-channel momentum before Higgs
CLHEP::HepLorentzVector const & qHp1, // t-channel momentum after Higgs
double mt, bool include_bottom, double mb
){
if (aptype==21&&bptype==21) // gg initial state
return ME_H_gg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else if (aptype==21&&bptype!=21) {
if (bptype > 0)
return ME_H_qg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4./9.;
else
return ME_H_qbarg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4./9.;
}
else if (bptype==21&&aptype!=21) {
if (aptype > 0)
return ME_H_qg(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb)*4./9.;
else
return ME_H_qbarg(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb)*4./9.;
}
else { // they are both quark
if (bptype>0) {
if (aptype>0)
return ME_H_qQ(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4.*4./(9.*9.);
else
return ME_H_qQbar(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4.*4./(9.*9.);
}
else {
if (aptype>0)
return ME_H_qQbar(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb)*4.*4./(9.*9.);
else
return ME_H_qbarQbar(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4.*4./(9.*9.);
}
}
throw std::logic_error("unknown particle types");
}
/** \brief Current matrix element squared with Higgs and unordered backward emission
* @param aptype Particle A PDG ID
* @param bptype Particle B PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param pg Unordered back Particle Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param qH t-channel momentum before Higgs
* @param qHp1 t-channel momentum after Higgs
* @returns ME Squared with Higgs and unordered backward emission
*
* @note This function assumes unordered gluon backwards from pa-p1 current.
* For unof, reverse call order
*/
double ME_Higgs_current_uno(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pg,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & qH, // t-channel momentum before Higgs
CLHEP::HepLorentzVector const & qHp1, // t-channel momentum after Higgs
double mt, bool include_bottom, double mb
){
if (bptype==21&&aptype!=21) {
if (aptype > 0)
return ME_H_unob_gQ(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb);
else
return ME_H_unob_gQbar(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb);
}
else { // they are both quark
if (aptype>0) {
if (bptype>0)
return ME_H_unob_qQ(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb);
else
return ME_H_unob_qbarQ(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb);
}
else {
if (bptype>0)
return ME_H_unob_qQbar(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb);
else
return ME_H_unob_qbarQbar(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb);
}
}
throw std::logic_error("unknown particle types");
}
CLHEP::HepLorentzVector to_HepLorentzVector(HEJ::Particle const & particle){
return {particle.p.px(), particle.p.py(), particle.p.pz(), particle.p.E()};
}
void validate(HEJ::MatrixElementConfig const & config) {
#ifndef HEJ_BUILD_WITH_QCDLOOP
if(!config.Higgs_coupling.use_impact_factors) {
throw std::invalid_argument{
"Invalid Higgs coupling settings.\n"
"HEJ without QCDloop support can only use impact factors.\n"
"Set use_impact_factors to true or recompile HEJ.\n"
};
}
#endif
if(config.Higgs_coupling.use_impact_factors
&& config.Higgs_coupling.mt != std::numeric_limits<double>::infinity()) {
throw std::invalid_argument{
"Conflicting settings: "
"impact factors may only be used in the infinite top mass limit"
};
}
}
} // namespace anonymous
namespace HEJ{
MatrixElement::MatrixElement(
std::function<double (double)> alpha_s,
MatrixElementConfig conf
):
alpha_s_{std::move(alpha_s)},
param_{std::move(conf)}
{
validate(param_);
}
double MatrixElement::tree_kin(
Event const & ev
) const {
if(! is_resummable(ev.type())) return 0.;
auto AWZH_boson = std::find_if(
begin(ev.outgoing()), end(ev.outgoing()),
[](Particle const & p){return is_AWZH_boson(p);}
);
if(AWZH_boson == end(ev.outgoing()))
return tree_kin_jets(ev);
switch(AWZH_boson->type){
case pid::Higgs:
return tree_kin_Higgs(ev);
case pid::Wp:
case pid::Wm:
return tree_kin_W(ev);
// TODO
case pid::photon:
case pid::Z:
default:
throw not_implemented("Emission of boson of unsupported type");
}
}
namespace{
constexpr int extremal_jet_idx = 1;
constexpr int no_extremal_jet_idx = 0;
bool treat_as_extremal(Particle const & parton){
return parton.p.user_index() == extremal_jet_idx;
}
template<class InputIterator>
double FKL_ladder_weight(
InputIterator begin_gluon, InputIterator end_gluon,
CLHEP::HepLorentzVector const & q0,
CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pn,
double lambda
){
double wt = 1;
auto qi = q0;
for(auto gluon_it = begin_gluon; gluon_it != end_gluon; ++gluon_it){
assert(gluon_it->type == pid::gluon);
const auto g = to_HepLorentzVector(*gluon_it);
const auto qip1 = qi - g;
if(treat_as_extremal(*gluon_it)){
wt *= C2Lipatovots(qip1, qi, pa, pb, lambda)*C_A;
} else{
wt *= C2Lipatovots(qip1, qi, pa, pb, p1, pn, lambda)*C_A;
}
qi = qip1;
}
return wt;
}
} // namespace anonymous
std::vector<Particle> MatrixElement::tag_extremal_jet_partons(
Event const & ev
) const{
auto out_partons = filter_partons(ev.outgoing());
if(out_partons.size() == ev.jets().size()){
// no additional emissions in extremal jets, don't need to tag anything
for(auto & parton: out_partons){
parton.p.set_user_index(no_extremal_jet_idx);
}
return out_partons;
}
// TODO: avoid reclustering
fastjet::ClusterSequence cs(to_PseudoJet(out_partons), ev.jet_def());
const auto jets = sorted_by_rapidity(cs.inclusive_jets(ev.min_jet_pt()));
assert(jets.size() >= 2);
auto most_backward = begin(jets);
auto most_forward = end(jets) - 1;
// skip jets caused by unordered emission or qqx
if(ev.type() == event_type::unob || ev.type() == event_type::qqxexb){
assert(jets.size() >= 3);
++most_backward;
}
else if(ev.type() == event_type::unof || ev.type() == event_type::qqxexf){
assert(jets.size() >= 3);
--most_forward;
}
const auto extremal_jet_indices = cs.particle_jet_indices(
{*most_backward, *most_forward}
);
assert(extremal_jet_indices.size() == out_partons.size());
for(size_t i = 0; i < out_partons.size(); ++i){
assert(HEJ::is_parton(out_partons[i]));
const int idx = (extremal_jet_indices[i]>=0)?
extremal_jet_idx:
no_extremal_jet_idx;
out_partons[i].p.set_user_index(idx);
}
return out_partons;
}
+ namespace {
+ template<class InIter, class partIter>
+ double tree_kin_jets_uno(InIter BeginIn, InIter EndIn, partIter BeginPart,
+ partIter EndPart, double lambda){
+
+ const auto pa = to_HepLorentzVector(*BeginIn);
+ const auto pb = to_HepLorentzVector(*(EndIn-1));
+
+ const auto pg = to_HepLorentzVector(*BeginPart);
+ const auto p1 = to_HepLorentzVector(*(BeginPart+1));
+ const auto pn = to_HepLorentzVector(*(EndPart-1));
+
+ const double current_factor = ME_uno_current(
+ (BeginIn)->type, (EndIn-1)->type, pg, pn, pb, p1, pa
+ )/(4.*(N_C*N_C - 1.));
+ const double ladder_factor = FKL_ladder_weight(
+ (BeginPart+2), (EndPart-1),
+ pa-p1-pg, pa, pb, p1, pn, lambda
+ );
+
+ return current_factor*ladder_factor;
+ }
+ }
+
double MatrixElement::tree_kin_jets(Event const & ev) const {
auto const & incoming = ev.incoming();
const auto partons = tag_extremal_jet_partons(ev);
- if(is_uno(ev.type())){
- throw not_implemented("unordered emission not implemented for pure jets");
- }
+ if (ev.type()==HEJ::event_type::unordered_backward){
+ return tree_kin_jets_uno(incoming.begin(), incoming.end(),
+ partons.begin(), partons.end(),
+ param_.regulator_lambda);
+ }
+ else if (ev.type()==HEJ::event_type::unordered_forward){
+ return tree_kin_jets_uno(incoming.rbegin(), incoming.rend(),
+ partons.rbegin(), partons.rend(),
+ param_.regulator_lambda);
+ }
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
const auto p1 = to_HepLorentzVector(partons.front());
const auto pn = to_HepLorentzVector(partons.back());
return ME_current(
incoming[0].type, incoming[1].type,
pn, pb, p1, pa
)/(4.*(N_C*N_C - 1.))*FKL_ladder_weight(
begin(partons) + 1, end(partons) - 1,
pa - p1, pa, pb, p1, pn,
param_.regulator_lambda
);
}
namespace{
double tree_kin_W_FKL(
int aptype, int bptype, HLV pa, HLV pb,
std::vector<Particle> const & partons,
HLV plbar, HLV pl,
double lambda
){
auto p1 = to_HepLorentzVector(partons[0]);
auto pn = to_HepLorentzVector(partons[partons.size() - 1]);
const auto begin_ladder = cbegin(partons) + 1;
const auto end_ladder = cend(partons) - 1;
bool wc = aptype==partons[0].type; //leg b emits w
auto q0 = pa - p1;
if(!wc)
q0 -= pl + plbar;
const double current_factor = ME_W_current(
aptype, bptype, pn, pb,
p1, pa, plbar, pl, wc
);
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder,
q0, pa, pb, p1, pn,
lambda
);
return current_factor*ladder_factor;
}
double tree_kin_W_unob(
int aptype, int bptype, HLV pa, HLV pb,
std::vector<Particle> const & partons,
HLV plbar, HLV pl,
double lambda
){
auto pg = to_HepLorentzVector(partons[0]);
auto p1 = to_HepLorentzVector(partons[1]);
auto pn = to_HepLorentzVector(partons[partons.size() - 1]);
const auto begin_ladder = cbegin(partons) + 2;
const auto end_ladder = cend(partons) - 1;
bool wc = aptype==partons[1].type; //leg b emits w
auto q0 = pa - p1 -pg;
if(!wc)
q0 -= pl + plbar;
const double current_factor = ME_W_unob_current(
aptype, bptype, pn, pb,
p1, pa, pg, plbar, pl, wc
);
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder,
q0, pa, pb, p1, pn,
lambda
);
return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor;
}
double tree_kin_W_unof(
int aptype, int bptype, HLV pa, HLV pb,
std::vector<Particle> const & partons,
HLV plbar, HLV pl,
double lambda
){
auto p1 = to_HepLorentzVector(partons[0]);
auto pn = to_HepLorentzVector(partons[partons.size() - 2]);
auto pg = to_HepLorentzVector(partons[partons.size() - 1]);
const auto begin_ladder = cbegin(partons) + 1;
const auto end_ladder = cend(partons) - 2;
bool wc = aptype==partons[0].type; //leg b emits w
auto q0 = pa - p1;
if(!wc)
q0 -= pl + plbar;
const double current_factor = ME_W_unof_current(
aptype, bptype, pn, pb,
p1, pa, pg, plbar, pl, wc
);
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder,
q0, pa, pb, p1, pn,
lambda
);
return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor;
}
double tree_kin_W_qqxb(
int aptype, int bptype, HLV pa, HLV pb,
std::vector<Particle> const & partons,
HLV plbar, HLV pl,
double lambda
){
HLV pq,pqbar;
if(is_quark(partons[0])){
pq = to_HepLorentzVector(partons[0]);
pqbar = to_HepLorentzVector(partons[1]);
}
else{
pq = to_HepLorentzVector(partons[1]);
pqbar = to_HepLorentzVector(partons[0]);
}
auto p1 = to_HepLorentzVector(partons[0]);
auto pn = to_HepLorentzVector(partons[partons.size() - 1]);
const auto begin_ladder = cbegin(partons) + 2;
const auto end_ladder = cend(partons) - 1;
bool wc = bptype!=partons.back().type; //leg b emits w
auto q0 = pa - pq - pqbar;
if(!wc)
q0 -= pl + plbar;
const double current_factor = ME_W_qqxb_current(
aptype, bptype, pa, pb,
pq, pqbar, pn, plbar, pl, wc
);
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder,
q0, pa, pb, p1, pn,
lambda
);
return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor;
}
double tree_kin_W_qqxf(
int aptype, int bptype, HLV pa, HLV pb,
std::vector<Particle> const & partons,
HLV plbar, HLV pl,
double lambda
){
HLV pq,pqbar;
if(is_quark(partons[partons.size() - 1])){
pq = to_HepLorentzVector(partons[partons.size() - 1]);
pqbar = to_HepLorentzVector(partons[partons.size() - 2]);
}
else{
pq = to_HepLorentzVector(partons[partons.size() - 2]);
pqbar = to_HepLorentzVector(partons[partons.size() - 1]);
}
auto p1 = to_HepLorentzVector(partons[0]);
auto pn = to_HepLorentzVector(partons[partons.size() - 1]);
const auto begin_ladder = cbegin(partons) + 1;
const auto end_ladder = cend(partons) - 2;
bool wc = aptype==partons.front().type; //leg b emits w
auto q0 = pa - p1;
if(!wc)
q0 -= pl + plbar;
const double current_factor = ME_W_qqxf_current(
aptype, bptype, pa, pb,
pq, pqbar, p1, plbar, pl, wc
);
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder,
q0, pa, pb, p1, pn,
lambda
);
return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor;
}
double tree_kin_W_qqxmid(
int aptype, int bptype, HLV pa, HLV pb,
std::vector<Particle> const & partons,
HLV plbar, HLV pl,
double lambda
){
HLV pq,pqbar;
const auto backmidquark = std::find_if(
begin(partons)+1, end(partons)-1,
[](Particle const & s){ return s.type != pid::gluon; }
);
assert(backmidquark!=end(partons)-1);
if (is_quark(backmidquark->type)){
pq = to_HepLorentzVector(*backmidquark);
pqbar = to_HepLorentzVector(*(backmidquark+1));
}
else {
pqbar = to_HepLorentzVector(*backmidquark);
pq = to_HepLorentzVector(*(backmidquark+1));
}
auto p1 = to_HepLorentzVector(partons[0]);
auto pn = to_HepLorentzVector(partons[partons.size() - 1]);
auto q0 = pa - p1;
// t-channel momentum after qqx
auto qqxt = q0;
bool wc, wqq;
if (backmidquark->type == -(backmidquark+1)->type){ // Central qqx does not emit
wqq=false;
if (aptype==partons[0].type) {
wc = true;
}
else{
wc = false;
q0-=pl+plbar;
}
}
else{
wqq = true;
wc = false;
qqxt-=pl+plbar;
}
const auto begin_ladder = cbegin(partons) + 1;
const auto end_ladder_1 = (backmidquark);
const auto begin_ladder_2 = (backmidquark+2);
const auto end_ladder = cend(partons) - 1;
for(auto parton_it = begin_ladder; parton_it < begin_ladder_2; ++parton_it){
qqxt -= to_HepLorentzVector(*parton_it);
}
int nabove = std::distance(begin_ladder, backmidquark);
int nbelow = std::distance(begin_ladder_2, end_ladder);
std::vector<HLV> partonsHLV;
partonsHLV.reserve(partons.size());
for (size_t i = 0; i != partons.size(); ++i) {
partonsHLV.push_back(to_HepLorentzVector(partons[i]));
}
const double current_factor = ME_W_qqxmid_current(
aptype, bptype, nabove, nbelow, pa, pb,
pq, pqbar, partonsHLV, plbar, pl, wqq, wc
);
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder_1,
q0, pa, pb, p1, pn,
lambda
)*FKL_ladder_weight(
begin_ladder_2, end_ladder,
qqxt, pa, pb, p1, pn,
lambda
);
return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor;
}
} // namespace anonymous
double MatrixElement::tree_kin_W(Event const & ev) const {
using namespace event_type;
auto const & incoming(ev.incoming());
auto const & decays(ev.decays());
HLV plbar, pl;
for (auto& x: decays) {
if (x.second.at(0).type < 0){
plbar = to_HepLorentzVector(x.second.at(0));
pl = to_HepLorentzVector(x.second.at(1));
}
else{
pl = to_HepLorentzVector(x.second.at(0));
plbar = to_HepLorentzVector(x.second.at(1));
}
}
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
const auto partons = tag_extremal_jet_partons(ev);
if(ev.type() == unordered_backward){
return tree_kin_W_unob(incoming[0].type, incoming[1].type,
pa, pb, partons, plbar, pl,
param_.regulator_lambda);
}
if(ev.type() == unordered_forward){
return tree_kin_W_unof(incoming[0].type, incoming[1].type,
pa, pb, partons, plbar, pl,
param_.regulator_lambda);
}
if(ev.type() == extremal_qqxb){
return tree_kin_W_qqxb(incoming[0].type, incoming[1].type,
pa, pb, partons, plbar, pl,
param_.regulator_lambda);
}
if(ev.type() == extremal_qqxf){
return tree_kin_W_qqxf(incoming[0].type, incoming[1].type,
pa, pb, partons, plbar, pl,
param_.regulator_lambda);
}
if(ev.type() == central_qqx){
return tree_kin_W_qqxmid(incoming[0].type, incoming[1].type,
pa, pb, partons, plbar, pl,
param_.regulator_lambda);
}
return tree_kin_W_FKL(incoming[0].type, incoming[1].type,
pa, pb, partons, plbar, pl,
param_.regulator_lambda);
}
double MatrixElement::tree_kin_Higgs(Event const & ev) const {
if(is_uno(ev.type())){
return tree_kin_Higgs_between(ev);
}
if(ev.outgoing().front().type == pid::Higgs){
return tree_kin_Higgs_first(ev);
}
if(ev.outgoing().back().type == pid::Higgs){
return tree_kin_Higgs_last(ev);
}
return tree_kin_Higgs_between(ev);
}
namespace {
// Colour acceleration multipliers, for gluons see eq. (7) in arXiv:0910.5113
#ifdef HEJ_BUILD_WITH_QCDLOOP
// TODO: code duplication with jets.cc
double K_g(double p1minus, double paminus) {
return 1./2.*(p1minus/paminus + paminus/p1minus)*(C_A - 1./C_A) + 1./C_A;
}
double K_g(
CLHEP::HepLorentzVector const & pout,
CLHEP::HepLorentzVector const & pin
) {
if(pin.z() > 0) return K_g(pout.plus(), pin.plus());
return K_g(pout.minus(), pin.minus());
}
double K(
ParticleID type,
CLHEP::HepLorentzVector const & pout,
CLHEP::HepLorentzVector const & pin
) {
if(type == ParticleID::gluon) return K_g(pout, pin);
return C_F;
}
#endif
// Colour factor in strict MRK limit
double K_MRK(ParticleID type) {
return (type == ParticleID::gluon)?C_A:C_F;
}
}
double MatrixElement::MH2_forwardH(
CLHEP::HepLorentzVector const & p1out,
CLHEP::HepLorentzVector const & p1in,
ParticleID type2,
CLHEP::HepLorentzVector const & p2out,
CLHEP::HepLorentzVector const & p2in,
CLHEP::HepLorentzVector const & pH,
double t1, double t2
) const{
ignore(p2out, p2in);
const double shat = p1in.invariantMass2(p2in);
// gluon case
#ifdef HEJ_BUILD_WITH_QCDLOOP
if(!param_.Higgs_coupling.use_impact_factors){
return K(type2, p2out, p2in)*C_A*1./(16*M_PI*M_PI)*t1/t2*ME_Houtside_gq(
p1out, p1in, p2out, p2in, pH,
param_.Higgs_coupling.mt, param_.Higgs_coupling.include_bottom,
param_.Higgs_coupling.mb
)/(4*(N_C*N_C - 1));
}
#endif
return K_MRK(type2)/C_A*9./2.*shat*shat*(
C2gHgp(p1in,p1out,pH) + C2gHgm(p1in,p1out,pH)
)/(t1*t2);
}
double MatrixElement::tree_kin_Higgs_first(Event const & ev) const {
auto const & incoming = ev.incoming();
auto const & outgoing = ev.outgoing();
assert(outgoing.front().type == pid::Higgs);
if(outgoing[1].type != pid::gluon) {
assert(incoming.front().type == outgoing[1].type);
return tree_kin_Higgs_between(ev);
}
const auto pH = to_HepLorentzVector(outgoing.front());
const auto partons = tag_extremal_jet_partons(
ev
);
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
const auto p1 = to_HepLorentzVector(partons.front());
const auto pn = to_HepLorentzVector(partons.back());
const auto q0 = pa - p1 - pH;
const double t1 = q0.m2();
const double t2 = (pn - pb).m2();
return MH2_forwardH(
p1, pa, incoming[1].type, pn, pb, pH,
t1, t2
)*FKL_ladder_weight(
begin(partons) + 1, end(partons) - 1,
q0, pa, pb, p1, pn,
param_.regulator_lambda
);
}
double MatrixElement::tree_kin_Higgs_last(Event const & ev) const {
auto const & incoming = ev.incoming();
auto const & outgoing = ev.outgoing();
assert(outgoing.back().type == pid::Higgs);
if(outgoing[outgoing.size()-2].type != pid::gluon) {
assert(incoming.back().type == outgoing[outgoing.size()-2].type);
return tree_kin_Higgs_between(ev);
}
const auto pH = to_HepLorentzVector(outgoing.back());
const auto partons = tag_extremal_jet_partons(
ev
);
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
auto p1 = to_HepLorentzVector(partons.front());
const auto pn = to_HepLorentzVector(partons.back());
auto q0 = pa - p1;
const double t1 = q0.m2();
const double t2 = (pn + pH - pb).m2();
return MH2_forwardH(
pn, pb, incoming[0].type, p1, pa, pH,
t2, t1
)*FKL_ladder_weight(
begin(partons) + 1, end(partons) - 1,
q0, pa, pb, p1, pn,
param_.regulator_lambda
);
}
namespace {
template<class InIter, class partIter>
double tree_kin_Higgs_uno(InIter BeginIn, InIter EndIn, partIter BeginPart,
partIter EndPart, const HLV & qH, const HLV & qHp1,
double mt, bool inc_bot, double mb){
const auto pa = to_HepLorentzVector(*BeginIn);
const auto pb = to_HepLorentzVector(*(EndIn-1));
const auto pg = to_HepLorentzVector(*BeginPart);
const auto p1 = to_HepLorentzVector(*(BeginPart+1));
const auto pn = to_HepLorentzVector(*(EndPart-1));
return ME_Higgs_current_uno(
(BeginIn)->type, (EndIn-1)->type, pg, pn, pb, p1, pa,
qH, qHp1, mt, inc_bot, mb
);
}
}
double MatrixElement::tree_kin_Higgs_between(Event const & ev) const {
using namespace event_type;
auto const & incoming = ev.incoming();
auto const & outgoing = ev.outgoing();
const auto the_Higgs = std::find_if(
begin(outgoing), end(outgoing),
[](Particle const & s){ return s.type == pid::Higgs; }
);
assert(the_Higgs != end(outgoing));
const auto pH = to_HepLorentzVector(*the_Higgs);
const auto partons = tag_extremal_jet_partons(ev);
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
auto p1 = to_HepLorentzVector(
partons[(ev.type() == unob)?1:0]
);
auto pn = to_HepLorentzVector(
partons[partons.size() - ((ev.type() == unof)?2:1)]
);
auto first_after_Higgs = begin(partons) + (the_Higgs-begin(outgoing));
assert(
(first_after_Higgs == end(partons) && (
(ev.type() == unob)
|| partons.back().type != pid::gluon
))
|| first_after_Higgs->rapidity() >= the_Higgs->rapidity()
);
assert(
(first_after_Higgs == begin(partons) && (
(ev.type() == unof)
|| partons.front().type != pid::gluon
))
|| (first_after_Higgs-1)->rapidity() <= the_Higgs->rapidity()
);
// always treat the Higgs as if it were in between the extremal FKL partons
if(first_after_Higgs == begin(partons)) ++first_after_Higgs;
else if(first_after_Higgs == end(partons)) --first_after_Higgs;
// t-channel momentum before Higgs
auto qH = pa;
for(auto parton_it = begin(partons); parton_it != first_after_Higgs; ++parton_it){
qH -= to_HepLorentzVector(*parton_it);
}
auto q0 = pa - p1;
auto begin_ladder = begin(partons) + 1;
auto end_ladder = end(partons) - 1;
double current_factor;
if(ev.type() == unob){
current_factor = HEJ::C_A*HEJ::C_A/2*tree_kin_Higgs_uno(
begin(incoming), end(incoming), begin(partons),
end(partons), qH, qH-pH, param_.Higgs_coupling.mt,
param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb
);
const auto p_unob = to_HepLorentzVector(partons.front());
q0 -= p_unob;
p1 += p_unob;
++begin_ladder;
}
else if(ev.type() == unof){
current_factor = HEJ::C_A*HEJ::C_A/2*tree_kin_Higgs_uno(
rbegin(incoming), rend(incoming), rbegin(partons),
rend(partons), qH-pH, qH, param_.Higgs_coupling.mt,
param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb
);
pn += to_HepLorentzVector(partons.back());
--end_ladder;
}
else{
current_factor = ME_Higgs_current(
incoming[0].type, incoming[1].type,
pn, pb, p1, pa, qH, qH - pH,
param_.Higgs_coupling.mt,
param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb
);
}
const double ladder_factor = FKL_ladder_weight(
begin_ladder, first_after_Higgs,
q0, pa, pb, p1, pn,
param_.regulator_lambda
)*FKL_ladder_weight(
first_after_Higgs, end_ladder,
qH - pH, pa, pb, p1, pn,
param_.regulator_lambda
);
return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor;
}
namespace {
double get_AWZH_coupling(Event const & ev, double alpha_s) {
const auto AWZH_boson = std::find_if(
begin(ev.outgoing()), end(ev.outgoing()),
[](auto const & p){return is_AWZH_boson(p);}
);
if(AWZH_boson == end(ev.outgoing())) return 1.;
switch(AWZH_boson->type){
case pid::Higgs:
return alpha_s*alpha_s;
case pid::Wp:
case pid::Wm:
return gw*gw*gw*gw/4.;
// TODO
case pid::photon:
case pid::Z:
default:
throw not_implemented("Emission of boson of unsupported type");
}
}
}
double MatrixElement::tree_param(Event const & ev, double mur) const {
assert(is_resummable(ev.type()));
const auto begin_partons = ev.begin_partons();
const auto end_partons = ev.end_partons();
const auto num_partons = std::distance(begin_partons, end_partons);
const double alpha_s = alpha_s_(mur);
const double gs2 = 4.*M_PI*alpha_s;
double res = std::pow(gs2, num_partons);
if(param_.log_correction){
// use alpha_s(q_perp), evolved to mur
assert(num_partons >= 2);
const auto first_emission = std::next(begin_partons);
const auto last_emission = std::prev(end_partons);
for(auto parton = first_emission; parton != last_emission; ++parton){
res *= 1. + alpha_s/(2.*M_PI)*beta0*log(mur/parton->perp());
}
}
return get_AWZH_coupling(ev, alpha_s)*res;
}
} // namespace HEJ
diff --git a/src/jets.cc b/src/jets.cc
index fb90832..c97ca19 100644
--- a/src/jets.cc
+++ b/src/jets.cc
@@ -1,276 +1,369 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "HEJ/jets.hh"
#include "HEJ/Constants.hh"
// Colour acceleration multiplier for gluons see eq. (7) in arXiv:0910.5113
// @TODO: this is not a current and should be moved somewhere else
double K_g(double p1minus, double paminus) {
return 1./2.*(p1minus/paminus + paminus/p1minus)*(HEJ::C_A - 1./HEJ::C_A) + 1./HEJ::C_A;
}
double K_g(
HLV const & pout,
HLV const & pin
) {
if(pin.z() > 0) return K_g(pout.plus(), pin.plus());
return K_g(pout.minus(), pin.minus());
}
CCurrent CCurrent::operator+(const CCurrent& other)
{
COM result_c0=c0 + other.c0;
COM result_c1=c1 + other.c1;
COM result_c2=c2 + other.c2;
COM result_c3=c3 + other.c3;
return CCurrent(result_c0,result_c1,result_c2,result_c3);
}
CCurrent CCurrent::operator-(const CCurrent& other)
{
COM result_c0=c0 - other.c0;
COM result_c1=c1 - other.c1;
COM result_c2=c2 - other.c2;
COM result_c3=c3 - other.c3;
return CCurrent(result_c0,result_c1,result_c2,result_c3);
}
CCurrent CCurrent::operator*(const double x)
{
COM result_c0=x*CCurrent::c0;
COM result_c1=x*CCurrent::c1;
COM result_c2=x*CCurrent::c2;
COM result_c3=x*CCurrent::c3;
return CCurrent(result_c0,result_c1,result_c2,result_c3);
}
CCurrent CCurrent::operator/(const double x)
{
COM result_c0=CCurrent::c0/x;
COM result_c1=CCurrent::c1/x;
COM result_c2=CCurrent::c2/x;
COM result_c3=CCurrent::c3/x;
return CCurrent(result_c0,result_c1,result_c2,result_c3);
}
CCurrent CCurrent::operator*(const COM x)
{
COM result_c0=x*CCurrent::c0;
COM result_c1=x*CCurrent::c1;
COM result_c2=x*CCurrent::c2;
COM result_c3=x*CCurrent::c3;
return CCurrent(result_c0,result_c1,result_c2,result_c3);
}
CCurrent CCurrent::operator/(const COM x)
{
COM result_c0=(CCurrent::c0)/x;
COM result_c1=(CCurrent::c1)/x;
COM result_c2=(CCurrent::c2)/x;
COM result_c3=(CCurrent::c3)/x;
return CCurrent(result_c0,result_c1,result_c2,result_c3);
}
std::ostream& operator <<(std::ostream& os, const CCurrent& cur)
{
os << "("<<cur.c0<< " ; "<<cur.c1<<" , "<<cur.c2<<" , "<<cur.c3<<")";
return os;
}
CCurrent operator * ( double x, CCurrent& m)
{
return m*x;
}
CCurrent operator * ( COM x, CCurrent& m)
{
return m*x;
}
CCurrent operator / ( double x, CCurrent& m)
{
return m/x;
}
CCurrent operator / ( COM x, CCurrent& m)
{
return m/x;
}
COM CCurrent::dot(HLV p1)
{
// Current goes (E,px,py,pz)
// Vector goes (px,py,pz,E)
return p1[3]*c0-p1[0]*c1-p1[1]*c2-p1[2]*c3;
}
COM CCurrent::dot(CCurrent p1)
{
return p1.c0*c0-p1.c1*c1-p1.c2*c2-p1.c3*c3;
}
//Current Functions
void joi(HLV pout, bool helout, HLV pin, bool helin, current &cur) {
cur[0]=0.;
cur[1]=0.;
cur[2]=0.;
cur[3]=0.;
const double sqpop = sqrt(pout.plus());
const double sqpom = sqrt(pout.minus());
const COM poperp = pout.x() + COM(0, 1) * pout.y();
if (helout != helin) {
throw std::invalid_argument{"Non-matching helicities"};
} else if (helout == false) { // negative helicity
if (pin.plus() > pin.minus()) { // if forward
const double sqpip = sqrt(pin.plus());
cur[0] = sqpop * sqpip;
cur[1] = sqpom * sqpip * poperp / abs(poperp);
cur[2] = -COM(0,1) * cur[1];
cur[3] = cur[0];
} else { // if backward
const double sqpim = sqrt(pin.minus());
cur[0] = -sqpom * sqpim * poperp / abs(poperp);
cur[1] = -sqpim * sqpop;
cur[2] = COM(0,1) * cur[1];
cur[3] = -cur[0];
}
} else { // positive helicity
if (pin.plus() > pin.minus()) { // if forward
const double sqpip = sqrt(pin.plus());
cur[0] = sqpop * sqpip;
cur[1] = sqpom * sqpip * conj(poperp) / abs(poperp);
cur[2] = COM(0,1) * cur[1];
cur[3] = cur[0];
} else { // if backward
const double sqpim = sqrt(pin.minus());
cur[0] = -sqpom * sqpim * conj(poperp) / abs(poperp);
cur[1] = -sqpim * sqpop;
cur[2] = -COM(0,1) * cur[1];
cur[3] = -cur[0];
}
}
}
CCurrent joi (HLV pout, bool helout, HLV pin, bool helin)
{
current cur;
joi(pout, helout, pin, helin, cur);
return CCurrent(cur[0],cur[1],cur[2],cur[3]);
}
void jio(HLV pin, bool helin, HLV pout, bool helout, current &cur) {
joi(pout, !helout, pin, !helin, cur);
}
CCurrent jio (HLV pin, bool helin, HLV pout, bool helout)
{
current cur;
jio(pin, helin, pout, helout, cur);
return CCurrent(cur[0],cur[1],cur[2],cur[3]);
}
void joo(HLV pi, bool heli, HLV pj, bool helj, current &cur) {
// Zero our current
cur[0] = 0.0;
cur[1] = 0.0;
cur[2] = 0.0;
cur[3] = 0.0;
if (heli!=helj) {
throw std::invalid_argument{"Non-matching helicities"};
} else if ( heli == true ) { // If positive helicity swap momenta
std::swap(pi,pj);
}
const double sqpjp = sqrt(pj.plus());
const double sqpjm = sqrt(pj.minus());
const double sqpip = sqrt(pi.plus());
const double sqpim = sqrt(pi.minus());
const COM piperp = pi.x() + COM(0,1) * pi.y();
const COM pjperp = pj.x() + COM(0,1) * pj.y();
const COM phasei = piperp / abs(piperp);
const COM phasej = pjperp / abs(pjperp);
cur[0] = sqpim * sqpjm * phasei * conj(phasej) + sqpip * sqpjp;
cur[1] = sqpim * sqpjp * phasei + sqpip * sqpjm * conj(phasej);
cur[2] = -COM(0, 1) * (sqpim * sqpjp * phasei - sqpip * sqpjm * conj(phasej));
cur[3] = -sqpim * sqpjm * phasei * conj(phasej) + sqpip * sqpjp;
}
CCurrent joo (HLV pi, bool heli, HLV pj, bool helj)
{
current cur;
joo(pi, heli, pj, helj, cur);
return CCurrent(cur[0],cur[1],cur[2],cur[3]);
}
namespace{
//@{
/**
* @brief Pure Jet FKL Contributions, function to handle all incoming types.
* @param p1out Outgoing Particle 1.
* @param p1in Incoming Particle 1.
* @param p2out Outgoing Particle 2
* @param p2in Incoming Particle 2
* @param aqlineb Bool. Is Backwards quark line an anti-quark line?
* @param aqlinef Bool. Is Forwards quark line an anti-quark line?
*
* Calculates j_\mu j^\mu.
* Handles all possible incoming states.
*/
double j_j(HLV p1out, HLV p1in, HLV p2out, HLV p2in, bool aqlineb, bool aqlinef){
HLV q1=p1in-p1out;
HLV q2=-(p2in-p2out);
current mj1m,mj1p,mj2m,mj2p;
// Note need to flip helicities in anti-quark case.
joi(p1out,!aqlineb, p1in,!aqlineb, mj1p);
joi(p1out, aqlineb, p1in, aqlineb, mj1m);
joi(p2out,!aqlinef, p2in,!aqlinef, mj2p);
joi(p2out, aqlinef, p2in, aqlinef, mj2m);
COM Mmp=cdot(mj1m,mj2p);
COM Mmm=cdot(mj1m,mj2m);
COM Mpp=cdot(mj1p,mj2p);
COM Mpm=cdot(mj1p,mj2m);
double sst=abs2(Mmm)+abs2(Mmp)+abs2(Mpp)+abs2(Mpm);
// Multiply by Cf^2
return HEJ::C_F*HEJ::C_F*(sst)/(q1.m2()*q2.m2());
}
} //anonymous namespace
double ME_qQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in){
return j_j(p1out, p1in, p2out, p2in, false, false);
}
double ME_qQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in){
return j_j(p1out, p1in, p2out, p2in, false, true);
}
double ME_qbarQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in){
return j_j(p1out, p1in, p2out, p2in, true, true);
}
double ME_qg(HLV p1out, HLV p1in, HLV p2out, HLV p2in){
return j_j(p1out, p1in, p2out, p2in, false, false)*K_g(p2out, p2in)/HEJ::C_F;
}
double ME_qbarg(HLV p1out, HLV p1in, HLV p2out, HLV p2in){
return j_j(p1out, p1in, p2out, p2in, true, false)*K_g(p2out, p2in)/(HEJ::C_F);
}
double ME_gg(HLV p1out, HLV p1in, HLV p2out, HLV p2in){
return j_j(p1out, p1in, p2out, p2in, false, false)*K_g(p1out, p1in)*K_g(p2out, p2in)/(HEJ::C_F*HEJ::C_F);
}
//@}
+
+namespace{
+ double juno_j(HLV pg, HLV p1out, HLV p1in, HLV p2out,
+ HLV p2in, bool aqlinea, bool aqlineb){
+ // This construction is taking rapidity order: pg > p1out >> p2out
+ HLV q1=p1in-p1out; // Top End
+ HLV q2=-(p2in-p2out); // Bottom End
+ HLV qg=p1in-p1out-pg; // Extra bit post-gluon
+
+ // Note <p1|eps|pa> current split into two by gauge choice.
+ // See James C's Thesis (p72). <p1|eps|pa> -> <p1|pg><pg|pa>
+ CCurrent mj1p=joi(p1out,!aqlinea,p1in,!aqlinea);
+ CCurrent mj1m=joi(p1out,aqlinea,p1in,aqlinea);
+ CCurrent jgap=joi(pg,!aqlinea,p1in,!aqlinea);
+ CCurrent jgam=joi(pg,aqlinea,p1in,aqlinea);
+
+ // Note for function joo(): <p1+|pg+> = <pg-|p1->.
+ CCurrent j2gp=joo(p1out,!aqlinea,pg,!aqlinea);
+ CCurrent j2gm=joo(p1out,aqlinea,pg,aqlinea);
+
+ CCurrent mj2p=joi(p2out,!aqlineb,p2in,!aqlineb);
+ CCurrent mj2m=joi(p2out,aqlineb,p2in,aqlineb);
+
+ // Dot products of these which occur again and again
+ COM Mmp=mj1m.dot(mj2p);
+ COM Mmm=mj1m.dot(mj2m);
+ COM Mpp=mj1p.dot(mj2p);
+ COM Mpm=mj1p.dot(mj2m);
+
+ CCurrent p1o(p1out),p2o(p2out),p2i(p2in),qsum(q1+qg),p1i(p1in);
+
+ CCurrent Lmm=(qsum*(Mmm)+(-2.*mj2m.dot(pg))*mj1m+2.*mj1m.dot(pg)*mj2m
+ +(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*Mmm/2.))/q1.m2();
+ CCurrent Lmp=(qsum*(Mmp) + (-2.*mj2p.dot(pg))*mj1m+2.*mj1m.dot(pg)*mj2p
+ +(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*Mmp/2.))/q1.m2();
+ CCurrent Lpm=(qsum*(Mpm) + (-2.*mj2m.dot(pg))*mj1p+2.*mj1p.dot(pg)*mj2m
+ +(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*Mpm/2.))/q1.m2();
+ CCurrent Lpp=(qsum*(Mpp) + (-2.*mj2p.dot(pg))*mj1p+2.*mj1p.dot(pg)*mj2p
+ +(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*Mpp/2.))/q1.m2();
+
+ CCurrent U1mm=(jgam.dot(mj2m)*j2gm+2.*p1o*Mmm)/(p1out+pg).m2();
+ CCurrent U1mp=(jgam.dot(mj2p)*j2gm+2.*p1o*Mmp)/(p1out+pg).m2();
+ CCurrent U1pm=(jgap.dot(mj2m)*j2gp+2.*p1o*Mpm)/(p1out+pg).m2();
+ CCurrent U1pp=(jgap.dot(mj2p)*j2gp+2.*p1o*Mpp)/(p1out+pg).m2();
+ CCurrent U2mm=((-1.)*j2gm.dot(mj2m)*jgam+2.*p1i*Mmm)/(p1in-pg).m2();
+ CCurrent U2mp=((-1.)*j2gm.dot(mj2p)*jgam+2.*p1i*Mmp)/(p1in-pg).m2();
+ CCurrent U2pm=((-1.)*j2gp.dot(mj2m)*jgap+2.*p1i*Mpm)/(p1in-pg).m2();
+ CCurrent U2pp=((-1.)*j2gp.dot(mj2p)*jgap+2.*p1i*Mpp)/(p1in-pg).m2();
+
+ constexpr double cf=HEJ::C_F;
+
+ double amm=cf*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*cf*cf/3.*vabs2(U1mm+U2mm);
+ double amp=cf*(2.*vre(Lmp-U1mp,Lmp+U2mp))+2.*cf*cf/3.*vabs2(U1mp+U2mp);
+ double apm=cf*(2.*vre(Lpm-U1pm,Lpm+U2pm))+2.*cf*cf/3.*vabs2(U1pm+U2pm);
+ double app=cf*(2.*vre(Lpp-U1pp,Lpp+U2pp))+2.*cf*cf/3.*vabs2(U1pp+U2pp);
+ double ampsq=-(amm+amp+apm+app);
+
+ //Divide by t-channels
+ ampsq/=q2.m2()*qg.m2();
+ ampsq/=16.;
+
+ // Factor of (Cf/Ca) for each quark to match j_j.
+ ampsq*=(HEJ::C_F*HEJ::C_F)/(HEJ::C_A*HEJ::C_A);
+
+ return ampsq;
+
+ }
+}
+
+//Unordered bits for pure jet
+double ME_unob_qQ (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in){
+ return juno_j(pg, p1out, p1in, p2out, p2in, false, false);
+}
+
+double ME_unob_qbarQ (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in){
+ return juno_j(pg, p1out, p1in, p2out, p2in, true, false);
+}
+
+double ME_unob_qQbar (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in){
+ return juno_j(pg, p1out, p1in, p2out, p2in, false, true);
+}
+
+double ME_unob_qbarQbar (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in){
+ return juno_j(pg, p1out, p1in, p2out, p2in, true, true);
+}
+
+double ME_unob_qg (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in){
+ return juno_j(pg, p1out, p1in, p2out, p2in, false, false)*K_g(p2out,p2in)/HEJ::C_F;
+}
+
+double ME_unob_qbarg (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in){
+ return juno_j(pg, p1out, p1in, p2out, p2in, true, false)*K_g(p2out,p2in)/HEJ::C_F;
+}
diff --git a/t/ME_data/ME_jets_tree.dat b/t/ME_data/ME_jets_tree.dat
index fc2be65..0f36db3 100644
--- a/t/ME_data/ME_jets_tree.dat
+++ b/t/ME_data/ME_jets_tree.dat
@@ -1,672 +1,742 @@
1089423.187
16.03009374
6634.608547
75414606.57
77554.50738
197575.0454
76983.13177
38.35809694
14.34964854
2.299306716e+10
769524899.8
20.89671172
76.68018946
8308.664153
2450736.574
52.24359459
497.7665557
5051.865522
94571.50832
1362099715
3690.443315
4732694.851
33060865.81
38851015.2
30289.18618
283948.732
14123.1659
40.56212546
6.779515122e+12
25724.24084
4.464903667
26.15970999
28358522.83
481820.3073
33.14640689
91.72649968
322690.8843
39.04279288
65541087.6
2256.832837
51.04918081
21658.09425
675.8508089
0.003214397574
413076.4025
1194042.703
15261824.18
1.190018997e+13
25852091.94
765101.1953
579.2342295
23845182.04
92.16369985
705.5684735
828764.9832
81.971147
1436684.808
70.01623522
230.5658106
8794.725516
890.1697598
97617236.05
31101.11974
10351668.97
187539232.1
93875.71518
43031260.49
4.067125353
6035.796288
16569782.71
46.2510317
52.88862266
0.5556956436
606771384.7
293585.0788
8415963.902
29.5069149
18.96420915
6.321441102e+14
3.592998225
129.9175932
19.85723932
20438567.28
2.633206648e+10
2051545.14
13.98077027
342.1399233
17.65245465
79250.23744
1030.69768
73407.42532
237883.1282
79.44231257
81668.43497
93415340.81
430.286772
2209.749603
568258.2296
173145.2499
1620285.613
337.6920671
153.562771
416.0616338
2187.267876
22.41576042
49.58992905
559801.7855
45.11194021
584508.9206
35.3041943
202.7925831
15.00277101
46.15144698
4.588470323e+13
597.3210575
39.62999028
5228016720
822065.5972
1.826581842e+10
6.479966482
486208349.5
62.22613221
21.46207828
792.2244925
98.46049107
192.195726
539903191.6
7888.570415
559946441.6
99.19582481
967.8349215
1186287.112
296.9970515
393.8848181
46.84865591
27.79939609
70102.51858
68.89626885
2702.415147
257.708207
155288432.1
1225359.249
90576.73064
9434.676956
88.90057162
4378317.363
15680420.87
6567.795332
94.80131393
44.21436446
61.79707787
44.41095665
7385.16094
227.1571585
94986.49609
67239897.71
96.37463492
2911.906879
49.889402
6585.792276
336579403.2
564.3980943
1554303.479
33573.8975
2.437369584e+10
226746.4443
75.83290452
27.75735686
94778193.37
59.97102623
64.15185493
7.774417686e+13
7438.779369
22.08782146
434613.7211
34293520.02
68775.95004
307606184.4
150615.387
899713.4452
9.223973371e+11
59817.02462
348.4932635
2609243018
24597079.8
98.70527236
83.33745442
12560815.61
1358729.011
6976819.18
7425.681281
488.6665953
13.67398052
401.1059245
1888.653087
4009.955386
175.1044182
864180.0454
37.7515354
47187493.74
109.4430381
359.7148155
2.340917207e+11
179.4280073
9827.615509
51.43706602
159.8742767
2066998.073
21197190.52
154.5419162
20.22828093
96.58078786
54.45061659
65.30848921
27.46863884
399638.3233
48.71138495
1579072.132
67717150.9
1034311422
2300.166812
542.347043
225.023696
15.39245502
5.045806172
408439.9118
2850069.206
412.6441358
461.2402375
22.40894016
37199025.96
16799601.72
1283263.103
19.12993911
4119931.525
138990.0255
122.7744597
18733.45006
56.21459033
343399063.9
40.69967214
1257.945403
2.549605794e+11
100938.8483
197079455.9
18898.20284
65974.22935
2007851.613
6824.687033
4076901842
7253481.757
3.002058086
1662154.497
3.070682383e+12
163.8653423
560506230.7
1510.132538
26693.22493
0.1285518109
494554.148
24594487.35
271745.3087
155.6774454
0.01468054527
3.34264067
1588023.004
11.02440854
2.276304783
2060679.506
1.060994634
7899476.433
11392931.55
0.002283073032
223.8249674
64.14820176
222013139.2
17285137.17
4904.948169
4.325556105e+10
2.025060646
64.72904144
0.7119503103
0.3108282312
127.8083553
231480500.9
3004499.276
35310522.85
9501.295688
34.48365857
0.1237148274
137.304491
22906.44193
2845.690168
1425.623527
8.059398488
2836.653576
545.9404471
6676238.836
453869217.1
3046455.67
3824.741211
1733.364341
163997517.5
46856740.63
18247.63995
55504734.26
1281.154525
0.01317003733
26796232.71
614.7809373
124.9917797
13.32414443
2.061438665
1437324478
7.653466138
215267780.4
304925.7977
12340855.35
1.292474046
4012100939
918.1699992
0.001829689464
266.1514907
19586.50145
9843.077092
34563.10338
2.665695419e+14
194168.4335
182.9208313
1.789617285
115011741.8
1.600475143
0.9039184064
54.16404691
3836249.251
0.1453799136
83576.30659
75918.72599
255875.4377
0.0940918961
8139.618634
+41.2113043
+41.2113043
+18.2434106
+18.2434106
+18.2434106
+18.2434106
497813992.4
9158737.332
81.5174511
+0.1756986828
+0.1756986828
+0.07788028858
+0.07788028858
+0.07788028858
+0.07788028858
85205.83347
1827.041926
5434428.481
34601.43212
645859.4651
8.409921219
+0.00878465753
+0.00878465753
+0.00878465753
+0.0200891802
+0.0200891802
102240.6791
0.3455555689
+5.186512555e-05
+5.186512555e-05
+2.095582296e-05
+2.095582296e-05
+2.095582296e-05
+2.095582296e-05
2.273623498e+11
2139040773
82682.82602
+1.697773045
+1.697773045
+1.697773045
+3.928017028
+3.928017028
164684299.3
6840.694229
6292.264217
1.266028837
23696.0817
9776772.75
146.8084842
40186.198
578656.2198
0.05760079494
1308.62247
2846240637
2269.666603
0.05728404415
97.64433608
816.0835471
+2.245768012
+2.245768012
+2.245768012
+5.059168555
+5.059168555
0.003537766648
715128749.1
581519.0926
28878.33931
606.7268882
8.540999968
1.843057112
14728.95688
0.09287199376
+0.0001302392694
+0.0001302392694
+5.236903901e-05
+5.236903901e-05
+5.236903901e-05
4.952792168
+0.03410256159
+0.03410256159
+0.01401286426
+0.01401286426
+0.01401286426
39703045.24
0.2253239938
0.4007086942
0.6679356978
8985.245646
1459.149121
222.574255
1410.034894
49.66459517
280057.7217
11990414.03
22.39502507
6652.821533
412.6309921
3976669.212
376.5293878
0.04993631096
325.2221445
0.001056663136
0.6094009495
0.01755460768
4169089.697
120922.3418
19.94386406
1.034905758
899064.7273
13370.52361
12.9511581
7.43997218
0.2195708041
22.69922572
1364456.758
37.31335252
193721994.5
5387171012
83382.68593
4356.219769
4.959390234
280348.136
4068.591503
497.1128967
1220.499616
37070808.46
5.029130002
13411.495
21606019.92
13.50033011
10356.03697
0.853804882
0.004144159842
66.5696697
215.1472766
1.80856617
47652781.08
187.4720271
0.02309225892
221825.3003
68502271.85
89524.33999
63.06068401
35193.88681
0.1324344626
0.07415048828
1.316536312
1914484.13
239.8833788
404.8578233
311149.7413
31.91203802
6765.191464
3683.120991
111.8561994
622181.6379
3.495906967
1.324392325
2932853.167
1221.022853
0.3964612919
772.5117573
1587.443569
16407.82538
1153370.819
697976.9971
5293.102236
0.6273574623
433333.6561
667.0010151
279.5253754
13.71499773
27581.98765
0.1265026607
0.1265655615
0.217545965
119235.7842
1207.825866
1.168584528
442927.4532
11025.74448
12920.20556
6.440028505
0.01867492246
0.009209754218
3.078810838
1.895501574
4354.355041
1376.801012
281.0956164
12.62981914
0.01896307384
1.469952339
1070.389041
0.3140297342
17576.52171
19.46940005
72639.28176
1.00064982
10.58896856
146.5709549
240963.5215
21441702.88
9436.867787
1807651761
253.6931177
33906.26088
73290.29485
58.88693751
0.1874080618
128856.1698
37.9170344
999.730091
79.29102322
4297075.439
2815203002
1.065757865
41.63356752
3132937.645
0.006570811459
114610718.1
19766422.36
104.7817608
978.6047854
1114.98705
2962094.183
4303175.509
40918.11534
6.458644406
0.01513394186
10888.55623
1429.371661
72470.14267
1239031.521
0.01533517898
9486.889743
9312.872953
0.001369262882
263.7126988
27397.40987
6.460887231e+11
20.47552795
11886.05681
359511.4015
16322071.14
5.495229
94.06819575
11159.01276
497280.7644
152.8232363
2819.864026
1840.922093
105.3416519
1.51196546
893.0699171
274.2733387
213.9840929
29.84435445
0.06575134529
68473.33981
0.04814411438
95307.18205
0.2103887883
87533.03201
170.1269759
12266.74779
4.906472106
1289.949588
206.2206679
1339736.494
3265180.524
817632.5321
777148.2976
45580.38902
112.3471097
40.46375548
113.0563114
0.5022425417
125.6469752
1250753.214
119.6687695
7890.066678
255592.6823
0.09544112407
80.43028795
168.8638851
231.3586602
9.003352253
3865310.035
0.0327080534
0.2259241849
0.3309640842
70628514.93
1.342293964
1558415395
1.535156072
1105643165
11.12340492
19553.68065
0.03262635444
43.60758633
149.96132
16.71798701
862303.2755
3779.707836
227.1354158
0.04181934127
0.007441956106
0.06547652512
562.1841603
0.1962099745
47427.77928
278196156.8
3272270.994
3.274695664
3.671659367
0.8773587215
3475.635402
150.1691185
79249376.57
100282855
1034000.854
34631.08123
0.2021471708
11338.25051
365.9254869
3614.640649
38892.33006
3206.818315
495.3831508
3131.024396
1216.930767
3.68554273
0.0007587392874
27267.82366
66215.88611
12270.6996
2181.519341
9.799873017
84.23360373
350624.9294
2324.233559
0.210264797
1438761.868
160.6723109
39.13994322
2317.668587
32196.71663
0.02293846799
177.7024259
294490.6299
963.9483044
683.7845476
12.22284439
4350.875864
132.6662523
800.0540763
0.6752102469
6384.094581
667.7145567
88.5979278
65957.42139
0.1009745291
7074.842515
512.2889019
11.48214352
2.877973951
1.31777398
7.553015302
0.01771209643
0.01420982393
+3.754178119e-07
+1.116826573e-06
+1.116826573e-06
+3.754178119e-07
+3.754178119e-07
+3.754178119e-07
+1.769703186e-06
+4.356179554e-06
+4.356179554e-06
+1.769703186e-06
+1.769703186e-06
+1.769703186e-06
5.350111784
5171856.193
0.6877419322
+1.22182262e-05
+1.22182262e-05
+1.22182262e-05
+2.818920729e-05
+2.818920729e-05
+2.708089303e-05
+2.708089303e-05
+2.708089303e-05
+6.339285468e-05
+6.339285468e-05
+222.6777557
+222.6777557
+222.6777557
+501.1839533
+501.1839533
diff --git a/t/ME_data/PSP_jets.lhe.gz b/t/ME_data/PSP_jets.lhe.gz
index b471729..3ac0130 100644
--- a/t/ME_data/PSP_jets.lhe.gz
+++ b/t/ME_data/PSP_jets.lhe.gz
@@ -1,3 +1,3 @@
version https://git-lfs.github.com/spec/v1
-oid sha256:98254005df93663441f19013925d3dd40575ffcbfd666387a05cfd6d1b64c489
-size 161645
+oid sha256:c8a0adcc7453343934c235a295515b2838b371f8aef4d52b6ffe61301774efe6
+size 163612