diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 7ac6ae8..e0315db 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -1,453 +1,471 @@ # --------------------------------- # - General Setup - # --------------------------------- stages: - build - test - FOG:build - FOG:test - clean_code - FOG:clear_code - long_test - publish workflow: rules: - when: always variables: GIT_SUBMODULE_STRATEGY: recursive # build directories HEJ_BUILD_DIR: tmp_HEJ/HEJ_build HEJ_INSTALL_DIR: tmp_HEJ/HEJ_installed FOG_BUILD_DIR: tmp_HEJ/FOG_build FOG_INSTALL_DIR: ${HEJ_INSTALL_DIR} # docker images DOCKER_BASIC: hejdock/hepenv # slc6 DOCKER_CLANG: 'hejdock/hepenv:clang9_ubuntu' # ubuntu with clang DOCKER_HEPMC3: hejdock/hepmc3env # centos7 DOCKER_HIGHFIVE: hejdock/highfiveenv # centos7 DOCKER_QCDLOOP: hejdock/qcdloopenv # slc6, includes rivet DOCKER_RIVET3: 'hejdock/rivetenv:rivet3_gcc9_ubuntu' # ubuntu DOCKER_RIVET3H2: 'hejdock/rivetenv:rivet3_hepmc2_gcc9_ubuntu' # ubuntu DOCKER_RIVET: hejdock/rivetenv # slc6 # default name of cmake CMAKE: cmake CTEST: ctest OVERWRITE_BOOST: "" # ----------- Macros ----------- after_script: - date .tags_template: tags: &tags_def - docker # default pipeline triggers .rules_template: rules: &rules_def - if: $CI_MERGE_REQUEST_ID when: on_success - if: $CI_COMMIT_TAG when: on_success - if: '$CI_COMMIT_BRANCH != null && $CI_COMMIT_BRANCH == "master"' when: on_success - if: '$CI_COMMIT_BRANCH != null && $CI_COMMIT_BRANCH =~ /^v\d+\.\d+$/' when: on_success - when: manual # non-blocking manual allow_failure: true # save complete history of failed tests .save_failure: artifacts: &artifacts_failed when: on_failure untracked: true expire_in: 3d # --------------------------------- # - Script Templates - # --------------------------------- # ----------- Build ----------- .HEJ_build: tags: *tags_def rules: *rules_def stage: build before_script: - date - source /cvmfs/pheno.egi.eu/HEJ/HEJ_env.sh || exit 1 # prepare build - t_HEJ_DIR=${PWD} - t_HEJ_INSTALL_DIR=${PWD}/${HEJ_INSTALL_DIR} - t_HEJ_BUILD_DIR=${PWD}/${HEJ_BUILD_DIR} # hack to allow setting "variable in variable" - t_CMAKE_FLAGS="" - if [[ -n ${OVERWRITE_BOOST} ]]; then declare -n t_boost_root=${OVERWRITE_BOOST}; t_CMAKE_FLAGS="-DBOOST_ROOT=${t_boost_root}"; fi - echo ${t_CMAKE_FLAGS} - mkdir -p ${t_HEJ_BUILD_DIR} - cd ${t_HEJ_BUILD_DIR} - ${CMAKE} ${t_HEJ_DIR} -DCMAKE_BUILD_TYPE=DEBUG -DCMAKE_INSTALL_PREFIX=${t_HEJ_INSTALL_DIR} ${t_CMAKE_FLAGS} script: - make -j $(nproc --ignore=1) - make install needs: [] # can run immediately artifacts: # save build and installed folder name: build expire_in: 1d paths: - ${HEJ_INSTALL_DIR} - ${HEJ_BUILD_DIR} # ----------- Test ----------- .HEJ_test: tags: *tags_def rules: *rules_def stage: test before_script: - date - source /cvmfs/pheno.egi.eu/HEJ/HEJ_env.sh || exit 1 # load HEJ - t_HEJ_DIR=${PWD} - t_HEJ_INSTALL_DIR=${PWD}/${HEJ_INSTALL_DIR} - export LD_LIBRARY_PATH=${t_HEJ_INSTALL_DIR}/lib:${LD_LIBRARY_PATH} - export PATH=${t_HEJ_INSTALL_DIR}/bin:${PATH} - cd ${HEJ_BUILD_DIR} - ${CMAKE} ${t_HEJ_DIR} # rerun cmake to create all test configure files script: - ${CTEST} --output-on-failure -j $(nproc --ignore=1) artifacts: *artifacts_failed .HEJrivet_test: extends: .HEJ_test script: - ${CTEST} --output-on-failure -j $(nproc --ignore=1) - bash -c '[ -f t/tst.yoda ]' && echo "found rivet output" - rivet-cmphistos t/tst.yoda - bash -c '[ -f MC_XS_XS.dat ]' && echo "yoda not empty" ## ----------- FOG build ----------- .FOG_build: tags: *tags_def rules: *rules_def stage: FOG:build before_script: - date - source /cvmfs/pheno.egi.eu/HEJ/HEJ_env.sh || exit 1 # load HEJ - t_HEJ_INSTALL_DIR=${PWD}/${HEJ_INSTALL_DIR} - export LD_LIBRARY_PATH=${t_HEJ_INSTALL_DIR}/lib:${LD_LIBRARY_PATH} - export PATH=${t_HEJ_INSTALL_DIR}/bin:${PATH} # prepare build - t_FOG_DIR=${PWD}/FixedOrderGen - t_FOG_INSTALL_DIR=${PWD}/${FOG_INSTALL_DIR} - t_FOG_BUILD_DIR=${PWD}/${FOG_BUILD_DIR} # hack to allow setting "variable in variable" - t_CMAKE_FLAGS="" - if [[ -n ${OVERWRITE_BOOST} ]]; then declare -n t_boost_root=${OVERWRITE_BOOST}; t_CMAKE_FLAGS="-DBOOST_ROOT=${t_boost_root}"; fi - mkdir -p ${t_FOG_BUILD_DIR} - cd ${t_FOG_BUILD_DIR} - ${CMAKE} ${t_FOG_DIR} -DCMAKE_BUILD_TYPE=DEBUG -DCMAKE_INSTALL_PREFIX=${t_FOG_INSTALL_DIR} ${t_CMAKE_FLAGS} script: - make -j $(nproc --ignore=1) - make install artifacts: # save build and installed folder name: FOG_build expire_in: 1d paths: - ${HEJ_INSTALL_DIR} - ${FOG_INSTALL_DIR} - ${FOG_BUILD_DIR} ## ----------- FOG test ----------- .FOG_test: tags: *tags_def rules: *rules_def stage: FOG:test before_script: - date - source /cvmfs/pheno.egi.eu/HEJ/HEJ_env.sh || exit 1 # load HEJ - t_FOG_DIR=${PWD}/FixedOrderGen - t_HEJ_INSTALL_DIR=${PWD}/${HEJ_INSTALL_DIR} - t_FOG_INSTALL_DIR=${PWD}/${FOG_INSTALL_DIR} - export LD_LIBRARY_PATH=${t_HEJ_INSTALL_DIR}/lib:${LD_LIBRARY_PATH} - export PATH=${t_HEJ_INSTALL_DIR}/bin:${t_FOG_INSTALL_DIR}/bin:${PATH} - t_FOG_BUILD_DIR=${PWD}/${FOG_BUILD_DIR} - cd ${t_FOG_BUILD_DIR} - ${CMAKE} ${t_FOG_DIR} # rerun cmake to create all test configure files script: - ${CTEST} --output-on-failure -j $(nproc --ignore=1) artifacts: *artifacts_failed # --------------------------------- # - Build & Test - # --------------------------------- # ----------- basic (always run) ----------- build:basic: image: ${DOCKER_BASIC} extends: .HEJ_build rules: - when: on_success test:basic: image: ${DOCKER_BASIC} extends: .HEJ_test needs: - job: build:basic artifacts: true rules: - when: on_success FOG:build:basic: image: ${DOCKER_BASIC} extends: .FOG_build needs: - job: build:basic artifacts: true rules: - when: on_success FOG:test:basic: image: ${DOCKER_BASIC} extends: .FOG_test needs: - job: FOG:build:basic artifacts: true rules: - when: on_success # ----------- HighFive/hdf5 ----------- build:HighFive: image: ${DOCKER_HIGHFIVE} extends: .HEJ_build variables: CMAKE: cmake3 test:HighFive: image: ${DOCKER_HIGHFIVE} extends: .HEJ_test variables: CMAKE: cmake3 CTEST: ctest3 needs: - job: build:HighFive artifacts: true # ----------- QCDloop ----------- build:qcdloop: image: ${DOCKER_QCDLOOP} extends: .HEJ_build test:qcdloop: image: ${DOCKER_QCDLOOP} extends: .HEJ_test needs: - job: build:qcdloop artifacts: true # ----------- rivet3 & HepMC3 ----------- build:rivet3: image: ${DOCKER_RIVET3} extends: .HEJ_build variables: OVERWRITE_BOOST: boost_ROOT_DIR test:rivet3: image: ${DOCKER_RIVET3} extends: .HEJrivet_test needs: - job: build:rivet3 artifacts: true # ----------- HepMC 3 ----------- # HepMC 3 already in rivet3 build:hepmc3: image: ${DOCKER_HEPMC3} extends: .HEJ_build variables: CMAKE: cmake3 test:hepmc3: image: ${DOCKER_HEPMC3} extends: .HEJ_test variables: CMAKE: cmake3 CTEST: ctest3 needs: - job: build:hepmc3 artifacts: true # ----------- rivet ----------- # rivet2 already in qcdloop build:rivet: image: ${DOCKER_RIVET} extends: .HEJ_build test:rivet: image: ${DOCKER_RIVET} extends: .HEJrivet_test needs: - job: build:rivet artifacts: true # ----------- rivet3 & HepMC2 ----------- # This shouldn't change too often build:rivet3h2: image: ${DOCKER_RIVET3H2} extends: .HEJ_build variables: OVERWRITE_BOOST: boost_ROOT_DIR test:rivet3h2: image: ${DOCKER_RIVET3H2} extends: .HEJrivet_test needs: - job: build:rivet3h2 artifacts: true # --------------------------------- # - Clean Code - # --------------------------------- No_tabs: stage: clean_code tags: *tags_def rules: *rules_def image: hejdock/git needs: [] # can run immediately script: - date - check_tabs artifacts: # don't save anything # ----------- Template no compiler warnings ----------- .Warning_build: extends: .HEJ_build stage: clean_code script: - ${CMAKE} ${t_HEJ_DIR} -DCMAKE_CXX_FLAGS="-Werror" -DCMAKE_BUILD_TYPE=RelWithDebInfo - make -j $(nproc --ignore=1) - make -j $(nproc --ignore=1) install artifacts: # don't save anything needs: [] # can run immediately .Warning_FOG: extends: .FOG_build stage: FOG:clear_code script: - ${CMAKE} ${t_FOG_DIR} -DCMAKE_CXX_FLAGS="-Werror" ${t_CMAKE_FLAGS} -DCMAKE_BUILD_TYPE=RelWithDebInfo - make -j $(nproc --ignore=1) - make -j $(nproc --ignore=1) install artifacts: # don't save anything # ----------- No gcc warnings ----------- No_Warning:basic: image: ${DOCKER_BASIC} extends: .Warning_build No_Warning:basic:FOG: image: ${DOCKER_BASIC} extends: .Warning_FOG needs: - job: build:basic artifacts: true # ----------- No clang warnings ----------- No_Warning:clang: image: ${DOCKER_CLANG} extends: .Warning_build variables: OVERWRITE_BOOST: boost_ROOT_DIR artifacts: # save installed folder name: build expire_in: 1d paths: - ${HEJ_INSTALL_DIR} No_Warning:clang:FOG: image: ${DOCKER_CLANG} extends: .Warning_FOG needs: - job: No_Warning:clang artifacts: true variables: OVERWRITE_BOOST: boost_ROOT_DIR # --------------------------------- # - Long tests - # --------------------------------- .HEJ_test_long: extends: .HEJ_build stage: long_test script: - ${CMAKE} ${t_HEJ_DIR} -DCMAKE_BUILD_TYPE=RelWithDebInfo ${t_CMAKE_FLAGS} -DTEST_ALL=TRUE - make -j $(nproc --ignore=1) install - ${CTEST} --output-on-failure -j $(nproc --ignore=1) needs: [] # can run immediately +.FOG_test_long: + extends: .FOG_build + stage: long_test + script: + - ${CMAKE} ${t_HEJ_DIR} -DCMAKE_BUILD_TYPE=RelWithDebInfo ${t_CMAKE_FLAGS} + -DTEST_ALL=TRUE + - make -j $(nproc --ignore=1) install + - ${CTEST} --output-on-failure -j $(nproc --ignore=1) + # ----------- QCDloop ----------- Long_test:qcdloop: image: ${DOCKER_QCDLOOP} extends: .HEJ_test_long +# ----------- Basic ----------- + +Long_test:FOG: + image: ${DOCKER_BASIC} + extends: .HEJ_test_long + needs: + - job: No_Warning:basic + artifacts: true + # --------------------------------- # - Publish - # --------------------------------- Publish_version: stage: publish image: hejdock/git dependencies: [] artifacts: # don't save anything rules: - if: '$CI_COMMIT_BRANCH != null && $CI_COMMIT_BRANCH =~ /^v\d+\.\d+$/' when: on_success before_script: - mkdir -p .ssh/ && echo "${SSH_KEY}" > .ssh/id_rsa && chmod 0600 .ssh/id_rsa - rm -rf ~/.ssh/id_rsa; mkdir -p ~/.ssh/ - ln -s $PWD/.ssh/id_rsa ~/.ssh/id_rsa && chmod 0600 ~/.ssh/id_rsa - ssh -T ${PUBLIC_GIT_PREFIX} -o "StrictHostKeyChecking no" || echo "added ssh" script: - git remote add public ${PUBLIC_GIT_PREFIX}${PUBLIC_GIT_POSTFIX} - git checkout $CI_COMMIT_REF_NAME - git branch - git pull - git push public - git push public --tags after_script: - rm -f /root/.ssh/id_rsa && rm -fr .ssh - git remote rm public diff --git a/CMakeLists.txt b/CMakeLists.txt index 73f6dc2..f91fc7c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,305 +1,305 @@ 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.6 LANGUAGES C CXX) ## User settings include(CMakeDependentOption) option(EXCLUDE_QCDloop "Do not include QCDloop" FALSE) option(EXCLUDE_HighFive "Do not include HighFive" FALSE) option(EXCLUDE_HepMC "Do not include HepMC version 2" FALSE) option(EXCLUDE_HepMC3 "Do not include HepMC version 3" FALSE) option(EXCLUDE_ROOT "Do not include ROOT" TRUE) # Only skip rivet if we both HepMC 2 & 3 are excluded, since we don't know which # HepMC version rivet needs CMAKE_DEPENDENT_OPTION(EXCLUDE_rivet "Do not include rivet" FALSE "NOT EXCLUDE_HepMC OR NOT EXCLUDE_HepMC3" TRUE) option(BUILD_EXAMPLES "Build examples" FALSE) option(TEST_ALL "Run additional (longer) tests" FALSE) option(TEST_COVERAGE "Generate test coverage with \"gcovr\"" FALSE) option(RUN_CLANG_TIDY "Run clang tidy" FALSE) # 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() ## Global flags for the compiler (all warnings) if (CMAKE_COMPILER_IS_GNUCC OR CMAKE_COMPILER_IS_GNUCXX) set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wextra -Wpedantic") elseif (MSVC) set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /W4 /WX /EHsc") endif() set(CMAKE_CXX_STANDARD_REQUIRED ON) set(CMAKE_CXX_STANDARD 14) ## Check if this is a git folder find_package(Git QUIET) if(GIT_FOUND AND EXISTS "${PROJECT_SOURCE_DIR}/.git") set(_is_git TRUE) endif() ## download FORM submodule # from https://cliutils.gitlab.io/modern-cmake/chapters/projects/submodule.html if(_is_git) option(GIT_SUBMODULE "Check submodules during build" ON) if(GIT_SUBMODULE AND NOT EXISTS "${PROJECT_SOURCE_DIR}/current_generator/form/Makefile.am") message(STATUS "Submodule update") execute_process(COMMAND ${GIT_EXECUTABLE} submodule update --init WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR} RESULT_VARIABLE GIT_SUBMOD_RESULT) if(NOT GIT_SUBMOD_RESULT EQUAL "0") message(FATAL_ERROR "git submodule update --init failed with ${GIT_SUBMOD_RESULT}, please checkout submodules") endif() endif() endif() if(NOT EXISTS "${PROJECT_SOURCE_DIR}/current_generator/form/Makefile.am") message(FATAL_ERROR "FORM was not downloaded! Please update git-submodules and try again.") endif() ## Create Version set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}/cmake/Modules/") # Get the latest abbreviated commit hash of the working branch if(_is_git) execute_process( COMMAND ${GIT_EXECUTABLE} 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_EXECUTABLE} rev-parse --abbrev-ref HEAD WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR} OUTPUT_VARIABLE PROJECT_GIT_BRANCH OUTPUT_STRIP_TRAILING_WHITESPACE ) else() set(PROJECT_GIT_REVISION "UNKNOWN") set(PROJECT_GIT_BRANCH "UNKNOWN") endif() ## 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) ## find dependencies find_package(CLHEP 2.3 REQUIRED) find_package(fastjet REQUIRED) find_package(LHAPDF 6 REQUIRED) ## 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) find_package(yaml-cpp) # requiring yaml does not work with fedora if(EXCLUDE_HepMC) message(STATUS "Skipping HepMC") else() find_package(HepMC 2 EXACT) 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 QUIET) # suppress CMake warning if HepMC3 not available ## print standard message find_package_handle_standard_args( HepMC3 FOUND_VAR HepMC3_FOUND REQUIRED_VARS HEPMC3_INCLUDE_DIR HEPMC3_LIBRARIES # VERSION_VAR # HEPMC3_VERSION ) endif() if(EXCLUDE_rivet) message(STATUS "Skipping rivet") else() find_package(rivet) if(rivet_FOUND) if(rivet_USE_HEPMC3 AND NOT HepMC3_FOUND) message(FATAL_ERROR "Rivet was compiled with HepMC version 3, " "but HepMC3 was not found!") elseif(NOT rivet_USE_HEPMC3 AND NOT HepMC_FOUND) message(FATAL_ERROR "Rivet was compiled with HepMC version 2, " "but HepMC2 was not found!") endif() endif() endif() if(EXCLUDE_QCDloop) message(STATUS "Skipping QCDloop") else() find_package(QCDloop 2) endif() if(NOT EXCLUDE_ROOT) # Don't print "Skipping ROOT" for default option (EXCLUDE_ROOT=TRUE) find_package(ROOT 6.14) # root targets are broken before 6.14 find_package_handle_standard_args( ROOT FOUND_VAR ROOT_FOUND REQUIRED_VARS ROOT_INCLUDE_DIRS ROOT_LIBRARIES VERSION_VAR ROOT_VERSION ) endif() if(EXCLUDE_HighFive) message(STATUS "Skipping HighFive") else() find_package(HighFive QUIET) find_package_handle_standard_args( HighFive CONFIG_MODE ) endif() include(RepairTargets) # more reliable cmake targets ## Save dependencies for templates set(HEJ_BUILD_WITH_HepMC2 ${HepMC_FOUND}) set(HEJ_BUILD_WITH_HepMC3 ${HepMC3_FOUND}) set(HEJ_BUILD_WITH_RIVET ${rivet_FOUND}) if(rivet_FOUND AND rivet_VERSION VERSION_LESS 3) set(HEJ_USE_RIVET2 TRUE) endif() set(HEJ_BUILD_WITH_QCDLOOP ${QCDloop_FOUND}) set(HEJ_BUILD_WITH_HDF5 ${HighFive_FOUND}) ## Template files configure_file( ${CMAKE_CURRENT_SOURCE_DIR}/cmake/Templates/Version.hh.in include/HEJ/Version.hh @ONLY ) configure_file( ${CMAKE_CURRENT_SOURCE_DIR}/cmake/Templates/ConfigFlags.hh.in include/HEJ/ConfigFlags.hh @ONLY ) configure_file( ${CMAKE_CURRENT_SOURCE_DIR}/cmake/Templates/HEJ-config.cc.in 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}) ## shortcut for HEJ specific includes set(HEJ_INCLUDE_DIR $ $ $ ) - +# Add test coverage if (TEST_COVERAGE AND (CMAKE_COMPILER_IS_GNUCC OR CMAKE_COMPILER_IS_GNUCXX) ) include(CodeCoverage) APPEND_COVERAGE_COMPILER_FLAGS() set(COVERAGE_GCOVR_EXCLUDES "${PROJECT_SOURCE_DIR}/include/LHEF/*" "${PROJECT_SOURCE_DIR}/include/cxxopts.hpp" "${PROJECT_SOURCE_DIR}/t/*" "${PROJECT_BINARY_DIR}/*") SETUP_TARGET_FOR_COVERAGE_GCOVR_HTML( NAME ctest_coverage EXECUTABLE ctest ) endif() add_subdirectory(current_generator) ## create main HEJ library add_subdirectory(src) ## define executable add_executable(HEJ_main src/bin/HEJ.cc) set_target_properties(HEJ_main PROPERTIES OUTPUT_NAME "HEJ") ## link libraries target_link_libraries(HEJ_main HEJ) add_executable(HEJ-config src/bin/HEJ-config.cc) target_include_directories(HEJ-config PRIVATE ${HEJ_INCLUDE_DIR}) 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_main HEJ-config DESTINATION ${INSTALL_BIN_DIR}) ## tests (prevent testing if current project is a subproject) if(CMAKE_PROJECT_NAME STREQUAL PROJECT_NAME) ## advanced version of enable_testing() ## adds "BUILD_TESTING" option to deactivate testing include(CTest) if(BUILD_TESTING) set(BUILD_EXAMPLES TRUE) # examples needed for testing set(include_tests TRUE) endif() endif() ## build examples if(BUILD_EXAMPLES) message(STATUS "Adding examples") add_subdirectory(examples) endif() ## build tests if(include_tests) message(STATUS "Adding tests") add_subdirectory(t) endif() add_subdirectory(doc) ## Clang-tidy if(RUN_CLANG_TIDY) find_program( CLANG_TIDY_EXE NAMES "clang-tidy" DOC "Path to clang-tidy executable" ) if(NOT CLANG_TIDY_EXE) message(FATAL_ERROR "clang-tidy not found, but requested. Please deactivate RUN_CLANG_TIDY" ) else() set_target_properties( HEJ PROPERTIES CXX_CLANG_TIDY "${CLANG_TIDY_EXE};-header-filter=${CMAKE_SOURCE_DIR}/include/HEJ;-fix;--fix-errors" ) endif() endif() diff --git a/Changes.md b/Changes.md index e6d2ade..8312713 100644 --- a/Changes.md +++ b/Changes.md @@ -1,140 +1,140 @@ # Changelog This is the log for changes to the HEJ program. Further changes to the HEJ API are documented in [`Changes-API.md`](Changes-API.md). If you are using HEJ as a library, please also read the changes there. ## Version 2.1 This release adds support for two new processes: * W boson with jets. * Jets with a charged lepton-antilepton pair via a virtual Z boson or photon. In addition, the complete set of first subleading processes (unordered gluon, central and extremal quark-antiquark pair) is implemented for pure jets and W + jets, see [arXiv:TODO](https://arxiv.org/abs/TODO). Unordered gluon emission is also supported for Higgs boson + jets and Z boson/photon + jets. This release include many changes to the code, which affect users of HEJ as a library (see [`Changes-API.md`](Changes-API.md)). ### 2.1.0 #### New Processes * 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 multiple jets together with a (off-shell) W bosons decaying into lepton & neutrino * Resummation for jets with a charged lepton-antilepton pair via a virtual Z boson or photon. Includes the `unordered` subleading process. * Resummation can now be performed on all subleading processes within pure jets also. This includes `unordered`, `extremal qqbar` and `central qqbar` processes. #### More Physics implementation * 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 * Use relative fraction for soft transverse momentum in extremal jets (`max ext soft pt fraction`). This supersedes `min extparton pt`, which is now deprecated and will be removed in future versions. #### Updates to Runcard * Allow multiplication and division of multiple scale functions e.g. `H_T/2*m_j1j2` * Grouped `event treatment` for subleading channels together in runcard - Rename `non-HEJ` processes to `non-resummable` * Read electro-weak constants from input - new mandatory setting `vev` to change vacuum expectation value - new mandatory settings `particle properties` to specify mass & width of bosons - - FOG: decays are now specified in `decays` setting (previously under + - `HEJFOG`: decays are now specified in `decays` setting (previously under `particle properties`) * Allow loading multiple analyses with `analyses`. The old `analysis` (with "i") is marked deprecated. * Optional setting to specify maximal number of Fixed Order events (`max events`, default is all) * Allow changing the regulator lambda in input (`regulator parameter`, only for advanced users) #### Changes to Input/Output * Added support to read & write `hdf5` event files suggested in [arXiv:1905.05120](https://arxiv.org/abs/1905.05120) (needs [HighFive](https://github.com/BlueBrain/HighFive)) * Support input with average weight equal to the cross section (`IDWTUP=1 or 4`) * Support unmodified Les Houches Event Files written by Sherpa with `cross section = sum(weights)/sum(trials)` * Analyses now get general run information (`LHEF::HEPRUP`) in the constructor. **This might break previously written, external analyses!** - external analyses should now be created with `make_analysis(YAML::Node const & config, LHEF::HEPRUP const & heprup)` * Support `rivet` version 3 with both `HepMC` version 2 and 3 - Multiple weights with `rivet 3` will only create one `.yoda` file (instead of one per weight/scale) * Added option to unweight only resummation events (`unweight: {type: resummation}`) * Added option for partially unweighting resummation events, similar to the fixed-order generator. * Improved unweighting algorithm. * Follow HepMC convention for particle Status codes: incoming = 11, decaying = 2, outgoing = 1 (unchanged) #### Miscellaneous * Print cross sections at end of run * Added example analysis & scale to `examples/`. Everything in `examples/` will be build when the flag `-DBUILD_EXAMPLES=TRUE` is set in `cmake`. * 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 * Require LHADPF version 6. Dropped support for all other versions. * Use `git-lfs` for raw data in test (`make test` now requires `git-lfs`) * Currents are now generated with [`FORM`](https://github.com/vermaseren/form) - `FORM` is included as a `git submodule`, use `git submodule update --init` to download `FORM` * Create [Sphinx](http://sphinx-doc.org/) and [Doxygen](http://doxygen.org/) documentation by `make sphinx` or `make doxygen` in your `build/` folder ## Version 2.0 First release of HEJ 2. Complete code rewrite compared to HEJ 1. Improved matching to Fixed Order ([arXiv:1805.04446](https://arxiv.org/abs/1805.04446)). Implemented processes: Higgs boson with jets (FKL and unordered gluon emission, with finite quark mass loop, [arXiv:1812.08072](https://arxiv.org/abs/1812.08072)), and pure jets (only FKL). See [arXiv:1902.08430](https://arxiv.org/abs/1902.08430) ### 2.0.6 * Fixed compiling rivet when YODA headers are _outside_ of rivet directory ### 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. ### 2.0.0 * First release diff --git a/FixedOrderGen/CMakeLists.txt b/FixedOrderGen/CMakeLists.txt index b17a9a1..b11a719 100644 --- a/FixedOrderGen/CMakeLists.txt +++ b/FixedOrderGen/CMakeLists.txt @@ -1,83 +1,89 @@ cmake_minimum_required(VERSION 3.1 FATAL_ERROR) set(CMAKE_LEGACY_CYGWIN_WIN32 0) set(CMAKE_EXPORT_COMPILE_COMMANDS ON) project("HEJ Fixed Order Generation" VERSION 2.0.6 LANGUAGES C CXX) +## User settings +include(CMakeDependentOption) + +option(TEST_ALL "Run additional (longer) tests" FALSE) +option(TEST_COVERAGE "Generate test coverage with \"gcovr\"" FALSE) + # 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() ## Global flags for the compiler (all warnings) if (CMAKE_COMPILER_IS_GNUCC OR CMAKE_COMPILER_IS_GNUCXX) set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wextra -Wpedantic") 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 ) CONFIGURE_FILE( ${CMAKE_CURRENT_SOURCE_DIR}/cmake/Templates/Version.hh.in include/Version.hh @ONLY ) ## Use cmake modules from HEJ src set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}/../cmake/Modules/") ## Find HEJ (only dependency) ## HEJ includes all sub dependencies (fastjet, lhapdf, ...) find_package(HEJ 2 REQUIRED) include(RepairTargets) # more reliable cmake targets -## define executable -file(GLOB HEJFOG_source ${CMAKE_CURRENT_SOURCE_DIR}/src/*.cc) -list(REMOVE_ITEM HEJFOG_source ${CMAKE_CURRENT_SOURCE_DIR}/src/main.cc) -add_library(hejfog_lib STATIC ${HEJFOG_source}) -target_include_directories(hejfog_lib - PUBLIC - ${PROJECT_SOURCE_DIR}/include - ${PROJECT_BINARY_DIR}/include -) -add_executable(HEJFOG ${CMAKE_CURRENT_SOURCE_DIR}/src/main.cc) - -## link libraries -set(libraries ${CMAKE_DL_LIBS}) +# Add test coverage +if (TEST_COVERAGE AND (CMAKE_COMPILER_IS_GNUCC OR CMAKE_COMPILER_IS_GNUCXX) ) + include(CodeCoverage) + APPEND_COVERAGE_COMPILER_FLAGS() + set(COVERAGE_GCOVR_EXCLUDES "${PROJECT_SOURCE_DIR}/t/*" + "${PROJECT_BINARY_DIR}/*") + SETUP_TARGET_FOR_COVERAGE_GCOVR_HTML( + NAME ctest_coverage + EXECUTABLE ctest + ) +endif() -target_link_libraries(hejfog_lib ${libraries} HEJ::HEJ) +## define executable +add_subdirectory(src) +add_executable(HEJFOG src/bin/main.cc) target_link_libraries(HEJFOG hejfog_lib) install(TARGETS HEJFOG DESTINATION bin) ## tests include(CTest) ## advanced version of enable_testing() ## adds "BUILD_TESTING" option to deactivate testing if(BUILD_TESTING) add_subdirectory(t) endif() diff --git a/FixedOrderGen/configFO.yml b/FixedOrderGen/configFO.yml index 6fca80b..2ffcb83 100644 --- a/FixedOrderGen/configFO.yml +++ b/FixedOrderGen/configFO.yml @@ -1,94 +1,95 @@ ## Number of generated events events: 200 jets: min pt: 20 # minimal jet pt, should be slightly below analysis cut peak pt: 30 # peak pt of jets, should be at analysis cut algorithm: antikt # jet clustering algorithm R: 0.4 # jet R parameter max rapidity: 5 # maximum jet rapidity ## Particle beam beam: energy: 6500 particles: [p, p] ## PDF ID according to LHAPDF pdf: 230000 ## Scattering process process: p p => h 4j ## Particle decays (multiple decays are allowed) decays: Higgs: {into: [photon, photon], branching ratio: 0.0023568762400521404} Wp: {into: [e+, nu_e], branching ratio: 1} Wm: {into: [e-, nu_e_bar]} # equivalent to "branching ratio: 1" ## Fraction of events with two extremal emissions in one direction ## that contain an subleading emission e.g. unordered emission subleading fraction: 0.05 ## Allow different subleading configurations. ## By default all processes are allowed. ## This does not check if the processes are implemented in HEJ! # # subleading channels: # - unordered -# - qqx +# - central qqx +# - extremal qqx ## Unweighting parameters ## remove to obtain weighted events # unweight: # sample size: 200 # should be similar to "events:", but not more than ~10000 # max deviation: 0 ## --- The following settings are equivalent to HEJ, --- ## ## --- see HEJ for reference. --- ## ## Central scale choice scales: max jet pperp ## Event output files event output: - HEJFO.lhe # - HEJFO_events.hepmc ## Analyses # # analyses: ## Rivet analysis # - rivet: MC_XS # rivet analysis name # output: HEJ # name of the yoda files, ".yoda" and scale suffix will be added ## Custom analysis # - plugin: /path/to/libmyanalysis.so # my analysis parameter: some value ## Selection of random number generator and seed random generator: name: mixmax # seed: 1 ## Vacuum expectation value vev: 246.2196508 ## Properties of the weak gauge bosons particle properties: Higgs: mass: 125 width: 0.004165 W: mass: 80.385 width: 2.085 Z: mass: 91.187 width: 2.495 ## Parameters for Higgs-gluon couplings ## this requires compilation with QCDloop # # Higgs coupling: # use impact factors: false # mt: 174 # include bottom: true # mb: 4.7 diff --git a/FixedOrderGen/include/EventGenerator.hh b/FixedOrderGen/include/EventGenerator.hh index 680deec..8848eeb 100644 --- a/FixedOrderGen/include/EventGenerator.hh +++ b/FixedOrderGen/include/EventGenerator.hh @@ -1,67 +1,68 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #pragma once #include #include "HEJ/EWConstants.hh" #include "HEJ/MatrixElement.hh" #include "HEJ/optional.hh" #include "HEJ/PDF.hh" #include "HEJ/ScaleFunction.hh" #include "Beam.hh" #include "Decay.hh" #include "JetParameters.hh" #include "Process.hh" #include "Status.hh" +#include "Subleading.hh" namespace HEJ { class Event; class HiggsCouplingSettings; class RNG; } //! Namespace for HEJ Fixed Order Generator namespace HEJFOG { class EventGenerator{ public: EventGenerator( Process process, Beam beam, HEJ::ScaleGenerator scale_gen, JetParameters jets, int pdf_id, - double subl_change, - unsigned int subl_channels, + double subl_chance, + Subleading subl_channels, ParticlesDecayMap particle_decays, HEJ::HiggsCouplingSettings Higgs_coupling, HEJ::EWConstants ew_parameters, std::shared_ptr ran ); HEJ::optional gen_event(); Status status() const { return status_; } private: HEJ::PDF pdf_; HEJ::MatrixElement ME_; HEJ::ScaleGenerator scale_gen_; Process process_; JetParameters jets_; Beam beam_; Status status_; - double subl_change_; - unsigned int subl_channels_; + double subl_chance_; + Subleading subl_channels_; ParticlesDecayMap particle_decays_; HEJ::EWConstants ew_parameters_; std::shared_ptr ran_; }; } diff --git a/FixedOrderGen/include/PhaseSpacePoint.hh b/FixedOrderGen/include/PhaseSpacePoint.hh index e88e82a..9873390 100644 --- a/FixedOrderGen/include/PhaseSpacePoint.hh +++ b/FixedOrderGen/include/PhaseSpacePoint.hh @@ -1,238 +1,298 @@ /** \file PhaseSpacePoint.hh * \brief Contains the PhaseSpacePoint Class * * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #pragma once #include #include #include #include #include +#include "boost/iterator/filter_iterator.hpp" + #include "HEJ/Event.hh" #include "HEJ/Particle.hh" #include "HEJ/PDG_codes.hh" #include "fastjet/PseudoJet.hh" #include "Decay.hh" #include "Status.hh" +#include "Subleading.hh" namespace HEJ { class EWConstants; class PDF; class RNG; } namespace HEJFOG { class JetParameters; class Process; //! A point in resummation phase space class PhaseSpacePoint{ public: + + //! Iterator over partons + using ConstPartonIterator = boost::filter_iterator< + bool (*)(HEJ::Particle const &), + std::vector::const_iterator + >; + //! Reverse Iterator over partons + using ConstReversePartonIterator = std::reverse_iterator< + ConstPartonIterator>; + //! Default PhaseSpacePoint Constructor PhaseSpacePoint() = delete; //! PhaseSpacePoint Constructor /** * @param proc The process to generate * @param jet_properties Jet defintion & cuts * @param pdf The pdf set (used for sampling) * @param E_beam Energie of the beam * @param subl_chance Chance to turn a potentially unordered * emission into an actual one * @param subl_channels Possible subleading channels. * see HEJFOG::Subleading * @param particle_properties Properties of producted boson * * Initially, only FKL phase space points are generated. subl_chance gives - * the change of turning one emissions into a subleading configuration, + * the chance of turning one emissions into a subleading configuration, * i.e. either unordered or central quark/anti-quark pair. Unordered * emissions require that the most extremal emission in any direction is * a quark or anti-quark and the next emission is a gluon. Quark/anti-quark * pairs are only generated for W processes. At most one subleading * emission will be generated in this way. */ PhaseSpacePoint( Process const & proc, JetParameters const & jet_properties, HEJ::PDF & pdf, double E_beam, double subl_chance, - unsigned int subl_channels, + Subleading subl_channels, ParticlesDecayMap const & particle_decays, HEJ::EWConstants const & ew_parameters, HEJ::RNG & ran ); //! Get Weight Function /** * @returns Weight of Event */ double weight() const{ return weight_; } Status status() const{ return status_; } //! Get Incoming Function /** * @returns Incoming Particles */ std::array const & incoming() const{ return incoming_; } //! Get Outgoing Function /** * @returns Outgoing Particles */ std::vector const & outgoing() const{ return outgoing_; } std::unordered_map> const & decays() const{ return decays_; } + //! Iterator to the first outgoing parton + ConstPartonIterator begin_partons() const; + //! Iterator to the first outgoing parton + ConstPartonIterator cbegin_partons() const; + + //! Iterator to the end of the outgoing partons + ConstPartonIterator end_partons() const; + //! Iterator to the end of the outgoing partons + ConstPartonIterator cend_partons() const; + + //! Reverse Iterator to the first outgoing parton + ConstReversePartonIterator rbegin_partons() const; + //! Reverse Iterator to the first outgoing parton + ConstReversePartonIterator crbegin_partons() const; + //! Reverse Iterator to the first outgoing parton + ConstReversePartonIterator rend_partons() const; + //! Reverse Iterator to the first outgoing parton + ConstReversePartonIterator crend_partons() const; + private: friend HEJ::Event::EventData to_EventData(PhaseSpacePoint psp); /** * @internal * @brief Generate LO parton momentum * * @param count Number of partons to generate * @param is_pure_jets If true ensures momentum conservation in x and y * @param jet_param Jet properties to fulfil * @param max_pt max allowed pt for a parton (typically E_CMS) * @param ran Random Number Generator * * @returns Momentum of partons * * Ensures that each parton is in its own jet. * Generation is independent of parton flavour. Output is sorted in rapidity. */ std::vector gen_LO_partons( int count, bool is_pure_jets, JetParameters const & jet_param, double max_pt, HEJ::RNG & ran ); HEJ::Particle gen_boson( HEJ::ParticleID bosonid, double mass, double width, HEJ::RNG & ran ); template fastjet::PseudoJet gen_last_momentum( ParticleMomenta const & other_momenta, double mass_square, double y ) const; bool jets_ok( std::vector const & Born_jets, std::vector const & partons ) const; /** * @internal * @brief Generate incoming partons according to the PDF * * @param uf Scale used in the PDF */ void reconstruct_incoming( - Process const & proc, unsigned int subl_channels, + Process const & proc, + double subl_chance, Subleading subl_channels, HEJ::PDF & pdf, double E_beam, double uf, HEJ::RNG & ran ); - /** - * @internal - * @brief Returns list of all allowed initial states partons - */ - std::array,2> filter_partons( - Process const & proc, unsigned int const subl_channels, - HEJ::RNG & ran - ); HEJ::ParticleID generate_incoming_id( std::size_t beam_idx, double x, double uf, HEJ::PDF & pdf, std::bitset<11> allowed_partons, HEJ::RNG & ran ); + //! @brief Returns list of all allowed initial states partons + std::array,2> allowed_incoming( + Process const & proc, + double subl_chance, Subleading subl_channels, + HEJ::RNG & ran + ); + + //! Ensure that incoming state compatible with A/W/Z production + std::array,2> incoming_AWZ( + Process const & proc, Subleading subl_channels, + std::array,2> allowed_partons, + HEJ::RNG & ran + ); + //! Ensure that incoming state compatible with extremal qqx + std::array,2> incoming_eqqx( + std::array,2> allowed_partons, HEJ::RNG & ran); + //! Ensure that incoming state compatible with unordered + std::array,2> incoming_uno( + std::array,2> allowed_partons, HEJ::RNG & ran); bool momentum_conserved(double ep) const; - HEJ::Particle const & most_backward_FKL( - std::vector const & partons - ) const; - HEJ::Particle const & most_forward_FKL( - std::vector const & partons - ) const; - HEJ::Particle & most_backward_FKL(std::vector & partons) const; - HEJ::Particle & most_forward_FKL(std::vector & partons) const; bool extremal_FKL_ok( std::vector const & partons ) const; double random_normal(double stddev, HEJ::RNG & ran); /** * @internal * @brief Turns a FKL configuration into a subleading one * - * @param chance Change to switch to subleading configuration + * @param chance Chance to switch to subleading configuration * @param channels Allowed channels for subleading process * @param proc Process to decide which subleading * configurations are allowed * * With a chance of "chance" the FKL configuration is either turned into * a unordered configuration or, for A/W/Z bosons, a configuration with * a central quark/anti-quark pair. */ - void maybe_turn_to_subl(double chance, unsigned int channels, + void maybe_turn_to_subl(double chance, Subleading channels, Process const & proc, HEJ::RNG & ran); + void turn_to_subl( + Subleading channels, + bool can_be_uno_backward, bool can_be_uno_forward, bool allow_strange, + HEJ::RNG & ran); void turn_to_uno(bool can_be_uno_backward, bool can_be_uno_forward, HEJ::RNG & ran); - void turn_to_qqx(bool allow_strange, HEJ::RNG & ran); + HEJ::ParticleID select_qqx_flavour(bool allow_strange, HEJ::RNG & ran); + void turn_to_cqqx(bool allow_strange, HEJ::RNG & ran); + void turn_to_eqqx(bool allow_strange, HEJ::RNG & ran); //! decay where we select the decay channel std::vector decay_boson( HEJ::Particle const & parent, std::vector const & decays, HEJ::RNG & ran ); //! generate decay products of a boson std::vector decay_boson( HEJ::Particle const & parent, std::vector const & decays, HEJ::RNG & ran ); /// @brief setup outgoing partons to ensure correct coupling to boson void couple_boson(HEJ::ParticleID boson, HEJ::RNG & ran); Decay select_decay_channel( std::vector const & decays, HEJ::RNG & ran ); double gen_hard_pt( int np, double ptmin, double ptmax, double y, HEJ::RNG & ran ); double gen_soft_pt(int np, double ptmax, HEJ::RNG & ran); double gen_parton_pt( int count, JetParameters const & jet_param, double ptmax, double y, HEJ::RNG & ran ); + //! Iterator over partons (non-const) + using PartonIterator = boost::filter_iterator< + bool (*)(HEJ::Particle const &), + std::vector::iterator + >; + //! Reverse Iterator over partons (non-const) + using ReversePartonIterator = std::reverse_iterator; + + //! Iterator to the first outgoing parton (non-const) + PartonIterator begin_partons(); + //! Iterator to the end of the outgoing partons (non-const) + PartonIterator end_partons(); + + //! Reverse Iterator to the first outgoing parton (non-const) + ReversePartonIterator rbegin_partons(); + //! Reverse Iterator to the first outgoing parton (non-const) + ReversePartonIterator rend_partons(); + double weight_; Status status_; std::array incoming_; std::vector outgoing_; //! Particle decays in the format {outgoing index, decay products} std::unordered_map> decays_; }; //! Extract HEJ::Event::EventData from PhaseSpacePoint HEJ::Event::EventData to_EventData(PhaseSpacePoint psp); } diff --git a/FixedOrderGen/include/Status.hh b/FixedOrderGen/include/Status.hh index 0b5e926..936421f 100644 --- a/FixedOrderGen/include/Status.hh +++ b/FixedOrderGen/include/Status.hh @@ -1,27 +1,27 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019 * \copyright GPLv2 or later */ #pragma once #include #include namespace HEJFOG { - enum Status{ + enum class Status{ good, not_enough_jets, too_much_energy }; inline std::string to_string(Status s){ switch(s){ - case good: return "good"; - case not_enough_jets: return "not enough jets"; - case too_much_energy: return "too much energy"; + case Status::good: return "good"; + case Status::not_enough_jets: return "not enough jets"; + case Status::too_much_energy: return "too much energy"; default:; } throw std::logic_error{"unreachable"}; } } diff --git a/FixedOrderGen/include/Subleading.hh b/FixedOrderGen/include/Subleading.hh index ef76017..c152cde 100644 --- a/FixedOrderGen/include/Subleading.hh +++ b/FixedOrderGen/include/Subleading.hh @@ -1,24 +1,34 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019 * \copyright GPLv2 or later */ #pragma once +#include + namespace HEJFOG { - namespace channels{ + namespace subleading { //!< @TODO confusing name with capital Subleading /** * Bit position of different subleading channels * e.g. (unsigned int) 1 => only unordered */ - enum Subleading: unsigned { - none = 0u, - all = ~0u, - uno = 1u, + enum Channels: unsigned { + uno, unordered = uno, - qqx = 2u + cqqx, + central_qqx = cqqx, + eqqx, + extremal_qqx = eqqx, + first = uno, + last = eqqx, }; } - using Subleading = channels::Subleading; + using Subleading = std::bitset; + + namespace subleading { + static constexpr Subleading ALL{~0u}; + static constexpr Subleading NONE{0u}; + } } diff --git a/FixedOrderGen/include/config.hh b/FixedOrderGen/include/config.hh index 212f059..69c62e4 100644 --- a/FixedOrderGen/include/config.hh +++ b/FixedOrderGen/include/config.hh @@ -1,49 +1,50 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #pragma once #include #include #include #include "HEJ/Config.hh" #include "HEJ/EWConstants.hh" #include "HEJ/Fraction.hh" #include "HEJ/HiggsCouplingSettings.hh" #include "HEJ/optional.hh" #include "HEJ/output_formats.hh" #include "yaml-cpp/yaml.h" #include "Beam.hh" #include "Decay.hh" #include "JetParameters.hh" #include "Process.hh" +#include "Subleading.hh" #include "UnweightSettings.hh" namespace HEJFOG { struct Config{ Process process; std::size_t events; JetParameters jets; Beam beam; int pdf_id; HEJ::Fraction subleading_fraction; - unsigned int subleading_channels; //! < see HEJFOG::Subleading + Subleading subleading_channels; //! < see HEJFOG::Subleading ParticlesDecayMap particle_decays; std::vector analyses_parameters; HEJ::ScaleConfig scales; std::vector output; HEJ::RNGConfig rng; HEJ::HiggsCouplingSettings Higgs_coupling; HEJ::EWConstants ew_parameters; HEJ::optional unweight; }; Config load_config(std::string const & config_file); } diff --git a/FixedOrderGen/src/CMakeLists.txt b/FixedOrderGen/src/CMakeLists.txt new file mode 100644 index 0000000..b03336d --- /dev/null +++ b/FixedOrderGen/src/CMakeLists.txt @@ -0,0 +1,13 @@ +file(GLOB src_files *.cc) + +add_library(hejfog_lib STATIC ${src_files}) +target_include_directories(hejfog_lib + PUBLIC + ${PROJECT_SOURCE_DIR}/include + ${PROJECT_BINARY_DIR}/include +) + +## link libraries +set(libraries ${CMAKE_DL_LIBS}) + +target_link_libraries(hejfog_lib ${libraries} HEJ::HEJ) diff --git a/FixedOrderGen/src/EventGenerator.cc b/FixedOrderGen/src/EventGenerator.cc index 84b03c7..935b1eb 100644 --- a/FixedOrderGen/src/EventGenerator.cc +++ b/FixedOrderGen/src/EventGenerator.cc @@ -1,99 +1,100 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #include "EventGenerator.hh" #include #include #include "HEJ/Config.hh" #include "HEJ/Event.hh" #include "HEJ/event_types.hh" #include "HEJ/EWConstants.hh" #include "HEJ/exceptions.hh" #include "HEJ/HiggsCouplingSettings.hh" #include "Process.hh" #include "Beam.hh" #include "JetParameters.hh" #include "PhaseSpacePoint.hh" namespace HEJFOG { EventGenerator::EventGenerator( Process process, Beam beam, HEJ::ScaleGenerator scale_gen, JetParameters jets, int pdf_id, - double subl_change, - unsigned int subl_channels, + double subl_chance, + Subleading subl_channels, ParticlesDecayMap particle_decays, HEJ::HiggsCouplingSettings Higgs_coupling, HEJ::EWConstants ew_parameters, std::shared_ptr ran ): pdf_{pdf_id, beam.particles[0], beam.particles[1]}, ME_{ [this](double mu){ return pdf_.Halphas(mu); }, HEJ::MatrixElementConfig{ false, std::move(Higgs_coupling), ew_parameters } }, scale_gen_{std::move(scale_gen)}, process_{std::move(process)}, jets_{std::move(jets)}, beam_{std::move(beam)}, - subl_change_{subl_change}, + subl_chance_{subl_chance}, subl_channels_{subl_channels}, particle_decays_{std::move(particle_decays)}, ew_parameters_{ew_parameters}, ran_{std::move(ran)} { } HEJ::optional EventGenerator::gen_event(){ HEJFOG::PhaseSpacePoint psp{ process_, jets_, pdf_, beam_.energy, - subl_change_, subl_channels_, + subl_chance_, subl_channels_, particle_decays_, ew_parameters_, *ran_ }; status_ = psp.status(); - if(status_ != good) return {}; + if(status_ != Status::good) return {}; HEJ::Event ev = scale_gen_( HEJ::Event{ to_EventData( std::move(psp) ).cluster( jets_.def, jets_.min_pt) } ); - if(!is_resummable(ev.type())) + if(!is_resummable(ev.type())){ throw HEJ::not_implemented("Tried to generate a event type, " "which is not yet implemented in HEJ."); + } ev.generate_colours(*ran_); const double shat = HEJ::shat(ev); const double xa = (ev.incoming()[0].E()-ev.incoming()[0].pz())/(2.*beam_.energy); const double xb = (ev.incoming()[1].E()+ev.incoming()[1].pz())/(2.*beam_.energy); // evaluate matrix element ev.parameters() *= ME_.tree(ev)/(shat*shat); // and PDFs ev.central().weight *= pdf_.pdfpt(0,xa,ev.central().muf, ev.incoming()[0].type); ev.central().weight *= pdf_.pdfpt(0,xb,ev.central().muf, ev.incoming()[1].type); for(std::size_t i = 0; i < ev.variations().size(); ++i){ auto & var = ev.variations(i); var.weight *= pdf_.pdfpt(0,xa,var.muf, ev.incoming()[0].type); var.weight *= pdf_.pdfpt(0,xb,var.muf, ev.incoming()[1].type); } return ev; } } diff --git a/FixedOrderGen/src/PhaseSpacePoint.cc b/FixedOrderGen/src/PhaseSpacePoint.cc index 7683358..66135a1 100644 --- a/FixedOrderGen/src/PhaseSpacePoint.cc +++ b/FixedOrderGen/src/PhaseSpacePoint.cc @@ -1,701 +1,974 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #include "PhaseSpacePoint.hh" #include #include #include #include #include #include #include #include #include #include #include "fastjet/ClusterSequence.hh" #include "HEJ/Constants.hh" #include "HEJ/EWConstants.hh" #include "HEJ/exceptions.hh" #include "HEJ/kinematics.hh" #include "HEJ/Particle.hh" #include "HEJ/PDF.hh" #include "HEJ/RNG.hh" #include "HEJ/utility.hh" #include "JetParameters.hh" #include "Process.hh" -#include "Subleading.hh" namespace { - using namespace HEJ; static_assert( std::numeric_limits::has_quiet_NaN, "no quiet NaN for double" ); constexpr double NaN = std::numeric_limits::quiet_NaN(); } // namespace anonymous namespace HEJFOG { - Event::EventData to_EventData(PhaseSpacePoint psp){ + HEJ::Event::EventData to_EventData(PhaseSpacePoint psp){ //! @TODO Same function already in HEJ - Event::EventData result; + HEJ::Event::EventData result; result.incoming = std::move(psp).incoming_; assert(result.incoming.size() == 2); result.outgoing = std::move(psp).outgoing_; // technically Event::EventData doesn't have to be sorted, // but PhaseSpacePoint should be anyway assert( std::is_sorted( begin(result.outgoing), end(result.outgoing), - rapidity_less{} + HEJ::rapidity_less{} ) ); assert(result.outgoing.size() >= 2); result.decays = std::move(psp).decays_; result.parameters.central = {NaN, NaN, psp.weight()}; return result; } + PhaseSpacePoint::ConstPartonIterator PhaseSpacePoint::begin_partons() const { + return cbegin_partons(); + } + PhaseSpacePoint::ConstPartonIterator PhaseSpacePoint::cbegin_partons() const { + return {HEJ::is_parton, cbegin(outgoing()), cend(outgoing())}; + } + + PhaseSpacePoint::ConstPartonIterator PhaseSpacePoint::end_partons() const { + return cend_partons(); + } + PhaseSpacePoint::ConstPartonIterator PhaseSpacePoint::cend_partons() const { + return {HEJ::is_parton, cend(outgoing()), cend(outgoing())}; + } + + PhaseSpacePoint::ConstReversePartonIterator PhaseSpacePoint::rbegin_partons() const { + return crbegin_partons(); + } + PhaseSpacePoint::ConstReversePartonIterator PhaseSpacePoint::crbegin_partons() const { + return std::reverse_iterator( cend_partons() ); + } + + PhaseSpacePoint::ConstReversePartonIterator PhaseSpacePoint::rend_partons() const { + return crend_partons(); + } + PhaseSpacePoint::ConstReversePartonIterator PhaseSpacePoint::crend_partons() const { + return std::reverse_iterator( cbegin_partons() ); + } + + PhaseSpacePoint::PartonIterator PhaseSpacePoint::begin_partons() { + return {HEJ::is_parton, begin(outgoing_), end(outgoing_)}; + } + + PhaseSpacePoint::PartonIterator PhaseSpacePoint::end_partons() { + return {HEJ::is_parton, end(outgoing_), end(outgoing_)}; + } + + PhaseSpacePoint::ReversePartonIterator PhaseSpacePoint::rbegin_partons() { + return std::reverse_iterator( end_partons() ); + } + + PhaseSpacePoint::ReversePartonIterator PhaseSpacePoint::rend_partons() { + return std::reverse_iterator( begin_partons() ); + } + namespace { bool can_swap_to_uno( - Particle const & p1, Particle const & p2 + HEJ::Particle const & p1, HEJ::Particle const & p2 ){ - return is_parton(p1) - && p1.type != pid::gluon - && p2.type == pid::gluon; + assert(is_parton(p1) && is_parton(p2)); + return p1.type != HEJ::pid::gluon + && p2.type == HEJ::pid::gluon; } - size_t count_gluons(std::vector::const_iterator first, - std::vector::const_iterator last + size_t count_gluons(PhaseSpacePoint::ConstPartonIterator first, + PhaseSpacePoint::ConstPartonIterator last ){ - return std::count_if(first, last, [](Particle const & p) - {return p.type == pid::gluon;}); + return std::count_if(first, last, [](HEJ::Particle const & p) + {return p.type == HEJ::pid::gluon;}); } /** assumes FKL configurations between first and last, * else there can be a quark in a non-extreme position * e.g. uno configuration gqg would pass */ - bool can_change_to_qqx( - std::vector::const_iterator first, - std::vector::const_iterator last + Subleading possible_qqx( + PhaseSpacePoint::ConstPartonIterator first, + PhaseSpacePoint::ConstReversePartonIterator last ){ - return 1 < count_gluons(first,last); + using namespace subleading; + assert( std::distance( first,last.base() )>2 ); + Subleading channels = ALL; + channels.reset(eqqx); + channels.reset(cqqx); + auto const ngluon = count_gluons(first,last.base()); + if(ngluon < 2) return channels; + if(first->type==HEJ::pid::gluon || last->type==HEJ::pid::gluon){ + channels.set(eqqx); + } + if(std::distance(first,last.base())>=4){ + channels.set(cqqx); + } + return channels; } + + template + bool uno_possible( + PartonIt first_parton, OutIt first_out + ){ + using namespace HEJ; + // Special case: Higgs can not be outside of uno + if(first_out->type == pid::Higgs + || std::next(first_out)->type==pid::Higgs){ + return false; + } + // decide what kind of subleading process is allowed + return can_swap_to_uno( *first_parton, *std::next(first_parton) ); + } + bool is_AWZ_proccess(Process const & proc){ - return proc.boson && is_AWZ_boson(*proc.boson); + return proc.boson && HEJ::is_AWZ_boson(*proc.boson); } - bool is_up_type(Particle const & part){ + bool is_up_type(HEJ::Particle const & part){ return is_anyquark(part) && !(std::abs(part.type)%2); } - bool is_down_type(Particle const & part){ + bool is_down_type(HEJ::Particle const & part){ return is_anyquark(part) && std::abs(part.type)%2; } - bool can_couple_to_W(Particle const & part, pid::ParticleID const W_id){ + bool can_couple_to_W( + HEJ::Particle const & part, HEJ::pid::ParticleID const W_id + ){ const int W_charge = W_id>0?1:-1; return std::abs(part.type)<5 && ( (W_charge*part.type > 0 && is_up_type(part)) || (W_charge*part.type < 0 && is_down_type(part)) ); } + + Subleading ensure_AWZ( + double & subl_chance, bool & allow_strange, + HEJ::ParticleID const boson, + PhaseSpacePoint::ConstPartonIterator first_parton, + PhaseSpacePoint::ConstPartonIterator last_parton + ){ + auto channels = subleading::ALL; + if(std::none_of(first_parton, last_parton, + [&boson](HEJ::Particle const & p){ + return can_couple_to_W(p, boson);})) { + // enforce qqx if A/W/Z can't couple somewhere else + // this is ensured to work through filter_partons in reconstruct_incoming + channels.reset(subleading::uno); + assert(channels.any()); + subl_chance = 1.; + // strange not allowed for W + if(std::abs(boson)== HEJ::pid::Wp) allow_strange = false; + } + return channels; + } + } + + void PhaseSpacePoint::turn_to_subl( + Subleading const channels, + bool const can_be_uno_backward, bool const can_be_uno_forward, + bool const allow_strange, + HEJ::RNG & ran + ){ + double const nchannels = channels.count(); + double const step = 1./nchannels; + weight_*=nchannels; + unsigned selected = subleading::first; + double rnd = nchannels>1?ran.flat():0.; + // @TODO optimise this sampling + for(; selected<=subleading::last; ++selected){ + assert(rnd>=0); + if(channels[selected]){ + if(rnd= 2); // decide what kind of subleading process is allowed - bool allow_uno = false; + bool const can_be_uno_backward = uno_possible(cbegin_partons(), + outgoing_.cbegin()); + bool const can_be_uno_forward = uno_possible(crbegin_partons(), + outgoing_.crbegin()); + if(channels[subleading::uno]){ + channels.set(subleading::uno, can_be_uno_backward || can_be_uno_forward); + } + channels &= possible_qqx(cbegin_partons(), crbegin_partons()); + bool allow_strange = true; - const size_t nout = outgoing_.size(); - const bool can_be_uno_backward = (channels&Subleading::uno) - && can_swap_to_uno(outgoing_[0], outgoing_[1]); - const bool can_be_uno_forward = (channels&Subleading::uno) - && can_swap_to_uno(outgoing_[nout-1], outgoing_[nout-2]); - allow_uno = can_be_uno_backward || can_be_uno_forward; - - bool allow_qqx = (channels&Subleading::qqx) - && can_change_to_qqx(outgoing_.cbegin(), outgoing_.cend()); if(is_AWZ_proccess(proc)) { - if(std::none_of(outgoing_.cbegin(), outgoing_.cend(), - [&proc](Particle const & p){ - return can_couple_to_W(p, *proc.boson);})) { - // enforce qqx if A/W/Z can't couple somewhere else - assert(allow_qqx); - allow_uno = false; - chance = 1.; - // strange not allowed for W - if(std::abs(*proc.boson)== pid::Wp) allow_strange = false; - } + channels &= ensure_AWZ(chance, allow_strange, *proc.boson, + cbegin_partons(), cend_partons()); } - if(!allow_uno && !allow_qqx) return; - if(ran.flat() < chance){ - weight_ /= chance; - if(allow_uno && !allow_qqx){ - turn_to_uno(can_be_uno_backward, can_be_uno_forward, ran); - } else if (!allow_uno && allow_qqx) { - turn_to_qqx(allow_strange, ran); - } else { - assert( allow_uno && allow_qqx); - if(ran.flat() < 0.5) turn_to_uno(can_be_uno_backward, can_be_uno_forward, ran); - else turn_to_qqx(allow_strange, ran); - weight_ *= 2.; - } - } else weight_ /= 1 - chance; + std::size_t const nchannels = channels.count(); + // no subleading + if(nchannels==0) return; + if(ran.flat() >= chance){ + weight_ /= 1 - chance; + return; + } + weight_ /= chance; + // select channel + return turn_to_subl( channels, can_be_uno_backward, can_be_uno_forward, + allow_strange, ran); } void PhaseSpacePoint::turn_to_uno( const bool can_be_uno_backward, const bool can_be_uno_forward, - RNG & ran + HEJ::RNG & ran ){ if(!can_be_uno_backward && !can_be_uno_forward) return; - const size_t nout = outgoing_.size(); if(can_be_uno_backward && can_be_uno_forward){ + weight_ *= 2.; if(ran.flat() < 0.5){ - std::swap(outgoing_[0].type, outgoing_[1].type); - } else { - std::swap(outgoing_[nout-1].type, outgoing_[nout-2].type); + return std::swap(begin_partons()->type, std::next(begin_partons())->type); } - weight_ *= 2.; - } else if(can_be_uno_backward){ - std::swap(outgoing_[0].type, outgoing_[1].type); - } else { - assert(can_be_uno_forward); - std::swap(outgoing_[nout-1].type, outgoing_[nout-2].type); + return std::swap(rbegin_partons()->type, std::next(rbegin_partons())->type); } + if(can_be_uno_backward){ + return std::swap(begin_partons()->type, std::next(begin_partons())->type); + } + assert(can_be_uno_forward); + std::swap(rbegin_partons()->type, std::next(rbegin_partons())->type); } - void PhaseSpacePoint::turn_to_qqx(const bool allow_strange, RNG & ran){ - /// find first and last gluon in FKL chain - auto first = std::find_if(outgoing_.begin(), outgoing_.end(), - [](Particle const & p){return p.type == pid::gluon;}); - std::vector FKL_gluons; - for(auto p = first; p!=outgoing_.end(); ++p){ - if(p->type == pid::gluon) FKL_gluons.push_back(&*p); - else if(is_anyquark(*p)) break; - } - const size_t ng = FKL_gluons.size(); - if(ng < 2) - throw std::logic_error("not enough gluons to create qqx"); - // select flavour of quark + //! select flavour of quark + HEJ::ParticleID PhaseSpacePoint::select_qqx_flavour( + const bool allow_strange, HEJ::RNG & ran + ){ const double r1 = 2.*ran.flat()-1.; - const double max_flavour = allow_strange?N_F:N_F-1; + const double max_flavour = allow_strange?HEJ::N_F:HEJ::N_F-1; weight_ *= max_flavour*2; - int flavour = pid::down + std::floor(std::abs(r1)*max_flavour); - flavour*=r1<0.?-1:1; + double const flavour = HEJ::pid::down + std::floor(std::abs(r1)*max_flavour); + return static_cast(flavour*(r1<0.?-1:1)); + } + + void PhaseSpacePoint::turn_to_cqqx(const bool allow_strange, HEJ::RNG & ran){ + // we assume all FKL partons to be gluons + auto first = std::next(begin_partons()); + auto last = std::next(rbegin_partons()); + auto const ng = std::distance(first, last.base()); + if(ng < 2) + throw std::logic_error("not enough gluons to create qqx"); + auto flavour = select_qqx_flavour(allow_strange, ran); // select gluon for switch - const size_t idx = std::floor((ng-1) * ran.flat()); - weight_ *= (ng-1); - FKL_gluons[idx]->type = ParticleID(flavour); - FKL_gluons[idx+1]->type = ParticleID(-flavour); + if(ng!=2){ + const double steps = 1./(ng-1); + weight_ /= steps; + for(auto rnd = ran.flat(); rnd>steps; ++first){ + rnd-=steps; + } + } + first->type = flavour; + std::next(first)->type = anti(flavour); + } + + void PhaseSpacePoint::turn_to_eqqx(const bool allow_strange, HEJ::RNG & ran){ + /// find first and last gluon in FKL chain + auto first = begin_partons(); + const bool can_forward = !is_anyquark(*first); + auto last = rbegin_partons(); + const bool can_backward = !is_anyquark(*last); + if(std::distance(first, last.base()) < 2) + throw std::logic_error("not enough gluons to create qqx"); + auto flavour = select_qqx_flavour(allow_strange, ran); + // select gluon for switch + if(can_forward && !can_backward){ + first->type = flavour; + std::next(first)->type = anti(flavour); + return; + } + if(!can_forward && can_backward){ + last->type = flavour; + std::next(last)->type = anti(flavour); + return; + } + assert(can_forward && can_backward); + weight_*=2.; + if(ran.flat()>0.5){ + first->type = flavour; + std::next(first)->type = anti(flavour); + return; + } + last->type = flavour; + std::next(last)->type = anti(flavour); } template fastjet::PseudoJet PhaseSpacePoint::gen_last_momentum( ParticleMomenta const & other_momenta, const double mass_square, const double y ) const { std::array pt{0.,0.}; for (auto const & p: other_momenta) { pt[0]-= p.px(); pt[1]-= p.py(); } const double mperp = std::sqrt(pt[0]*pt[0]+pt[1]*pt[1]+mass_square); const double pz=mperp*std::sinh(y); const double E=mperp*std::cosh(y); return {pt[0], pt[1], pz, E}; } namespace { //! adds a particle to target (in correct rapidity ordering) //! @returns positon of insertion - auto insert_particle(std::vector & target, - Particle && particle + auto insert_particle(std::vector & target, + HEJ::Particle && particle ){ const auto pos = std::upper_bound( - begin(target),end(target),particle,rapidity_less{} + begin(target),end(target),particle,HEJ::rapidity_less{} ); target.insert(pos, std::move(particle)); return pos; } } PhaseSpacePoint::PhaseSpacePoint( Process const & proc, JetParameters const & jet_param, - PDF & pdf, double E_beam, + HEJ::PDF & pdf, double E_beam, double const subl_chance, - unsigned int const subl_channels, + Subleading subl_channels, ParticlesDecayMap const & particle_decays, - EWConstants const & ew_parameters, - RNG & ran + HEJ::EWConstants const & ew_parameters, + HEJ::RNG & ran ){ assert(proc.njets >= 2); - status_ = good; + status_ = Status::good; weight_ = 1; + // ensure that all setting are consistent + if(subl_chance == 0.) + subl_channels.reset(); + const int nout = proc.njets + (proc.boson?1:0) + proc.boson_decay.size(); outgoing_.reserve(nout); // generate parton momenta const bool is_pure_jets = (nout == proc.njets); auto partons = gen_LO_partons( proc.njets, is_pure_jets, jet_param, E_beam, ran ); // pre fill flavour with gluons for(auto&& p_out: partons) { - outgoing_.emplace_back(Particle{pid::gluon, std::move(p_out), {}}); + outgoing_.emplace_back(HEJ::Particle{HEJ::pid::gluon, std::move(p_out), {}}); } - if(status_ != good) return; + if(status_ != Status::good) return; if(proc.boson){ // decay boson auto const & boson_prop = ew_parameters.prop(*proc.boson) ; auto boson{ gen_boson(*proc.boson, boson_prop.mass, boson_prop.width, ran) }; const auto pos{insert_particle(outgoing_, std::move(boson))}; const size_t boson_idx = std::distance(begin(outgoing_), pos); auto const & boson_decay = particle_decays.find(*proc.boson); if( !proc.boson_decay.empty() ){ // decay given in proc decays_.emplace( boson_idx, decay_boson(outgoing_[boson_idx], proc.boson_decay, ran) ); } else if( boson_decay != particle_decays.end() && !boson_decay->second.empty() ){ // decay given explicitly decays_.emplace( boson_idx, decay_boson(outgoing_[boson_idx], boson_decay->second, ran) ); } } // normalisation of momentum-conserving delta function weight_ *= std::pow(2*M_PI, 4); /** @TODO * uf (jet_param.min_pt) doesn't correspond to our final scale choice. * The HEJ scale generators currently expect a full event as input, * so fixing this is not completely trivial */ - reconstruct_incoming(proc, subl_channels, pdf, E_beam, jet_param.min_pt, ran); - if(status_ != good) return; + reconstruct_incoming(proc, subl_chance, subl_channels, pdf, E_beam, jet_param.min_pt, ran); + if(status_ != Status::good) return; // set outgoing states - most_backward_FKL(outgoing_).type = incoming_[0].type; - most_forward_FKL(outgoing_).type = incoming_[1].type; + begin_partons()->type = incoming_[0].type; + rbegin_partons()->type = incoming_[1].type; maybe_turn_to_subl(subl_chance, subl_channels, proc, ran); if(proc.boson) couple_boson(*proc.boson, ran); } // pt generation, see eq:pt_sampling in developer manual double PhaseSpacePoint::gen_hard_pt( const int np , const double ptmin, const double ptmax, const double /* y */, - RNG & ran - ) { + HEJ::RNG & ran + ){ // heuristic parameter for pt sampling, see eq:pt_par in developer manual const double ptpar = ptmin + np/5.; const double arctan = std::atan((ptmax - ptmin)/ptpar); const double xpt = ran.flat(); const double pt = ptmin + ptpar*std::tan(xpt*arctan); const double cosine = std::cos(xpt*arctan); weight_ *= pt*ptpar*arctan/(cosine*cosine); return pt; } - double PhaseSpacePoint::gen_soft_pt(int np, double max_pt, RNG & ran) { + double PhaseSpacePoint::gen_soft_pt(int np, double max_pt, HEJ::RNG & ran) { constexpr double ptpar = 4.; const double r = ran.flat(); const double pt = max_pt + ptpar/np*std::log(r); weight_ *= pt*ptpar/(np*r); return pt; } double PhaseSpacePoint::gen_parton_pt( int count, JetParameters const & jet_param, double max_pt, double y, - RNG & ran + HEJ::RNG & ran ) { constexpr double p_small_pt = 0.02; if(! jet_param.peak_pt) { return gen_hard_pt(count, jet_param.min_pt, max_pt, y, ran); } const double r = ran.flat(); if(r > p_small_pt) { weight_ /= 1. - p_small_pt; return gen_hard_pt(count, *jet_param.peak_pt, max_pt, y, ran); } weight_ /= p_small_pt; const double pt = gen_soft_pt(count, *jet_param.peak_pt, ran); if(pt < jet_param.min_pt) { weight_=0.0; - status_ = not_enough_jets; + status_ = Status::not_enough_jets; return jet_param.min_pt; } return pt; } std::vector PhaseSpacePoint::gen_LO_partons( int np, bool is_pure_jets, JetParameters const & jet_param, double max_pt, - RNG & ran + HEJ::RNG & ran ){ if (np<2) throw std::invalid_argument{"Not enough partons in gen_LO_partons"}; weight_ /= std::pow(16.*std::pow(M_PI,3),np); weight_ /= std::tgamma(np+1); //remove rapidity ordering std::vector partons; partons.reserve(np); for(int i = 0; i < np; ++i){ const double y = -jet_param.max_y + 2*jet_param.max_y*ran.flat(); weight_ *= 2*jet_param.max_y; const bool is_last_parton = i+1 == np; if(is_pure_jets && is_last_parton) { constexpr double parton_mass_sq = 0.; partons.emplace_back(gen_last_momentum(partons, parton_mass_sq, y)); break; } const double phi = 2*M_PI*ran.flat(); weight_ *= 2.0*M_PI; const double pt = gen_parton_pt(np, jet_param, max_pt, y, ran); if(weight_ == 0.0) return {}; partons.emplace_back(fastjet::PtYPhiM(pt, y, phi)); assert(jet_param.min_pt <= partons[i].pt()); assert(partons[i].pt() <= max_pt+1e-5); } // Need to check that at LO, the number of jets = number of partons; fastjet::ClusterSequence cs(partons, jet_param.def); auto cluster_jets=cs.inclusive_jets(jet_param.min_pt); if (cluster_jets.size()!=unsigned(np)){ weight_=0.0; - status_ = not_enough_jets; + status_ = Status::not_enough_jets; return {}; } - std::sort(begin(partons), end(partons), rapidity_less{}); + std::sort(begin(partons), end(partons), HEJ::rapidity_less{}); return partons; } - Particle PhaseSpacePoint::gen_boson( - ParticleID bosonid, double mass, double width, - RNG & ran + HEJ::Particle PhaseSpacePoint::gen_boson( + HEJ::ParticleID bosonid, double mass, double width, + HEJ::RNG & ran ){ // Usual phase space measure weight_ /= 16.*std::pow(M_PI, 3); // Generate a y Gaussian distributed around 0 /// @TODO check magic numbers for different boson Higgs /// @TODO better sampling for W const double stddev_y = 1.6; const double y = random_normal(stddev_y, ran); const double r1 = ran.flat(); const double s_boson = mass*( mass + width*std::tan(M_PI/2.*r1 + (r1-1.)*std::atan(mass/width)) ); // off-shell s_boson sampling, compensates for Breit-Wigner /// @TODO use a flag instead - if(std::abs(bosonid) == pid::Wp){ + if(std::abs(bosonid) == HEJ::pid::Wp){ weight_/=M_PI*M_PI*16.; //Corrects B-W factors, see git issue 132 weight_*= mass*width*( M_PI+2.*std::atan(mass/width) ) / ( 1. + std::cos( M_PI*r1 + 2.*(r1-1.)*std::atan(mass/width) ) ); } auto p = gen_last_momentum(outgoing_, s_boson, y); - return Particle{bosonid, std::move(p), {}}; - } - - Particle const & PhaseSpacePoint::most_backward_FKL( - std::vector const & partons - ) const{ - if(!is_parton(partons[0])) return partons[1]; - return partons[0]; - } - - Particle const & PhaseSpacePoint::most_forward_FKL( - std::vector const & partons - ) const{ - const size_t last_idx = partons.size() - 1; - if(!is_parton(partons[last_idx])) return partons[last_idx-1]; - return partons[last_idx]; - } - - Particle & PhaseSpacePoint::most_backward_FKL( - std::vector & partons - ) const{ - if(!is_parton(partons[0])) return partons[1]; - return partons[0]; - } - - Particle & PhaseSpacePoint::most_forward_FKL( - std::vector & partons - ) const{ - const size_t last_idx = partons.size() - 1; - if(!is_parton(partons[last_idx])) return partons[last_idx-1]; - return partons[last_idx]; + return HEJ::Particle{bosonid, std::move(p), {}}; } namespace { /// partons are ordered: even = anti, 0 = gluon - ParticleID index_to_pid(size_t i){ - if(!i) return pid::gluon; - return static_cast( i%2 ? (i+1)/2 : -i/2 ); + HEJ::ParticleID index_to_pid(size_t i){ + if(!i) return HEJ::pid::gluon; + return static_cast( i%2 ? (i+1)/2 : -i/2 ); } /// partons are ordered: even = anti, 0 = gluon - size_t pid_to_index(ParticleID id){ - if(id==pid::gluon) return 0; + size_t pid_to_index(HEJ::ParticleID id){ + if(id==HEJ::pid::gluon) return 0; return id>0 ? id*2-1 : std::abs(id)*2; } - std::bitset<11> init_allowed(ParticleID const id){ - if(std::abs(id) == pid::proton) + using part_mask = std::bitset<11>; //!< Selection mask for partons + + part_mask init_allowed(HEJ::ParticleID const id){ + if(std::abs(id) == HEJ::pid::proton) return ~0; - std::bitset<11> out{0}; - if(is_parton(id)) + part_mask out{0}; + if(HEJ::is_parton(id)) out[pid_to_index(id)] = 1; return out; } /// decides which "index" (see index_to_pid) are allowed for process - std::bitset<11> allowed_quarks(ParticleID const boson){ - std::bitset<11> allowed = ~0; - if(std::abs(boson) == pid::Wp){ - // special case W: - // Wp: anti-down or up-type quark, no b/t - // Wm: down or anti-up-type quark, no b/t - allowed = boson>0? 0b00011001101 - :0b00100110011; + part_mask allowed_quarks(HEJ::ParticleID const boson){ + if(std::abs(boson) != HEJ::pid::Wp){ + return ~1; // not a gluon + } + // special case W: + // Wp: anti-down or up-type quark, no b/t + // Wm: down or anti-up-type quark, no b/t + return boson>0?0b00011001100 + :0b00100110010; + } + } + + std::array PhaseSpacePoint::incoming_AWZ( + Process const & proc, Subleading const subl_channels, + std::array allowed_partons, + HEJ::RNG & ran + ){ + assert(proc.boson); + auto couple_parton = allowed_quarks(*proc.boson); + if( // coupling possible through input + allowed_partons[0] == (couple_parton&allowed_partons[0]) + || allowed_partons[1] == (couple_parton&allowed_partons[1]) + // cqqx possible + || (proc.njets >= 4 && subl_channels[subleading::cqqx]) + ){ + return allowed_partons; + } + // eqqx only possible if one incoming is a gluon + if(proc.njets >= 3 && subl_channels[subleading::eqqx]){ + couple_parton.set(pid_to_index(HEJ::ParticleID::gluon)); + } + // only first can couple + if( (allowed_partons[0]&couple_parton).any() + &&(allowed_partons[1]&couple_parton).none() + ){ + return {allowed_partons[0]&couple_parton, allowed_partons[1]}; + } + // only second can couple + if( (allowed_partons[0]&couple_parton).none() + && (allowed_partons[1]&couple_parton).any() + ){ + return {allowed_partons[0], allowed_partons[1]&couple_parton}; + } + // both can couple + if( (allowed_partons[0]&couple_parton).any() + && (allowed_partons[1]&couple_parton).any() + ){ + double rnd = ran.flat(); + weight_*=3.; + if(rnd<1./3.){ + return { + allowed_partons[0] & couple_parton, + allowed_partons[1] & ~couple_parton + }; + } + if(rnd<2./3.){ + return { + allowed_partons[0] & ~couple_parton, + allowed_partons[1] & couple_parton + }; + } + return { + allowed_partons[0] & couple_parton, + allowed_partons[1] & couple_parton + }; + } + throw std::invalid_argument{"Incoming state not allowed."}; + } + + std::array PhaseSpacePoint::incoming_eqqx( + std::array allowed_partons, HEJ::RNG & ran + ){ + auto const gluon_idx = pid_to_index(HEJ::pid::gluon); + auto & first_beam = allowed_partons[0]; + auto & second_beam = allowed_partons[1]; + if(first_beam[gluon_idx] && !second_beam[gluon_idx]){ + first_beam.reset(); + first_beam.set(gluon_idx); + return allowed_partons; + } + if(!first_beam[gluon_idx] && second_beam[gluon_idx]) { + second_beam.reset(); + second_beam.set(gluon_idx); + return allowed_partons; + } + if(first_beam[gluon_idx] && second_beam[gluon_idx]) { + // both beams can be gluons + // if one can't be a quark everything is good + auto first_quarks = first_beam; + first_quarks.reset(gluon_idx); + auto second_quarks = second_beam; + second_quarks.reset(gluon_idx); + if(first_quarks.none() || second_quarks.none()){ + return allowed_partons; + } + // else choose one to be a gluon + double rnd = ran.flat(); + weight_*=3.; + if(rnd<1./3.){ + allowed_partons[0].reset(); + allowed_partons[0].set(gluon_idx); + allowed_partons[1].reset(gluon_idx); + } else if(rnd<2./3.){ + allowed_partons[1].reset(); + allowed_partons[1].set(gluon_idx); + allowed_partons[0].reset(gluon_idx); + } else { + allowed_partons[0].reset(); + allowed_partons[0].set(gluon_idx); + allowed_partons[1].reset(); + allowed_partons[1].set(gluon_idx); } - return allowed; + return allowed_partons; } + throw std::invalid_argument{ + "Incoming state not allowed for pure extremal qqx."}; + } + + std::array PhaseSpacePoint::incoming_uno( + std::array allowed_partons, HEJ::RNG & ran + ){ + auto const gluon_idx = pid_to_index(HEJ::pid::gluon); + auto & first_beam = allowed_partons[0]; + auto first_quarks = first_beam; + first_quarks.reset(gluon_idx); + auto & second_beam = allowed_partons[1]; + auto second_quarks = second_beam; + second_quarks.reset(gluon_idx); + if(first_quarks.any() && second_quarks.none()){ + first_beam.reset(gluon_idx); + return allowed_partons; + } + if(first_quarks.none() && second_quarks.any()) { + second_beam.reset(gluon_idx); + return allowed_partons; + } + if(first_quarks.any() && second_quarks.any()) { + // both beams can be quarks + // if one can't be gluon everything is good + if(!first_beam[gluon_idx] || !second_beam[gluon_idx]){ + return allowed_partons; + } + // else choose one to be a quark + double rnd = ran.flat(); + weight_*=3.; + if(rnd<1./3.){ + allowed_partons[0].reset(gluon_idx); + allowed_partons[1].reset(); + allowed_partons[1].set(gluon_idx); + } else if(rnd<2./3.){ + allowed_partons[1].reset(gluon_idx); + allowed_partons[0].reset(); + allowed_partons[0].set(gluon_idx); + } else { + allowed_partons[0].reset(gluon_idx); + allowed_partons[1].reset(gluon_idx); + } + return allowed_partons; + } + throw std::invalid_argument{ + "Incoming state not allowed for pure unordered."}; } /** + * @brief Returns list of all allowed initial states partons + * * checks which partons are allowed as initial state: * 1. only allow what is given in the Runcard (p -> all) * 2. A/W/Z require something to couple to * a) no qqx => no incoming gluon * b) 2j => no incoming gluon - * c) 3j => can couple OR is gluon => 2 gluons become qqx later + * c) >3j => can couple OR is gluon => 2 gluons become qqx later + * 3. pure eqqx requires at least one gluon + * 4. pure uno requires at least one quark */ - std::array,2> PhaseSpacePoint::filter_partons( - Process const & proc, unsigned int const subl_channels, RNG & ran + std::array PhaseSpacePoint::allowed_incoming( + Process const & proc, + double const subl_chance, Subleading const subl_channels, + HEJ::RNG & ran ){ - std::array,2> allowed_partons{ + // all possible incoming states + std::array allowed_partons{ init_allowed(proc.incoming[0]), init_allowed(proc.incoming[1]) }; - bool const allow_qqx = subl_channels&Subleading::qqx; // special case A/W/Z - if(is_AWZ_proccess(proc) && ((proc.njets < 4) || !allow_qqx)){ - // all possible incoming states - auto allowed(allowed_quarks(*proc.boson)); - if(proc.njets == 2 || !allow_qqx) allowed[0]=0; - - // possible states per leg - std::array,2> const maybe_partons{ - allowed_partons[0]&allowed, allowed_partons[1]&allowed}; - - if(maybe_partons[0].any() && maybe_partons[1].any()){ - // two options to get allowed initial state => choose one at random - const size_t idx = ran.flat() < 0.5; - allowed_partons[idx] = maybe_partons[idx]; - // else choose the possible - } else if(maybe_partons[0].any()) { - allowed_partons[0] = maybe_partons[0]; - } else if(maybe_partons[1].any()) { - allowed_partons[1] = maybe_partons[1]; - } else{ - throw std::invalid_argument{"Incoming state not allowed."}; - } + if(proc.boson && is_AWZ_proccess(proc)){ + allowed_partons = incoming_AWZ(proc, subl_channels, allowed_partons, ran); + } + // special case: pure subleading + if(subl_chance!=1.){ + return allowed_partons; + } + auto other_channels = subl_channels; + // pure eqqx + other_channels.reset(subleading::eqqx); + if(other_channels.none()){ + return incoming_eqqx(allowed_partons, ran); + } + other_channels = subl_channels; + // pure uno + other_channels.reset(subleading::uno); + if(other_channels.none()){ + return incoming_uno(allowed_partons, ran); } return allowed_partons; } + void PhaseSpacePoint::reconstruct_incoming( - Process const & proc, unsigned int const subl_channels, - PDF & pdf, double E_beam, + Process const & proc, + double const subl_chance, Subleading const subl_channels, + HEJ::PDF & pdf, double E_beam, double uf, - RNG & ran + HEJ::RNG & ran ){ std::tie(incoming_[0].p, incoming_[1].p) = incoming_momenta(outgoing_); // calculate xa, xb const double sqrts=2*E_beam; const double xa=(incoming_[0].E()-incoming_[0].pz())/sqrts; const double xb=(incoming_[1].E()+incoming_[1].pz())/sqrts; // abort if phase space point is outside of collider energy reach if (xa>1. || xb>1.){ weight_=0; - status_ = too_much_energy; + status_ = Status::too_much_energy; return; } auto const & ids = proc.incoming; - std::array,2> allowed_partons( - filter_partons(proc, subl_channels, ran)); + std::array allowed_partons + = allowed_incoming(proc, subl_chance, subl_channels, ran); for(size_t i = 0; i < 2; ++i){ - if(ids[i] == pid::proton || ids[i] == pid::p_bar){ + if(ids[i] == HEJ::pid::proton || ids[i] == HEJ::pid::p_bar){ // pick ids according to pdfs incoming_[i].type = generate_incoming_id(i, i?xb:xa, uf, pdf, allowed_partons[i], ran); } else { assert(allowed_partons[i][pid_to_index(ids[i])]); incoming_[i].type = ids[i]; } } assert(momentum_conserved(1e-7)); } - ParticleID PhaseSpacePoint::generate_incoming_id( + HEJ::ParticleID PhaseSpacePoint::generate_incoming_id( size_t const beam_idx, double const x, double const uf, - PDF & pdf, std::bitset<11> allowed_partons, RNG & ran + HEJ::PDF & pdf, part_mask allowed_partons, HEJ::RNG & ran ){ std::array pdf_wt; pdf_wt[0] = allowed_partons[0]? - std::fabs(pdf.pdfpt(beam_idx,x,uf,pid::gluon)):0.; + std::fabs(pdf.pdfpt(beam_idx,x,uf,HEJ::pid::gluon)):0.; double pdftot = pdf_wt[0]; for(size_t i = 1; i < pdf_wt.size(); ++i){ pdf_wt[i] = allowed_partons[i]? 4./9.*std::fabs(pdf.pdfpt(beam_idx,x,uf,index_to_pid(i))):0; pdftot += pdf_wt[i]; } const double r1 = pdftot * ran.flat(); double sum = 0; for(size_t i=0; i < pdf_wt.size(); ++i){ if (r1 < (sum+=pdf_wt[i])){ weight_*= pdftot/pdf_wt[i]; return index_to_pid(i); } } std::cerr << "Error in choosing incoming parton: "< allowed_parts; - for(auto & part: outgoing_){ + std::vector allowed_parts; + for(auto part_it=begin_partons(); part_it!=end_partons(); ++part_it){ // Wp -> up OR anti-down, Wm -> anti-up OR down, no bottom - if ( can_couple_to_W(part, boson) ) - allowed_parts.push_back(&part); + if ( can_couple_to_W(*part_it, boson) ) + allowed_parts.push_back(part_it); } if(allowed_parts.size() == 0){ throw std::logic_error{"Found no parton for coupling with boson"}; } // select one and flip it size_t idx = 0; if(allowed_parts.size() > 1){ /// @TODO more efficient sampling /// old code: probability[i] = exp(parton[i].y - W.y) idx = std::floor(ran.flat()*allowed_parts.size()); weight_ *= allowed_parts.size(); } const int W_charge = boson>0?1:-1; allowed_parts[idx]->type = - static_cast( allowed_parts[idx]->type - W_charge ); + static_cast( allowed_parts[idx]->type - W_charge ); } - double PhaseSpacePoint::random_normal( - double stddev, - RNG & ran - ){ + double PhaseSpacePoint::random_normal( double stddev, HEJ::RNG & ran ){ const double r1 = ran.flat(); const double r2 = ran.flat(); const double lninvr1 = -std::log(r1); const double result = stddev*std::sqrt(2.*lninvr1)*std::cos(2.*M_PI*r2); weight_ *= exp(result*result/(2*stddev*stddev))*std::sqrt(2.*M_PI)*stddev; return result; } bool PhaseSpacePoint::momentum_conserved(double ep) const{ fastjet::PseudoJet diff; for(auto const & in: incoming()) diff += in.p; for(auto const & out: outgoing()) diff -= out.p; - return nearby_ep(diff, fastjet::PseudoJet{}, ep); + return HEJ::nearby_ep(diff, fastjet::PseudoJet{}, ep); } Decay PhaseSpacePoint::select_decay_channel( std::vector const & decays, - RNG & ran + HEJ::RNG & ran ){ double br_total = 0.; for(auto const & decay: decays) br_total += decay.branching_ratio; // adjust weight // this is given by (channel branching ratio)/(chance to pick channel) // where (chance to pick channel) = // (channel branching ratio)/(total branching ratio) weight_ *= br_total; if(decays.size()==1) return decays.front(); const double r1 = br_total*ran.flat(); double br_sum = 0.; for(auto const & decay: decays){ br_sum += decay.branching_ratio; if(r1 < br_sum) return decay; } throw std::logic_error{"unreachable"}; } - std::vector PhaseSpacePoint::decay_boson( - Particle const & parent, + std::vector PhaseSpacePoint::decay_boson( + HEJ::Particle const & parent, std::vector const & decays, - RNG & ran + HEJ::RNG & ran ){ const auto channel = select_decay_channel(decays, ran); return decay_boson(parent, channel.products, ran); } - std::vector PhaseSpacePoint::decay_boson( - Particle const & parent, - std::vector const & decays, - RNG & ran + std::vector PhaseSpacePoint::decay_boson( + HEJ::Particle const & parent, + std::vector const & decays, + HEJ::RNG & ran ){ if(decays.size() != 2){ - throw not_implemented{ + throw HEJ::not_implemented{ "only decays into two particles are implemented" }; } - std::vector decay_products(decays.size()); + std::vector decay_products(decays.size()); for(size_t i = 0; i < decays.size(); ++i){ decay_products[i].type = decays[i]; } // choose polar and azimuth angle in parent rest frame const double E = parent.m()/2; const double theta = 2.*M_PI*ran.flat(); const double cos_phi = 2.*ran.flat()-1.; // Jacobian Factors for W in line 418 const double sin_phi = std::sqrt(1. - cos_phi*cos_phi); // Know 0 < phi < pi const double px = E*std::cos(theta)*sin_phi; const double py = E*std::sin(theta)*sin_phi; const double pz = E*cos_phi; decay_products[0].p.reset(px, py, pz, E); decay_products[1].p.reset(-px, -py, -pz, E); for(auto & particle: decay_products) particle.p.boost(parent.p); return decay_products; } } diff --git a/FixedOrderGen/src/main.cc b/FixedOrderGen/src/bin/main.cc similarity index 97% rename from FixedOrderGen/src/main.cc rename to FixedOrderGen/src/bin/main.cc index c9b7fec..579428f 100644 --- a/FixedOrderGen/src/main.cc +++ b/FixedOrderGen/src/bin/main.cc @@ -1,293 +1,293 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include "fastjet/ClusterSequence.hh" #include "fastjet/PseudoJet.hh" #include "HEJ/Analysis.hh" #include "HEJ/CombinedEventWriter.hh" #include "HEJ/CrossSectionAccumulator.hh" #include "HEJ/Event.hh" #include "HEJ/exceptions.hh" #include "HEJ/get_analysis.hh" #include "HEJ/make_RNG.hh" #include "HEJ/ProgressBar.hh" #include "HEJ/ScaleFunction.hh" #include "HEJ/Unweighter.hh" #include "LHEF/LHEF.h" #include "config.hh" #include "EventGenerator.hh" #include "Status.hh" #include "Version.hh" namespace YAML { class Node; } namespace { constexpr auto banner = " __ ___ __ ______ __ " " __ \n / / / (_)___ _/ /_ / ____/___ " " ___ _________ ___ __ / /__ / /______ \n " " / /_/ / / __ `/ __ \\ / __/ / __ \\/ _ \\/ ___/ __ `/ / / / __ / / _" " \\/ __/ ___/ \n / __ / / /_/ / / / / / /___/ /" " / / __/ / / /_/ / /_/ / / /_/ / __/ /_(__ ) " " \n /_/ /_/_/\\__, /_/ /_/ /_____/_/ /_/\\___/_/ \\__, /\\__, / \\___" "_/\\___/\\__/____/ \n ____///__/ " "__ ____ ///__//____/ ______ __ " " \n / ____(_) _____ ____/ / / __ \\_________/ /__ _____ / " "____/__ ____ ___ _________ _/ /_____ _____\n / /_ / / |/_/ _ \\/ __" " / / / / / ___/ __ / _ \\/ ___/ / / __/ _ \\/ __ \\/ _ \\/ ___/ __ `/ " "__/ __ \\/ ___/\n / __/ / /> > get_analyses( std::vector const & parameters, LHEF::HEPRUP const & heprup ){ try{ return HEJ::get_analyses(parameters, heprup); } catch(std::exception const & exc){ std::cerr << "Failed to load analysis: " << exc.what() << '\n'; std::exit(EXIT_FAILURE); } } template auto make_lowpt_filter(Iterator begin, Iterator end, HEJ::optional peak_pt){ return boost::make_filter_iterator( [peak_pt](HEJ::Event const & ev){ assert(! ev.jets().empty()); double min_pt = peak_pt?(*peak_pt):0.; const auto softest_jet = fastjet::sorted_by_pt(ev.jets()).back(); return softest_jet.pt() > min_pt; }, begin, end ); } int main(int argn, char** argv) { using namespace std::string_literals; if (argn < 2) { std::cerr << "\n# Usage:\n." << argv[0] << " config_file\n"; return EXIT_FAILURE; } std::cout << banner; std::cout << "Version " << HEJFOG::Version::String() << ", revision " << HEJFOG::Version::revision() << std::endl; fastjet::ClusterSequence::print_banner(); using clock = std::chrono::system_clock; const auto start_time = clock::now(); // read configuration auto config = load_config(argv[1]); std::shared_ptr ran{ HEJ::make_RNG(config.rng.name, config.rng.seed)}; assert(ran != nullptr); HEJ::ScaleGenerator scale_gen{ config.scales.base, config.scales.factors, config.scales.max_ratio }; HEJFOG::EventGenerator generator{ config.process, config.beam, std::move(scale_gen), config.jets, config.pdf_id, config.subleading_fraction, config.subleading_channels, config.particle_decays, config.Higgs_coupling, config.ew_parameters, ran }; // prepare process information for output LHEF::HEPRUP heprup; heprup.IDBMUP=std::pair(config.beam.particles[0], config.beam.particles[1]); heprup.EBMUP=std::make_pair(config.beam.energy, config.beam.energy); heprup.PDFGUP=std::make_pair(0,0); heprup.PDFSUP=std::make_pair(config.pdf_id,config.pdf_id); heprup.NPRUP=1; heprup.XSECUP=std::vector(1.); heprup.XERRUP=std::vector(1.); heprup.LPRUP=std::vector{1}; heprup.generators.emplace_back(LHEF::XMLTag{}); heprup.generators.back().name = HEJFOG::Version::package_name(); heprup.generators.back().version = HEJFOG::Version::String(); HEJ::CombinedEventWriter writer{config.output, heprup}; std::vector> analyses = get_analyses( config.analyses_parameters, heprup ); assert(analyses.empty() || analyses.front() != nullptr); // warm-up phase to train unweighter HEJ::optional unweighter{}; std::map status_counter; std::vector events; std::uint64_t trials = 0; if(config.unweight) { std::cout << "Calibrating unweighting ...\n"; const auto warmup_start = clock::now(); const std::size_t warmup_events = config.unweight->sample_size; HEJ::ProgressBar warmup_progress{std::cout, warmup_events}; for(; events.size() < warmup_events; ++trials){ auto ev = generator.gen_event(); ++status_counter[generator.status()]; - assert( (generator.status() == HEJFOG::good) == bool(ev) ); - if(generator.status() != HEJFOG::good) continue; + assert( (generator.status() == HEJFOG::Status::good) == bool(ev) ); + if(generator.status() != HEJFOG::Status::good) continue; const bool pass_cuts = analyses.empty() || std::any_of( begin(analyses), end(analyses), [&ev](auto const & analysis) { return analysis->pass_cuts(*ev, *ev); } ); if(pass_cuts) { events.emplace_back(std::move(*ev)); ++warmup_progress; } } std::cout << std::endl; unweighter = HEJ::Unweighter(); unweighter->set_cut_to_peakwt( make_lowpt_filter(events.cbegin(), events.cend(), config.jets.peak_pt), make_lowpt_filter(events.cend(), events.cend(), config.jets.peak_pt), config.unweight->max_dev ); std::vector unweighted_events; for(auto && ev: events) { auto unweighted = unweighter->unweight(std::move(ev), *ran); if(unweighted) { unweighted_events.emplace_back(std::move(*unweighted)); } } events = std::move(unweighted_events); if(events.empty()) { std::cerr << "Failed to generate events. Please increase \"unweight: sample size\"" " or reduce \"unweight: max deviation\"\n"; return EXIT_FAILURE; } const auto warmup_end = clock::now(); const double completion = static_cast(events.size())/config.events; const std::chrono::duration remaining_time = (warmup_end- warmup_start)*(1./completion - 1); const auto finish = clock::to_time_t( std::chrono::time_point_cast(warmup_end + remaining_time) ); std::cout << "Generated " << events.size() << "/" << config.events << " events (" << static_cast(std::round(100*completion)) << "%)\n" << "Estimated remaining generation time: " << remaining_time.count() << " seconds (" << std::put_time(std::localtime(&finish), "%c") << ")\n\n"; } // end unweighting warm-up // main generation loop // event weight is wrong, need to divide by "total number of trials" afterwards HEJ::ProgressBar progress{std::cout, config.events}; progress.increment(events.size()); events.reserve(config.events); for(; events.size() < config.events; ++trials){ auto ev = generator.gen_event(); ++status_counter[generator.status()]; - assert( (generator.status() == HEJFOG::good) == bool(ev) ); - if(generator.status() != HEJFOG::good) continue; + assert( (generator.status() == HEJFOG::Status::good) == bool(ev) ); + if(generator.status() != HEJFOG::Status::good) continue; const bool pass_cuts = analyses.empty() || std::any_of( begin(analyses), end(analyses), [&ev](auto const & analysis) { return analysis->pass_cuts(*ev, *ev); } ); if(pass_cuts) { if(unweighter) { auto unweighted = unweighter->unweight(std::move(*ev), *ran); if(! unweighted) continue; ev = std::move(unweighted); } events.emplace_back(std::move(*ev)); ++progress; } } std::cout << std::endl; // final run though events with correct weight HEJ::CrossSectionAccumulator xs; for(auto & ev: events){ ev.parameters() *= invGeV2_to_pb/trials; for(auto const & analysis: analyses) { if(analysis->pass_cuts(ev, ev)) { analysis->fill(ev, ev); } } writer.write(ev); xs.fill(ev); } for(auto const & analysis: analyses) { analysis->finalise(); } writer.finish(); // Print final informations const std::chrono::duration run_time = (clock::now() - start_time); std::cout << "\nTask Runtime: " << run_time.count() << " seconds for " << events.size() << " Events (" << events.size()/run_time.count() << " evts/s)\n" << std::endl; std::cout << xs << "\n"; for(auto && entry: status_counter){ const double fraction = static_cast(entry.second)/trials; const int percent = std::round(100*fraction); std::cout << "status " << std::left << std::setw(16) << (to_string(entry.first) + ":") << " ["; for(int i = 0; i < percent/2; ++i) std::cout << '#'; for(int i = percent/2; i < 50; ++i) std::cout << ' '; std::cout << "] " << percent << "%" << std::endl; } return EXIT_SUCCESS; } diff --git a/FixedOrderGen/src/config.cc b/FixedOrderGen/src/config.cc index 4d869b4..91f75cb 100644 --- a/FixedOrderGen/src/config.cc +++ b/FixedOrderGen/src/config.cc @@ -1,439 +1,440 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #include "config.hh" #include #include #include #include #include #include #include #include "HEJ/exceptions.hh" #include "HEJ/PDG_codes.hh" #include "HEJ/YAMLreader.hh" -#include "Subleading.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", "vev", "event output", "analyses", "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", "W", "Z"}){ for(auto && particle_opt: {"mass", "width"}){ supported["particle properties"][particle_type][particle_opt] = ""; } } for(auto && particle_type: {"Higgs", "Wp", "W+", "Wm", "W-", "Z"}){ for(auto && particle_opt: {"into", "branching ratio"}){ supported["decays"][particle_type][particle_opt] = ""; } } 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 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 split( std::string const & str, std::string const & delims ){ std::vector 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 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 "+name(ids[0])+" and "+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 namespace HEJFOG::channels; + using namespace HEJFOG::subleading; set_from_yaml(name, yaml); if(name == "none") - return none; + return NONE; if(name == "all") - return all; + return ALL; + Subleading channel = NONE; if(name == "unordered" || name == "uno") - return uno; - if(name == "qqx") - return qqx; + return channel.set(uno); + if(name == "central qqx" || name == "cqqx") + return channel.set(cqqx); + if(name == "extremal qqx" || name == "eqqx") + return channel.set(eqqx); throw HEJ::unknown_option("Unknown subleading channel '"+name+"'"); } - unsigned int get_subleading_channels(YAML::Node const & node){ + Subleading get_subleading_channels(YAML::Node const & node){ using YAML::NodeType; - using namespace HEJFOG::channels; + using namespace HEJFOG::subleading; // all channels allowed by default - if(!node) return all; + if(!node) return ALL; switch(node.Type()){ case NodeType::Undefined: - return all; + return ALL; case NodeType::Null: - return none; + 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; + Subleading channels; 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, std::string const & entry, std::string const & boson ){ Decay decay; set_from_yaml(decay.products, node, entry, boson, "into"); decay.branching_ratio=1; set_from_yaml_if_defined(decay.branching_ratio, node, entry, boson, "branching ratio"); return decay; } std::vector get_decays(YAML::Node const & node, std::string const & entry, std::string const & boson ){ using YAML::NodeType; if(!node[entry][boson]) return {}; switch(node[entry][boson].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, entry, boson)}; case NodeType::Sequence: std::vector result; for(auto && decay_str: node[entry][boson]){ result.emplace_back(get_decay(decay_str, entry, boson)); } return result; } throw std::logic_error{"unreachable"}; } ParticlesDecayMap get_all_decays(YAML::Node const & node, std::string const & entry ){ if(!node[entry]) return {}; if(!node[entry].IsMap()) throw HEJ::invalid_type{entry + " have to be a map"}; ParticlesDecayMap result; for(auto const & sub_node: node[entry]) { const auto boson = sub_node.first.as(); const auto id = HEJ::to_ParticleID(boson); result.emplace(id, get_decays(node, entry, boson)); } return result; } HEJ::ParticleProperties get_particle_properties( YAML::Node const & node, std::string const & entry, std::string const & boson ){ HEJ::ParticleProperties result; set_from_yaml(result.mass, node, entry, boson, "mass"); set_from_yaml(result.width, node, entry, boson, "width"); return result; } HEJ::EWConstants get_ew_parameters(YAML::Node const & node){ HEJ::EWConstants result; double vev; set_from_yaml(vev, node, "vev"); result.set_vevWZH(vev, get_particle_properties(node, "particle properties", "W"), get_particle_properties(node, "particle properties", "Z"), get_particle_properties(node, "particle properties", "Higgs") ); 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){ auto const & 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_channels = Subleading::none; + config.subleading_channels.reset(); else config.subleading_channels = get_subleading_channels(yaml["subleading channels"]); config.ew_parameters = get_ew_parameters(yaml); config.particle_decays = get_all_decays(yaml, "decays"); if(config.process.boson // check that Ws always decay && std::abs(*config.process.boson) == HEJ::ParticleID::Wp && config.process.boson_decay.empty() ){ auto const & decay = config.particle_decays.find(*config.process.boson); if(decay == config.particle_decays.end() || decay->second.empty()) throw std::invalid_argument{ "Decay for "+name(*config.process.boson)+" is required"}; } set_from_yaml_if_defined(config.analyses_parameters, yaml, "analyses"); if(yaml["analysis"]){ std::cerr << "WARNING: Configuration entry 'analysis' is deprecated. " " Use 'analyses' instead.\n"; YAML::Node ana; set_from_yaml(ana, yaml, "analysis"); if(!ana.IsNull()){ config.analyses_parameters.push_back(ana); } } 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/FixedOrderGen/t/CMakeLists.txt b/FixedOrderGen/t/CMakeLists.txt index eb00af1..ef2d76f 100644 --- a/FixedOrderGen/t/CMakeLists.txt +++ b/FixedOrderGen/t/CMakeLists.txt @@ -1,108 +1,205 @@ set(tst_dir "${CMAKE_CURRENT_SOURCE_DIR}") +set(runcard_dir ${tst_dir}/runcards) + +function(add_long_test) + set(options NONE) + set(oneValueArgs NAME) + set(multiValueArgs COMMAND) + cmake_parse_arguments(TEST "${options}" "${oneValueArgs}" "${multiValueArgs}" ${ARGN}) + add_test( + NAME ${TEST_NAME}_short + COMMAND ${TEST_COMMAND} short + ) + if(TEST_ALL) + add_test( + NAME ${TEST_NAME} + COMMAND ${TEST_COMMAND} + ) + endif() +endfunction() + foreach(tst W_reconstruct_enu W_2j_classify W_nj_classify) add_executable(test_${tst} ${tst_dir}/${tst}.cc) target_link_libraries(test_${tst} hejfog_lib) - add_test(NAME ${tst} COMMAND test_${tst} WORKING_DIRECTORY ${tst_dir}) + add_test(NAME ${tst} + COMMAND test_${tst} + WORKING_DIRECTORY ${runcard_dir} + ) endforeach() # this only tests if the runcard actually works, not if the result is correct add_test( NAME main_2j - COMMAND HEJFOG ${tst_dir}/config_2j.yml - ) + COMMAND HEJFOG ${runcard_dir}/2j.yml +) add_test( NAME main_h2j - COMMAND HEJFOG ${tst_dir}/config_h2j.yml - ) + COMMAND HEJFOG ${runcard_dir}/h2j.yml +) add_test( NAME main_h2j_decay - COMMAND HEJFOG ${tst_dir}/config_h2j_decay.yml - ) + COMMAND HEJFOG ${runcard_dir}/h2j_decay.yml +) add_test( NAME peakpt - COMMAND HEJFOG ${tst_dir}/config_2j_peak.yml - ) + COMMAND HEJFOG ${runcard_dir}/2j_peak.yml +) +# special case where coupling not trivial +add_test( + NAME main_W4j_uno+eqqx + COMMAND HEJFOG ${runcard_dir}/Wp4j_uno+eqqx.yml +) # check that uno emission doesn't change FKL xs add_executable(FKL_uno FKL_uno.cc) target_link_libraries(FKL_uno hejfog_lib) add_test( NAME FKL_uno # calculated with HEJ revision 9570e3809613272ac4b8bf3236279ba23cf64d20 - COMMAND FKL_uno ${tst_dir}/config_h3j_uno.yml 0.0243548 0.000119862 + COMMAND FKL_uno ${runcard_dir}/h3j_uno.yml 0.0243548 0.000119862 ) # xs tests add_executable(xs_gen xs_gen.cc) target_link_libraries(xs_gen hejfog_lib) ## Higgs add_test( + NAME xs_hqQ + # calculated with Sherpa see #132 + COMMAND xs_gen ${runcard_dir}/hqQ.yml 1.612e-02 1.26303e-04 +) +add_test( NAME xs_h2j # calculated with HEJ revision 9570e3809613272ac4b8bf3236279ba23cf64d20 - COMMAND xs_gen ${tst_dir}/config_h2j.yml 2.04928 0.00956022 + COMMAND xs_gen ${runcard_dir}/h2j.yml 2.04928 0.00956022 ) add_test( NAME xs_h3j # calculated with HEJ revision bd4388fe55cbc3f5a7b6139096456c551294aa31 - COMMAND xs_gen ${tst_dir}/config_h3j.yml 1.07807 0.0132409 + COMMAND xs_gen ${runcard_dir}/h3j.yml 1.07807 0.0132409 +) +add_long_test( + NAME xs_h3j_uno + # calculated with HEJ revision 9570e3809613272ac4b8bf3236279ba23cf64d20 + COMMAND xs_gen ${runcard_dir}/h3j_uno.yml 0.00347538 3.85875e-05 ) add_test( NAME xs_h5j # calculated with HEJ revision dbde2ffbb3b383ae6709b2424d8f0f9d5658270b - COMMAND xs_gen ${tst_dir}/config_h5j.yml 0.0112504 0.000199633 -) -add_test( - NAME xs_h3j_uno - # calculated with HEJ revision 9570e3809613272ac4b8bf3236279ba23cf64d20 - COMMAND xs_gen ${tst_dir}/config_h3j_uno.yml 0.00347538 3.85875e-05 + COMMAND xs_gen ${runcard_dir}/h5j.yml 0.0112504 0.000199633 ) add_test( NAME xs_h2j_decay # calculated with HEJ revision 9570e3809613272ac4b8bf3236279ba23cf64d20 - COMMAND xs_gen ${tst_dir}/config_h2j_decay.yml 0.00466994 2.20995e-05 + COMMAND xs_gen ${runcard_dir}/h2j_decay.yml 0.00466994 2.20995e-05 ) ## pure jets add_test( + NAME xs_qQ + # calculated with Sherpa see #132 + COMMAND xs_gen ${runcard_dir}/qQ.yml 7.354e+05 1.905e+03 +) +add_test( NAME xs_2j # calculated with "combined" HEJ svn r3480 - COMMAND xs_gen ${tst_dir}/config_2j.yml 86.42031848e06 590570 + COMMAND xs_gen ${runcard_dir}/2j.yml 86.42031848e06 590570 ) add_test( NAME xs_3j_uno - # calculated with HEJ revision b40e0212b144cdd267df19333b574cc7b91847f1 - COMMAND xs_gen ${tst_dir}/config_3j_uno.yml 520697 8948.18 + # calculated with HEJ revision 9401196fba75b5d16bc33f2a309175fecaca00f1 + COMMAND xs_gen ${runcard_dir}/3j_uno.yml 900 14.3131 ) add_test( NAME xs_3j_qqx - # calculated with HEJ revision b40e0212b144cdd267df19333b574cc7b91847f1 - COMMAND xs_gen ${tst_dir}/config_3j_qqx.yml 122233 1613.04 + # calculated with HEJ revision 9401196fba75b5d16bc33f2a309175fecaca00f1 + COMMAND xs_gen ${runcard_dir}/3j_qqx.yml 62040 1005 ) add_test( NAME xs_4j_qqx - # calculated with HEJ revision b40e0212b144cdd267df19333b574cc7b91847f1 - COMMAND xs_gen ${tst_dir}/config_4j_qqx.yml 29805.4 581.473 + # calculated with HEJ revision 9401196fba75b5d16bc33f2a309175fecaca00f1 + COMMAND xs_gen ${runcard_dir}/4j_qqx.yml 28936 550 ) add_test( NAME xs_4j # calculated with HEJ revision 13207b5f67a5f40a2141aa7ee515b022bd4efb65 - COMMAND xs_gen ${tst_dir}/config_4j.yml 915072 15402.4 + COMMAND xs_gen ${runcard_dir}/4j.yml 915072 15402.4 ) ## W add_test( - NAME xs_W2j - # calculated with HEJ revision 987bb30a7985a24a7961f98cfbbc74d3a7992970 - COMMAND xs_gen ${tst_dir}/config_Wm2j.yml 231.404 3.43798 + NAME xs_WqQ + # calculated with Sherpa see #132 + COMMAND xs_gen ${runcard_dir}/WpqQ.yml 3.086e+00 4.511e-02 ) add_test( + NAME xs_W2j + # calculated with HEJ revision 449f2f6b597020e9c9e35699568edc05c827fc11+1 + COMMAND xs_gen ${runcard_dir}/Wm2j.yml 4.177443e+02 7.446928e+00 +) +add_long_test( + NAME xs_W3j_eqqx + # calculated with HEJ revision 449f2f6b597020e9c9e35699568edc05c827fc11+1 + COMMAND xs_gen ${runcard_dir}/Wp3j_qqx.yml 8.710751e+00 1.245725e-01 +) +add_long_test( NAME xs_W3j_uno - # calculated with HEJ revision 987bb30a7985a24a7961f98cfbbc74d3a7992970 - COMMAND xs_gen ${tst_dir}/config_Wp3j_uno.yml 0.139372 0.00248345 + # calculated with HEJ revision 449f2f6b597020e9c9e35699568edc05c827fc11+1 + COMMAND xs_gen ${runcard_dir}/Wp3j_uno.yml 3.000955e-01 5.831799e-03 ) +add_long_test( + NAME xs_W4j_qqx + # calculated with HEJ revision be065dc9a21e5965ce57583f6c0a3d953664b82b + COMMAND xs_gen ${runcard_dir}/Wp4j_qqx.yml 9.274718e-02 4.875742e-03 +) + +# Test that sum of partons == proton +add_executable(PSP_channel PSP_channel.cc) +target_link_libraries(PSP_channel hejfog_lib) +# pure jets add_test( - NAME xs_W3j_eqqx - # calculated with HEJ revision 987bb30a7985a24a7961f98cfbbc74d3a7992970 - COMMAND xs_gen ${tst_dir}/config_Wp3j_qqx.yml 8.45323 0.103449 + NAME channel_2j + COMMAND PSP_channel ${runcard_dir}/2j.yml ) add_test( - NAME xs_W5j_qqx - # calculated with HEJ revision 987bb30a7985a24a7961f98cfbbc74d3a7992970 - COMMAND xs_gen ${tst_dir}/config_Wp4j_qqx.yml 0.0851908 0.00300447 + NAME channel_3j_qqx + COMMAND PSP_channel ${runcard_dir}/3j_qqx.yml +) +add_test( + NAME channel_3j_uno + COMMAND PSP_channel ${runcard_dir}/3j_uno.yml +) +add_long_test( + NAME channel_4j_qqx + COMMAND PSP_channel ${runcard_dir}/4j_qqx.yml +) +# W+jets +# pure jets +add_test( + NAME channel_W2j + COMMAND PSP_channel ${runcard_dir}/Wm2j.yml +) +add_test( + NAME channel_W3j_uno + COMMAND PSP_channel ${runcard_dir}/Wp3j_uno.yml +) +add_test( + NAME channel_W3j_eqqx + COMMAND PSP_channel ${runcard_dir}/Wp3j_qqx.yml +) +add_long_test( + NAME channel_W4j_qqx + COMMAND PSP_channel ${runcard_dir}/Wp4j_qqx.yml +) + +# Test that each subleading channel is generated consistently +add_executable(PSP_subl PSP_subl.cc) +target_link_libraries(PSP_subl hejfog_lib) +add_long_test( + NAME subl_5j + COMMAND PSP_subl ${runcard_dir}/5j.yml +) +add_long_test( + NAME subl_h5j + COMMAND PSP_subl ${runcard_dir}/h5j.yml +) +add_long_test( + NAME subl_W5j + COMMAND PSP_subl ${runcard_dir}/Wm5j.yml ) diff --git a/FixedOrderGen/t/FKL_uno.cc b/FixedOrderGen/t/FKL_uno.cc index 1e766f6..f8449f3 100644 --- a/FixedOrderGen/t/FKL_uno.cc +++ b/FixedOrderGen/t/FKL_uno.cc @@ -1,90 +1,90 @@ /** * check that adding uno emissions doesn't change the FKL cross section * * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #include #include #include #include #include #include "HEJ/Event.hh" #include "HEJ/event_types.hh" #include "HEJ/PDG_codes.hh" #include "HEJ/Ranlux64.hh" #include "HEJ/ScaleFunction.hh" #include "config.hh" #include "EventGenerator.hh" #include "Subleading.hh" #include "Status.hh" //! throw error if condition not fulfilled #define ASSERT(x) if(!(x)) { \ - throw std::logic_error("ASSERTion '" #x "' failed."); \ + throw std::logic_error("Assertion '" #x "' failed."); \ } namespace { constexpr double invGeV2_to_pb = 389379292.; } int main(int argn, char const *argv[]){ using namespace HEJFOG; if(argn != 4){ std::cerr << "Usage: " << argv[0] << " config.yml xs xs_err"; return EXIT_FAILURE; } const double xs_ref = std::stod(argv[2]); const double err_ref = std::stod(argv[3]); (void) err_ref; auto config = load_config(argv[1]); config.subleading_fraction = 0.37; std::shared_ptr ran{std::make_shared()}; HEJFOG::EventGenerator generator{ config.process, config.beam, HEJ::ScaleGenerator{ config.scales.base, config.scales.factors, config.scales.max_ratio }, config.jets, config.pdf_id, config.subleading_fraction, config.subleading_channels, config.particle_decays, config.Higgs_coupling, config.ew_parameters, ran }; double xs = 0., xs_err = 0.; std::size_t uno_found = 0; for (std::size_t trials = 0; trials < config.events; ++trials){ auto ev = generator.gen_event(); - if(generator.status() != good) continue; + if(generator.status() != Status::good) continue; ASSERT(ev); if(ev->type() != HEJ::event_type::FKL){ ++uno_found; continue; } xs += ev->central().weight; xs_err += ev->central().weight*ev->central().weight; } xs *= invGeV2_to_pb/config.events; xs_err = std::sqrt(xs_err); xs_err *= invGeV2_to_pb/config.events; std::cout << xs_ref << " ~ " << xs << " +- " << xs_err << '\n'; std::cout << uno_found << " events with unordered emission" << std::endl; ASSERT(uno_found > 0); ASSERT(std::abs(xs - xs_ref) < 3*xs_err); ASSERT(xs_err < 0.05*xs); return EXIT_SUCCESS; } diff --git a/FixedOrderGen/t/PSP_channel.cc b/FixedOrderGen/t/PSP_channel.cc new file mode 100644 index 0000000..dfde2b0 --- /dev/null +++ b/FixedOrderGen/t/PSP_channel.cc @@ -0,0 +1,222 @@ +/** + * check that the sum of all possible quarks is the same as a proton + * + * \authors The HEJ collaboration (see AUTHORS for details) + * \date 2020 + * \copyright GPLv2 or later + */ +#include +#include +#include +#include + +#include "HEJ/Mixmax.hh" +#include "HEJ/Event.hh" + +#include "config.hh" +#include "EventGenerator.hh" + +//! throw error if condition not fulfilled +#define ASSERT(x) if(!(x)) { \ + throw std::logic_error("Assertion '" #x "' failed."); \ + } + +//! throw error if condition not fulfilled +#define ASSERT_THROW(x, exception) try { \ + x; \ + std::cerr << "'" #x "' did not throw an exception.\n"; \ + throw; \ + } catch(exception const & ex){ \ + std::cout << "throw triggered: " << ex.what() << std::endl; \ + } \ + catch (...) { \ + std::cerr << "Unexpected exception thrown for '" #x "'.\n"; \ + throw; \ + } + +namespace { + + HEJFOG::EventGenerator make_generator( + HEJFOG::Config const & config, + std::shared_ptr ran + ){ + return HEJFOG::EventGenerator{ + config.process, + config.beam, + HEJ::ScaleGenerator{ + config.scales.base, + config.scales.factors, + config.scales.max_ratio + }, + config.jets, + config.pdf_id, + config.subleading_fraction, + config.subleading_channels, + config.particle_decays, + config.Higgs_coupling, + config.ew_parameters, + ran + }; + } + + const static std::array PARTONS{ + HEJ::ParticleID::g, + HEJ::ParticleID::d, HEJ::ParticleID::u, + HEJ::ParticleID::s, HEJ::ParticleID::c, + HEJ::ParticleID::b, + HEJ::ParticleID::d_bar, HEJ::ParticleID::u_bar, + HEJ::ParticleID::s_bar, HEJ::ParticleID::c_bar, + HEJ::ParticleID::b_bar + }; + + bool only_channel( + HEJFOG::Subleading channels, + HEJFOG::subleading::Channels selected + ){ + return channels[selected] && channels.reset(selected).none(); + } + + void generate_till_throw(HEJFOG::EventGenerator & gen ){ + for(std::size_t i=0; i<100; ++i){ + auto ev = gen.gen_event(); + if(gen.status()==HEJFOG::Status::good){ + std::cerr << "Wrongfully generated valid event!" << *ev <<"\n"; + throw; + } + } + std::cerr << "Unable to generate a proper throw!\n"; + throw; + } + + bool can_couple_to_W(HEJ::ParticleID const id, HEJ::ParticleID const boson){ + using namespace HEJ; + if(!is_anyquark(id)){ + return false; + } + if(std::abs(id)==HEJ::ParticleID::b || std::abs(id)==HEJ::ParticleID::top) + return false; + // Wp: + if(boson>0){ + // anti-down + if(id%2==0){ + return is_quark(id); + } + // or up-type quark + return is_antiquark(id); + } + // Wm: + // down + if(id%2==0){ + return is_antiquark(id); + } + // or anti-up-type quark + return is_quark(id); + } +} + +int main(int argc, char const *argv[]) +{ + if(argc != 2 && argc != 3){ + std::cerr << "Usage: " << argv[0] << " config.yaml\n"; + return EXIT_FAILURE; + } + const bool short_only = argc==3; + std::cout <(); + + if(short_only) + config.events/=10.; + + double tot_weight = 0.; + double tot_err = 0.; + config.process.incoming[0] = config.process.incoming[1] = HEJ::ParticleID::proton; + { + HEJFOG::EventGenerator gen = make_generator(config, ran); + for(std::size_t i=0; icentral().weight/config.events; + tot_weight += wgt; + tot_err += wgt*wgt; + } + } + ASSERT(tot_weight!=0.); + } + tot_err = std::sqrt(tot_err); + + + config.events /= PARTONS.size(); + double tot_channels = 0.; + double err_channels = 0.; + for(auto b1: PARTONS){ + for(auto b2: PARTONS){ + double wgt_channel = 0.; + double err_channel = 0.; + config.process.incoming[0] = b1; + config.process.incoming[1] = b2; + HEJFOG::EventGenerator gen = make_generator(config, ran); + std::cout << name(b1) << "+" << name(b2) << " " <3 + && config.subleading_fraction>0. + && config.subleading_channels[HEJFOG::subleading::cqqx]){ + // this will force central qqx + } else { + auto const boson = *config.process.boson; + if(!can_couple_to_W(b1, boson) + && !can_couple_to_W(b2, boson) + ){ + std::cout << "bad " << name(boson) << std::endl; + ASSERT_THROW(generate_till_throw(gen), std::invalid_argument); + continue; + } + } + } + for(std::size_t i=0; icentral().weight/config.events; + wgt_channel += wgt; + err_channel += wgt*wgt; + } + } + ASSERT(wgt_channel!=0.); + tot_channels += wgt_channel; + err_channels += err_channel; + err_channel = std::sqrt(err_channel); + std::cout << "=> " << wgt_channel << " +/- " << err_channel << std::endl; + } + } + err_channels = std::sqrt(err_channels); + std::cout << tot_weight << " +/- " << tot_err + << " vs. " << tot_channels << " +/- " << err_channels << std::endl; + + ASSERT(std::abs(tot_weight - tot_channels) + < 2.*sqrt(err_channels*err_channels+tot_err*tot_err)); + + return EXIT_SUCCESS; +} diff --git a/FixedOrderGen/t/PSP_subl.cc b/FixedOrderGen/t/PSP_subl.cc new file mode 100644 index 0000000..e6de902 --- /dev/null +++ b/FixedOrderGen/t/PSP_subl.cc @@ -0,0 +1,218 @@ +/** + * check that each subprocess xs is correctly generated when only selecting one + * + * \authors The HEJ collaboration (see AUTHORS for details) + * \date 2020 + * \copyright GPLv2 or later + */ +#include +#include +#include +#include + +#include "HEJ/Mixmax.hh" +#include "HEJ/Event.hh" +#include "HEJ/CrossSectionAccumulator.hh" + +#include "config.hh" +#include "EventGenerator.hh" + +//! throw error if condition not fulfilled +#define ASSERT(x) if(!(x)) { \ + throw std::logic_error("Assertion '" #x "' failed."); \ + } + +//! throw error if condition not fulfilled +#define ASSERT_THROW(x, exception) try { \ + x; \ + std::cerr << "'" #x "' did not throw an exception.\n"; \ + throw; \ + } catch(exception const & ex){ \ + std::cout << "throw triggered: " << ex.what() << std::endl; \ + } \ + catch (...) { \ + std::cerr << "Unexpected exception thrown for '" #x "'.\n"; \ + throw; \ + } + +namespace { + const static std::array ALL_TYPES{ + HEJ::event_type::FKL, + HEJ::event_type::unordered_backward, + HEJ::event_type::unordered_forward, + HEJ::event_type::extremal_qqxb, + HEJ::event_type::extremal_qqxf, + HEJ::event_type::central_qqx + }; + + HEJFOG::EventGenerator make_generator( + HEJFOG::Config const & config, + std::shared_ptr ran + ){ + return HEJFOG::EventGenerator{ + config.process, + config.beam, + HEJ::ScaleGenerator{ + config.scales.base, + config.scales.factors, + config.scales.max_ratio + }, + config.jets, + config.pdf_id, + config.subleading_fraction, + config.subleading_channels, + config.particle_decays, + config.Higgs_coupling, + config.ew_parameters, + ran + }; + } + + bool valid_type( + HEJFOG::Subleading const channels, HEJ::event_type::EventType type + ){ + using namespace HEJ::event_type; + using c = HEJFOG::subleading::Channels; + if(channels.none()) + return type == FKL; + if(channels.count()!=1) + return false; + if(channels[c::uno]) + return type == unob || type == unof; + if(channels[c::eqqx]) + return type == qqxexb || type == qqxexf; + if(channels[c::cqqx]) + return type == qqxmid; + throw HEJ::unknown_option{"wrong channel"}; + } + + void generate_till_throw(HEJFOG::EventGenerator & gen ){ + for(std::size_t i=0; i<100; ++i){ + auto ev = gen.gen_event(); + if(gen.status()==HEJFOG::Status::good){ + std::cerr << "Wrongfully generated valid event!" << *ev <<"\n"; + throw; + } + } + std::cerr << "Unable to generate a proper throw!\n"; + throw; + } +} + +int main(int argc, char const *argv[]) +{ + if(argc != 2 && argc != 3){ + std::cerr << "Usage: " << argv[0] << " config.yaml\n"; + return EXIT_FAILURE; + } + const bool short_only = argc==3; + std::cout <(); + + if(short_only) + config.events/=10.; + + // mixes sample + HEJ::CrossSectionAccumulator xs_tot; + config.process.incoming[0] = config.process.incoming[1] = HEJ::ParticleID::proton; + config.subleading_fraction=0.8; + if(config.process.boson && *config.process.boson == HEJ::ParticleID::Higgs){ + config.subleading_channels.reset(); + config.subleading_channels.set(HEJFOG::subleading::uno); + } else { + config.subleading_channels.set(); + } + { + HEJFOG::EventGenerator gen = make_generator(config, ran); + for(std::size_t i=0; icentral().weight /= config.events; + xs_tot.fill(*ev); + } + } + ASSERT(xs_tot.total().value!=0.); + } + + // config.events /= HEJFOG::subleading::last+1; + + // pure FKL + HEJ::CrossSectionAccumulator xs_subl; + { + config.subleading_fraction = 0.; + HEJFOG::EventGenerator gen = make_generator(config, ran); + for(std::size_t i=0; icentral().weight /= config.events; + xs_subl.fill(*ev); + ASSERT(valid_type(0,ev->type())); + } + } + ASSERT(xs_subl.total().value!=0.); + std::cout << "=>\n" << xs_subl << std::endl; + } + + // pure subleading + config.subleading_fraction = 1.; + for(unsigned channel = HEJFOG::subleading::first; + channel <= HEJFOG::subleading::last; ++channel + ){ + config.subleading_channels.reset(); + config.subleading_channels.set(channel); + HEJFOG::EventGenerator gen = make_generator(config, ran); + HEJ::CrossSectionAccumulator xs_channel; + std::cout << config.subleading_channels << " " <type()==HEJ::event_type::FKL) + continue; + ev->central().weight /= config.events; + xs_subl.fill(*ev); + xs_channel.fill(*ev); + ASSERT(valid_type(config.subleading_channels,ev->type())); + } + } + ASSERT(xs_subl.total().value!=0.); + std::cout << "=>\n" << xs_channel << std::endl; + } + std::cout << "Total:\n" << xs_tot << " vs.\n" << xs_subl << std::endl; + + for(auto type: ALL_TYPES){ + double diff = 0.; + double err = 1.; + try { + auto const & tot = xs_tot[type]; + auto const & subl = xs_subl[type]; + diff = tot.value - subl.value; + err = sqrt(tot.error+subl.error); + } catch (std::out_of_range const &){ + std::cout << name(type) << " not set" << std::endl; + ASSERT_THROW(xs_tot[type], std::out_of_range); + ASSERT_THROW(xs_subl[type], std::out_of_range); + continue; + } + if(std::abs(diff) > 2.*err){ + std::cerr << "Large difference in " << name(type) + << " (" << (diff/err) << " sigma)\n"; + return EXIT_FAILURE; + } + } + + return EXIT_SUCCESS; +} diff --git a/FixedOrderGen/t/W_2j_classify.cc b/FixedOrderGen/t/W_2j_classify.cc index 622e288..414cb31 100644 --- a/FixedOrderGen/t/W_2j_classify.cc +++ b/FixedOrderGen/t/W_2j_classify.cc @@ -1,168 +1,165 @@ /** * \brief check that the PSP generates only "valid" W + 2 jets events * * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ -#ifdef NDEBUG -#undef NDEBUG -#endif #include #include #include #include #include "fastjet/JetDefinition.hh" #include "HEJ/EWConstants.hh" #include "HEJ/exceptions.hh" #include "HEJ/Mixmax.hh" #include "HEJ/Particle.hh" #include "HEJ/PDF.hh" #include "HEJ/PDG_codes.hh" #include "Decay.hh" #include "JetParameters.hh" #include "PhaseSpacePoint.hh" #include "Process.hh" #include "Status.hh" #include "Subleading.hh" namespace { using namespace HEJFOG; using namespace HEJ; void print_psp(PhaseSpacePoint const & psp){ std::cerr << "Process:\n" << psp.incoming()[0].type << " + "<< psp.incoming()[1].type << " -> "; for(auto const & out: psp.outgoing()){ std::cerr << out.type << " "; } std::cerr << "\n"; } void bail_out(PhaseSpacePoint const & psp, std::string msg){ print_psp(psp); throw std::logic_error{msg}; } bool is_up_type(Particle const & part){ return HEJ::is_anyquark(part) && !(std::abs(part.type)%2); } bool is_down_type(Particle const & part){ return HEJ::is_anyquark(part) && std::abs(part.type)%2; } bool check_W2j(PhaseSpacePoint const & psp, ParticleID const W_type){ bool found_quark = false; bool found_anti = false; std::vector out_partons; std::vector Wp; for(auto const & p: psp.outgoing()){ if(p.type == W_type) Wp.push_back(p); else if(is_parton(p)) out_partons.push_back(p); else bail_out(psp, "Found particle with is not " +std::to_string(int(W_type))+" or parton"); } if(Wp.size() != 1 || out_partons.size() != 2){ bail_out(psp, "Found wrong number of outgoing partons"); } for(std::size_t j=0; j<2; ++j){ auto const & in = psp.incoming()[j]; auto const & out = out_partons[j]; if(is_quark(in) || is_antiquark(in)) { found_quark = true; if(in.type != out.type) { // switch in quark type -> Wp couples to it if(found_anti){ // already found qq for coupling to W bail_out(psp, "Found second up/down pair"); } else if(std::abs(in.type)>4 || std::abs(out.type)>4){ bail_out(psp, "Found bottom/top pair"); } found_anti = true; if( is_up_type(in)) { // "up" in if(W_type > 0){ // -> only allowed u -> Wp + d if(in.type < 0 || is_up_type(out) || out.type < 0) bail_out(psp, "u -/> Wp + d"); } else { // -> only allowed ux -> Wm + dx if(in.type > 0 || is_up_type(out) || out.type > 0) bail_out(psp, "ux -/> Wm + dx"); } } else { // "down" in if(W_type > 0){ // -> only allowed dx -> Wp + ux if(in.type > 0 || is_down_type(out) || out.type > 0) bail_out(psp, "dx -/> Wp + ux"); } else { // -> only allowed d -> Wm + u if(in.type < 0 || is_down_type(out) || out.type < 0) bail_out(psp, "d -/> Wm + u"); } } } } } if(!found_quark) { bail_out(psp, "Found no initial quarks"); } else if(!found_anti){ bail_out(psp, "Found no up/down pair"); } return true; } } int main(){ constexpr std::size_t n_psp_base = 1337; const JetParameters jet_para{ fastjet::JetDefinition(fastjet::JetAlgorithm::antikt_algorithm, 0.4), 30, 5, 30}; PDF pdf(11000, pid::proton, pid::proton); constexpr double E_cms = 13000.; - constexpr double subl_change = 0.5; - constexpr auto subl_channels = Subleading::all; + constexpr double subl_chance = 0.5; + const Subleading subl_channels = subleading::ALL; const ParticlesDecayMap boson_decays{ {pid::Wp, {Decay{ {pid::e_bar, pid::nu_e}, 1.} }}, {pid::Wm, {Decay{ {pid::e, pid::nu_e_bar}, 1.} }} }; const EWConstants ew_constants{246.2196508, ParticleProperties{80.385, 2.085}, ParticleProperties{91.187, 2.495}, ParticleProperties{125, 0.004165} }; HEJ::Mixmax ran{}; // Wp2j Process proc {{pid::proton,pid::proton}, 2, pid::Wp, {}}; std::size_t n_psp = n_psp_base; for( std::size_t i = 0; i try again ++n_psp; } } std::cout << "Wp+2j: Took " << n_psp << " to generate " << n_psp_base << " successfully PSP (" << 1.*n_psp/n_psp_base << " trials/PSP)" << std::endl; // Wm2j proc = Process{{pid::proton,pid::proton}, 2, pid::Wm, {}}; n_psp = n_psp_base; for( std::size_t i = 0; i try again ++n_psp; } } std::cout << "Wm+2j: Took " << n_psp << " to generate " << n_psp_base << " successfully PSP (" << 1.*n_psp/n_psp_base << " trials/PSP)" << std::endl; std::cout << "All processes passed." << std::endl; return EXIT_SUCCESS; } diff --git a/FixedOrderGen/t/W_nj_classify.cc b/FixedOrderGen/t/W_nj_classify.cc index db1609e..2cabb8e 100644 --- a/FixedOrderGen/t/W_nj_classify.cc +++ b/FixedOrderGen/t/W_nj_classify.cc @@ -1,207 +1,227 @@ /** * \brief check that the PSP generates the all W+jet subleading processes * * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ -#ifdef NDEBUG -#undef NDEBUG -#endif #include #include #include #include #include #include #include #include #include "HEJ/Event.hh" #include "HEJ/event_types.hh" #include "HEJ/EWConstants.hh" #include "HEJ/exceptions.hh" #include "HEJ/Mixmax.hh" #include "HEJ/PDF.hh" #include "HEJ/PDG_codes.hh" #include "fastjet/JetDefinition.hh" #include "Decay.hh" #include "JetParameters.hh" #include "PhaseSpacePoint.hh" #include "Process.hh" #include "Status.hh" #include "Subleading.hh" namespace { using namespace HEJFOG; using namespace HEJ; void print_psp(PhaseSpacePoint const & psp){ std::cerr << "Process:\n" << psp.incoming()[0].type << " + "<< psp.incoming()[1].type << " -> "; for(auto const & out: psp.outgoing()){ std::cerr << out.type << " "; } std::cerr << "\n"; } void bail_out(PhaseSpacePoint const & psp, std::string msg){ print_psp(psp); throw std::logic_error{msg}; } + + #if !defined(__clang__) && defined(__GNUC__) && (__GNUC__ < 6) + // gcc version < 6 explicitly needs hash function for enum + // see https://gcc.gnu.org/bugzilla/show_bug.cgi?id=60970 + using type_map = std::unordered_map>; + #else + using type_map = std::unordered_map; + #endif } int main(){ constexpr std::size_t n_psp_base = 10375; const JetParameters jet_para{ fastjet::JetDefinition(fastjet::JetAlgorithm::antikt_algorithm, 0.4), 30, 5, 30}; PDF pdf(11000, pid::proton, pid::proton); constexpr double E_cms = 13000.; - constexpr double subl_change = 0.8; + constexpr double subl_chance = 0.8; const ParticlesDecayMap boson_decays{ {pid::Wp, {Decay{ {pid::e_bar, pid::nu_e}, 1.} }}, {pid::Wm, {Decay{ {pid::e, pid::nu_e_bar}, 1.} }} }; const EWConstants ew_constants{246.2196508, ParticleProperties{80.385, 2.085}, ParticleProperties{91.187, 2.495}, ParticleProperties{125, 0.004165} }; HEJ::Mixmax ran{}; - auto subl_channels = Subleading::all; + Subleading subl_channels = subleading::ALL; std::vector allowed_types{event_type::FKL, event_type::unob, event_type::unof, event_type::qqxexb, event_type::qqxexf}; std::cout << "Wp3j" << std::endl; // Wp3j Process proc {{pid::proton,pid::proton}, 3, pid::Wp, {}}; std::size_t n_psp = n_psp_base; - #if !defined(__clang__) && defined(__GNUC__) && (__GNUC__ < 6) - // gcc version < 6 explicitly needs hash function for enum - // see https://gcc.gnu.org/bugzilla/show_bug.cgi?id=60970 - std::unordered_map> type_counter; - #else - std::unordered_map type_counter; - #endif + type_map type_counter; for( std::size_t i = 0; i try again ++n_psp; } } std::cout << "Wp+3j: Took " << n_psp << " to generate " << n_psp_base << " successfully PSP (" << 1.*n_psp/n_psp_base << " trials/PSP)" << std::endl; std::cout << "States by classification:\n"; for(auto const & entry: type_counter){ const double fraction = static_cast(entry.second)/n_psp_base; const int percent = std::round(100*fraction); std::cout << std::left << std::setw(25) << (name(entry.first) + std::string(":")) << entry.second << " (" << percent << "%)\n"; } for(auto const & t: allowed_types){ if(type_counter[t] < 0.05 * n_psp_base){ std::cerr << "Less than 5% of the events are of type " << name(t) << std::endl; return EXIT_FAILURE; } } // Wm3j - only uno proc = Process{{pid::proton,pid::proton}, 3, pid::Wm, {}}; n_psp = n_psp_base; - subl_channels = Subleading::uno; + subl_channels.reset(); + subl_channels.set(subleading::uno); allowed_types = {event_type::FKL, event_type::unob, event_type::unof}; type_counter.clear(); for( std::size_t i = 0; i try again ++n_psp; } } std::cout << "Wm+3j (only uno): Took " << n_psp << " to generate " << n_psp_base << " successfully PSP (" << 1.*n_psp/n_psp_base << " trials/PSP)" << std::endl; std::cout << "States by classification:\n"; for(auto const & entry: type_counter){ const double fraction = static_cast(entry.second)/n_psp_base; const int percent = std::round(100*fraction); std::cout << std::left << std::setw(25) << (name(entry.first) + std::string(":")) << entry.second << " (" << percent << "%)\n"; } for(auto const & t: allowed_types){ if(type_counter[t] < 0.05 * n_psp_base){ std::cerr << "Less than 5% of the events are of type " << name(t) << std::endl; return EXIT_FAILURE; } } // Wm4j proc = Process{{pid::proton,pid::proton}, 4, pid::Wm, {}}; n_psp = n_psp_base; - subl_channels = Subleading::all; + subl_channels.set(); allowed_types = {event_type::FKL, event_type::unob, event_type::unof, event_type::qqxexb, event_type::qqxexf, event_type::qqxmid}; type_counter.clear(); + std::array wpos_counter; // position of the W boson (back, central, forward) for( std::size_t i = 0; itype==pid::Wm){ + ++wpos_counter[0][ev.type()]; + } else if(ev.outgoing().rbegin()->type==pid::Wm){ + ++wpos_counter[2][ev.type()]; + } else { + ++wpos_counter[1][ev.type()]; + } } else { // bad process -> try again ++n_psp; } } std::cout << "Wm+4j: Took " << n_psp << " to generate " << n_psp_base << " successfully PSP (" << 1.*n_psp/n_psp_base << " trials/PSP)" << std::endl; std::cout << "States by classification:\n"; for(auto const & entry: type_counter){ const double fraction = static_cast(entry.second)/n_psp_base; const int percent = std::round(100*fraction); std::cout << std::left << std::setw(25) << (name(entry.first) + std::string(":")) << entry.second << " (" << percent << "%)\n"; } for(auto const & t: allowed_types){ if(type_counter[t] < 0.03 * n_psp_base){ std::cerr << "Less than 3% of the events are of type " << name(t) << std::endl; return EXIT_FAILURE; } } + std::cout << "Stats by Wpos:\n"; + for(std::size_t i=0; i #include #include #include #include #include "HEJ/Event.hh" #include "HEJ/Mixmax.hh" #include "HEJ/ScaleFunction.hh" #include "config.hh" #include "EventGenerator.hh" #include "Status.hh" +//! throw error if condition not fulfilled +#define ASSERT(x) if(!(x)) { \ + throw std::logic_error("Assertion '" #x "' failed."); \ + } + namespace { using namespace HEJFOG; using namespace HEJ; constexpr std::size_t num_events = 1000; constexpr double invGeV2_to_pb = 389379292.; } double get_xs(std::string config_name){ auto config { load_config(config_name) }; config.events = num_events; std::shared_ptr ran{std::make_shared(11)}; HEJFOG::EventGenerator generator{ config.process, config.beam, HEJ::ScaleGenerator{ config.scales.base, config.scales.factors, config.scales.max_ratio }, config.jets, config.pdf_id, config.subleading_fraction, config.subleading_channels, config.particle_decays, config.Higgs_coupling, config.ew_parameters, ran }; double xs = 0.; for (std::size_t trials = 0; trials < config.events; ++trials){ auto ev = generator.gen_event(); - if(generator.status() != good) continue; - assert(ev); + if(generator.status() != Status::good) continue; + ASSERT(ev); ev->central().weight *= invGeV2_to_pb; ev->central().weight /= config.events; xs += ev->central().weight; } return xs; } int main(){ - double xs_W{ get_xs("config_Wp_2j.yml")}; - double xs_enu{ get_xs("config_Wp_2j_decay.yml")}; + double xs_W{ get_xs("Wp_2j.yml")}; + double xs_enu{ get_xs("Wp_2j_decay.yml")}; if(std::abs(xs_W/xs_enu-1.)>1e-6){ std::cerr << "Reconstructing the W in the runcard gave a different results (" << xs_W << " vs. "<< xs_enu << " -> " << std::abs(xs_W/xs_enu-1.)*100 << "%)\n"; return EXIT_FAILURE; } return EXIT_SUCCESS; } diff --git a/FixedOrderGen/t/config_2j.yml b/FixedOrderGen/t/runcards/2j.yml similarity index 90% rename from FixedOrderGen/t/config_2j.yml rename to FixedOrderGen/t/runcards/2j.yml index 31ebe8b..137650c 100644 --- a/FixedOrderGen/t/config_2j.yml +++ b/FixedOrderGen/t/runcards/2j.yml @@ -1,42 +1,43 @@ events: 100000 jets: min pt: 30 R: 0.4 algorithm: antikt max rapidity: 5 beam: energy: 6500 particles: [p, p] pdf: 11000 process: p p => 2j subleading fraction: 0.1 subleading channels: - - qqx - - uno + - extremal qqx + - central qqx + - unordered scales: 125 random generator: name: mixmax particle properties: Higgs: mass: 125 width: 0 W: mass: 80.385 width: 2.085 Z: mass: 91.187 width: 2.495 vev: 246.2196508 unweight: sample size: 100 max deviation: -10 diff --git a/FixedOrderGen/t/config_2j_peak.yml b/FixedOrderGen/t/runcards/2j_peak.yml similarity index 94% rename from FixedOrderGen/t/config_2j_peak.yml rename to FixedOrderGen/t/runcards/2j_peak.yml index 346e13f..deeec6e 100644 --- a/FixedOrderGen/t/config_2j_peak.yml +++ b/FixedOrderGen/t/runcards/2j_peak.yml @@ -1,41 +1,41 @@ events: 1000 jets: min pt: 20 R: 0.4 algorithm: antikt max rapidity: 5 peak pt: 30 unweight: sample size: 100 max deviation: 1 beam: energy: 6500 particles: [p, p] pdf: 11000 process: p p => 2j subleading fraction: 0.1 -subleading channels: none +subleading channels: all scales: 125 random generator: name: mixmax particle properties: Higgs: mass: 125 width: 0 W: mass: 80.385 width: 2.085 Z: mass: 91.187 width: 2.495 vev: 246.2196508 diff --git a/FixedOrderGen/t/config_3j_qqx.yml b/FixedOrderGen/t/runcards/3j_qqx.yml similarity index 94% rename from FixedOrderGen/t/config_3j_qqx.yml rename to FixedOrderGen/t/runcards/3j_qqx.yml index 610a4a8..2873ade 100644 --- a/FixedOrderGen/t/config_3j_qqx.yml +++ b/FixedOrderGen/t/runcards/3j_qqx.yml @@ -1,41 +1,41 @@ -events: 200000 +events: 300000 jets: min pt: 50 R: 0.4 algorithm: antikt max rapidity: 4 beam: energy: 6500 particles: [p, p] pdf: 11000 process: p p => 3j subleading fraction: 1 subleading channels: - - qqx + - eqqx scales: 91 random generator: name: mixmax particle properties: Higgs: mass: 125 width: 0 W: mass: 80.385 width: 2.085 Z: mass: 91.187 width: 2.495 vev: 246.2196508 unweight: sample size: 100 max deviation: -10 diff --git a/FixedOrderGen/t/config_3j_uno.yml b/FixedOrderGen/t/runcards/3j_uno.yml similarity index 100% rename from FixedOrderGen/t/config_3j_uno.yml rename to FixedOrderGen/t/runcards/3j_uno.yml diff --git a/FixedOrderGen/t/config_4j.yml b/FixedOrderGen/t/runcards/4j.yml similarity index 93% copy from FixedOrderGen/t/config_4j.yml copy to FixedOrderGen/t/runcards/4j.yml index 164f5de..8bd1134 100644 --- a/FixedOrderGen/t/config_4j.yml +++ b/FixedOrderGen/t/runcards/4j.yml @@ -1,42 +1,43 @@ -events: 200000 +events: 250000 jets: min pt: 30 R: 0.4 algorithm: antikt max rapidity: 5 beam: energy: 6500 particles: [p, p] pdf: 11000 process: p p => 4j subleading fraction: 0.1 subleading channels: - - qqx + - cqqx + - eqqx - uno scales: 125 random generator: name: mixmax particle properties: Higgs: mass: 125 width: 0 W: mass: 80.385 width: 2.085 Z: mass: 91.187 width: 2.495 vev: 246.2196508 unweight: sample size: 100 max deviation: -10 diff --git a/FixedOrderGen/t/config_4j_qqx.yml b/FixedOrderGen/t/runcards/4j_qqx.yml similarity index 96% rename from FixedOrderGen/t/config_4j_qqx.yml rename to FixedOrderGen/t/runcards/4j_qqx.yml index 1c10656..f16787e 100644 --- a/FixedOrderGen/t/config_4j_qqx.yml +++ b/FixedOrderGen/t/runcards/4j_qqx.yml @@ -1,41 +1,42 @@ events: 1100000 jets: min pt: 30 R: 0.4 algorithm: antikt max rapidity: 5 beam: energy: 3500 particles: [p, p] pdf: 11000 process: p p => 4j subleading fraction: 1 subleading channels: - - qqx + - cqqx + - eqqx scales: 125 random generator: name: mixmax particle properties: Higgs: mass: 125 width: 0 W: mass: 80.385 width: 2.085 Z: mass: 91.187 width: 2.495 vev: 246.2196508 unweight: sample size: 100 max deviation: -10 diff --git a/FixedOrderGen/t/config_4j.yml b/FixedOrderGen/t/runcards/5j.yml similarity index 77% rename from FixedOrderGen/t/config_4j.yml rename to FixedOrderGen/t/runcards/5j.yml index 164f5de..666455e 100644 --- a/FixedOrderGen/t/config_4j.yml +++ b/FixedOrderGen/t/runcards/5j.yml @@ -1,42 +1,39 @@ -events: 200000 +events: 500000 jets: min pt: 30 R: 0.4 algorithm: antikt - max rapidity: 5 + max rapidity: 6 beam: - energy: 6500 + energy: 7000 particles: [p, p] pdf: 11000 -process: p p => 4j +process: p p => 5j subleading fraction: 0.1 -subleading channels: - - qqx - - uno scales: 125 random generator: name: mixmax particle properties: Higgs: mass: 125 width: 0 W: mass: 80.385 width: 2.085 Z: mass: 91.187 width: 2.495 vev: 246.2196508 unweight: sample size: 100 max deviation: -10 diff --git a/FixedOrderGen/t/config_Wm2j.yml b/FixedOrderGen/t/runcards/Wm2j.yml similarity index 100% rename from FixedOrderGen/t/config_Wm2j.yml rename to FixedOrderGen/t/runcards/Wm2j.yml diff --git a/FixedOrderGen/t/config_Wp3j_uno.yml b/FixedOrderGen/t/runcards/Wm5j.yml similarity index 68% copy from FixedOrderGen/t/config_Wp3j_uno.yml copy to FixedOrderGen/t/runcards/Wm5j.yml index 9b39135..56f13e7 100644 --- a/FixedOrderGen/t/config_Wp3j_uno.yml +++ b/FixedOrderGen/t/runcards/Wm5j.yml @@ -1,40 +1,38 @@ events: 2000000 jets: - min pt: 100 - R: 0.5 + min pt: 50 + R: 0.4 algorithm: antikt max rapidity: 5 beam: - energy: 8000 + energy: 15000 particles: [p, p] pdf: 11000 -process: u u => W+ 3j +process: p p => W- 5j subleading fraction: 1. -subleading channels: - - unordered -scales: 125 +scales: 91 particle properties: Higgs: mass: 125 width: 0 W: mass: 80.385 width: 2.085 Z: mass: 91.187 width: 2.495 decays: - Wp: {into: [e+, nu_e], branching ratio: 1} + Wm: {into: [e-, electron_antineutrino], branching ratio: 1} random generator: name: mixmax vev: 246.2196508 diff --git a/FixedOrderGen/t/config_Wp3j_qqx.yml b/FixedOrderGen/t/runcards/Wp3j_qqx.yml similarity index 98% rename from FixedOrderGen/t/config_Wp3j_qqx.yml rename to FixedOrderGen/t/runcards/Wp3j_qqx.yml index 8b7dfe1..997713c 100644 --- a/FixedOrderGen/t/config_Wp3j_qqx.yml +++ b/FixedOrderGen/t/runcards/Wp3j_qqx.yml @@ -1,40 +1,40 @@ events: 1000000 jets: min pt: 20 R: 0.4 algorithm: antikt max rapidity: 5 beam: energy: 5000 particles: [p, p] pdf: 11000 process: u u => W+ 3j subleading fraction: 1. subleading channels: - - qqx + - cqqx scales: 125 particle properties: Higgs: mass: 125 width: 0 W: mass: 80.385 width: 2.085 Z: mass: 91.187 width: 2.495 decays: Wp: {into: [e+, nu_e], branching ratio: 1} random generator: name: mixmax vev: 246.2196508 diff --git a/FixedOrderGen/t/config_Wp3j_uno.yml b/FixedOrderGen/t/runcards/Wp3j_uno.yml similarity index 100% rename from FixedOrderGen/t/config_Wp3j_uno.yml rename to FixedOrderGen/t/runcards/Wp3j_uno.yml diff --git a/FixedOrderGen/t/config_Wp4j_qqx.yml b/FixedOrderGen/t/runcards/Wp4j_qqx.yml similarity index 96% copy from FixedOrderGen/t/config_Wp4j_qqx.yml copy to FixedOrderGen/t/runcards/Wp4j_qqx.yml index 66a842d..db3fb3f 100644 --- a/FixedOrderGen/t/config_Wp4j_qqx.yml +++ b/FixedOrderGen/t/runcards/Wp4j_qqx.yml @@ -1,40 +1,41 @@ events: 1000000 jets: min pt: 40 R: 0.4 algorithm: antikt max rapidity: 6 beam: energy: 6000 particles: [p, p] pdf: 11000 process: u u => W+ 4j subleading fraction: 1. subleading channels: - - qqx + - cqqx + - eqqx scales: 125 particle properties: Higgs: mass: 125 width: 0 W: mass: 80.385 width: 2.085 Z: mass: 91.187 width: 2.495 decays: Wp: {into: [e+, nu_e], branching ratio: 1} random generator: name: mixmax vev: 246.2196508 diff --git a/FixedOrderGen/t/config_Wp4j_qqx.yml b/FixedOrderGen/t/runcards/Wp4j_uno+eqqx.yml similarity index 81% rename from FixedOrderGen/t/config_Wp4j_qqx.yml rename to FixedOrderGen/t/runcards/Wp4j_uno+eqqx.yml index 66a842d..fc223db 100644 --- a/FixedOrderGen/t/config_Wp4j_qqx.yml +++ b/FixedOrderGen/t/runcards/Wp4j_uno+eqqx.yml @@ -1,40 +1,42 @@ -events: 1000000 +events: 10000 jets: min pt: 40 R: 0.4 algorithm: antikt max rapidity: 6 beam: energy: 6000 particles: [p, p] pdf: 11000 -process: u u => W+ 4j +# code should enforce coupling at proton leg +process: d p => W+ 4j subleading fraction: 1. subleading channels: - - qqx + - eqqx + - uno scales: 125 particle properties: Higgs: mass: 125 width: 0 W: mass: 80.385 width: 2.085 Z: mass: 91.187 width: 2.495 decays: Wp: {into: [e+, nu_e], branching ratio: 1} random generator: name: mixmax vev: 246.2196508 diff --git a/FixedOrderGen/t/config_Wp_2j.yml b/FixedOrderGen/t/runcards/Wp_2j.yml similarity index 100% copy from FixedOrderGen/t/config_Wp_2j.yml copy to FixedOrderGen/t/runcards/Wp_2j.yml diff --git a/FixedOrderGen/t/config_Wp_2j_decay.yml b/FixedOrderGen/t/runcards/Wp_2j_decay.yml similarity index 100% copy from FixedOrderGen/t/config_Wp_2j_decay.yml copy to FixedOrderGen/t/runcards/Wp_2j_decay.yml diff --git a/FixedOrderGen/t/config_Wp_2j.yml b/FixedOrderGen/t/runcards/WpqQ.yml similarity index 67% copy from FixedOrderGen/t/config_Wp_2j.yml copy to FixedOrderGen/t/runcards/WpqQ.yml index 14b9216..f57da2d 100644 --- a/FixedOrderGen/t/config_Wp_2j.yml +++ b/FixedOrderGen/t/runcards/WpqQ.yml @@ -1,35 +1,35 @@ -events: 200000 +events: 400000 jets: - min pt: 30 - R: 0.4 + min pt: 20 algorithm: antikt - max rapidity: 5 + R: 0.4 + max rapidity: 6 beam: energy: 6500 particles: [p, p] pdf: 11000 -process: p p => e+ nu_e 2j +process: u s => e+ nu_e 2j + +subleading fraction: .4 -subleading fraction: 0.37 +scales: 125. + +random generator: + name: mixmax -scales: 125 +vev: 246.2196508 particle properties: Higgs: mass: 125 width: 0.004165 W: mass: 80.385 width: 2.085 Z: - mass: 91.187 + mass: 91.1876 width: 2.495 - -random generator: - name: mixmax - -vev: 246.2196508 diff --git a/FixedOrderGen/t/config_h2j.yml b/FixedOrderGen/t/runcards/h2j.yml similarity index 100% rename from FixedOrderGen/t/config_h2j.yml rename to FixedOrderGen/t/runcards/h2j.yml diff --git a/FixedOrderGen/t/config_h2j_decay.yml b/FixedOrderGen/t/runcards/h2j_decay.yml similarity index 100% rename from FixedOrderGen/t/config_h2j_decay.yml rename to FixedOrderGen/t/runcards/h2j_decay.yml diff --git a/FixedOrderGen/t/config_h3j.yml b/FixedOrderGen/t/runcards/h3j.yml similarity index 93% rename from FixedOrderGen/t/config_h3j.yml rename to FixedOrderGen/t/runcards/h3j.yml index 47c3b1e..bf7c437 100644 --- a/FixedOrderGen/t/config_h3j.yml +++ b/FixedOrderGen/t/runcards/h3j.yml @@ -1,36 +1,36 @@ events: 200000 jets: min pt: 30 R: 0.4 algorithm: antikt max rapidity: 5 beam: energy: 6500 particles: [p, p] pdf: 11000 process: p p => h 3j subleading fraction: 0.37 -subleading channels: none +subleading channels: scales: 125 particle properties: Higgs: mass: 125 width: 0 W: mass: 80.385 width: 2.085 Z: mass: 91.187 width: 2.495 random generator: name: mixmax vev: 246.2196508 diff --git a/FixedOrderGen/t/config_h3j_uno.yml b/FixedOrderGen/t/runcards/h3j_uno.yml similarity index 100% rename from FixedOrderGen/t/config_h3j_uno.yml rename to FixedOrderGen/t/runcards/h3j_uno.yml diff --git a/FixedOrderGen/t/config_h5j.yml b/FixedOrderGen/t/runcards/h5j.yml similarity index 100% rename from FixedOrderGen/t/config_h5j.yml rename to FixedOrderGen/t/runcards/h5j.yml diff --git a/FixedOrderGen/t/config_Wp_2j_decay.yml b/FixedOrderGen/t/runcards/hqQ.yml similarity index 60% rename from FixedOrderGen/t/config_Wp_2j_decay.yml rename to FixedOrderGen/t/runcards/hqQ.yml index 37ce351..a3849ea 100644 --- a/FixedOrderGen/t/config_Wp_2j_decay.yml +++ b/FixedOrderGen/t/runcards/hqQ.yml @@ -1,38 +1,40 @@ -events: 200000 +events: 400000 jets: - min pt: 30 - R: 0.4 + min pt: 20 algorithm: antikt - max rapidity: 5 + R: 0.4 + max rapidity: 6 beam: energy: 6500 particles: [p, p] pdf: 11000 -process: p p => Wp 2j +process: u d => h 2j + +subleading fraction: .4 -subleading fraction: 0.37 +scales: 125. -scales: 125 +random generator: + name: mixmax + +decays: + Higgs: + branching ratio: 1 + into: [photon, photon] + +vev: 246.2196508 particle properties: Higgs: mass: 125 width: 0.004165 W: mass: 80.385 width: 2.085 Z: - mass: 91.187 + mass: 91.1876 width: 2.495 - -decays: - Wp: {into: [e+, nu_e], branching ratio: 1} - -random generator: - name: mixmax - -vev: 246.2196508 diff --git a/FixedOrderGen/t/config_Wp_2j.yml b/FixedOrderGen/t/runcards/qQ.yml similarity index 64% rename from FixedOrderGen/t/config_Wp_2j.yml rename to FixedOrderGen/t/runcards/qQ.yml index 14b9216..483f8f1 100644 --- a/FixedOrderGen/t/config_Wp_2j.yml +++ b/FixedOrderGen/t/runcards/qQ.yml @@ -1,35 +1,35 @@ -events: 200000 +events: 400000 jets: - min pt: 30 - R: 0.4 + min pt: 20 algorithm: antikt - max rapidity: 5 + R: 0.4 + max rapidity: 6 beam: - energy: 6500 + energy: 7000 particles: [p, p] pdf: 11000 -process: p p => e+ nu_e 2j +process: u d => 2j + +subleading fraction: .4 -subleading fraction: 0.37 +scales: 125. + +random generator: + name: mixmax -scales: 125 +vev: 246.2196508 particle properties: Higgs: mass: 125 width: 0.004165 W: mass: 80.385 width: 2.085 Z: - mass: 91.187 + mass: 91.1876 width: 2.495 - -random generator: - name: mixmax - -vev: 246.2196508 diff --git a/FixedOrderGen/t/xs_gen.cc b/FixedOrderGen/t/xs_gen.cc index b674a45..1b02581 100644 --- a/FixedOrderGen/t/xs_gen.cc +++ b/FixedOrderGen/t/xs_gen.cc @@ -1,106 +1,114 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2020 * \copyright GPLv2 or later */ #include #include #include #include #include #include "HEJ/Event.hh" #include "HEJ/Mixmax.hh" #include "HEJ/PDG_codes.hh" #include "HEJ/ScaleFunction.hh" #include "HEJ/utility.hh" #include "config.hh" #include "EventGenerator.hh" #include "Status.hh" //! throw error if condition not fulfilled #define ASSERT(x) if(!(x)) { \ throw std::logic_error("Assertion '" #x "' failed."); \ } namespace { constexpr double invGeV2_to_pb = 389379292.; } int main(int argn, char const *argv[]){ using namespace HEJFOG; - if(argn != 4){ + if(argn != 4 && argn != 5){ std::cerr << "Usage: " << argv[0] << " config.yml xs xs_err"; return EXIT_FAILURE; } + const bool short_only = argn==5; + const double xs_ref = std::stod(argv[2]); const double err_ref = std::stod(argv[3]); - const auto config = load_config(argv[1]); + auto config = load_config(argv[1]); std::shared_ptr ran{std::make_shared()}; HEJFOG::EventGenerator generator{ config.process, config.beam, HEJ::ScaleGenerator{ config.scales.base, config.scales.factors, config.scales.max_ratio }, config.jets, config.pdf_id, config.subleading_fraction, config.subleading_channels, config.particle_decays, config.Higgs_coupling, config.ew_parameters, ran }; + if(short_only) + config.events/=10.; double xs = 0., xs_err = 0.; for (std::size_t trials = 0; trials < config.events; ++trials){ auto ev = generator.gen_event(); - if(generator.status() != good) continue; + if(generator.status() != Status::good) continue; ASSERT(ev); if(config.process.boson){ const auto boson = std::find_if( begin(ev->outgoing()), end(ev->outgoing()), [&](HEJ::Particle const & p){ return p.type == *config.process.boson; } ); ASSERT(boson != end(ev->outgoing())); if(!config.process.boson_decay.empty()){ ASSERT(ev->decays().size() == 1); const auto decay = ev->decays().begin(); ASSERT(ev->outgoing().size() > decay->first); ASSERT(decay->first == static_cast( std::distance(ev->outgoing().begin(), boson))); ASSERT(decay->second.size() == 2); auto const & decay_part = decay->second; ASSERT(decay_part[0].type == config.process.boson_decay[0]); ASSERT(decay_part[1].type == config.process.boson_decay[1]); ASSERT(HEJ::nearby_ep(decay_part[0].p + decay_part[1].p, boson->p, 1e-6)); } } xs += ev->central().weight; xs_err += ev->central().weight*ev->central().weight; } xs_err = std::sqrt(xs_err); xs *= invGeV2_to_pb/config.events; xs_err *= invGeV2_to_pb/config.events; - std::cout << xs_ref << " +- " << err_ref << " (" << err_ref/xs_ref*100. << "%)" - << " ~ " << xs << " +- " << xs_err - << " (" << xs_err/xs*100. << "%)" << std::endl; + std::cout < " << std::abs(xs - xs_ref)/sqrt(xs_err*xs_err+err_ref*err_ref) << " sigma" << std::endl; ASSERT(std::abs(xs - xs_ref) < 2.*sqrt(xs_err*xs_err+err_ref*err_ref)); - ASSERT(std::abs(err_ref - xs_err) < 0.1*xs_err); + if(!short_only){ + ASSERT(std::abs(err_ref - xs_err) < 0.1*xs_err); + } return EXIT_SUCCESS; } diff --git a/doc/sphinx/HEJFOG.rst b/doc/sphinx/HEJFOG.rst index 72bc54f..b485fe4 100644 --- a/doc/sphinx/HEJFOG.rst +++ b/doc/sphinx/HEJFOG.rst @@ -1,296 +1,299 @@ 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. +provided that the ``CT10nlo`` PDF sets 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 `_ 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 might further :ref:`decay` - 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:`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 `decays`_. .. _`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 `_ 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 `_ 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`. 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:`unordered` or :code:`central qqx`. Only has an effect if :code:`subleading fraction` is non-zero. Subleading processes requested here also have to be available from HEJ, for a list of implemented processes see :ref:`event treatment`. 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. - - :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:`extremal qqx`: Production of a quark-antiquark pair as extremal + partons in rapidity. If the :code:`central 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 pair in place of two gluons adjacent in rapidity at the most forward or backwards quarks. + - :code:`central qqx`: Production of a central quark-antiquark pair. This + setting is similar to :code:`extremal qqx`, but :code:`subleading fraction` + gives the probability of emitting a quark-antiquark pair in place of two + gluons adjacent in rapidity at any point *inside* the :code:`FKL` gluon + chain. .. _`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. .. _`decays`: **decays** Optional setting specifying the decays of the particle. Only the decay into two particles is implemented. Each decay has the form :code:`boson: {into: [p1,p2], branching ratio: r}` where :code:`boson` is the corresponding boson, :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 specified or set to :code:`1`. .. _`FOG scales`: **scales** Specifies the renormalisation and factorisation scales for the output events. For details, see the corresponding entry in the HEJ 2 :ref:`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. .. _`FOG 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:`settings` for details. .. _`RanLux init`: .. _`FOG random generator`: **random generator** Sets parameters for random number generation. See the HEJ 2 :ref:`settings` for details. .. _`FOG analyses`: **analyses** Specifies the name and settings for one or more analyses library. This can be useful to specify cuts at the fixed-order level. See the corresponding entry in the HEJ 2 :ref:`settings` for details. .. _`FOG vev`: **vev** Higgs vacuum expectation value in GeV. All electro-weak constants are derived from this together with the :ref:`particle properties `. .. _`FOG particle properties`: **particle properties** Specifies various properties of the different particles (Higgs, W or Z). These have to be specified for all bosons. See the corresponding entry in the HEJ 2 :ref:`settings` for details. .. _`FOG 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:`settings` for details. diff --git a/include/HEJ/CrossSectionAccumulator.hh b/include/HEJ/CrossSectionAccumulator.hh index a8589ea..42bdd63 100644 --- a/include/HEJ/CrossSectionAccumulator.hh +++ b/include/HEJ/CrossSectionAccumulator.hh @@ -1,99 +1,107 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #pragma once #include #include #include #include #include "HEJ/event_types.hh" namespace HEJ { class Event; //! Collection of Cross Section with its uncertainty template struct XSWithError { T value = T{}; //!< Cross Section - T error = T{}; //!< Error + T error = T{}; //!< Squared error (Variance) }; /** * @brief Sum of Cross Section for different subproccess */ class CrossSectionAccumulator { public: //! Fill with single event //! @note for multiple resummation events use fill_correlated() instead void fill(Event const & ev); //! Fill by weight and type //! @note for multiple resummation events use fill_correlated() instead void fill(double wt, event_type::EventType type); //! Fill by weight, error and type //! @note The error will be _added_ to the current error void fill(double wt, double err, event_type::EventType type); /** * @brief Fill with multiple correlated weights * @details This should be used to fill multiple reweighted events, * coming from the same fixed order point. * Total error for \f$N\f$ fixed order points each giving \f$M_i\f$ * resummation events is: * \f[ * \delta^2=\sum_i \left(\sum_j w_{i,j}\right)^2 * +\sum_{i,j} \left(w_{i,j}\right)^2, * \f] * @note This is equivalent to fill() for only one reweighted event * coming from each fixed order point (\f$M_i=1\f$) */ void fill_correlated(std::vector const & evts); //! iterator implementation of fill_correlated() template void fill_correlated(ConstIt begin, ConstIt end); //! explicit version of fill_correlated() by giving sum(wt) and sum(wt^2) void fill_correlated(double sum, double sum2, event_type::EventType type); //! begin of Cross Section and error for subprocesses auto begin() const { return std::begin(xs_); } //! end of Cross Section and error for subprocesses auto end() const { return std::end(xs_); } //! total Cross Section and error XSWithError const & total() const { return total_; } + //! Cross Section and error of specific type + XSWithError const & at(event_type::EventType type) const { + return xs_.at(type); + } + //! Cross Section and error of specific type + XSWithError const & operator[](event_type::EventType type) const { + return at(type); + } private: std::map> xs_; XSWithError total_; }; //! Print CrossSectionAccumulator to stream std::ostream& operator<<(std::ostream& os, const CrossSectionAccumulator& xs); // ------------ Implementation ------------ template void CrossSectionAccumulator::fill_correlated(ConstIt begin, ConstIt end){ if(std::distance(begin, end) < 2){ // only one event fill(*begin); return; } double sum = 0.; double sum2 = 0.; const auto type = begin->type(); for(; begin != end; ++begin){ double const wt = begin->central().weight; sum += wt; sum2 += wt*wt; } fill_correlated(sum, sum2, type); } } // namespace HEJ diff --git a/include/HEJ/PDF.hh b/include/HEJ/PDF.hh index adbc025..44a09b2 100644 --- a/include/HEJ/PDF.hh +++ b/include/HEJ/PDF.hh @@ -1,81 +1,81 @@ /** \file * * \brief Contains all the necessary classes and functions for interaction with PDFs. * * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #pragma once #include #include #include #include "HEJ/PDG_codes.hh" namespace LHAPDF { class PDF; } namespace HEJ { //! Class for interaction with a PDF set class PDF { public: /** * \brief PDF Constructor * @param id Particle ID according to PDG * @param beam1 Particle ID of particle in beam 1 * @param beam2 Particle ID of particle in beam 2 */ PDF(int id, ParticleID beam1, ParticleID beam2); /** * \brief Calculate the pdf value x*f(x, q) * @param beam_idx Beam number (0 or 1) * @param x Momentum fraction * @param q Energy scale * @param id PDG particle id * @returns x*f(x, q) * * Returns 0 if x or q are outside the range covered by the PDF set */ double pdfpt(std::size_t beam_idx, double x, double q, ParticleID id) const; /** * \brief Value of the strong coupling \f$\alpha_s(q)\f$ at the given scale * @param q Renormalisation scale * @returns Value of the strong coupling constant */ double Halphas(double q) const; //! Check if the energy scale is within the range covered by the PDF set /** * @param q Energy Scale * @returns true if q is within the covered range, false otherwise */ bool inRangeQ(double q) const; //! Check if the momentum fraction is within the range covered by the PDF set /** * @param x Momentum Fraction * @returns true if x is within the covered range, false otherwise */ bool inRangeX(double x) const; //! PDF id of the current set int id() const; PDF(PDF const & other) = delete; //!< Copy not possible PDF & operator=(PDF const & other) = delete; //!< Copy not possible - PDF(PDF && other) = default; //!< Moving allowed - PDF & operator=(PDF && other) = default; //!< Moving allowed + PDF(PDF && other); //!< Moving allowed + PDF & operator=(PDF && other); //!< Moving allowed ~PDF(); private: //! @internal unique_ptr does not allow copy std::unique_ptr pdf_; std::array beamtype_; }; } // namespace HEJ diff --git a/include/HEJ/PDG_codes.hh b/include/HEJ/PDG_codes.hh index 4f47e3c..9c4571c 100644 --- a/include/HEJ/PDG_codes.hh +++ b/include/HEJ/PDG_codes.hh @@ -1,238 +1,241 @@ /** \file PDG_codes.hh * \brief Contains the Particle IDs of all relevant SM particles. * * Large enumeration included which has multiple entries for potential * alternative names of different particles. There are also functions * which can be used to determine if a particle is a parton or if * it is a non-gluon boson. * * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #pragma once #include #include namespace HEJ { //! particle ids according to PDG namespace pid { //! The possible particle identities. We use PDG IDs as standard. enum ParticleID: int{ //! Unspecified type, should never be used!, debug only unspecified = 0, d = 1, /*!< Down Quark */ down = d, /*!< Down Quark */ u = 2, /*!< Up Quark */ up = u, /*!< Up Quark */ s = 3, /*!< Strange Quark */ strange = s, /*!< Strange Quark */ c = 4, /*!< Charm Quark */ charm = c, /*!< Charm Quark */ b = 5, /*!< Bottom Quark */ bottom = b, /*!< Bottom Quark */ t = 6, /*!< Top Quark */ top = t, /*!< Top Quark */ e = 11, /*!< Electron */ electron = e, /*!< Electron */ nu_e = 12, /*!< Electron Neutrino */ electron_neutrino = nu_e, /*!< Electron neutrino */ mu = 13, /*!< Muon */ muon = mu, /*!< Muon */ nu_mu = 14, /*!< Muon Neutrino */ muon_neutrino = nu_mu, /*!< Muon Neutrino */ tau = 15, /*!< Tau */ nu_tau = 16, /*!< Tau Neutrino */ tau_neutrino = nu_tau, /*!< Tau Neutrino */ d_bar = -d, /*!< Anti-Down Quark */ antidown = d_bar, /*!< Anti-Down Quark */ u_bar = -u, /*!< Anti-Up quark */ antiup = -u, /*!< Anti-Up quark */ s_bar = -s, /*!< Anti-Strange Quark */ antistrange = -s, /*!< Anti-Strange Quark */ c_bar = -c, /*!< Anti-Charm Quark */ anticharm = -c, /*!< Anti-Charm Quark */ b_bar = -b, /*!< Anti-Bottom Quark */ antibottom = -b, /*!< Anti-Bottom Quark */ t_bar = -t, /*!< Anti-Top Quark */ antitop = -t, /*!< Anti-Top Quark */ e_bar = -e, /*!< Positron */ positron = e_bar, /*!< Positron */ antielectron = positron, /*!< Positron */ nu_e_bar = -nu_e, /*!< Electron Anti-Neutrino */ electron_antineutrino = nu_e_bar,/*!< Electron Anti-Neutrino */ mu_bar = -mu, /*!< Anti-Muon */ antimuon = -mu, /*!< Anti-Muon */ nu_mu_bar = -nu_mu, /*!< Muon Anti-Neutrino */ muon_antineutrino = nu_mu_bar, /*!< Muon Anti-Neutrino */ tau_bar = -tau, /*!< Anti-Tau */ antitau = tau_bar, /*!< Anti-Tau */ nu_tau_bar = -nu_tau, /*!< Tau Anti-Neutrino */ tau_antineutrino = nu_tau_bar, /*!< Tau Anti-Neutrino */ gluon = 21, /*!< Gluon */ g = gluon, /*!< Gluon */ photon = 22, /*!< Photon */ gamma = photon, /*!< Photon */ Z = 23, /*!< Z Boson */ Z_photon_mix = 81, /*!< Z/photon superposition */ Z_gamma_mix = Z_photon_mix, /*!< Z/photon superposition */ Wp = 24, /*!< W- Boson */ Wm = -Wp, /*!< W+ Boson */ h = 25, /*!< Higgs Boson */ Higgs = h, /*!< Higgs Boson */ higgs = h, /*!< Higgs Boson */ p = 2212, /*!< Proton */ proton = p, /*!< Proton */ p_bar = -p, /*!< Anti-Proton */ antiproton = p_bar, /*!< Anti-Proton */ }; //! Get the of the particle with the given PDG ID std::string name(ParticleID id); + //! return the negative flavour of the given PDG ID + ParticleID anti(ParticleID id); + } // namespace pid using ParticleID = pid::ParticleID; //! Convert a particle name to the corresponding PDG particle ID ParticleID to_ParticleID(std::string const & name); /** * \brief Function to determine if particle is a quark * @param id PDG ID of particle * @returns true if the particle is a quark, false otherwise */ inline constexpr bool is_quark(ParticleID id){ return (id >= pid::down && id <= pid::top); } /** * \brief Function to determine if particle is an antiquark * @param id PDG ID of particle * @returns true if the particle is an antiquark, false otherwise */ inline constexpr bool is_antiquark(ParticleID id){ return (id <= pid::d_bar && id >= pid::t_bar); } /** * \brief Function to determine if particle is an (anti-)quark * @param id PDG ID of particle * @returns true if the particle is a quark or antiquark, false otherwise */ inline constexpr bool is_anyquark(ParticleID id){ return is_quark(id) || is_antiquark(id); } /** * \brief Function to determine if particle is a gluon * @param id PDG ID of particle * @returns true if the particle is a gluon, false otherwise */ inline constexpr bool is_gluon(ParticleID id){ return id == pid::gluon; } /** * \brief Function to determine if particle is a parton * @param id PDG ID of particle * @returns true if the particle is a parton, false otherwise */ inline constexpr bool is_parton(ParticleID id){ return is_gluon(id) || (is_anyquark(id) && std::abs(id) != pid::top); } /** * \brief function to determine if the particle is a photon, W or Z * @param id PDG ID of particle * @returns true if the partice is an A,W,Z, or H, false otherwise */ inline constexpr bool is_AWZ_boson(ParticleID id){ return id == pid::Wm || (id >= pid::photon && id <= pid::Wp) || id == pid::Z_photon_mix; } /** * \brief function to determine if the particle is a photon, W, Z, or Higgs * boson * @param id PDG ID of particle * @returns true if the partice is an A,W,Z, or H, false otherwise */ inline constexpr bool is_AWZH_boson(ParticleID id){ return is_AWZ_boson(id) || (id == pid::Higgs); } /** * \brief Function to determine if particle is a lepton * @param id PDG ID of particle * @returns true if the particle is a lepton, false otherwise */ inline constexpr bool is_lepton(ParticleID id){ return (id >= pid::electron && id <= pid::tau_neutrino); } /** * \brief Function to determine if particle is an antilepton * @param id PDG ID of particle * @returns true if the particle is an antilepton, false otherwise */ inline constexpr bool is_antilepton(ParticleID id){ return (id <= pid::positron && id >= pid::nu_tau_bar); } /** * \brief Function to determine if particle is an (anti-)lepton * @param id PDG ID of particle * @returns true if the particle is a lepton or antilepton, * false otherwise */ inline constexpr bool is_anylepton(ParticleID id){ return ( is_lepton(id) || is_antilepton(id)); } /** * \brief Function to determine if particle is a neutrino * @param id PDG ID of particle * @returns true if the particle is a neutrino, false otherwise */ inline constexpr bool is_neutrino(ParticleID id){ return (id == pid::nu_e || id == pid::tau_neutrino || id == pid::muon_neutrino); } /** * \brief Function to determine if particle is an antineutrino * @param id PDG ID of particle * @returns true if the particle is an antineutrino, false otherwise */ inline constexpr bool is_antineutrino(ParticleID id){ return (id == pid::nu_e_bar || id == pid::nu_tau_bar || id == pid::nu_mu_bar); } /** * \brief Function to determine if particle is an (anti-)neutrino * @param id PDG ID of particle * @returns true if the particle is a neutrino or antineutrino, * false otherwise */ inline constexpr bool is_anyneutrino(ParticleID id){ return ( is_neutrino(id) || is_antineutrino(id)); } } // namespace HEJ diff --git a/src/Event.cc b/src/Event.cc index ba02818..998daef 100644 --- a/src/Event.cc +++ b/src/Event.cc @@ -1,1179 +1,1171 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #include "HEJ/Event.hh" #include #include #include #include #include #include #include #include #include #include #include "fastjet/ClusterSequence.hh" #include "fastjet/JetDefinition.hh" #include "fastjet/PseudoJet.hh" #include "LHEF/LHEF.h" #include "HEJ/Constants.hh" #include "HEJ/PDG_codes.hh" #include "HEJ/RNG.hh" #include "HEJ/exceptions.hh" #include "HEJ/optional.hh" namespace HEJ { namespace { using std::size_t; //! LHE status codes namespace lhe_status { enum Status: int { in = -1, decay = 2, out = 1, }; } using LHE_Status = lhe_status::Status; //! true if leptonic W decay bool valid_W_decay( int const w_type, // sign of W std::vector const & decays ){ if(decays.size() != 2) // no 1->2 decay return false; const int pidsum = decays[0].type + decays[1].type; if( std::abs(pidsum) != 1 || pidsum != w_type ) // correct charge return false; // leptonic decay (only check first, second follows from pidsum) if( w_type == 1 ) // W+ return is_antilepton(decays[0]) || is_neutrino(decays[0]); // W- return is_lepton(decays[0]) || is_antineutrino(decays[0]); } //! true for Z decay to charged leptons bool valid_Z_decay(std::vector const & decays){ if(decays.size() != 2) // no 1->2 decay return false; const int pidsum = decays[0].type + decays[1].type; if( std::abs(pidsum) != 0 ) // correct charge return false; // leptonic decay (only check first, second follows from pidsum) return is_anylepton(decays[0]) && !is_anyneutrino(decays[0]); } /// @name helper functions to determine event type //@{ /** * \brief check if final state valid for HEJ * * check if there is at most one photon, W, H, Z in the final state * and all the rest are quarks or gluons */ bool final_state_ok(Event const & ev){ std::vector const & outgoing = ev.outgoing(); if(ev.decays().size() > 1) // at most one decay return false; bool has_AWZH_boson = false; for( size_t i=0; ifirst != i ) return false; if( !valid_W_decay(out.type>0?+1:-1, ev.decays().cbegin()->second) ) return false; } // valid decay for Z if(out.type == ParticleID::Z_photon_mix){ // exactly 1 decay if( ev.decays().size() != 1 || ev.decays().cbegin()->first != i ) return false; if( !valid_Z_decay(ev.decays().cbegin()->second) ) return false; } } else if(! is_parton(out.type)) return false; } return true; } /** * returns all EventTypes implemented in HEJ */ size_t implemented_types(std::vector const & bosons){ using namespace event_type; if(bosons.empty()) return FKL | unob | unof | qqxexb | qqxexf | qqxmid; if(bosons.size()>1) return non_resummable; // multi boson switch (bosons[0].type) { case ParticleID::Wp: case ParticleID::Wm: return FKL | unob | unof | qqxexb | qqxexf | qqxmid; case ParticleID::Z_photon_mix: return FKL | unob | unof; case ParticleID::h: return FKL | unob | unof; default: return non_resummable; } } /** * \brief function which determines if type change is consistent with Wp emission. * @param in incoming Particle id * @param out outgoing Particle id * @param qqx Current both incoming/both outgoing? * * \see is_Wm_Change */ bool is_Wp_Change(ParticleID in, ParticleID out, bool qqx){ - if(!qqx && (in==-1 || in== 2 || in==-3 || in== 4)) return out== (in-1); - if( qqx && (in== 1 || in==-2 || in== 3 || in==-4)) return out==-(in+1); + using namespace pid; + if(!qqx && (in==d_bar || in==u || in==s_bar || in==c)) + return out == (in-1); + if( qqx && (in==d || in==u_bar || in==s || in==c_bar)) + return out == -(in+1); return false; } /** * \brief function which determines if type change is consistent with Wm emission. * @param in incoming Particle id * @param out outgoing Particle id * @param qqx Current both incoming/both outgoing? * * Ensures that change type of quark line is possible by a flavour changing * Wm emission. Allows checking of qqx currents also. */ bool is_Wm_Change(ParticleID in, ParticleID out, bool qqx){ - if(!qqx && (in== 1 || in==-2 || in== 3 || in==-4)) return out== (in+1); - if( qqx && (in==-1 || in== 2 || in==-3 || in== 4)) return out==-(in-1); + using namespace pid; + if(!qqx && (in==d || in==u_bar || in==s || in==c_bar)) + return out == (in+1); + if( qqx && (in==d_bar || in==u || in==s_bar || in==c)) + return out == -(in-1); return false; } /** * \brief checks if particle type remains same from incoming to outgoing * @param in incoming Particle * @param out outgoing Particle * @param qqx Current both incoming/outgoing? */ bool no_flavour_change(ParticleID in, ParticleID out, bool qqx){ const int qqxCurrent = qqx?-1:1; if(std::abs(in)<=pid::top || in==pid::gluon) return (in==out*qqxCurrent); return false; } bool has_2_jets(Event const & event){ return event.jets().size() >= 2; } /** * \brief check if we have a valid Impact factor * @param in incoming Particle * @param out outgoing Particle * @param qqx Current both incoming/outgoing? * @param W_change returns +1 if Wp, -1 if Wm, else 0 */ bool is_valid_impact_factor( ParticleID in, ParticleID out, bool qqx, int & W_change ){ if( no_flavour_change(in, out, qqx) ){ return true; } if( is_Wp_Change(in, out, qqx) ) { W_change+=1; return true; } if( is_Wm_Change(in, out, qqx) ) { W_change-=1; return true; } return false; } //! Returns all possible classifications from the impact factors // the beginning points are changed s.t. after the the classification they // point to the beginning of the (potential) FKL chain // sets W_change: + if Wp change // 0 if no change // - if Wm change // This function can be used with forward & backwards iterators template size_t possible_impact_factors( ParticleID incoming_id, // incoming OutIterator & begin_out, OutIterator const & end_out, // outgoing int & W_change, std::vector const & boson, bool const backward // backward? ){ using namespace event_type; assert(boson.size() < 2); // keep track of all states that we don't test size_t not_tested = qqxmid; if(backward) not_tested |= unof | qqxexf; else not_tested |= unob | qqxexb; // Is this LL current? if( is_valid_impact_factor(incoming_id, begin_out->type, false, W_change) ){ ++begin_out; return not_tested | FKL; } // or NLL current? // -> needs two partons in two different jets if( std::distance(begin_out, end_out)>=2 ){ auto next = std::next(begin_out); // Is this unordered emisson? if( incoming_id!=pid::gluon && begin_out->type==pid::gluon ){ if( is_valid_impact_factor( incoming_id, next->type, false, W_change ) ){ // veto Higgs inside uno assert(next!=end_out); if( !boson.empty() && boson.front().type == ParticleID::h ){ if( (backward && boson.front().rapidity() < next->rapidity()) ||(!backward && boson.front().rapidity() > next->rapidity())) return non_resummable; } begin_out = std::next(next); return not_tested | (backward?unob:unof); } } // Is this QQbar? else if( incoming_id==pid::gluon ){ if( is_valid_impact_factor( begin_out->type, next->type, true, W_change ) ){ // veto Higgs inside qqx assert(next!=end_out); if( !boson.empty() && boson.front().type == ParticleID::h ){ if( (backward && boson.front().rapidity() < next->rapidity()) ||(!backward && boson.front().rapidity() > next->rapidity())) return non_resummable; } begin_out = std::next(next); return not_tested | (backward?qqxexb:qqxexf); } } } return non_resummable; } //! Returns all possible classifications from central emissions // the beginning points are changed s.t. after the the classification they // point to the end of the emission chain // sets W_change: + if Wp change // 0 if no change // - if Wm change template size_t possible_central( OutIterator & begin_out, OutIterator const & end_out, int & W_change, std::vector const & boson ){ using namespace event_type; assert(boson.size() < 2); // if we already passed the central chain, // then it is not a valid all-order state if(std::distance(begin_out, end_out) < 0) return non_resummable; // keep track of all states that we don't test size_t possible = unob | unof | qqxexb | qqxexf; // Find the first non-gluon/non-FKL while( (begin_out->type==pid::gluon) && (begin_out!=end_out) ){ ++begin_out; } // end of chain -> FKL if( begin_out==end_out ){ return possible | FKL; } // is this a qqbar-pair? // needs two partons in two separate jets auto next = std::next(begin_out); if( is_valid_impact_factor( begin_out->type, next->type, true, W_change ) ){ // veto Higgs inside qqx if( !boson.empty() && boson.front().type == ParticleID::h && boson.front().rapidity() > begin_out->rapidity() && boson.front().rapidity() < next->rapidity() ){ return non_resummable; } begin_out = std::next(next); // remaining chain should be pure gluon/FKL for(; begin_out!=end_out; ++begin_out){ if(begin_out->type != pid::gluon) return non_resummable; } return possible | qqxmid; } return non_resummable; } /** * \brief Checks for all event types * @param ev Event * @returns Event Type * */ event_type::EventType classify(Event const & ev){ using namespace event_type; if(! has_2_jets(ev)) return no_2_jets; // currently we can't handle multiple boson states in the ME. So they are // considered "bad_final_state" even though the "classify" could work with // them. if(! final_state_ok(ev)) return bad_final_state; // initialise variables auto const & in = ev.incoming(); // range for current checks auto begin_out{ev.cbegin_partons()}; auto end_out{ev.crbegin_partons()}; assert(std::distance(begin(in), end(in)) == 2); assert(std::distance(begin_out, end_out.base()) >= 2); assert(std::is_sorted(begin_out, end_out.base(), rapidity_less{})); auto const boson{ filter_AWZH_bosons(ev.outgoing()) }; // we only allow one boson through final_state_ok assert(boson.size()<=1); // keep track of potential W couplings, at the end the sum should be 0 int remaining_Wp = 0; int remaining_Wm = 0; if(!boson.empty() && std::abs(boson.front().type) == ParticleID::Wp ){ if(boson.front().type>0) ++remaining_Wp; else ++remaining_Wm; } int W_change = 0; size_t final_type = ~(no_2_jets | bad_final_state); // check forward impact factor final_type &= possible_impact_factors( in.front().type, begin_out, end_out.base(), W_change, boson, true ); if( final_type == non_resummable ) return non_resummable; if(W_change>0) remaining_Wp-=W_change; else if(W_change<0) remaining_Wm+=W_change; W_change = 0; // check backward impact factor final_type &= possible_impact_factors( in.back().type, end_out, std::make_reverse_iterator(begin_out), W_change, boson, false ); if( final_type == non_resummable ) return non_resummable; if(W_change>0) remaining_Wp-=W_change; else if(W_change<0) remaining_Wm+=W_change; W_change = 0; // check central emissions final_type &= possible_central( begin_out, end_out.base(), W_change, boson ); if( final_type == non_resummable ) return non_resummable; if(W_change>0) remaining_Wp-=W_change; else if(W_change<0) remaining_Wm+=W_change; // Check whether the right number of Ws are present if( remaining_Wp != 0 || remaining_Wm != 0 ) return non_resummable; // result has to be unique if( (final_type & (final_type-1)) != 0) return non_resummable; // check that each sub processes is implemented // (has to be done at the end) if( (final_type & ~implemented_types(boson)) != 0 ) return non_resummable; return static_cast(final_type); } //@} Particle extract_particle(LHEF::HEPEUP const & hepeup, size_t i){ auto id = static_cast(hepeup.IDUP[i]); auto colour = is_parton(id)?hepeup.ICOLUP[i]:optional(); return { id, { hepeup.PUP[i][0], hepeup.PUP[i][1], hepeup.PUP[i][2], hepeup.PUP[i][3] }, colour }; } bool is_decay_product(std::pair const & mothers){ if(mothers.first == 0) return false; return mothers.second == 0 || mothers.first == mothers.second; } } // namespace Event::EventData::EventData(LHEF::HEPEUP const & hepeup){ parameters.central = EventParameters{ hepeup.scales.mur, hepeup.scales.muf, hepeup.XWGTUP }; size_t in_idx = 0; for (int i = 0; i < hepeup.NUP; ++i) { // skip decay products // we will add them later on, but we have to ensure that // the decayed particle is added before if(is_decay_product(hepeup.MOTHUP[i])) continue; auto particle = extract_particle(hepeup, i); // needed to identify mother particles for decay products particle.p.set_user_index(i+1); if(hepeup.ISTUP[i] == LHE_Status::in){ if(in_idx > incoming.size()) { throw std::invalid_argument{ "Event has too many incoming particles" }; } incoming[in_idx++] = std::move(particle); } else outgoing.emplace_back(std::move(particle)); } // add decay products for (int i = 0; i < hepeup.NUP; ++i) { if(!is_decay_product(hepeup.MOTHUP[i])) continue; const int mother_id = hepeup.MOTHUP[i].first; const auto mother = std::find_if( begin(outgoing), end(outgoing), [mother_id](Particle const & particle){ return particle.p.user_index() == mother_id; } ); if(mother == end(outgoing)){ throw std::invalid_argument{"invalid decay product parent"}; } const int mother_idx = std::distance(begin(outgoing), mother); assert(mother_idx >= 0); decays[mother_idx].emplace_back(extract_particle(hepeup, i)); } } Event::Event( UnclusteredEvent const & ev, fastjet::JetDefinition const & jet_def, double const min_jet_pt ): Event( Event::EventData{ ev.incoming, ev.outgoing, ev.decays, Parameters{ev.central, ev.variations} }.cluster(jet_def, min_jet_pt) ) {} //! @TODO remove in HEJ 2.2.0 UnclusteredEvent::UnclusteredEvent(LHEF::HEPEUP const & hepeup){ Event::EventData const evData{hepeup}; incoming = evData.incoming; outgoing = evData.outgoing; decays = evData.decays; central = evData.parameters.central; variations = evData.parameters.variations; } void Event::EventData::sort(){ // sort particles std::sort( begin(incoming), end(incoming), [](Particle const & o1, Particle const & o2){return o1.p.pz() idx(old_outgoing.size()); std::iota(idx.begin(), idx.end(), 0); std::sort(idx.begin(), idx.end(), [&old_outgoing](size_t i, size_t j){ return old_outgoing[i].rapidity() < old_outgoing[j].rapidity(); }); outgoing.clear(); outgoing.reserve(old_outgoing.size()); for(size_t i: idx) { outgoing.emplace_back(std::move(old_outgoing[i])); } // find decays again if(!decays.empty()){ auto old_decays = std::move(decays); decays.clear(); for(size_t i=0; isecond)); } assert(old_decays.size() == decays.size()); } } namespace { Particle reconstruct_boson(std::vector const & leptons) { Particle decayed_boson; decayed_boson.p = leptons[0].p + leptons[1].p; const int pidsum = leptons[0].type + leptons[1].type; if(pidsum == +1) { assert(is_antilepton(leptons[0])); if(is_antineutrino(leptons[0])) { throw not_implemented{"lepton-flavour violating final state"}; } assert(is_neutrino(leptons[1])); // charged antilepton + neutrino means we had a W+ decayed_boson.type = pid::Wp; } else if(pidsum == -1) { assert(is_antilepton(leptons[0])); if(is_neutrino(leptons[1])) { throw not_implemented{"lepton-flavour violating final state"}; } assert(is_antineutrino(leptons[0])); // charged lepton + antineutrino means we had a W- decayed_boson.type = pid::Wm; } else if(pidsum == 0) { assert(is_anylepton(leptons[0])); if(is_anyneutrino(leptons[0])) { throw not_implemented{"final state with two neutrinos"}; } // charged lepton-antilepton pair means we had a Z/photon decayed_boson.type = pid::Z_photon_mix; } else { throw not_implemented{ "final state with leptons " + name(leptons[0].type) + " and " + name(leptons[1].type) }; } return decayed_boson; } } // namespace void Event::EventData::reconstruct_intermediate() { const auto begin_leptons = std::partition( begin(outgoing), end(outgoing), [](Particle const & p) {return !is_anylepton(p);} ); // We can only reconstruct FS with 2 leptons if(std::distance(begin_leptons, end(outgoing)) != 2) return; std::vector leptons(begin_leptons, end(outgoing)); std::sort( begin(leptons), end(leptons), [](Particle const & p0, Particle const & p1) { assert(is_anylepton(p0) && is_anylepton(p1)); return p0.type < p1.type; } ); // `reconstruct_boson` can throw, it should therefore be called before // changing `outgoing` to allow the user to recover the original EventData auto boson = reconstruct_boson(leptons); outgoing.erase(begin_leptons, end(outgoing)); outgoing.emplace_back(std::move(boson)); decays.emplace(outgoing.size()-1, std::move(leptons)); } Event Event::EventData::cluster( fastjet::JetDefinition const & jet_def, double const min_jet_pt ){ sort(); return Event{ std::move(incoming), std::move(outgoing), std::move(decays), std::move(parameters), jet_def, min_jet_pt }; } Event::Event( std::array && incoming, std::vector && outgoing, std::unordered_map> && decays, Parameters && parameters, fastjet::JetDefinition const & jet_def, double const min_jet_pt ): incoming_{std::move(incoming)}, outgoing_{std::move(outgoing)}, decays_{std::move(decays)}, parameters_{std::move(parameters)}, cs_{ to_PseudoJet( filter_partons(outgoing_) ), jet_def }, min_jet_pt_{min_jet_pt} { jets_ = sorted_by_rapidity(cs_.inclusive_jets(min_jet_pt_)); assert(std::is_sorted(begin(outgoing_), end(outgoing_), rapidity_less{})); type_ = classify(*this); } namespace { //! check that Particles have a reasonable colour bool correct_colour(Particle const & part){ ParticleID id{ part.type }; if(!is_parton(id)) return !part.colour; if(!part.colour) return false; Colour const & col{ *part.colour }; if(is_quark(id)) return col.first != 0 && col.second == 0; if(is_antiquark(id)) return col.first == 0 && col.second != 0; assert(id==ParticleID::gluon); return col.first != 0 && col.second != 0 && col.first != col.second; } //! Connect parton to t-channel colour line & update the line //! returns false if connection not possible template bool try_connect_t(OutIterator const & it_part, Colour & line_colour){ if( line_colour.first == it_part->colour->second ){ line_colour.first = it_part->colour->first; return true; } if( line_colour.second == it_part->colour->first ){ line_colour.second = it_part->colour->second; return true; } return false; } //! Connect parton to u-channel colour line & update the line //! returns false if connection not possible template bool try_connect_u(OutIterator & it_part, Colour & line_colour){ auto it_next = std::next(it_part); if( try_connect_t(it_next, line_colour) && try_connect_t(it_part, line_colour) ){ it_part=it_next; return true; } return false; } } // namespace bool Event::is_leading_colour() const { if( !correct_colour(incoming()[0]) || !correct_colour(incoming()[1]) ) return false; Colour line_colour = *incoming()[0].colour; std::swap(line_colour.first, line_colour.second); // reasonable colour if(!std::all_of(outgoing().cbegin(), outgoing().cend(), correct_colour)) return false; for(auto it_part = cbegin_partons(); it_part!=cend_partons(); ++it_part){ switch (type()) { case event_type::FKL: if( !try_connect_t(it_part, line_colour) ) return false; break; case event_type::unob: case event_type::qqxexb: { if( !try_connect_t(it_part, line_colour) // u-channel only allowed at impact factor && (std::distance(cbegin_partons(), it_part)!=0 || !try_connect_u(it_part, line_colour))) return false; break; } case event_type::unof: case event_type::qqxexf: { if( !try_connect_t(it_part, line_colour) // u-channel only allowed at impact factor && (std::distance(it_part, cend_partons())!=2 || !try_connect_u(it_part, line_colour))) return false; break; } case event_type::qqxmid:{ auto it_next = std::next(it_part); if( !try_connect_t(it_part, line_colour) // u-channel only allowed at qqx/qxq pair && ( ( !(is_quark(*it_part) && is_antiquark(*it_next)) && !(is_antiquark(*it_part) && is_quark(*it_next))) || !try_connect_u(it_part, line_colour)) ) return false; break; } default: throw std::logic_error{"unreachable"}; } // no colour singlet exchange/disconnected diagram if(line_colour.first == line_colour.second) return false; } return (incoming()[1].colour->first == line_colour.first) && (incoming()[1].colour->second == line_colour.second); } namespace { //! connect incoming Particle to colour flow void connect_incoming(Particle & in, int & colour, int & anti_colour){ in.colour = std::make_pair(anti_colour, colour); // gluon if(in.type == pid::gluon) return; if(in.type > 0){ // quark assert(is_quark(in)); in.colour->second = 0; colour*=-1; return; } // anti-quark assert(is_antiquark(in)); in.colour->first = 0; anti_colour*=-1; } //! connect outgoing Particle to t-channel colour flow template void connect_tchannel( OutIterator & it_part, int & colour, int & anti_colour, RNG & ran ){ assert(colour>0 || anti_colour>0); if(it_part->type == ParticleID::gluon){ // gluon if(colour>0 && anti_colour>0){ // on g line => connect to colour OR anti-colour (random) if(ran.flat() < 0.5){ it_part->colour = std::make_pair(colour+2,colour); colour+=2; } else { it_part->colour = std::make_pair(anti_colour, anti_colour+2); anti_colour+=2; } } else if(colour > 0){ // on q line => connect to available colour it_part->colour = std::make_pair(colour+2, colour); colour+=2; } else { assert(colour<0 && anti_colour>0); // on qx line => connect to available anti-colour it_part->colour = std::make_pair(anti_colour, anti_colour+2); anti_colour+=2; } } else if(is_quark(*it_part)) { // quark assert(anti_colour>0); if(colour>0){ // on g line => connect and remove anti-colour it_part->colour = std::make_pair(anti_colour, 0); anti_colour+=2; anti_colour*=-1; } else { // on qx line => new colour colour*=-1; it_part->colour = std::make_pair(colour, 0); } } else if(is_antiquark(*it_part)) { // anti-quark assert(colour>0); if(anti_colour>0){ // on g line => connect and remove colour it_part->colour = std::make_pair(0, colour); colour+=2; colour*=-1; } else { // on q line => new anti-colour anti_colour*=-1; it_part->colour = std::make_pair(0, anti_colour); } } else { // not a parton assert(!is_parton(*it_part)); it_part->colour = {}; } } //! connect to t- or u-channel colour flow template void connect_utchannel( OutIterator & it_part, int & colour, int & anti_colour, RNG & ran ){ OutIterator it_first = it_part++; if(ran.flat()<.5) {// t-channel connect_tchannel(it_first, colour, anti_colour, ran); connect_tchannel(it_part, colour, anti_colour, ran); } else { // u-channel connect_tchannel(it_part, colour, anti_colour, ran); connect_tchannel(it_first, colour, anti_colour, ran); } } } // namespace bool Event::generate_colours(RNG & ran){ // generate only for HEJ events if(!event_type::is_resummable(type())) return false; assert(std::is_sorted( begin(outgoing()), end(outgoing()), rapidity_less{})); assert(incoming()[0].pz() < incoming()[1].pz()); // positive (anti-)colour -> can connect // negative (anti-)colour -> not available/used up by (anti-)quark int colour = COLOUR_OFFSET; int anti_colour = colour+1; // initialise first connect_incoming(incoming_[0], colour, anti_colour); // reset outgoing colours std::for_each(outgoing_.begin(), outgoing_.end(), [](Particle & part){ part.colour = {};}); for(auto it_part = begin_partons(); it_part!=end_partons(); ++it_part){ switch (type()) { // subleading can connect to t- or u-channel case event_type::unob: case event_type::qqxexb: { if( std::distance(begin_partons(), it_part)==0) connect_utchannel(it_part, colour, anti_colour, ran); else connect_tchannel(it_part, colour, anti_colour, ran); break; } case event_type::unof: case event_type::qqxexf: { if( std::distance(it_part, end_partons())==2) connect_utchannel(it_part, colour, anti_colour, ran); else connect_tchannel(it_part, colour, anti_colour, ran); break; } case event_type::qqxmid:{ auto it_next = std::next(it_part); if( std::distance(begin_partons(), it_part)>0 && std::distance(it_part, end_partons())>2 && ( (is_quark(*it_part) && is_antiquark(*it_next)) || (is_antiquark(*it_part) && is_quark(*it_next)) ) ) connect_utchannel(it_part, colour, anti_colour, ran); else connect_tchannel(it_part, colour, anti_colour, ran); break; } default: // rest has to be t-channel connect_tchannel(it_part, colour, anti_colour, ran); } } // Connect last connect_incoming(incoming_[1], anti_colour, colour); assert(is_leading_colour()); return true; } // generate_colours namespace { bool valid_parton( std::vector const & jets, Particle const & parton, int const idx, double const max_ext_soft_pt_fraction, double const min_extparton_pt ){ // TODO code overlap with PhaseSpacePoint::pass_extremal_cuts if(min_extparton_pt > parton.pt()) return false; if(idx<0) return false; assert(static_cast(jets.size())>=idx); auto const & jet{ jets[idx] }; return (parton.p - jet).pt()/jet.pt() <= max_ext_soft_pt_fraction; } } // namespace // this should work with multiple types bool Event::valid_hej_state(double const max_frac, double const min_pt ) const { using namespace event_type; if(!is_resummable(type())) return false; auto const & jet_idx{ particle_jet_indices() }; auto idx_begin{ jet_idx.cbegin() }; auto idx_end{ jet_idx.crbegin() }; auto part_begin{ cbegin_partons() }; auto part_end{ crbegin_partons() }; // always seperate extremal jets if( !valid_parton(jets(), *part_begin, *idx_begin, max_frac, min_pt) ) return false; ++part_begin; ++idx_begin; if( !valid_parton(jets(), *part_end, *idx_end, max_frac, min_pt) ) return false; ++part_end; ++idx_end; // unob -> second parton in own jet if( type() & (unob | qqxexb) ){ if( !valid_parton(jets(), *part_begin, *idx_begin, max_frac, min_pt) ) return false; ++part_begin; ++idx_begin; } if( type() & (unof | qqxexf) ){ if( !valid_parton(jets(), *part_end, *idx_end, max_frac, min_pt) ) return false; ++part_end; // ++idx_end; // last check, we don't need idx_end afterwards } if( type() & qqxmid ){ // find qqx pair auto begin_qqx{ std::find_if( part_begin, part_end.base(), [](Particle const & part) -> bool { return part.type != ParticleID::gluon; } )}; assert(begin_qqx != part_end.base()); long int qqx_pos{ std::distance(part_begin, begin_qqx) }; assert(qqx_pos >= 0); idx_begin+=qqx_pos; - if( !( valid_parton(jets(),*begin_qqx, *idx_begin, max_frac,min_pt) - && valid_parton(jets(),*(++begin_qqx),*(++idx_begin),max_frac,min_pt) + if( !( valid_parton(jets(), *begin_qqx, *idx_begin, + max_frac,min_pt) + && valid_parton(jets(), *std::next(begin_qqx), *std::next(idx_begin), + max_frac,min_pt) )) return false; } return true; } Event::ConstPartonIterator Event::begin_partons() const { return cbegin_partons(); } Event::ConstPartonIterator Event::cbegin_partons() const { - return boost::make_filter_iterator( - static_cast(is_parton), - cbegin(outgoing()), - cend(outgoing()) - ); + return {HEJ::is_parton, cbegin(outgoing()), cend(outgoing())}; } Event::ConstPartonIterator Event::end_partons() const { return cend_partons(); } Event::ConstPartonIterator Event::cend_partons() const { - return boost::make_filter_iterator( - static_cast(is_parton), - cend(outgoing()), - cend(outgoing()) - ); + return {HEJ::is_parton, cend(outgoing()), cend(outgoing())}; } Event::ConstReversePartonIterator Event::rbegin_partons() const { return crbegin_partons(); } Event::ConstReversePartonIterator Event::crbegin_partons() const { return std::reverse_iterator( cend_partons() ); } Event::ConstReversePartonIterator Event::rend_partons() const { return crend_partons(); } Event::ConstReversePartonIterator Event::crend_partons() const { return std::reverse_iterator( cbegin_partons() ); } Event::PartonIterator Event::begin_partons() { - return boost::make_filter_iterator( - static_cast(is_parton), - begin(outgoing_), - end(outgoing_) - ); + return {HEJ::is_parton, begin(outgoing_), end(outgoing_)}; } Event::PartonIterator Event::end_partons() { - return boost::make_filter_iterator( - static_cast(is_parton), - end(outgoing_), - end(outgoing_) - ); + return {HEJ::is_parton, end(outgoing_), end(outgoing_)}; } Event::ReversePartonIterator Event::rbegin_partons() { return std::reverse_iterator( end_partons() ); } Event::ReversePartonIterator Event::rend_partons() { return std::reverse_iterator( begin_partons() ); } namespace { void print_momentum(std::ostream & os, fastjet::PseudoJet const & part){ constexpr int prec = 6; const std::streamsize orig_prec = os.precision(); os < const & col){ constexpr int width = 3; if(!col) os << "(no color)"; // American spelling for better alignment else os << "(" <first << ", " <second << ")"; } } // namespace std::ostream& operator<<(std::ostream & os, Event const & ev){ constexpr int prec = 4; constexpr int wtype = 3; // width for types const std::streamsize orig_prec = os.precision(); os < rapidity=" < rapidity=" < rapidity=" < incoming{ event.incoming() }; // First incoming should be positive pz according to LHE standard // (or at least most (everyone?) do it this way, and Pythia assumes it) if(incoming[0].pz() < incoming[1].pz()) std::swap(incoming[0], incoming[1]); for(Particle const & in: incoming){ result.IDUP.emplace_back(in.type); result.ISTUP.emplace_back(LHE_Status::in); result.PUP.push_back({in.p[0], in.p[1], in.p[2], in.p[3], in.p.m()}); result.MOTHUP.emplace_back(0, 0); assert(in.colour); result.ICOLUP.emplace_back(*in.colour); } // outgoing for(size_t i = 0; i < event.outgoing().size(); ++i){ Particle const & out = event.outgoing()[i]; result.IDUP.emplace_back(out.type); const int status = event.decays().count(i) != 0u ?LHE_Status::decay :LHE_Status::out; result.ISTUP.emplace_back(status); result.PUP.push_back({out.p[0], out.p[1], out.p[2], out.p[3], out.p.m()}); result.MOTHUP.emplace_back(1, 2); if(out.colour) result.ICOLUP.emplace_back(*out.colour); else{ result.ICOLUP.emplace_back(std::make_pair(0,0)); } } // decays for(auto const & decay: event.decays()){ for(auto const & out: decay.second){ result.IDUP.emplace_back(out.type); result.ISTUP.emplace_back(LHE_Status::out); result.PUP.push_back({out.p[0], out.p[1], out.p[2], out.p[3], out.p.m()}); const size_t mother_idx = 1 + event.incoming().size() + decay.first; result.MOTHUP.emplace_back(mother_idx, mother_idx); result.ICOLUP.emplace_back(0,0); } } assert(result.ICOLUP.size() == num_particles); static constexpr double unknown_spin = 9.; //per Les Houches accord result.VTIMUP = std::vector(num_particles, unknown_spin); result.SPINUP = result.VTIMUP; return result; } } // namespace HEJ diff --git a/src/PDF.cc b/src/PDF.cc index f30a540..4cdd837 100644 --- a/src/PDF.cc +++ b/src/PDF.cc @@ -1,74 +1,76 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #include "HEJ/PDF.hh" #include #include #include #include #include #include "LHAPDF/LHAPDF.h" namespace HEJ { namespace { int to_beam(ParticleID id){ if(std::abs(id) == pid::proton){ return (id > 0)?1:-1; } throw std::invalid_argument( "unknown beam type: " + std::to_string(id) ); } } // namespace PDF::PDF(int id, ParticleID beam1, ParticleID beam2): pdf_{LHAPDF::mkPDF(id)}, beamtype_{{to_beam(beam1), to_beam(beam2)}} {} double PDF::pdfpt(size_t beam_idx, double x, double q, ParticleID id) const{ if(!(inRangeQ(q) && inRangeX(x))) return 0.; if(id == pid::gluon){ return pdf_->xfxQ(pid::gluon,x,q); } if(std::abs(id) <= pid::top){ return pdf_->xfxQ(id*beamtype_[beam_idx],x,q); } std::cerr << "particle type unknown: "<< id << std::endl; return 0.0; } double PDF::Halphas(double q) const{ double as = pdf_->alphasQ(q); if (std::isnan(as) || as > 0.5) { as = 0.5; } return as; } int PDF::id() const{ return pdf_->lhapdfID(); } bool PDF::inRangeQ(double q) const{ return pdf_->inRangeQ(q); } bool PDF::inRangeX(double x) const{ return pdf_->inRangeX(x); } // can't be in header since we forward declare LHAPDF::PDF PDF::~PDF() = default; + PDF::PDF(PDF && other) = default; + PDF & PDF::operator=(PDF && other) = default; } // namespace HEJ diff --git a/src/PDG_codes.cc b/src/PDG_codes.cc index 59a58d9..7c3ffe7 100644 --- a/src/PDG_codes.cc +++ b/src/PDG_codes.cc @@ -1,113 +1,117 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #include "HEJ/PDG_codes.hh" #include #include "HEJ/exceptions.hh" namespace HEJ { ParticleID to_ParticleID(std::string const & name){ using namespace HEJ::pid; static const std::map known = { {"d", d}, {"down", down}, {"1",static_cast(1)}, {"u", u}, {"up", up}, {"2",static_cast(2)}, {"s", s}, {"strange", strange}, {"3",static_cast(3)}, {"c", c}, {"charm", charm}, {"4",static_cast(4)}, {"b", b}, {"bottom", bottom}, {"5",static_cast(5)}, {"t", t}, {"top", top}, {"6",static_cast(6)}, {"e", e}, {"electron", electron}, {"e-", e}, {"11",static_cast(11)}, {"nu_e", nu_e}, {"electron_neutrino", electron_neutrino}, {"electron-neutrino", electron_neutrino}, {"12",static_cast(12)}, {"mu", mu}, {"muon", muon}, {"mu-", mu}, {"13",static_cast(13)}, {"nu_mu", nu_mu}, {"muon_neutrino", muon_neutrino}, {"muon-neutrino", muon_neutrino}, {"14",static_cast(14)}, {"tau", tau}, {"tau-", tau}, {"15",static_cast(15)}, {"nu_tau", nu_tau}, {"tau_neutrino", tau_neutrino}, {"tau-neutrino", tau_neutrino}, {"16",static_cast(16)}, {"d_bar", d_bar}, {"antidown", antidown}, {"-1",static_cast(-1)}, {"u_bar", u_bar}, {"antiup", antiup}, {"-2",static_cast(-2)}, {"s_bar", s_bar}, {"antistrange", antistrange}, {"-3",static_cast(-3)}, {"c_bar", c_bar}, {"anticharm", anticharm}, {"-4",static_cast(-4)}, {"b_bar", b_bar}, {"antibottom", antibottom}, {"-5",static_cast(-5)}, {"t_bar", t_bar}, {"antitop", antitop}, {"-6",static_cast(-6)}, {"e_bar", e_bar}, {"antielectron", antielectron}, {"positron", positron}, {"e+", e_bar}, {"-11",static_cast(-11)}, {"nu_e_bar", nu_e_bar}, {"electron_antineutrino", electron_antineutrino}, {"electron-antineutrino", electron_antineutrino}, {"-12",static_cast(-12)}, {"mu_bar", mu_bar}, {"mu+", mu_bar}, {"antimuon", antimuon}, {"-13",static_cast(-13)}, {"nu_mu_bar", nu_mu_bar}, {"muon_antineutrino", muon_antineutrino}, {"muon-antineutrino", muon_antineutrino}, {"-14",static_cast(-14)}, {"tau_bar", tau_bar}, {"tau+", tau_bar}, {"antitau", antitau}, {"-15",static_cast(-15)}, {"nu_tau_bar", nu_tau_bar}, {"tau_antineutrino", tau_antineutrino}, {"tau-antineutrino", tau_antineutrino}, {"-16",static_cast(-16)}, {"gluon", gluon}, {"g", g}, {"21",static_cast(21)}, {"photon", photon}, {"gamma", gamma}, {"22",static_cast(22)}, {"Z", Z}, {"23",static_cast(23)}, {"Z_photon_mix", Z_photon_mix}, {"Z_gamma_mix", Z_gamma_mix}, {"Z/photon superposition",Z_photon_mix}, {"81",static_cast(81)}, {"Wp", Wp}, {"W+", Wp}, {"24",static_cast(24)}, {"Wm", Wm}, {"W-", Wm}, {"-24",static_cast(-24)}, {"h", h}, {"H", h}, {"Higgs", Higgs}, {"higgs", higgs}, {"25",static_cast(25)}, {"p", p}, {"proton", proton}, {"2212",static_cast(2212)}, {"p_bar", p_bar}, {"antiproton", antiproton}, {"-2212",static_cast(-2212)} }; const auto res = known.find(name); if(res == known.end()){ throw std::invalid_argument("Unknown particle " + name); } return res->second; } namespace pid { std::string name(ParticleID id) { using namespace HEJ::pid; switch (id) { case unspecified: return "unspecified"; case down: return "down"; case up: return "up"; case strange: return "strange"; case charm: return "charm"; case bottom: return "bottom"; case top: return "top"; case electron: return "electron"; case muon: return "muon"; case tau: return "tau"; case electron_neutrino: return "electron-neutrino"; case muon_neutrino: return "muon-neutrino"; case tau_neutrino: return "tau-neutrino"; case antidown: return "antidown"; case antiup: return "antiup"; case antistrange: return "antistrange"; case anticharm: return "anticharm"; case antibottom: return "antibottom"; case antitop: return "antitop"; case positron: return "positron"; case antimuon: return "antimuon"; case antitau: return "antitau"; case electron_antineutrino: return "electron-antineutrino"; case muon_antineutrino: return "muon-antineutrino"; case tau_antineutrino: return "tau-antineutrino"; case gluon: return "gluon"; case photon: return "photon"; case Z: return "Z"; case Z_photon_mix: return "Z/photon superposition"; case Wp: return "W+"; case Wm: return "W-"; case Higgs: return "Higgs"; case proton: return "proton"; case antiproton: return "antiproton"; } throw std::logic_error{"unreachable"}; } + + ParticleID anti(ParticleID const id){ + return static_cast(-id); + } } // namespace pid } // namespace HEJ