Page MenuHomeHEPForge
Files F8723468

No OneTemporary
Actions

View Options

diff --git a/CMakeLists.txt b/CMakeLists.txt
index a84ac3c..26370a8 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,216 +1,210 @@
cmake_minimum_required(VERSION 3.1 FATAL_ERROR)
set(CMAKE_LEGACY_CYGWIN_WIN32 0)
set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
project("Reversed HEJ" VERSION 0.0.1 LANGUAGES C CXX)
## 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 -Werror")
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
)
set(INSTALL_INCLUDE_DIR "include/RHEJ/")
set(INSTALL_BIN_DIR "bin/")
set(INSTALL_LIB_DIR "lib/")
CONFIGURE_FILE( ${CMAKE_CURRENT_SOURCE_DIR}/cmake/Templates/Version.hh.in
${PROJECT_BINARY_DIR}/include/RHEJ/Version.hh @ONLY )
configure_file( ${CMAKE_CURRENT_SOURCE_DIR}/cmake/Templates/rHEJ-config.cc.in
${PROJECT_BINARY_DIR}/src/bin/rHEJ-config.cc @ONLY )
## Add directories and find dependences
include_directories(${CMAKE_CURRENT_SOURCE_DIR}/include ${PROJECT_BINARY_DIR}/include)
find_package(fastjet REQUIRED)
include_directories(${fastjet_INCLUDE_PATH})
find_package(clhep 2.3 REQUIRED)
include_directories(${clhep_INCLUDE_PATH})
find_package(lhapdf REQUIRED)
include_directories(${lhapdf_INCLUDE_PATH})
find_package(gsl REQUIRED)
include_directories(${gsl_INCLUDE_PATH})
find_package(Boost REQUIRED COMPONENTS iostreams)
include_directories(${Boost_INCLUDE_DIRS})
find_package(yaml-cpp)
include_directories(${YAML_CPP_INCLUDE_DIR})
find_package(HepMC 2)
if(${HepMC_FOUND})
message (STATUS "HepMC installation found: ${HepMC_INCLUDE_DIRS}")
set(
CMAKE_CXX_FLAGS
"${CMAKE_CXX_FLAGS} -DRHEJ_BUILD_WITH_HepMC_VERSION=${HepMC_VERSION_MAJOR}"
)
include_directories(${HepMC_INCLUDE_DIRS})
find_package(rivet)
if(${rivet_FOUND})
include_directories(${rivet_INCLUDE_PATH})
set(
CMAKE_CXX_FLAGS
"${CMAKE_CXX_FLAGS} -DRHEJ_BUILD_WITH_RIVET"
)
endif()
endif()
find_package(QCDloop 2)
if(${QCDloop_FOUND})
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DRHEJ_BUILD_WITH_QCDLOOP")
include_directories(SYSTEM ${QCDloop_INCLUDE_DIRS})
endif()
add_subdirectory(src)
## define executable
add_executable(rHEJ src/bin/rHEJ.cc)
## link libraries
target_link_libraries(rHEJ rhej)
add_executable(rHEJ-config src/bin/rHEJ-config.cc)
file(GLOB rhej_headers ${CMAKE_CURRENT_SOURCE_DIR}/include/RHEJ/*.hh ${PROJECT_BINARY_DIR}/include/RHEJ/*.hh)
file(GLOB lhef_headers ${CMAKE_CURRENT_SOURCE_DIR}/include/LHEF/*.h)
install(FILES ${rhej_headers} DESTINATION ${INSTALL_INCLUDE_DIR})
install(FILES ${lhef_headers} DESTINATION include/LHEF/)
install(TARGETS rHEJ rHEJ-config DESTINATION ${INSTALL_BIN_DIR})
## tests
enable_testing()
set(tst_dir "${CMAKE_CURRENT_SOURCE_DIR}/t")
-add_executable(test_Matrix ${tst_dir}/test_Matrix.cc)
add_executable(test_classify ${tst_dir}/test_classify.cc)
add_executable(test_psp ${tst_dir}/test_psp.cc)
add_executable(test_ME_generic ${tst_dir}/test_ME_generic.cc)
add_executable(check_res ${tst_dir}/check_res.cc)
add_executable(check_lhe ${tst_dir}/check_lhe.cc)
add_library(scales SHARED ${tst_dir}/scales.cc)
add_executable(test_scale_import ${tst_dir}/test_scale_import)
add_executable(test_descriptions ${tst_dir}/test_descriptions)
-target_link_libraries(test_Matrix rhej)
target_link_libraries(test_classify rhej)
target_link_libraries(test_psp rhej)
target_link_libraries(test_ME_generic rhej)
target_link_libraries(check_res rhej)
target_link_libraries(check_lhe rhej)
target_link_libraries(test_scale_import rhej)
target_link_libraries(test_descriptions rhej)
## add tests
add_test(
- NAME t_matrix
- COMMAND test_Matrix
- )
-add_test(
NAME t_classify
COMMAND test_classify ${tst_dir}/classify.lhe.gz
)
add_test(
NAME t_psp
COMMAND test_psp ${tst_dir}/psp_gen.lhe.gz
)
set(tst_ME_data_dir "${tst_dir}/ME_data")
add_test(
NAME t_ME_j
COMMAND test_ME_generic ${tst_ME_data_dir}/config_mtinf.yml ${tst_ME_data_dir}/ME_jets.dat ${tst_ME_data_dir}/PSP_jets.lhe.gz
)
add_test(
NAME t_ME_h
COMMAND test_ME_generic ${tst_ME_data_dir}/config_mtinf.yml ${tst_ME_data_dir}/ME_h_mtinf.dat ${tst_ME_data_dir}/PSP_h.lhe.gz
)
if(${QCDloop_FOUND})
add_test(
NAME t_ME_h_mt
COMMAND test_ME_generic ${tst_ME_data_dir}/config_mt.yml ${tst_ME_data_dir}/ME_h_mt.dat ${tst_ME_data_dir}/PSP_h.lhe.gz
)
add_test(
NAME t_ME_h_mtmb
COMMAND test_ME_generic ${tst_ME_data_dir}/config_mtmb.yml ${tst_ME_data_dir}/ME_h_mtmb.dat ${tst_ME_data_dir}/PSP_h.lhe.gz
)
endif()
add_test(
NAME t_2j
COMMAND check_res ${tst_dir}/2j.lhe.gz 3.49391e+07 419684
)
add_test(
NAME t_3j
COMMAND check_res ${tst_dir}/3j.lhe.gz 2.37902e+06 25746.6
)
add_test(
NAME t_4j
COMMAND check_res ${tst_dir}/4j.lhe.gz 603713 72822.6
)
add_test(
NAME t_h_3j
COMMAND check_res ${tst_dir}/h_3j.lhe.gz 0.821622 0.0220334
)
add_test(
NAME t_h_3j_uno
COMMAND check_res ${tst_dir}/h_3j_uno.lhe.gz 0.0261968 0.000341549 uno
)
if(${HepMC_FOUND})
file(READ "${tst_dir}/jet_config.yml" config)
file(WRITE "${tst_dir}/jet_config_withHepMC.yml" "${config} - tst.hepmc")
if(${rivet_FOUND})
file(READ "${tst_dir}/jet_config_withHepMC.yml" config)
file(WRITE "${tst_dir}/jet_config_withRivet.yml" "${config}\n\nanalysis:\n rivet: MC_XS\n output: tst")
add_test(
NAME t_main
COMMAND rHEJ ${tst_dir}/jet_config_withRivet.yml ${tst_dir}/2j.lhe.gz
)
else()
add_test(
NAME t_main
COMMAND rHEJ ${tst_dir}/jet_config_withHepMC.yml ${tst_dir}/2j.lhe.gz
)
endif()
if(${HepMC_VERSION_MAJOR} GREATER 2)
add_executable(check_hepmc ${tst_dir}/check_hepmc.cc)
target_link_libraries(check_hepmc rhej)
add_test(
NAME t_hepmc
COMMAND check_hepmc tst.hepmc
)
endif()
else()
add_test(
NAME t_main
COMMAND rHEJ ${tst_dir}/jet_config.yml ${tst_dir}/2j.lhe.gz
)
endif()
add_test(
NAME t_lhe
COMMAND check_lhe tst.lhe
)
add_test(
NAME t_scale_import
COMMAND test_scale_import ${tst_dir}/jet_config_with_import.yml
)
add_test(
NAME t_descriptions
COMMAND test_descriptions
)
diff --git a/include/RHEJ/Matrix.hh b/include/RHEJ/Matrix.hh
deleted file mode 100644
index 306da42..0000000
--- a/include/RHEJ/Matrix.hh
+++ /dev/null
@@ -1,57 +0,0 @@
-/** \file
- * \brief This header file contains the Matrix class
- */
-
-#pragma once
-
-#include <cstddef> // needed for size_t
-#include <string>
-
-#include <gsl/gsl_matrix.h>
-
-namespace RHEJ{
-
-
- //! Class for matrices
- class Matrix{
- public:
- Matrix() = default;
- //! Constructor
- /**
- * @param rows Number of rows
- * @param columns Number of columns
- */
- Matrix(size_t rows, size_t columns);
- //! Access element in given row and column
- double & operator()(size_t row, size_t column);
- //! Access element in given row and column
- double const & operator()(size_t row, size_t column) const;
-
- //! Get number of rows
- size_t rows() const;
- //! Get number of columns
- size_t columns() const;
-
- Matrix(Matrix const & other);
- Matrix(Matrix && other) noexcept;
-
- Matrix & operator=(Matrix const & other);
- Matrix & operator=(Matrix && other) noexcept;
-
- ~Matrix();
- private:
- gsl_matrix* m;
-
- friend double det(Matrix m);
- friend void swap(Matrix & a, Matrix & b) noexcept;
- };
-
- //! Swap two matrices
- void swap(Matrix & a, Matrix & b) noexcept;
-
- //! Compute the determinant
- double det(Matrix m);
-
- //! Convert a matrix to a string
- std::string to_string(Matrix const & m);
-}
diff --git a/src/Jacobian.cc b/src/Jacobian.cc
index 835265a..1fa7d7f 100644
--- a/src/Jacobian.cc
+++ b/src/Jacobian.cc
@@ -1,44 +1,71 @@
#include "RHEJ/Jacobian.hh"
-#include "RHEJ/Matrix.hh"
+#include <boost/numeric/ublas/matrix.hpp>
+#include <boost/numeric/ublas/lu.hpp>
namespace RHEJ{
namespace{
enum coordinates : size_t {
x1, x2
};
+
+
+ namespace ublas = boost::numeric::ublas;
+
+ template<class Matrix>
+ double det(ublas::matrix_expression<Matrix> const& m) {
+
+ ublas::permutation_matrix<size_t> pivots{m().size1()};
+ Matrix mLu{m()};
+
+ const auto is_singular = lu_factorize(mLu, pivots);
+
+ if(is_singular) return 0.;
+
+ double det = 1.0;
+ for (std::size_t i = 0; i < pivots.size(); ++i){
+ if (pivots(i) != i) det = -det;
+
+ det *= mLu(i,i);
+ }
+
+ return det;
+ }
+
+ using ublas::matrix;
}
double Jacobian(
std::vector<fastjet::PseudoJet> const & jets,
fastjet::PseudoJet const & q
) {
+
static constexpr int num_coordinates = 2;
- auto Jacobian = Matrix{
+ auto Jacobian = matrix<double>{
num_coordinates*jets.size(),
num_coordinates*jets.size()
};
double P_perp = 0.;
for(auto const & J: jets) P_perp += J.perp();
for(size_t l = 0; l < jets.size(); ++l){
const double Jl_perp = jets[l].perp();
for(size_t lp = 0; lp < jets.size(); ++lp){
const int delta_l = l == lp;
const double Jlp_perp = jets[lp].perp();
for(size_t x = x1; x <= x2; ++x){
for(size_t xp = x1; xp <= x2; ++xp){
const int delta_x = x == xp;
Jacobian(2*l + x, 2*lp + xp) =
+ delta_l*delta_x
- q[x]*jets[lp][xp]/(P_perp*Jlp_perp)*(
+ delta_l - Jl_perp/P_perp
);
}
}
}
}
return det(Jacobian);
}
}
diff --git a/src/Matrix.cc b/src/Matrix.cc
deleted file mode 100644
index 2833975..0000000
--- a/src/Matrix.cc
+++ /dev/null
@@ -1,91 +0,0 @@
-#include "RHEJ/Matrix.hh"
-
-#include <gsl/gsl_linalg.h>
-
-namespace RHEJ{
-
- Matrix::Matrix(size_t rows, size_t columns):
- m(gsl_matrix_alloc(rows, columns))
- {}
-
- size_t Matrix::rows() const{
- return m->size1;
- }
-
- size_t Matrix::columns() const{
- return m->size2;
- }
-
- Matrix::Matrix(Matrix const & other):
- m(gsl_matrix_alloc(other.rows(), other.columns()))
- {
- gsl_matrix_memcpy(m, other.m);
- }
-
- Matrix::Matrix(Matrix && other) noexcept:
- m(other.m)
- {
- other.m = gsl_matrix_alloc(0, 0);
- }
-
- double & Matrix::operator()(size_t row, size_t column){
- return *gsl_matrix_ptr(m, row, column);
- }
-
- double const & Matrix::operator()(size_t row, size_t column) const{
- return *gsl_matrix_const_ptr(m, row, column);
- }
-
- Matrix & Matrix::operator=(Matrix const & other){
- Matrix o{other};
- swap(*this, o);
- return *this;
- }
-
- Matrix & Matrix::operator=(Matrix && other) noexcept{
- std::swap(m, other.m);
- return *this;
- }
-
- Matrix::~Matrix(){
- gsl_matrix_free(m);
- }
-
- void swap(Matrix & a, Matrix & b) noexcept{
- std::swap(a.m, b.m);
- }
-
- namespace{
- struct permutation{
- gsl_permutation* p;
- permutation(size_t size):
- p{gsl_permutation_alloc(size)}
- {}
- permutation(permutation const & other) = delete;
- void operator=(permutation const & other) = delete;
-
- ~permutation(){
- gsl_permutation_free(p);
- }
- };
- }
-
- double det(Matrix m){
- int signum;
- permutation p(m.rows());
- gsl_linalg_LU_decomp(m.m , p.p , &signum);
- return gsl_linalg_LU_det(m.m, signum);
- }
-
- std::string to_string(Matrix const & m){
- std::string result;
- for(size_t i = 0; i < m.rows(); ++i){
- for(size_t j = 0; j < m.columns(); ++j){
- result += std::to_string(m(i, j));
- result += " ";
- }
- result.back() = '\n';
- }
- return result;
- }
-}
diff --git a/t/test_Matrix.cc b/t/test_Matrix.cc
deleted file mode 100644
index 76fff71..0000000
--- a/t/test_Matrix.cc
+++ /dev/null
@@ -1,47 +0,0 @@
-#ifdef NDEBUG
-#undef NDEBUG
-#endif
-
-#include <cassert>
-#include <cmath>
-#include <iostream>
-
-#include "RHEJ/Matrix.hh"
-
-int main(){
- constexpr size_t dim = 4;
- constexpr double ep = 1e-16;
-
- auto m = RHEJ::Matrix{dim, dim};
- assert(m.rows() == dim);
- assert(m.columns() == dim);
- m(0,0) = 0.35388638718713294;
- assert(m(0,0) == 0.35388638718713294);
- m(0,1) = 0.441472069363144;
- m(0,2) = 0.7821385102801595;
- m(0,3) = 0.4506533375385655;
- m(1,0) = 0.5954699211986723;
- m(1,1) = 0.7337301649423917;
- m(1,2) = 0.24953309317812855;
- m(1,3) = 0.5747608628993972;
- m(2,0) = 0.5241499726390944;
- m(2,1) = 0.5202946315779862;
- m(2,2) = 0.5178464229701134;
- m(2,3) = 0.22486555579427625;
- m(3,0) = 0.28294299023923397;
- m(3,1) = 0.23061982229765166;
- m(3,2) = 0.9472220498156341;
- m(3,3) = 0.68515373407338;
- auto mcpy = m;
- assert(mcpy.rows() == m.rows());
- assert(mcpy.columns() == m.columns());
- assert(mcpy(0,0) == 0.35388638718713294);
- m(0,0) = 0;
- assert(m(0,0) == 0);
- assert(mcpy(0,0) == 0.35388638718713294);
- m = std::move(mcpy);
- assert(m.rows() == dim);
- assert(m.columns() == dim);
- assert(m(0,0) == 0.35388638718713294);
- assert(std::abs(det(m) + 0.028722124827027316) < ep);
-}

File Metadata

Mime Type
text/x-diff
Expires
Mon, Jan 20, 8:37 PM (16 h, 26 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4242264
Default Alt Text
(13 KB)

Event Timeline