Page MenuHomeHEPForge
Files F7877704

No OneTemporary
Actions

View Options

diff --git a/FixedOrderGen/CMakeLists.txt b/FixedOrderGen/CMakeLists.txt
index e1d2529..ef16662 100644
--- a/FixedOrderGen/CMakeLists.txt
+++ b/FixedOrderGen/CMakeLists.txt
@@ -1,113 +1,113 @@
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.5 LANGUAGES C CXX)
# Set a default build type if none was specified
set(default_build_type "RelWithDebInfo")
if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES)
message(STATUS "Setting build type to '${default_build_type}' as none was specified.")
set(CMAKE_BUILD_TYPE "${default_build_type}" CACHE
STRING "Choose the type of build." FORCE)
# Set the possible values of build type for cmake-gui
set_property(CACHE CMAKE_BUILD_TYPE PROPERTY STRINGS
"Debug" "Release" "MinSizeRel" "RelWithDebInfo")
endif()
## Flags for the compiler. No warning allowed.
if (CMAKE_COMPILER_IS_GNUCC OR CMAKE_COMPILER_IS_GNUCXX)
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wextra")
elseif (MSVC)
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /W4 /WX /EHsc")
endif()
set(CMAKE_CXX_STANDARD_REQUIRED ON)
set(CMAKE_CXX_STANDARD 14)
## Create Version
set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}/cmake/Modules/")
# Get the latest abbreviated commit hash of the working branch
execute_process(
COMMAND git rev-parse HEAD
WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
OUTPUT_VARIABLE PROJECT_GIT_REVISION
OUTPUT_STRIP_TRAILING_WHITESPACE
)
# Get the current working branch
execute_process(
COMMAND git rev-parse --abbrev-ref HEAD
WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
OUTPUT_VARIABLE PROJECT_GIT_BRANCH
OUTPUT_STRIP_TRAILING_WHITESPACE
)
CONFIGURE_FILE( ${CMAKE_CURRENT_SOURCE_DIR}/cmake/Templates/Version.hh.in
${PROJECT_BINARY_DIR}/include/Version.hh @ONLY )
## Add directories and find dependences
include_directories(${CMAKE_CURRENT_SOURCE_DIR}/include ${PROJECT_BINARY_DIR}/include)
set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}/../cmake/Modules/")
find_package(HEJ 2 REQUIRED)
include_directories(${HEJ_INCLUDE_DIR})
find_package(fastjet REQUIRED)
include_directories(${fastjet_INCLUDE_DIRS})
find_package(CLHEP 2.3 REQUIRED)
include_directories(${CLHEP_INCLUDE_DIRS})
find_package(LHAPDF REQUIRED)
include_directories(${LHAPDF_INCLUDE_DIRS})
## Amend unintuitive behaviour of FindBoost.cmake
## Priority of BOOST_ROOT over BOOSTROOT matches FindBoost.cmake
## at least for cmake 3.12
if(DEFINED BOOST_ROOT)
message("BOOST_ROOT set - only looking for Boost in ${BOOST_ROOT}")
set(Boost_NO_BOOST_CMAKE ON)
elseif(DEFINED BOOSTROOT)
message("BOOSTROOT set - only looking for Boost in ${BOOSTROOT}")
set(Boost_NO_BOOST_CMAKE ON)
endif()
find_package(Boost REQUIRED COMPONENTS iostreams)
include_directories(${Boost_INCLUDE_DIRS})
find_package(yaml-cpp) # requiring yaml does not work with fedora
include_directories(${YAML_CPP_INCLUDE_DIR})
## 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 STATIC ${HEJFOG_source})
add_executable(HEJFOG ${CMAKE_CURRENT_SOURCE_DIR}/src/main.cc)
## link libraries
set(libraries ${CMAKE_DL_LIBS} ${LHAPDF_LIBRARIES} ${CLHEP_LIBRARIES}
${fastjet_LIBRARIES} ${Boost_LIBRARIES} ${YAML_CPP_LIBRARIES} ${HEJ_LIBRARIES})
target_link_libraries(hejfog ${libraries})
target_link_libraries(HEJFOG hejfog)
install(TARGETS HEJFOG DESTINATION bin)
## tests
enable_testing()
set(tst_dir "${CMAKE_CURRENT_SOURCE_DIR}/t")
-foreach(tst h_2j h_3j h_5j h_3j_uno1 h_3j_uno2 h_2j_decay 2j 4j)# W_2j_classify W_nj_classify)
+foreach(tst h_2j h_3j h_5j h_3j_uno1 h_3j_uno2 h_2j_decay 2j 4j W_2j_classify W_nj_classify)
add_executable(test_${tst} ${tst_dir}/${tst}.cc)
target_link_libraries(test_${tst} hejfog)
add_test(NAME t_${tst} COMMAND test_${tst} WORKING_DIRECTORY ${tst_dir})
endforeach()
add_test(
NAME t_main_2j
COMMAND HEJFOG ${tst_dir}/config_2j.yml
)
add_test(
NAME t_main_h2j
COMMAND HEJFOG ${tst_dir}/config_h_2j.yml
)
add_test(
NAME t_main_h2j_decay
COMMAND HEJFOG ${tst_dir}/config_h_2j_decay.yml
)
diff --git a/FixedOrderGen/t/W_2j_classify.cc b/FixedOrderGen/t/W_2j_classify.cc
index 5dfc5c8..bf1d22b 100644
--- a/FixedOrderGen/t/W_2j_classify.cc
+++ b/FixedOrderGen/t/W_2j_classify.cc
@@ -1,152 +1,152 @@
/**
* \brief check that the PSP generates only "valid" W + 2 jets events
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#ifdef NDEBUG
#undef NDEBUG
#endif
#include "JetParameters.hh"
#include "ParticleProperties.hh"
#include "PhaseSpacePoint.hh"
#include "Process.hh"
#include "Subleading.hh"
#include "HEJ/Mixmax.hh"
#include "HEJ/PDF.hh"
#include "HEJ/utility.hh"
using namespace HEJFOG;
using namespace HEJ;
namespace {
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) && !(abs(part.type)%2);
}
bool is_down_type(Particle const & part){
return HEJ::is_anyquark(part) && abs(part.type)%2;
}
bool check_W2j(PhaseSpacePoint const & psp, ParticleID const W_type){
bool found_quark = false;
bool found_anti = false;
std::vector<Particle> out_partons;
std::vector<Particle> 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(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(abs(in.type)>4 || 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 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;
const ParticlesPropMap boson_prop{
{pid::Wp, {91.1876, 2.085, {Decay{ {pid::e_bar, pid::nu_e}, 1.}} }},
{pid::Wm, {91.1876, 2.085, {Decay{ {pid::e, pid::nu_e_bar}, 1.}} }}
};
HEJ::Mixmax ran{};
// Wp2j
- Process proc {{pid::proton,pid::proton}, 2, pid::Wp};
+ Process proc {{pid::proton,pid::proton}, 2, pid::Wp, {}};
size_t n_psp = n_psp_base;
for( size_t i = 0; i<n_psp; ++i){
const PhaseSpacePoint psp{proc,jet_para,pdf,E_cms, subl_change,subl_channels,
boson_prop, ran};
if(psp.status()==good){
check_W2j(psp, *proc.boson);
} else { // bad process -> 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};
+ proc = Process{{pid::proton,pid::proton}, 2, pid::Wm, {}};
n_psp = n_psp_base;
for( size_t i = 0; i<n_psp; ++i){
const PhaseSpacePoint psp{proc,jet_para,pdf,E_cms, subl_change,subl_channels,
boson_prop, ran};
if(psp.status()==good){
check_W2j(psp, *proc.boson);
} else { // bad process -> 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 9e832ab..6844405 100644
--- a/FixedOrderGen/t/W_nj_classify.cc
+++ b/FixedOrderGen/t/W_nj_classify.cc
@@ -1,184 +1,184 @@
/**
* \brief check that the PSP generates the all W+jet subleading processes
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#ifdef NDEBUG
#undef NDEBUG
#endif
#include <algorithm>
#include "JetParameters.hh"
#include "ParticleProperties.hh"
#include "PhaseSpacePoint.hh"
#include "Process.hh"
#include "Subleading.hh"
#include "HEJ/Event.hh"
#include "HEJ/Mixmax.hh"
#include "HEJ/PDF.hh"
#include "HEJ/utility.hh"
using namespace HEJFOG;
using namespace HEJ;
namespace {
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};
}
}
int main(){
constexpr 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;
const ParticlesPropMap boson_prop{
{pid::Wp, {91.1876, 2.085, {Decay{ {pid::e_bar, pid::nu_e}, 1.}} }},
{pid::Wm, {91.1876, 2.085, {Decay{ {pid::e, pid::nu_e_bar}, 1.}} }}
};
HEJ::Mixmax ran{};
auto subl_channels = Subleading::all;
std::vector<event_type::EventType> 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};
+ Process proc {{pid::proton,pid::proton}, 3, pid::Wp, {}};
size_t n_psp = n_psp_base;
std::unordered_map<event_type::EventType, size_t> type_counter;
for( size_t i = 0; i<n_psp; ++i){
const PhaseSpacePoint psp{proc,jet_para,pdf,E_cms, subl_change,subl_channels,
boson_prop, ran};
if(psp.status()==good){
const Event ev{ to_EventData(psp).cluster(jet_para.def, jet_para.min_pt) };
++type_counter[ev.type()];
if( std::find(allowed_types.cbegin(), allowed_types.cend(), ev.type())
== allowed_types.cend()) {
bail_out(psp, "Found not allowed event of type "
+std::string(event_type::names[ev.type()]));
}
} else { // bad process -> 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<double>(entry.second)/n_psp_base;
const int percent = std::round(100*fraction);
std::cout << std::left << std::setw(25)
<< (event_type::names[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 " << event_type::names[t] << std::endl;
return EXIT_FAILURE;
}
}
// Wm3j - only uno
- proc = Process{{pid::proton,pid::proton}, 3, pid::Wm};
+ proc = Process{{pid::proton,pid::proton}, 3, pid::Wm, {}};
n_psp = n_psp_base;
subl_channels = Subleading::uno;
allowed_types = {event_type::FKL, event_type::unob, event_type::unof};
type_counter.clear();
for( size_t i = 0; i<n_psp; ++i){
const PhaseSpacePoint psp{proc,jet_para,pdf,E_cms, subl_change,subl_channels,
boson_prop, ran};
if(psp.status()==good){
const Event ev{ to_EventData(psp).cluster(jet_para.def, jet_para.min_pt) };
++type_counter[ev.type()];
if( std::find(allowed_types.cbegin(), allowed_types.cend(), ev.type())
== allowed_types.cend()) {
bail_out(psp, "Found not allowed event of type "
+std::string(event_type::names[ev.type()]));
}
} else { // bad process -> 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<double>(entry.second)/n_psp_base;
const int percent = std::round(100*fraction);
std::cout << std::left << std::setw(25)
<< (event_type::names[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 " << event_type::names[t] << std::endl;
return EXIT_FAILURE;
}
}
// Wm4j
- proc = Process{{pid::proton,pid::proton}, 4, pid::Wm};
+ proc = Process{{pid::proton,pid::proton}, 4, pid::Wm, {}};
n_psp = n_psp_base;
subl_channels = Subleading::all;
allowed_types = {event_type::FKL,
event_type::unob, event_type::unof, event_type::qqxexb, event_type::qqxexf,
event_type::qqxmid};
type_counter.clear();
for( size_t i = 0; i<n_psp; ++i){
const PhaseSpacePoint psp{proc,jet_para,pdf,E_cms, subl_change,subl_channels,
boson_prop, ran};
if(psp.status()==good){
const Event ev{ to_EventData(psp).cluster(jet_para.def, jet_para.min_pt)};
++type_counter[ev.type()];
if( std::find(allowed_types.cbegin(), allowed_types.cend(), ev.type())
== allowed_types.cend()) {
bail_out(psp, "Found not allowed event of type "
+std::string(event_type::names[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<double>(entry.second)/n_psp_base;
const int percent = std::round(100*fraction);
std::cout << std::left << std::setw(25)
<< (event_type::names[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 " << event_type::names[t] << std::endl;
return EXIT_FAILURE;
}
}
std::cout << "All processes passed." << std::endl;
return EXIT_SUCCESS;
}

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 4:17 PM (1 d, 15 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3796898
Default Alt Text
(16 KB)

Event Timeline