Page Menu
Home
HEPForge
Search
Configure Global Search
Log In
Files
F8724096
No One
Temporary
Actions
View File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Flag For Later
Size
199 KB
Subscribers
None
View Options
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 48df8dd..a84ac3c 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,212 +1,216 @@
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
- COMMAND test_ME_generic ${tst_ME_data_dir}/config_mtinf.yml ${tst_ME_data_dir}/ME_mtinf.dat ${tst_ME_data_dir}/PSP.lhe.gz
+ 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_mt
- COMMAND test_ME_generic ${tst_ME_data_dir}/config_mt.yml ${tst_ME_data_dir}/ME_mt.dat ${tst_ME_data_dir}/PSP.lhe.gz
+ 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_mtmb
- COMMAND test_ME_generic ${tst_ME_data_dir}/config_mtmb.yml ${tst_ME_data_dir}/ME_mtmb.dat ${tst_ME_data_dir}/PSP.lhe.gz
+ 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
- ) # TODO replace this by fixed pure Jet result, current value not validated
+ )
add_test(
NAME t_3j
COMMAND check_res ${tst_dir}/3j.lhe.gz 2.37902e+06 25746.6
- ) # TODO replace this by fixed pure Jet result, current value not validated
+ )
add_test(
NAME t_4j
COMMAND check_res ${tst_dir}/4j.lhe.gz 603713 72822.6
- ) # TODO replace this by fixed pure Jet result, current value not validated
+ )
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/FixedOrderGen/CMakeLists.txt b/FixedOrderGen/CMakeLists.txt
index 9ce16dd..e1962f9 100644
--- a/FixedOrderGen/CMakeLists.txt
+++ b/FixedOrderGen/CMakeLists.txt
@@ -1,86 +1,86 @@
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 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
)
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
${CMAKE_CURRENT_SOURCE_DIR}/../include)
set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}/../cmake/Modules/")
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(Boost REQUIRED COMPONENTS iostreams)
include_directories(${Boost_INCLUDE_DIRS})
find_package(yaml-cpp)
include_directories(${YAML_CPP_INCLUDE_DIR})
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()
## define executable
file(GLOB FOgen_source ${CMAKE_CURRENT_SOURCE_DIR}/src/*.cc)
list(REMOVE_ITEM FOgen_source ${CMAKE_CURRENT_SOURCE_DIR}/src/main.cc)
add_library(hejfog STATIC ${FOgen_source})
add_executable(FOgen ${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})
if(${QCDloop_FOUND})
list(APPEND libraries ${QCDloop_LIBRARIES} quadmath)
endif()
# add libraries for reversed HEJ <by hand>
list(APPEND libraries rhej)
target_link_libraries(hejfog ${libraries})
target_link_libraries(FOgen hejfog)
install(TARGETS FOgen DESTINATION bin)
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)
+foreach(tst h_2j h_3j h_5j h_3j_uno1 h_3j_uno2 h_2j_decay 2j 4j)
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()
diff --git a/FixedOrderGen/include/PhaseSpacePoint.hh b/FixedOrderGen/include/PhaseSpacePoint.hh
index fadd10b..1471f4c 100644
--- a/FixedOrderGen/include/PhaseSpacePoint.hh
+++ b/FixedOrderGen/include/PhaseSpacePoint.hh
@@ -1,160 +1,164 @@
/** \file PhaseSpacePoint.hh
* \brief Contains the PhaseSpacePoint Class
*/
#pragma once
#include <vector>
#include "RHEJ/utility.hh"
#include "RHEJ/Event.hh"
#include "RHEJ/PDG_codes.hh"
#include "RHEJ/PDF.hh"
#include "RHEJ/RNG.hh"
#include "Status.hh"
#include "JetParameters.hh"
#include "HiggsProperties.hh"
namespace HEJFOG{
class Process;
using RHEJ::Particle;
//! A point in resummation phase space
class PhaseSpacePoint{
public:
//! Default PhaseSpacePoint Constructor
PhaseSpacePoint() = default;
//! PhaseSpacePoint Constructor
/**
* @param proc The process to generate
* @param jet_def The Jet Definition Used
* @param jetptmin Minimum Jet Transverse Momentum
* @param rapmax Maximum parton rapidity
* @param pdf The pdf set (used for sampling)
* @param uno_chance Chance to turn a potentially unordered
* emission into an actual one
*
* Initially, only FKL phase space points are generated. If the most
* extremal emission in any direction is a quark or anti-quark and the
* next emission is a gluon, uno_chance is the chance to turn this into
* an unordered emission event by swapping the two flavours. At most one
* unordered emission will be generated in this way.
*/
PhaseSpacePoint(
Process const & proc,
JetParameters const & jet_properties,
RHEJ::PDF & pdf, double E_beam,
double uno_chance,
HiggsProperties const & higgs_properties,
RHEJ::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<Particle, 2> const & incoming() const{
return incoming_;
}
//! Get Outgoing Function
/**
* @returns Outgoing Particles
*/
std::vector<Particle> const & outgoing() const{
return outgoing_;
}
std::unordered_map<int, std::vector<Particle>> const & decays() const{
return decays_;
}
private:
std::vector<fastjet::PseudoJet> gen_LO_partons(
- int count,
+ int count, bool is_pure_jets,
JetParameters const & jet_param,
double max_pt,
RHEJ::RNG & ran
);
Particle gen_boson(
RHEJ::ParticleID bosonid, double mass, double width,
RHEJ::RNG & ran
);
- fastjet::PseudoJet gen_last_momentum(double mass_square, double y) const;
+ template<class ParticleMomenta>
+ fastjet::PseudoJet gen_last_momentum(
+ ParticleMomenta const & other_momenta,
+ double mass_square, double y
+ ) const;
bool jets_ok(
std::vector<fastjet::PseudoJet> const & Born_jets,
std::vector<fastjet::PseudoJet> const & partons
) const;
void reconstruct_incoming(
std::array<RHEJ::ParticleID, 2> const & ids,
RHEJ::PDF & pdf, double E_beam,
double uf,
RHEJ::RNG & ran
);
RHEJ::ParticleID generate_incoming_id(
double x, double uf, RHEJ::PDF & pdf,
RHEJ::RNG & ran
);
bool momentum_conserved(double ep) const;
RHEJ::Particle const & most_backward_FKL(
std::vector<RHEJ::Particle> const & partons
) const;
RHEJ::Particle const & most_forward_FKL(
std::vector<RHEJ::Particle> const & partons
) const;
RHEJ::Particle & most_backward_FKL(std::vector<RHEJ::Particle> & partons) const;
RHEJ::Particle & most_forward_FKL(std::vector<RHEJ::Particle> & partons) const;
bool extremal_FKL_ok(
std::vector<fastjet::PseudoJet> const & partons
) const;
double random_normal(double stddev, RHEJ::RNG & ran);
void maybe_turn_to_uno(double chance, RHEJ::RNG & ran);
std::vector<Particle> decay_boson(
RHEJ::Particle const & parent,
std::vector<Decay> const & decays,
RHEJ::RNG & ran
);
Decay select_decay_channel(
std::vector<Decay> const & decays,
RHEJ::RNG & ran
);
double gen_hard_pt(
int np, double ptmin, double ptmax, double y,
RHEJ::RNG & ran
);
double gen_soft_pt(int np, double ptmax, RHEJ::RNG & ran);
double gen_parton_pt(
int count, JetParameters const & jet_param, double ptmax, double y,
RHEJ::RNG & ran
);
double weight_;
Status status_;
std::array<Particle, 2> incoming_;
std::vector<Particle> outgoing_;
//! Particle decays in the format {outgoing index, decay products}
std::unordered_map<int, std::vector<Particle>> decays_;
};
RHEJ::UnclusteredEvent to_UnclusteredEvent(PhaseSpacePoint const & psp);
}
diff --git a/FixedOrderGen/src/PhaseSpacePoint.cc b/FixedOrderGen/src/PhaseSpacePoint.cc
index 0e31d49..a2af7e4 100644
--- a/FixedOrderGen/src/PhaseSpacePoint.cc
+++ b/FixedOrderGen/src/PhaseSpacePoint.cc
@@ -1,476 +1,489 @@
#include "PhaseSpacePoint.hh"
#include <random>
#include "RHEJ/kinematics.hh"
#include "RHEJ/utility.hh"
#include "RHEJ/exceptions.hh"
#include "Process.hh"
#include <CLHEP/Random/Randomize.h>
#include <CLHEP/Random/RanluxEngine.h>
using namespace RHEJ;
namespace HEJFOG{
static_assert(
std::numeric_limits<double>::has_quiet_NaN,
"no quiet NaN for double"
);
constexpr double NaN = std::numeric_limits<double>::quiet_NaN();
RHEJ::UnclusteredEvent to_UnclusteredEvent(PhaseSpacePoint const & psp){
RHEJ::UnclusteredEvent result;
result.incoming = psp.incoming();
std::sort(
begin(result.incoming), end(result.incoming),
[](Particle o1, Particle o2){return o1.p.pz()<o2.p.pz();}
);
assert(result.incoming.size() == 2);
result.outgoing = psp.outgoing();
assert(
std::is_sorted(
begin(result.outgoing), end(result.outgoing),
RHEJ::rapidity_less{}
)
);
assert(result.outgoing.size() >= 2);
result.decays = psp.decays();
result.central.mur = NaN;
result.central.muf = NaN;
result.central.weight = psp.weight();
return result;
}
namespace{
bool can_swap_to_uno(
RHEJ::Particle const & p1, RHEJ::Particle const & p2
){
return is_parton(p1)
&& p1.type != RHEJ::pid::gluon
&& p2.type == RHEJ::pid::gluon;
}
}
void PhaseSpacePoint::maybe_turn_to_uno(
double chance,
RHEJ::RNG & ran
){
assert(outgoing_.size() >= 2);
const size_t nout = outgoing_.size();
const bool can_be_uno_backward = can_swap_to_uno(
outgoing_[0], outgoing_[1]
);
const bool can_be_uno_forward = can_swap_to_uno(
outgoing_[nout-1], outgoing_[nout-2]
);
if(!can_be_uno_backward && !can_be_uno_forward) return;
if(ran.flat() < chance){
weight_ /= chance;
if(can_be_uno_backward && can_be_uno_forward){
if(ran.flat() < 0.5){
std::swap(outgoing_[0].type, outgoing_[1].type);
}
else{
std::swap(outgoing_[nout-1].type, outgoing_[nout-2].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);
}
}
else weight_ /= 1 - chance;
}
+ template<class ParticleMomenta>
+ fastjet::PseudoJet PhaseSpacePoint::gen_last_momentum(
+ ParticleMomenta const & other_momenta,
+ const double mass_square, const double y
+ ) const {
+ std::array<double,2> pt{0.,0.};
+ for (auto const & p: other_momenta) {
+ pt[0]-= p.px();
+ pt[1]-= p.py();
+ }
+
+ const double mperp = sqrt(pt[0]*pt[0]+pt[1]*pt[1]+mass_square);
+ const double pz=mperp*sinh(y);
+ const double E=mperp*cosh(y);
+
+ return {pt[0], pt[1], pz, E};
+ }
+
PhaseSpacePoint::PhaseSpacePoint(
Process const & proc,
JetParameters const & jet_param,
RHEJ::PDF & pdf, double E_beam,
double uno_chance,
HiggsProperties const & h,
RHEJ::RNG & ran
)
{
assert(proc.njets >= 2);
+ status_ = good;
+ weight_ = 1;
+
const int nout = proc.njets + (proc.boson?1:0);
- status_ = good;
- weight_ = 1;
- weight_ /= std::tgamma(nout);
-
- outgoing_.reserve(nout);
-
- for(auto&& p_out: gen_LO_partons(proc.njets, jet_param, E_beam, ran)){
- outgoing_.emplace_back(Particle{pid::gluon, std::move(p_out)});
- }
- if(status_ != good) return;
-
- if(proc.boson && *proc.boson == pid::Higgs){
- // The Higgs
- auto Hparticle=gen_boson(pid::higgs, h.mass, h.width, ran);
- auto pos=std::upper_bound(
- begin(outgoing_),end(outgoing_),Hparticle,rapidity_less{}
- );
- outgoing_.insert(pos, Hparticle);
- if(! h.decays.empty()){
- const int boson_idx = std::distance(begin(outgoing_), pos);
- decays_.emplace(
- boson_idx,
- decay_boson(outgoing_[boson_idx], h.decays, ran)
- );
- }
- }
- // normalisation of momentum-conserving delta function
- weight_ *= pow(2*M_PI, 4);
-
- reconstruct_incoming(proc.incoming, pdf, E_beam, jet_param.min_pt, ran);
- if(status_ != good) return;
- // set outgoing states
- most_backward_FKL(outgoing_).type = incoming_[0].type;
- most_forward_FKL(outgoing_).type = incoming_[1].type;
-
- if(proc.njets > 2) maybe_turn_to_uno(uno_chance, ran);
+ outgoing_.reserve(nout);
+
+ const bool is_pure_jets = !proc.boson;
+ auto partons = gen_LO_partons(
+ proc.njets, is_pure_jets, jet_param, E_beam, ran
+ );
+ for(auto&& p_out: partons) {
+ outgoing_.emplace_back(Particle{pid::gluon, std::move(p_out)});
+ }
+ if(status_ != good) return;
+
+ if(proc.boson && *proc.boson == pid::Higgs){
+ // The Higgs
+ auto Hparticle=gen_boson(pid::higgs, h.mass, h.width, ran);
+ auto pos=std::upper_bound(
+ begin(outgoing_),end(outgoing_),Hparticle,rapidity_less{}
+ );
+ outgoing_.insert(pos, Hparticle);
+ if(! h.decays.empty()){
+ const int boson_idx = std::distance(begin(outgoing_), pos);
+ decays_.emplace(
+ boson_idx,
+ decay_boson(outgoing_[boson_idx], h.decays, ran)
+ );
+ }
+ }
+ // normalisation of momentum-conserving delta function
+ weight_ *= pow(2*M_PI, 4);
+
+ reconstruct_incoming(proc.incoming, pdf, E_beam, jet_param.min_pt, ran);
+ if(status_ != good) return;
+ // set outgoing states
+ most_backward_FKL(outgoing_).type = incoming_[0].type;
+ most_forward_FKL(outgoing_).type = incoming_[1].type;
+
+ if(proc.njets > 2) maybe_turn_to_uno(uno_chance, ran);
}
double PhaseSpacePoint::gen_hard_pt(
int np , double ptmin, double ptmax, double y,
RHEJ::RNG & ran
) {
// heuristic parameters for pt sampling
const double ptpar = ptmin + np/5.;
const double arg_small_y = atan((ptmax - ptmin)/ptpar);
const double y_cut = 3.;
const double r1 = ran.flat();
if(y < y_cut){
const double pt = ptmin + ptpar*tan(r1*arg_small_y);
const double temp = cos(r1*arg_small_y);
weight_ *= pt*ptpar*arg_small_y/(temp*temp);
return pt;
}
const double ptpar2 = ptpar/(1 + 5*(y-y_cut));
const double temp = 1. - std::exp((ptmin-ptmax)/ptpar2);
const double pt = ptmin - ptpar2*std::log(1-r1*temp);
weight_ *= pt*ptpar2*temp/(1-r1*temp);
return pt;
}
double PhaseSpacePoint::gen_soft_pt(int np, double max_pt, RHEJ::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,
RHEJ::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;
return jet_param.min_pt;
}
return pt;
}
std::vector<fastjet::PseudoJet> PhaseSpacePoint::gen_LO_partons(
- int np ,
+ int np, bool is_pure_jets,
JetParameters const & jet_param,
double max_pt,
RHEJ::RNG & ran
){
if (np<2) throw std::invalid_argument{"Not enough partons in gen_LO_partons"};
weight_ /= pow(16.*pow(M_PI,3),np);
+ weight_ /= std::tgamma(np+1); //remove rapidity ordering
std::vector<fastjet::PseudoJet> 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;
return {};
}
std::sort(begin(partons), end(partons), rapidity_less{});
return partons;
}
Particle PhaseSpacePoint::gen_boson(
RHEJ::ParticleID bosonid, double mass, double width,
RHEJ::RNG & ran
){
// Usual phase space measure
weight_ /= 16.*pow(M_PI, 3);
// Generate a y Gaussian distributed around 0
// TODO: magic number only for Higgs
const double y = random_normal(1.6, ran);
const double r1 = ran.flat();
const double sH = mass*(
mass + width*tan(M_PI/2.*r1 + (r1-1.)*atan(mass/width))
);
- auto p = gen_last_momentum(sH, y);
+ auto p = gen_last_momentum(outgoing_, sH, y);
return Particle{bosonid, std::move(p)};
}
- fastjet::PseudoJet PhaseSpacePoint::gen_last_momentum(
- const double mass_square, const double y
- ) const {
- std::array<double,2> pt{0.,0.};
- for (auto const & parton:outgoing_) {
- pt[0]-=parton.px();
- pt[1]-=parton.py();
- }
-
- const double mperp = sqrt(pt[0]*pt[0]+pt[1]*pt[1]+mass_square);
- const double pz=mperp*sinh(y);
- const double E=mperp*cosh(y);
-
- return {pt[0], pt[1], pz, E};
- }
-
Particle const & PhaseSpacePoint::most_backward_FKL(
std::vector<Particle> const & partons
) const{
if(!RHEJ::is_parton(partons[0])) return partons[1];
return partons[0];
}
Particle const & PhaseSpacePoint::most_forward_FKL(
std::vector<Particle> const & partons
) const{
const size_t last_idx = partons.size() - 1;
if(!RHEJ::is_parton(partons[last_idx])) return partons[last_idx-1];
return partons[last_idx];
}
Particle & PhaseSpacePoint::most_backward_FKL(
std::vector<Particle> & partons
) const{
if(!RHEJ::is_parton(partons[0])) return partons[1];
return partons[0];
}
Particle & PhaseSpacePoint::most_forward_FKL(
std::vector<Particle> & partons
) const{
const size_t last_idx = partons.size() - 1;
if(!RHEJ::is_parton(partons[last_idx])) return partons[last_idx-1];
return partons[last_idx];
}
void PhaseSpacePoint::reconstruct_incoming(
std::array<RHEJ::ParticleID, 2> const & ids,
RHEJ::PDF & pdf, double E_beam,
double uf,
RHEJ::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].p.e()-incoming_[0].p.pz())/sqrts;
const double xb=(incoming_[1].p.e()+incoming_[1].p.pz())/sqrts;
// abort if phase space point is outside of collider energy reach
if (xa>1. || xb>1.){
weight_=0;
status_ = too_much_energy;
return;
}
// pick pdfs
/** TODO:
* ufa, ufb don't correspond to our final scale choice.
* The reversed HEJ scale generators currently expect a full event as input,
* so fixing this is not completely trivial
*/
if(ids[0] == pid::proton || ids[0] == pid::p_bar){
const double ufa=uf;
incoming_[0].type = generate_incoming_id(xa, ufa, pdf, ran);
}
else {
incoming_[0].type = ids[0];
}
if(ids[1] == pid::proton || ids[1] == pid::p_bar){
const double ufb=uf;
incoming_[1].type = generate_incoming_id(xb, ufb, pdf, ran);
}
else {
incoming_[1].type = ids[1];
}
assert(momentum_conserved(1e-7));
}
RHEJ::ParticleID PhaseSpacePoint::generate_incoming_id(
double x, double uf, RHEJ::PDF & pdf,
RHEJ::RNG & ran
){
const double pdfg=fabs(pdf.pdfpt(0,x,uf,pid::gluon));
const double pdfu=fabs(pdf.pdfpt(0,x,uf,pid::up));
const double pdfd=fabs(pdf.pdfpt(0,x,uf,pid::down));
const double pdfux=fabs(pdf.pdfpt(0,x,uf,pid::u_bar));
const double pdfdx=fabs(pdf.pdfpt(0,x,uf,pid::d_bar));
const double pdfc=fabs(pdf.pdfpt(0,x,uf,pid::charm));
const double pdfs=fabs(pdf.pdfpt(0,x,uf,pid::strange));
const double pdfsx=fabs(pdf.pdfpt(0,x,uf,pid::s_bar));
const double pdfb=fabs(pdf.pdfpt(0,x,uf,pid::b));
const double pdftot=pdfg+4.0/9.0*(pdfu + pdfd + pdfux + pdfdx +pdfs +pdfsx + 2.0*(pdfc +pdfb ));
const double r1=pdftot*ran.flat();
double sum;
if (r1<(sum=pdfg)) {
weight_*=pdftot/pdfg;
return pid::gluon;
}
if (r1<(sum+=(4./9.)*pdfu)) {
weight_*=pdftot/((4./9.)*pdfu);
return pid::up;
}
else if (r1<(sum+=(4./9.)*pdfd)) {
weight_*=pdftot/((4./9.)*pdfd);
return pid::down;
}
else if (r1<(sum+=(4./9.)*pdfux)) {
weight_*=pdftot/((4./9.)*pdfux);
return pid::u_bar;
}
else if (r1<(sum+=(4./9.)*pdfdx)) {
weight_*=pdftot/((4./9.)*pdfdx);
return pid::d_bar;
}
else if (r1<(sum+=(4./9.)*pdfc)) {
weight_*=pdftot/((4./9.)*pdfc);
return pid::c;
}
else if (r1<(sum+=(4./9.)*pdfc)){
weight_*=pdftot/((4./9.)*pdfc);
return pid::c_bar;
}
else if (r1<(sum+=(4./9.)*pdfs)) {
weight_*=pdftot/((4./9.)*pdfs);
return pid::s;
}
else if (r1<(sum+=(4./9.)*pdfsx)) {
weight_*=pdftot/((4./9.)*pdfsx);
return pid::s_bar;
}
else if (r1<(sum+=(4./9.)*pdfb)) {
weight_*=pdftot/((4./9.)*pdfb);
return pid::b;
}
else if (r1<=(sum+=(4./9.)*pdfb)) {
weight_*=pdftot/((4./9.)*pdfb);
return pid::b_bar;
}
std::cout << "Error in choosing incoming parton: "<<x<<" "<<uf<<" "<<sum<<" "<<pdftot<<" "<<r1;
std::cout << " "<<pdfg+4./9.*(pdfu+pdfux+pdfd+pdfdx+pdfs+pdfsx+2.*(pdfc+pdfb))<<std::endl;
throw std::logic_error{"Failed to choose parton flavour"};
}
double PhaseSpacePoint::random_normal(
double stddev,
RHEJ::RNG & ran
){
const double r1 = ran.flat();
const double r2 = ran.flat();
const double lninvr1 = -log(r1);
const double result = stddev*sqrt(2.*lninvr1)*cos(2.*M_PI*r2);
weight_ *= exp(result*result/(2*stddev*stddev))*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);
}
Decay PhaseSpacePoint::select_decay_channel(
std::vector<Decay> const & decays,
RHEJ::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;
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<Particle> PhaseSpacePoint::decay_boson(
RHEJ::Particle const & parent,
std::vector<Decay> const & decays,
RHEJ::RNG & ran
){
const auto channel = select_decay_channel(decays, ran);
if(channel.products.size() != 2){
throw RHEJ::not_implemented{
"only decays into two particles are implemented"
};
}
std::vector<Particle> decay_products(channel.products.size());
for(size_t i = 0; i < channel.products.size(); ++i){
decay_products[i].type = channel.products[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.;
const double sin_phi = sqrt(1. - cos_phi*cos_phi); // Know 0 < phi < pi
const double px = E*cos(theta)*sin_phi;
const double py = E*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/main.cc
index b4640d4..e90f0e5 100644
--- a/FixedOrderGen/src/main.cc
+++ b/FixedOrderGen/src/main.cc
@@ -1,245 +1,247 @@
/**
* Name: main.cc
* Authors: Jeppe R. Andersen
*/
#include <fstream>
#include <algorithm>
#include <memory>
#include <chrono>
#include <iostream>
#include <map>
#include "yaml-cpp/yaml.h"
#include "config.hh"
#include "LHEF/LHEF.h"
#include "RHEJ/CombinedEventWriter.hh"
#include "RHEJ/get_analysis.hh"
#include "RHEJ/utility.hh"
//#include "RHEJ/EventReweighter.hh"
#include "RHEJ/stream.hh"
#include "RHEJ/LesHouchesWriter.hh"
#include "RHEJ/ProgressBar.hh"
#include "RHEJ/make_RNG.hh"
#include "EventGenerator.hh"
#include "PhaseSpacePoint.hh"
#include "Unweighter.hh"
#include "CrossSectionAccumulator.hh"
#include "Version.hh"
namespace{
constexpr auto banner =
" __ ___ __ ______ __ __ \n"
" / / / (_)___ _/ /_ / ____/___ ___ _________ ___ __ / /__ / /______ \n"
" / /_/ / / __ `/ __ \\ / __/ / __ \\/ _ \\/ ___/ __ `/ / / / __ / / _ \\/ __/ ___/ \n"
" / __ / / /_/ / / / / / /___/ / / / __/ / / /_/ / /_/ / / /_/ / __/ /_(__ ) \n"
" /_/ /_/_/\\__, /_/ /_/ /_____/_/ /_/\\___/_/ \\__, /\\__, / \\____/\\___/\\__/____/ \n"
" ____///__/ __ ____ ///__//____/ ______ __ \n"
" / ____(_) _____ ____/ / / __ \\_________/ /__ _____ / ____/__ ____ ___ _________ _/ /_____ _____\n"
" / /_ / / |/_/ _ \\/ __ / / / / / ___/ __ / _ \\/ ___/ / / __/ _ \\/ __ \\/ _ \\/ ___/ __ `/ __/ __ \\/ ___/\n"
" / __/ / /> </ __/ /_/ / / /_/ / / / /_/ / __/ / / /_/ / __/ / / / __/ / / /_/ / /_/ /_/ / / \n"
" /_/ /_/_/|_|\\___/\\__,_/ \\____/_/ \\__,_/\\___/_/ \\____/\\___/_/ /_/\\___/_/ \\__,_/\\__/\\____/_/ \n"
;
constexpr double invGeV2_to_pb = 389379292.;
constexpr long long max_warmup_events = 10000;
}
HEJFOG::Config load_config(char const * filename){
try{
return HEJFOG::load_config(filename);
}
catch(std::exception const & exc){
std::cerr << "Error: " << exc.what() << '\n';
std::exit(EXIT_FAILURE);
}
}
std::unique_ptr<RHEJ::Analysis> get_analysis(
YAML::Node const & parameters
){
try{
return RHEJ::get_analysis(parameters);
}
catch(std::exception const & exc){
std::cerr << "Failed to load analysis: " << exc.what() << '\n';
std::exit(EXIT_FAILURE);
}
}
void rescale_weights(RHEJ::Event & ev, double factor) {
ev.central().weight *= factor;
for(auto & var: ev.variations()){
var.weight *= factor;
}
}
int main(int argn, char** argv) {
using namespace std::string_literals;
if (argn < 2) {
std::cerr << "\n# Usage:\n.FOgen 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::unique_ptr<RHEJ::Analysis> analysis = get_analysis(
config.analysis_parameters
);
assert(analysis != nullptr);
auto ran = RHEJ::make_RNG(config.rng.name, config.rng.seed);
assert(ran != nullptr);
RHEJ::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.unordered_fraction,
config.Higgs_properties,
config.Higgs_coupling,
*ran
};
- //TODO: fix Les Houches init block (HEPRUP)
- LHEF::HEPRUP lhefheprup;
- lhefheprup.IDBMUP=std::pair<long,long>(config.beam.particles[0], config.beam.particles[1]);
- lhefheprup.EBMUP=std::make_pair(config.beam.energy, config.beam.energy);
- lhefheprup.PDFGUP=std::make_pair(0,0);
- lhefheprup.PDFSUP=std::make_pair(config.pdf_id,config.pdf_id);
- lhefheprup.NPRUP=1;
- lhefheprup.XSECUP=std::vector<double>(1.);
- lhefheprup.XERRUP=std::vector<double>(1.);
- lhefheprup.LPRUP=std::vector<int>{1};
-
- RHEJ::CombinedEventWriter writer{config.output, lhefheprup};
+ LHEF::HEPRUP heprup;
+ heprup.IDBMUP=std::pair<long,long>(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<double>(1.);
+ heprup.XERRUP=std::vector<double>(1.);
+ heprup.LPRUP=std::vector<int>{1};
+ heprup.generators.emplace_back(LHEF::XMLTag{});
+ heprup.generators.back().name = HEJFOG::Version::package_name();
+ heprup.generators.back().version = HEJFOG::Version::String();
+
+ RHEJ::CombinedEventWriter writer{config.output, heprup};
RHEJ::optional<HEJFOG::Unweighter> unweighter{};
std::map<HEJFOG::Status, int> status_counter;
std::vector<RHEJ::Event> events;
int trials = 0;
// warm-up phase to train unweighter
if(config.unweight) {
std::cout << "Calibrating unweighting ...\n";
const auto warmup_start = clock::now();
const size_t warmup_events = config.unweight->sample_size;
RHEJ::ProgressBar<size_t> warmup_progress{std::cout, warmup_events};
for(; events.size() < warmup_events; ++trials){
auto ev = generator.gen_event();
++status_counter[generator.status()];
if(generator.status() == HEJFOG::good && analysis->pass_cuts(ev, ev)) {
events.emplace_back(std::move(ev));
++warmup_progress;
}
}
std::cout << std::endl;
unweighter = HEJFOG::Unweighter{
begin(events), end(events), config.unweight->max_dev, *ran,
config.jets.peak_pt?(*config.jets.peak_pt):0.
};
std::vector<RHEJ::Event> unweighted_events;
for(auto && ev: events) {
auto unweighted = unweighter->unweight(std::move(ev));
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<double>(events.size())/config.events;
const std::chrono::duration<double> remaining_time =
(warmup_end- warmup_start)*(1./completion - 1);
const auto finish = clock::to_time_t(
std::chrono::time_point_cast<std::chrono::seconds>(warmup_end + remaining_time)
);
std::cout
<< "Generated " << events.size() << "/" << config.events << " events ("
<< static_cast<int>(std::round(100*completion)) << "%)\n"
<< "Estimated remaining generation time: "
<< remaining_time.count() << " seconds ("
<< std::put_time(std::localtime(&finish), "%c") << ")\n\n";
}
RHEJ::ProgressBar<long long> progress{std::cout, config.events};
progress.increment(events.size());
events.reserve(config.events);
for(; events.size() < static_cast<size_t>(config.events); ++trials){
auto ev = generator.gen_event();
++status_counter[generator.status()];
if(generator.status() == HEJFOG::good && analysis->pass_cuts(ev, ev)) {
if(unweighter) {
auto unweighted = unweighter->unweight(std::move(ev));
if(! unweighted) continue;
ev = std::move(*unweighted);
}
events.emplace_back(std::move(ev));
++progress;
}
}
std::cout << std::endl;
HEJFOG::CrossSectionAccumulator xs;
for(auto & ev: events){
rescale_weights(ev, invGeV2_to_pb/trials);
analysis->fill(ev, ev);
writer.write(ev);
xs.fill(ev);
}
analysis->finalise();
const std::chrono::duration<double> run_time = (clock::now() - start_time);
std::cout << "\nTask Runtime: " << run_time.count() << " seconds.\n\n";
std::cout.precision(10);
std::cout.setf(std::ios::fixed);
std::cout
<< " " << std::left << std::setw(20)
<< "Cross section: " << xs.total().value
<< " +- " << std::sqrt(xs.total().error) << " (pb)\n";
for(auto const & xs_type: xs) {
std::cout
<< " " << std::left << std::setw(20)
<< (RHEJ::event_type::names[xs_type.first] + ": "s)
<< xs_type.second.value << " +- "
<< std::sqrt(xs_type.second.error) << " (pb)\n";
}
std::cout << '\n';
for(auto && entry: status_counter){
const double fraction = static_cast<double>(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 << "%\n";
}
}
diff --git a/FixedOrderGen/t/2j.cc b/FixedOrderGen/t/2j.cc
new file mode 100644
index 0000000..1e078c3
--- /dev/null
+++ b/FixedOrderGen/t/2j.cc
@@ -0,0 +1,57 @@
+#ifdef NDEBUG
+#undef NDEBUG
+#endif
+
+#include <algorithm>
+#include <cmath>
+#include <cassert>
+#include <iostream>
+
+#include "config.hh"
+#include "EventGenerator.hh"
+#include "RHEJ/Mixmax.hh"
+
+#include "RHEJ/PDF.hh"
+#include "RHEJ/MatrixElement.hh"
+
+using namespace HEJFOG;
+
+int main(){
+ constexpr double invGeV2_to_pb = 389379292.;
+ constexpr double xs_ref = 86.42031848*1e6; //calculated with "combined" HEJ svn r3480
+
+ auto config = load_config("config_2j.yml");
+
+ RHEJ::Mixmax ran{};
+ HEJFOG::EventGenerator generator{
+ config.process,
+ config.beam,
+ RHEJ::ScaleGenerator{
+ config.scales.base,
+ config.scales.factors,
+ config.scales.max_ratio
+ },
+ config.jets,
+ config.pdf_id,
+ config.unordered_fraction,
+ config.Higgs_properties,
+ config.Higgs_coupling,
+ ran
+ };
+
+ double xs = 0., xs_err = 0.;
+ for (int trials = 0; trials < config.events; ++trials){
+ auto ev = generator.gen_event();
+ if(generator.status() != good) continue;
+ ev.central().weight *= invGeV2_to_pb;
+ ev.central().weight /= config.events;
+
+ xs += ev.central().weight;
+ xs_err += ev.central().weight*ev.central().weight;
+ }
+ xs_err = std::sqrt(xs_err);
+ std::cout << xs_ref << " ~ " << xs << " +- " << xs_err << '\n';
+
+ assert(std::abs(xs - xs_ref) < 3*xs_err);
+ assert(xs_err < 0.01*xs);
+}
diff --git a/FixedOrderGen/t/4j.cc b/FixedOrderGen/t/4j.cc
new file mode 100644
index 0000000..5ff903c
--- /dev/null
+++ b/FixedOrderGen/t/4j.cc
@@ -0,0 +1,58 @@
+#ifdef NDEBUG
+#undef NDEBUG
+#endif
+
+#include <algorithm>
+#include <cmath>
+#include <cassert>
+#include <iostream>
+
+#include "config.hh"
+#include "EventGenerator.hh"
+#include "RHEJ/Mixmax.hh"
+
+#include "RHEJ/PDF.hh"
+#include "RHEJ/MatrixElement.hh"
+
+using namespace HEJFOG;
+
+int main(){
+ constexpr double invGeV2_to_pb = 389379292.;
+ constexpr double xs_ref = 0.81063619*1e6; //calculated with "combined" HEJ svn r3480
+
+ auto config = load_config("config_2j.yml");
+ config.process.njets = 4;
+
+ RHEJ::Mixmax ran{};
+ HEJFOG::EventGenerator generator{
+ config.process,
+ config.beam,
+ RHEJ::ScaleGenerator{
+ config.scales.base,
+ config.scales.factors,
+ config.scales.max_ratio
+ },
+ config.jets,
+ config.pdf_id,
+ config.unordered_fraction,
+ config.Higgs_properties,
+ config.Higgs_coupling,
+ ran
+ };
+
+ double xs = 0., xs_err = 0.;
+ for (int trials = 0; trials < config.events; ++trials){
+ auto ev = generator.gen_event();
+ if(generator.status() != good) continue;
+ ev.central().weight *= invGeV2_to_pb;
+ ev.central().weight /= config.events;
+
+ xs += ev.central().weight;
+ xs_err += ev.central().weight*ev.central().weight;
+ }
+ xs_err = std::sqrt(xs_err);
+ std::cout << xs_ref << " ~ " << xs << " +- " << xs_err << '\n';
+
+ assert(std::abs(xs - xs_ref) < 3*xs_err);
+ assert(xs_err < 0.03*xs);
+}
diff --git a/FixedOrderGen/t/config_2j.yml b/FixedOrderGen/t/config_2j.yml
new file mode 100644
index 0000000..67625dc
--- /dev/null
+++ b/FixedOrderGen/t/config_2j.yml
@@ -0,0 +1,26 @@
+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
+
+unordered fraction: 0
+
+scales: 125
+
+random generator:
+ name: mixmax
+
+Higgs properties:
+ mass: 125
+ width: 0
diff --git a/doc/developer_manual/developer_manual.tex b/doc/developer_manual/developer_manual.tex
index 8efba55..34a8c93 100644
--- a/doc/developer_manual/developer_manual.tex
+++ b/doc/developer_manual/developer_manual.tex
@@ -1,1113 +1,1111 @@
\documentclass[a4paper,11pt]{article}
\usepackage{fourier}
\usepackage[T1]{fontenc}
\usepackage{microtype}
\usepackage{geometry}
\usepackage{enumitem}
\setlist[description]{leftmargin=\parindent,labelindent=\parindent}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage[utf8x]{inputenc}
\usepackage{graphicx}
\usepackage{xcolor}
\usepackage{todonotes}
\usepackage{listings}
\usepackage{xspace}
\usepackage{tikz}
\usetikzlibrary{arrows.meta}
\usetikzlibrary{shapes}
\usetikzlibrary{calc}
\usepackage[colorlinks,linkcolor={blue!50!black}]{hyperref}
\graphicspath{{build/figures/}{figures/}}
\emergencystretch \hsize
\newcommand{\HEJ}{{\tt HEJ}\xspace}
\newcommand{\HIGHEJ}{\emph{High Energy Jets}\xspace}
\newcommand{\cmake}{\href{https://cmake.org/}{cmake}\xspace}
\newcommand{\html}{\href{https://www.w3.org/html/}{html}\xspace}
\newcommand{\YAML}{\href{http://yaml.org/}{YAML}\xspace}
\definecolor{darkgreen}{rgb}{0,0.4,0}
\lstset{ %
backgroundcolor=\color{lightgray}, % choose the background color; you must add \usepackage{color} or \usepackage{xcolor}
basicstyle=\footnotesize\usefont{T1}{DejaVuSansMono-TLF}{m}{n}, % the size of the fonts that are used for the code
breakatwhitespace=false, % sets if automatic breaks should only happen at whitespace
breaklines=false, % sets automatic line breaking
captionpos=t, % sets the caption-position to bottom
commentstyle=\color{red}, % comment style
deletekeywords={...}, % if you want to delete keywords from the given language
escapeinside={\%*}{*)}, % if you want to add LaTeX within your code
extendedchars=true, % lets you use non-ASCII characters; for 8-bits encodings only, does not work with UTF-8
frame=false, % adds a frame around the code
keepspaces=true, % keeps spaces in text, useful for keeping indentation of code (possibly needs columns=flexible)
keywordstyle=\color{blue}, % keyword style
otherkeywords={}, % if you want to add more keywords to the set
numbers=none, % where to put the line-numbers; possible values are (none, left, right)
numbersep=5pt, % how far the line-numbers are from the code
rulecolor=\color{black}, % if not set, the frame-color may be changed on line-breaks within not-black text (e.g. comments (green here))
showspaces=false, % show spaces everywhere adding particular underscores; it overrides 'showstringspaces'
showstringspaces=false, % underline spaces within strings only
showtabs=false, % show tabs within strings adding particular underscores
stepnumber=2, % the step between two line-numbers. If it's 1, each line will be numbered
stringstyle=\color{gray}, % string literal style
tabsize=2, % sets default tabsize to 2 spaces
title=\lstname,
emph={},
emphstyle=\color{darkgreen}
}
\begin{document}
\tikzstyle{mynode}=[rectangle split,rectangle split parts=2, draw,rectangle split part fill={lightgray, none}]
\title{Reversed HEJ developer manual}
\author{}
\maketitle
\tableofcontents
\newpage
\section{Overview}
\label{sec:overview}
Reversed HEJ is a C++ program and library implementing an algorithm to
apply \HIGHEJ resummation~\cite{Andersen:2008ue,Andersen:2008gc} to
pre-generated fixed-order events. This document is intended to give an
overview over the concepts and structure of this implementation.
\subsection{Project structure}
\label{sec:project}
Reversed HEJ is developed under the \href{https://git-scm.com/}{git}
version control system. The main repository is on the IPPP
\href{https://gitlab.com/}{gitlab} server under
\url{https://gitlab.dur.scotgrid.ac.uk/hej/reversed_hej}. To get a local
copy, get an account on the gitlab server and use
\begin{lstlisting}[language=sh,caption={}]
git clone git@gitlab.dur.scotgrid.ac.uk:hej/reversed_hej.git
\end{lstlisting}
This should create a directory \texttt{reversed\_hej} with the following
contents:
\begin{description}
\item[doc:] Contains additional documentation, see section~\ref{sec:doc}.
\item[include:] Contains the C++ header files.
\item[src:] Contains the C++ source files.
\item[t:] Contains the source code for the automated tests.
\item[CMakeLists.txt:] Configuration file for the \cmake build
system. See section~\ref{sec:cmake}.
\item[cmake:] Auxiliary files for \cmake. This includes modules for
finding installed software in \texttt{cmake/Modules} and templates for
code generation during the build process in \texttt{cmake/Templates}.
\item[config.yml:] Sample configuration file for running reversed HEJ.
\item[events.lhe:] Fixed-order event sample that can be used as input
for reversed HEJ.
\item[FixedOrderGen:] Contains the code for the fixed-order generator,
see section~\ref{sec:HEJFOG}.
\end{description}
In the following all paths are given relative to the
\texttt{reversed\_hej} directory.
\subsection{Documentation}
\label{sec:doc}
The \texttt{doc} directory contains user documentation in
\texttt{doc/sphinx} and the configuration to generate source code
documentation in \texttt{doc/doxygen}.
The user documentation explains how to install and run reversed HEJ. The
format is
\href{http://docutils.sourceforge.net/rst.html}{reStructuredText}, which
is mostly human-readable. Other formats, like \html, can be generated with the
\href{http://www.sphinx-doc.org/en/master/}{sphinx} generator with
\begin{lstlisting}[language=sh,caption={}]
make html
\end{lstlisting}
To document the source code we use
\href{https://www.stack.nl/~dimitri/doxygen/}{doxygen}. To generate
\html documentation, use the command
\begin{lstlisting}[language=sh,caption={}]
doxygen Doxyfile
\end{lstlisting}
in the \texttt{doc/doxygen} directory.
\subsection{Build system}
\label{sec:cmake}
For the most part, reversed HEJ is a library providing classes and
functions that can be used to add resummation to fixed-order events. In
addition, there is a relatively small executable program leveraging this
library to read in events from an input file and produce resummation
events. Both the library and the program are built and installed with
the help of \cmake.
Debug information can be turned on by using
\begin{lstlisting}[language=sh,caption={}]
cmake base/directory -DCMAKE_BUILD_TYPE=Debug
make install
\end{lstlisting}
This facilitates the use of debuggers like \href{https://www.gnu.org/software/gdb/}{gdb}.
The main \cmake configuration file is \texttt{CMakeLists.txt}. It defines the
compiler flags, software prerequisites, header and source files used to
build reversed HEJ, and the automated tests.
\texttt{cmake/Modules} contains module files that help with the
detection of the software prerequisites and \texttt{cmake/Templates}
template files for the automatic generation of header and
source files. For example, this allows to only keep the version
information in one central location (\texttt{CMakeLists.txt}) and
automatically generate a header file from the template \texttt{Version.hh.in} to propagate this to the C++ code.
\subsection{General coding guidelines}
\label{sec:notes}
The goal is to make the reversed HEJ code well-structured and
readable. Here are a number of guidelines to this end.
\begin{description}
\item[Observe the boy scout rule.] Always leave the code cleaner
than how you found it. Ugly hacks can be useful for testing, but
shouldn't make their way into the main branch.
\item[Ask if something is unclear.] Often there is a good reason why
code is written the way it is. Sometimes that reason is only obvious to
the original author (use \lstinline!git blame! to find them), in which
case they should be poked to add a comment. Sometimes there is no good
reason, but nobody has had the time to come up with something better,
yet. In some places the code might just be bad.
\item[Don't break tests.] There are a number of tests in the \texttt{t}
directory, which can be run with \lstinline!make test!. Ideally, all
tests should run successfully in each git revision. If your latest
commit broke a test and you haven't pushed to the central repository
yet, you can fix it with \lstinline!git commit --amend!. If an earlier
local commit broke a test, you can use \lstinline!git rebase -i! if
you feel confident.
\item[Test your new code.] When you add some new functionality, also add an
automated test. This can be useful even if you don't know the
``correct'' result because it prevents the code from changing its behaviour
silently in the future. \href{http://www.valgrind.org/}{valgrind} is a
very useful tool to detect potential memory leaks.
\item[Stick to the coding style.] It is somewhat easier to read code
that has a uniform coding and indentation style. We don't have a
strict style, but it helps if your code looks similar to what is
already there.
\end{description}
\section{Program flow}
\label{sec:flow}
A run of the reversed HEJ program has three stages: initialisation,
event processing, and cleanup. The following sections outline these
stages and their relations to the various classes and functions in the
code. Unless denoted otherwise, all classes and functions are part of
the \lstinline!RHEJ! namespace. The code for the reversed HEJ program is
in \texttt{src/main.cc}, all other code comprises the reversed HEJ
library. Classes and free functions are usually implemented in header
and source files with a corresponding name, i.e. the code for
\lstinline!MyClass! can usually be found in
\texttt{include/RHEJ/MyClass.hh} and \texttt{src/MyClass.cc}.
\subsection{Initialisation}
\label{sec:init}
The first step is to load and parse the \YAML configuration file. The
entry point for this is the \lstinline!load_config! function and the
related code can be found in \texttt{include/RHEJ/YAMLreader.hh},
\texttt{include/RHEJ/config.hh} and the corresponding \texttt{.cc} files
in the \texttt{src} directory. The implementation is based on the
\href{https://github.com/jbeder/yaml-cpp}{yaml-cpp} library.
The \lstinline!load_config! function returns a \lstinline!Config! object
containing all settings. To detect potential mistakes as early as
possible, we throw an exception whenever one of the following errors
occurs:
\begin{itemize}
\item There is an unknown option in the \YAML file.
\item A setting is invalid, for example a string is given where a number
would be expected.
\item An option value is not set.
\end{itemize}
The third rule is sometimes relaxed for ``advanced'' settings with an
obvious default, like for importing custom scales or analyses.
The information stored in the \lstinline!Config! object is then used to
initialise various objects required for the event processing stage
described in section~\ref{sec:processing}. First, the
\lstinline!get_analysis! function creates an object that inherits from
the \lstinline!Analysis! interface.\footnote{In the context of C++ the
proper technical expression is ``pure abstract class''.} Using an
interface allows us to decide the concrete type of the analysis at run
time instead of having to make a compile-time decision. Depending on the
settings, \lstinline!get_analysis! creates either a user-defined
analysis loaded from an external library (see the user documentation
\todo{link}) or the default \lstinline!EmptyAnalysis!, which does
nothing.
Together with a number of further objects, whose roles are described in
section~\ref{sec:processing}, we also initialise the global random
number generator. We again use an interface to defer deciding the
concrete type until the program is actually run. Currently, we support the
\href{https://mixmax.hepforge.org/}{MIXMAX}
(\texttt{include/RHEJ/Mixmax.hh}) and Ranlux64
(\texttt{include/RHEJ/Ranlux64.hh}) random number generators, both are provided
by \href{http://proj-clhep.web.cern.ch/}{CLHEP}.
We also set up a \lstinline!LHEF::Reader! object (see
\href{http://home.thep.lu.se/~leif/LHEF/}{\texttt{include/LHEF/LHEF.h}}) for
reading events from a file in the Les
Houches event file format~\cite{Alwall:2006yp}. A small wrapper around
the
\href{https://www.boost.org/doc/libs/1_67_0/libs/iostreams/doc/index.html}{boost
iostreams} library allows us to also read event files compressed with
\href{https://www.gnu.org/software/gzip/}{gzip}. The wrapper code is in
\texttt{include/RHEJ/stream.hh} and the \texttt{src/stream.cc}.
\subsection{Event processing}
\label{sec:processing}
In the second stage events are continously read from the event
file. After jet clustering, a number of corresponding resummation events
are generated for each input event and fed into the analysis and a
number of output files. The roles of various classes and functions are
illustrated in the following flow chart:
\begin{center}
\begin{tikzpicture}[node distance=2cm and 5mm]
\node (reader) [mynode]
{\lstinline!LHEF::Reader::readEvent!\nodepart{second}{read event}};
\node
(cluster) [mynode,below=of reader]
{\lstinline!Event! constructor\nodepart{second}{cluster jets}};
\node
(resum) [mynode,below=of cluster]
{\lstinline!EventReweighter::reweight!\nodepart{second}{perform resummation}};
\node
(cut) [mynode,below=of resum]
{\lstinline!Analysis::pass_cuts!\nodepart{second}{apply cuts}};
\node
(fill) [mynode,below left=of cut]
{\lstinline!Analysis::fill!\nodepart{second}{analyse event}};
\node
(write) [mynode,below right=of cut]
{\lstinline!CombinedEventWriter::write!\nodepart{second}{write out event}};
\node
(control) [below=of cut] {};
\draw[-{Latex[length=3mm, width=1.5mm]}]
(reader.south) -- node[left] {\lstinline!LHEF::HEPEUP!} (cluster.north);
\draw[-{Latex[length=3mm, width=1.5mm]}]
(cluster.south) -- node[left] {\lstinline!Event!} (resum.north);
\draw[-{Latex[length=3mm, width=1.5mm]}]
(resum.south) -- (cut.north);
\draw[-{Latex[length=3mm, width=1.5mm]}]
($(resum.south)+(7mm, 0cm)$) -- ($(cut.north)+(7mm, 0cm)$);
\draw[-{Latex[length=3mm, width=1.5mm]}]
($(resum.south)-(7mm, 0cm)$) -- node[left] {\lstinline!Event!} ($(cut.north)-(7mm, 0cm)$);
\draw[-{Latex[length=3mm, width=1.5mm]}]
($(cut.south)-(3mm,0mm)$) .. controls ($(control)-(3mm,0mm)$) ..node[left] {\lstinline!Event!} (fill.east);
\draw[-{Latex[length=3mm, width=1.5mm]}]
($(cut.south)-(3mm,0mm)$) .. controls ($(control)-(3mm,0mm)$) .. (write.west);
\draw[-{Latex[length=3mm, width=1.5mm]}]
($(cut.south)+(3mm,0mm)$) .. controls ($(control)+(3mm,0mm)$) .. (fill.east);
\draw[-{Latex[length=3mm, width=1.5mm]}]
($(cut.south)+(3mm,0mm)$) .. controls ($(control)+(3mm,0mm)$) ..node[right] {\lstinline!Event!} (write.west);
\end{tikzpicture}
\end{center}
The resummation is performed by the \lstinline!EventReweighter! class,
which is described in more detail in section~\ref{sec:resum}. The
\lstinline!CombinedEventWriter! writes events to zero or more output
files. To this end, it contains a number of objects implementing the
\lstinline!EventWriter! interface. These event writers typically write
the events to a file in a given format. We currently have the
\lstinline!LesHouchesWriter! for event files in the Les Houches Event
File format and the \lstinline!HepMCWriter! for the
\href{https://hepmc.web.cern.ch/hepmc/}{HepMC} format (Version 2 and 3).
\subsection{Resummation}
\label{sec:resum}
In the \lstinline!EventReweighter::reweight! member function, we first
classify the input fixed-order event (FKL, unordered, non-HEJ, \dots)
and decide according to the user settings whether to discard, keep, or
resum the event. If we perform resummation for the given event, we
generate a number of trial \lstinline!PhaseSpacePoint! objects. Phase
space generation is discussed in more detail in
section~\ref{sec:pspgen}. We then perform jet clustering according to
the settings for the resummation jets on each
\lstinline!PhaseSpacePoint!, update the factorisation and
renormalisation scale in the resulting \lstinline!Event! and reweight it
according to the ratio of pdf factors and \HEJ matrix elements between
resummation and original fixed-order event:
\begin{center}
\begin{tikzpicture}[node distance=2cm and 5mm]
\node (in) {};
\node (treat) [diamond,draw,below=of in,minimum size=3.5cm,
label={[anchor=west, inner sep=8pt]west:discard},
label={[anchor=east, inner sep=14pt]east:keep},
label={[anchor=south, inner sep=20pt]south:reweight}
] {};
\draw (treat.north west) -- (treat.south east);
\draw (treat.north east) -- (treat.south west);
\node
(psp) [mynode,below=of treat]
{\lstinline!PhaseSpacePoint! constructor};
\node
(cluster) [mynode,below=of psp]
{\lstinline!Event! constructor\nodepart{second}{cluster jets}};
\node
(gen_scales) [mynode,below=of cluster]
{\lstinline!ScaleGenerator::operator()!\nodepart{second}{update scales}};
\node
(rescale) [mynode,below=of gen_scales]
{\lstinline!PDF::pdfpt!,
\lstinline!MatrixElement!\nodepart{second}{reweight}};
\node (out) [below of=rescale] {};
\draw[-{Latex[length=3mm, width=1.5mm]}]
(in.south) -- node[left] {\lstinline!Event!} (treat.north);
\draw[-{Latex[length=3mm, width=1.5mm]}]
(treat.south) -- node[left] {\lstinline!Event!} (psp.north);
\draw[-{Latex[length=3mm, width=1.5mm]}]
(psp.south) -- (cluster.north);
\draw[-{Latex[length=3mm, width=1.5mm]}]
($(psp.south)+(7mm, 0cm)$) -- ($(cluster.north)+(7mm, 0cm)$);
\draw[-{Latex[length=3mm, width=1.5mm]}]
($(psp.south)-(7mm, 0cm)$) -- node[left]
{\lstinline!PhaseSpacePoint!} ($(cluster.north)-(7mm, 0cm)$);
\draw[-{Latex[length=3mm, width=1.5mm]}]
(cluster.south) -- (gen_scales.north);
\draw[-{Latex[length=3mm, width=1.5mm]}]
($(cluster.south)+(7mm, 0cm)$) -- ($(gen_scales.north)+(7mm, 0cm)$);
\draw[-{Latex[length=3mm, width=1.5mm]}]
($(cluster.south)-(7mm, 0cm)$) -- node[left]
{\lstinline!Event!} ($(gen_scales.north)-(7mm, 0cm)$);
\draw[-{Latex[length=3mm, width=1.5mm]}]
(gen_scales.south) -- (rescale.north);
\draw[-{Latex[length=3mm, width=1.5mm]}]
($(gen_scales.south)+(7mm, 0cm)$) -- ($(rescale.north)+(7mm, 0cm)$);
\draw[-{Latex[length=3mm, width=1.5mm]}]
($(gen_scales.south)-(7mm, 0cm)$) -- node[left]
{\lstinline!Event!} ($(rescale.north)-(7mm, 0cm)$);
\draw[-{Latex[length=3mm, width=1.5mm]}]
(rescale.south) -- (out.north);
\draw[-{Latex[length=3mm, width=1.5mm]}]
($(rescale.south)+(7mm, 0cm)$) -- ($(out.north)+(7mm, 0cm)$);
\draw[-{Latex[length=3mm, width=1.5mm]}]
($(rescale.south)-(7mm, 0cm)$) -- node[left]
{\lstinline!Event!} ($(out.north)-(7mm, 0cm)$);
\node (helper) at ($(treat.east) + (15mm,0cm)$) {};
\draw[-{Latex[length=3mm, width=1.5mm]}]
(treat.east) -- ($(treat.east) + (15mm,0cm)$)
-- node[left] {\lstinline!Event!} (helper |- gen_scales.east) -- (gen_scales.east)
;
\end{tikzpicture}
\end{center}
\subsection{Phase space point generation}
\label{sec:pspgen}
The resummed and matched \HEJ cross section for pure jet production of
FKL configurations is given by (cf. eq. (3) of~\cite{Andersen:2018tnm})
\begin{align}
\label{eq:resumdijetFKLmatched2}
% \begin{split}
\sigma&_{2j}^\mathrm{resum, match}=\sum_{f_1, f_2}\ \sum_m
\prod_{j=1}^m\left(
\int_{p_{j\perp}^B=0}^{p_{j\perp}^B=\infty}
\frac{\mathrm{d}^2\mathbf{p}_{j\perp}^B}{(2\pi)^3}\ \int
\frac{\mathrm{d} y_j^B}{2} \right) \
(2\pi)^4\ \delta^{(2)}\!\!\left(\sum_{k=1}^{m}
\mathbf{p}_{k\perp}^B\right)\nonumber\\
&\times\ x_a^B\ f_{a, f_1}(x_a^B, Q_a^B)\ x_b^B\ f_{b, f_2}(x_b^B, Q_b^B)\
\frac{\overline{\left|\mathcal{M}_\text{LO}^{f_1f_2\to f_1g\cdots
gf_2}\big(\big\{p^B_j\big\}\big)\right|}^2}{(\hat {s}^B)^2}\nonumber\\
& \times (2\pi)^{-4+3m}\ 2^m \nonumber\\
&\times\ \sum_{n=2}^\infty\
\int_{p_{1\perp}=p_{\perp,\mathrm{min}} }^{p_{1\perp}=\infty}
\frac{\mathrm{d}^2\mathbf{p}_{1\perp}}{(2\pi)^3}\
\int_{p_{n\perp}=p_{\perp,\mathrm{min}}}^{p_{n\perp}=\infty}
\frac{\mathrm{d}^2\mathbf{p}_{n\perp}}{(2\pi)^3}\
\prod_{i=2}^{n-1}\int_{p_{i\perp}=\lambda}^{p_{i\perp}=\infty}
\frac{\mathrm{d}^2\mathbf{p}_{i\perp}}{(2\pi)^3}\ (2\pi)^4\ \delta^{(2)}\!\!\left(\sum_{k=1}^n
\mathbf{p}_{k\perp}\right )\\
&\times \ \mathbf{T}_y \prod_{i=1}^n
\left(\int \frac{\mathrm{d} y_i}{2}\right)\
\mathcal{O}_{mj}^e\
\left(\prod_{l=1}^{m-1}\delta^{(2)}(\mathbf{p}_{\mathcal{J}_{l}\perp}^B -
\mathbf{j}_{l\perp})\right)\
\left(\prod_{l=1}^m\delta(y^B_{\mathcal{J}_l}-y_{\mathcal{J}_l})\right)
\ \mathcal{O}_{2j}(\{p_i\})\nonumber\\
&\times \frac{(\hat{s}^B)^2}{\hat{s}^2}\ \frac{x_a f_{a,f_1}(x_a, Q_a)\ x_b f_{b,f_2}(x_b, Q_b)}{x_a^B\ f_{a,f_1}(x_a^B, Q_a^B)\ x_b^B\ f_{b,f_2}(x_b^B, Q_b^B)}\ \frac{\overline{\left|\mathcal{M}_{\mathrm{HEJ}}^{f_1 f_2\to f_1 g\cdots
gf_2}(\{ p_i\})\right|}^2}{\overline{\left|\mathcal{M}_\text{LO, HEJ}^{f_1f_2\to f_1g\cdots
gf_2}\big(\big\{p^B_j\big\}\big)\right|}^{2}} \,.\nonumber
% \end{split}
\end{align}
The first two lines correspond to the generation of the fixed-order
input events with incoming partons $f_1, f_2$ and outgoing momenta
$p_j^B$, where $\mathbf{p}_{j\perp}^B$ and $y_j^B$ denote the respective
transverse momentum and rapidity. Note that, at leading order, these
coincide with the fixed-order jet momenta $p_{\mathcal{J}_j}^B$.
$f_{a,f_1}(x_a, Q_a),f_{b,f_2}(x_b, Q_b)$ are the pdf factors for the incoming partons with
momentum fractions $x_a$ and $x_b$. The square of the partonic
centre-of-mass energy is denoted by $\hat{s}^B$ and
$\mathcal{M}_\text{LO}^{f_1f_2\to f_1g\cdots gf_2}$ is the
leading-order matrix element.
The third line is a factor accounting for the different multiplicities
between fixed-order and resummation events. Lines four and five are
the integration over the resummation phase space described in this
section. $p_i$ are the momenta of the outgoing partons in resummation
phase space. $\mathbf{T}_y$ denotes rapidity
ordering and $\mathcal{O}_{mj}^e$ projects out the exclusive $m$-jet
component. The relation between resummation and fixed-order momenta is
fixed by the $\delta$ functions. The first sets each transverse fixed-order jet
momentum to some function $\mathbf{j_{l\perp}}$ of the resummation
momenta. The exact form is described in section~\ref{sec:ptj_res}. The second
$\delta$ forces the rapidities of resummation and fixed-order jets to be
the same. Finally, the last line is the reweighting of pdf and matrix
element factors already shown in section~\ref{sec:resum}.
There are two kinds of cut-off in the integration over the resummation
partons. $\lambda$ is a technical cut-off connected to the cancellation
of infrared divergencies between real and virtual corrections. Its
numerical value is set in
\texttt{include/RHEJ/Constants.h}. $p_{\perp,\mathrm{min}}$ regulates
and \emph{uncancelled} divergence in the extremal parton momenta. Its
size is set by the user configuration \todo{link}.
It is straightforward to generalise eq.~(\ref{eq:resumdijetFKLmatched2})
to unordered configurations and processes with additional colourless
emissions, for example a Higgs or electroweak boson. In the latter case only
the fixed-order integration and the matrix elements change.
\subsubsection{Gluon Multiplicity}
\label{sec:psp_ng}
The first step in evaluating the resummation phase space in
eq.~(\ref{eq:resumdijetFKLmatched2}) is to randomly pick terms in the
sum over the number of emissions. This sampling of the gluon
multiplicity is done in the \lstinline!PhaseSpacePoint::sample_ng!
function in \texttt{src/PhaseSpacePoint.cc}.
The typical number of extra emissions depends strongly on the rapidity
span of the underlying fixed-order event. Let us, for example, consider
a fixed-order FKL-type multi-jet configuration with rapidities
$y_{j_f},\,y_{j_b}$ of the most forward and backward jets,
respectively. By eq.~(\ref{eq:resumdijetFKLmatched2}), the jet
multiplicity and the rapidity of each jet are conserved when adding
resummation. This implies that additional hard radiation is restricted
to rapidities $y$ within a region $y_{j_b} \lesssim y \lesssim
y_{j_f}$. Within \HEJ, we require the most forward and most backward
emissions to be hard \todo{specify how hard} in order to avoid divergences, so this constraint
in fact applies to \emph{all} additional radiation.
To simplify the remaining discussion, let us remove the FKL rapidity
ordering
\begin{equation}
\label{eq:remove_y_order}
\mathbf{T}_y \prod_{i=1}^n\int \frac{\mathrm{d}y_i}{2} =
\frac{1}{n!}\prod_{i=1}^n\int
\frac{\mathrm{d}y_i}{2}\,,
\end{equation}
where all rapidity integrals now cover a region which is approximately
bounded by $y_{j_b}$ and $y_{j_f}$. Each of the $m$ jets has to contain at least
one parton; selecting random emissions we can rewrite the phase space
integrals as
\begin{equation}
\label{eq:select_jets}
\frac{1}{n!}\prod_{i=1}^n\int [\mathrm{d}p_i] =
\left(\prod_{i=1}^{m}\int [\mathrm{d}p_i]\ {\cal J}_i(p_i)\right)
\frac{1}{n_g!}\prod_{i=m+1}^{m+n_g}\int [\mathrm{d}p_i]
\end{equation}
with jet selection functions
\begin{equation}
\label{eq:def_jet_selection}
{\cal J}_i(p) =
\begin{cases}
1 &p\text{ clustered into jet }i\\
0 & \text{otherwise}
\end{cases}
\end{equation}
and $n_g \equiv n - m$. Here and in the following we use the short-hand
notation $[\mathrm{d}p_i]$ to denote the phase-space measure for parton
$i$. As is evident from eq.~\eqref{eq:select_jets}, adding an extra emission
$n_g+1$ introduces a suppression factor $\tfrac{1}{n_g+1}$. However, the
additional phase space integral also results in an enhancement proportional
to $\Delta y_{j_f j_b} = y_{j_f} - y_{j_b}$. This is a result of the
rapidity-independence of the MRK limit of the integrand, consisting of the
matrix elements divided by the flux factor. Indeed, we observe that the
typical number of gluon emissions is to a good approximation proportional to
the rapidity separation and the phase space integral is dominated by events
with $n_g \approx \Delta y_{j_f j_b}$.
For the actual phase space sampling, we assume a Poisson distribution
and extract the mean number of gluon emissions in different rapidity
bins and fit the results to a linear function in $\Delta y_{j_f j_b}$,
finding a coefficient of $0.975$ for the inclusive production of a Higgs
boson with two jets. Here are the observed and fitted average gluon
multiplicities as a function of $\Delta y_{j_f j_b}$:
\begin{center}
\includegraphics[width=.75\textwidth]{ng_mean}
\end{center}
As shown for two rapidity slices the assumption of a Poisson
distribution is also a good approximation:
\begin{center}
\includegraphics[width=.49\textwidth]{{ng_1.5}.pdf}\hfill
\includegraphics[width=.49\textwidth]{{ng_5.5}.pdf}
\end{center}
\subsubsection{Number of Gluons inside Jets}
\label{sec:psp_ng_jet}
For each of the $n_g$ gluon emissions we can split the phase-space
integral into a (disconnected) region inside the jets and a remainder:
\begin{equation}
\label{eq:psp_split}
\int [\mathrm{d}p_i] = \int [\mathrm{d}p_i]\,
\theta\bigg(\sum_{j=1}^{m}{\cal J}_j(p_i)\bigg) + \int [\mathrm{d}p_i]\,
\bigg[1-\theta\bigg(\sum_{j=1}^{m}{\cal J}_j(p_i)\bigg)\bigg]\,.
\end{equation}
The next step is to decide how many of the gluons will form part of a
jet. This is done in the \lstinline!PhaseSpacePoint::sample_ng_jets!
function.
We choose an importance sampling which is flat in the plane
spanned by the azimuthal angle $\phi$ and the rapidity $y$. This is
observed in BFKL and valid in the limit of Multi-Regge-Kinematics
(MRK). Furthermore, we assume anti-$k_t$ jets, which cover an area of
$\pi R^2$.
In principle, the total accessible area in the $y$-$\phi$ plane is given
by $2\pi \Delta y_{fb}$, where $\Delta y_{fb}\geq \Delta y_{j_f j_b}$ is
the a priori unknown rapidity separation between the most forward and
backward partons. In most cases the extremal jets consist of single
partons, so that $\Delta y_{fb} = \Delta y_{j_f j_b}$. For the less common
case of two partons forming a jet we observe a maximum distance of $R$
between the constituents and the jet centre. In rare cases jets have
more than two constituents. Empirically, they are always within a
distance of $\tfrac{5}{3}R$ to the centre of the jet, so
$\Delta y_{fb} \leq \Delta y_{j_f j_b} + \tfrac{10}{3} R$. In practice, the
extremal partons are required to carry a large fraction of the jet
transverse momentum and will therefore be much closer to the jet axis.
In summary, for sufficiently large rapidity separations we can use the
approximation $\Delta y_{fb} \approx \Delta y_{j_f j_b}$. This scenario
is depicted here:
\begin{center}
\includegraphics[width=0.5\linewidth]{ps_large_y}
\end{center}
If there is no overlap between jets, the probability $p_{\cal J, >}$ for
an extra gluon to end up inside a jet is then given by
\begin{equation}
\label{eq:p_J_large}
p_{\cal J, >} = \frac{(m - 1)\*R^2}{2\Delta y_{j_f j_b}}\,.
\end{equation}
For a very small rapidity separation, eq.~\eqref{eq:p_J_large}
obviously overestimates the true probability. The maximum phase space
covered by jets in the limit of a vanishing rapidity distance between
all partons is $2mR \Delta y_{fb}$:
\begin{center}
\includegraphics[width=0.5\linewidth]{ps_small_y}
\end{center}
We therefore estimate the probability for a parton to end up inside a jet as
\begin{equation}
\label{eq:p_J}
p_{\cal J} = \min\bigg(\frac{(m - 1)\*R^2}{2\Delta y_{j_f j_b}}, \frac{mR}{\pi}\bigg)\,.
\end{equation}
Here we compare this estimate with the actually observed
fraction of additional emissions into jets as a function of the rapidity
separation:
\begin{center}
\includegraphics[width=0.75\linewidth]{pJ}
\end{center}
\subsubsection{Gluons outside Jets}
\label{sec:gluons_nonjet}
Using our estimate for the probability of a gluon to be a jet
constituent, we choose a number $n_{g,{\cal J}}$ of gluons inside
jets, which also fixes the number $n_g - n_{g,{\cal J}}$ of gluons
outside jets. As explained later on, we need to generate the momenta of
the gluons outside jets first. This is done in
\lstinline!PhaseSpacePoint::gen_non_jet!.
The azimuthal angle $\phi$ is generated flat within $0\leq \phi \leq 2
\pi$. The allowed rapidity interval is set by the most forward and
backward partons, which are necessarily inside jets. Since these parton
rapidities are not known at this point, we also have to postpone the
rapidity generation for the gluons outside jets. For the scalar
transverse momentum $p_\perp = |\mathbf{p}_\perp|$ of a gluon outside
jets we use the parametrisation
\begin{equation}
\label{eq:p_nonjet}
p_\perp = \lambda + \tilde{p}_\perp\*\tan(\tau\*r)\,, \qquad
\tau = \arctan\bigg(\frac{p_{\perp{\cal J}_\text{min}} - \lambda}{\tilde{p}_\perp}\bigg)\,.
\end{equation}
For $r \in [0,1)$, $p_\perp$ is always less than the minimum momentum
$p_{\perp{\cal J}_\text{min}}$ required for a jet. $\tilde{p}_\perp$ is
a free parameter, a good empirical value is $\tilde{p}_\perp = [1.3 +
0.2\*(n_g - n_{g,\cal J})]\,$GeV
\subsubsection{Resummation jet momenta}
\label{sec:ptj_res}
On the one hand, each jet momentum is given by the sum of its
constituent momenta. On the other hand, the resummation jet momenta are
fixed by the constraints in line five of the master
equation~\eqref{eq:resumdijetFKLmatched2}. We therefore have to
calculate the resummation jet momenta from these constraints before
generating the momenta of the gluons inside jets. This is done in
\lstinline!PhaseSpacePoint::reshuffle! and in the free
\lstinline!resummation_jet_momenta! function.
The resummation jet momenta are determined by the $\delta$ functions in
line five of eq.~(\ref{eq:resumdijetFKLmatched2}). The rapidities are
fixed to the rapidities of the jets in the input fixed-order events, so
that the FKL ordering is guaranteed to be preserved. For the transverse
momentum components we currently use the traditional \HEJ reshuffling
relation
\begin{equation}
\label{eq:ptreassign}
\mathbf{p}^B_{\mathcal{J}_{l\perp}} = \mathbf{j}_{l\perp} \equiv \mathbf{p}_{\mathcal{J}_{l}\perp} + \mathbf{q}_\perp \,
\frac{|\mathbf{p}_{\mathcal{J}_{l}\perp}|}{P_\perp},
\end{equation}
where $\mathbf{q}_\perp = \sum_{j=1}^n \mathbf{p}_{i\perp}
\bigg[1-\theta\bigg(\sum_{j=1}^{m}{\cal J}_j(p_i)\bigg)\bigg] $ is the
total transverse momentum of all partons \emph{outside} jets and
$P_\perp = \sum_{j=1}^m |\mathbf{p}_{\mathcal{J}_{j}\perp}|$. Since the
total transverse momentum of an event vanishes, we can also use
$\mathbf{q}_\perp = - \sum_{j=1}^m
\mathbf{p}_{\mathcal{J}_{j}\perp}$. Eq.~(\ref{eq:ptreassign}) is a
non-linear system of equations in the resummation jet momenta
$\mathbf{p}_{\mathcal{J}_{l}\perp}$. To solve it we use
\href{https://www.gnu.org/software/gsl/}{GSL} routines.
The reshuffling relation~\eqref{eq:ptreassign} allows the transverse
momenta $p^B_{\mathcal{J}_{l\perp}}$ of the fixed-order jets to be
somewhat below the minimum transverse momentum of resummation jets. It
is crucial that this difference does not become too large, as the
fixed-order cross section diverges for vanishing transverse momenta. In
the production of a Higgs boson with resummation jets above $30\,$GeV we observe
that the contribution from fixed-order events with jets softer than
about $20\,$GeV can be safely neglected. This is shown in the following
plot of the differential cross section over the transverse momentum of
the softest fixed-order jet:
\begin{center}
\includegraphics[width=.75\textwidth]{ptBMin}
\end{center}
Finally, to perform the integration over the resummation phase space, we
need the $\delta$ functions
$\delta^{(2)}(\mathbf{p}_{\mathcal{J}_{l}\perp}^B -
\mathbf{j}_{l\perp})$ form eq.~(\ref{eq:resumdijetFKLmatched2}) to be linear in the resummation parton momenta
instead of the fixed-order momenta. This introduces an extra Jacobian $\Delta$:
\begin{equation}
\label{eq:delta_rewrite}
\prod_{l=1}^m \delta^{(2)}(\mathbf{p}_{\mathcal{J}_{l\perp}}^B -
\mathbf{j}_{l\perp}) = \frac{1}{\Delta} \prod_{l=1}^m \delta^{(2)}(\mathbf{p}_{\mathcal{J}_{l\perp}} -
\mathbf{j}_{l\perp}^B)\,,
\end{equation}
where $\mathbf{j}_{l\perp}^B$ only depends on the Born momenta. We have
extended the product to run to $m$ instead of $m-1$ by eliminating the
last $\delta$ function $\delta^{(2)}\!\!\left(\sum_{k=1}^n \mathbf{p}_{k\perp}\right )$.
The Jacobian $\Delta$ is the determinant of a $2m \times 2m$ matrix with $l, l' = 1,\dots,m$
and $X, X' = x,y$.
\begin{equation}
\label{eq:jacobian}
\Delta = \bigg| \delta_{l l'} \delta_{X X'} + \frac{q_X\, p_{{\cal
J}_{l'}X'}}{|\mathbf{p}_{{\cal J}_{l'} \perp}| P_\perp}\bigg(\delta_{l l'}
- \frac{|\mathbf{p}_{{\cal J}_l \perp}|}{P_\perp}\bigg)\bigg|\,.
\end{equation}
\subsubsection{Gluons inside Jets}
\label{sec:gluons_jet}
After the steps outlined in section~\ref{sec:psp_ng_jet}, we have a
total number of $m + n_{g,{\cal J}}$ constituents. In
\lstinline!PhaseSpacePoint::distribute_jet_partons! we distribute them
randomly among the jets such that each jet has at least one
constituent. We then generate their momenta in
\lstinline!PhaseSpacePoint::split! using the \lstinline!Splitter! class.
The phase space integral for a jet ${\cal J}$ is given by
\begin{equation}
\label{eq:ps_jetparton} \prod_{i\text{ in }{\cal J}} \bigg(\int
\mathrm{d}\mathbf{p}_{i\perp}\ \int \mathrm{d} y_i
\bigg)\delta^{(2)}\Big(\sum_{i\text{ in }{\cal J}} \mathbf{p}_{i\perp} -
\mathbf{j}_{\perp}^B\Big)\delta(y_{\mathcal{J}}-y^B_{\mathcal{J}})\,.
\end{equation}
For jets with a single constituent, the parton momentum is obiously equal to the
jet momentum. In the case of two constituents, we observe that the
partons are always inside the jet cone with radius $R$ and often very
close to the jet centre. The following plots show the typical relative
distance $\Delta R/R$ for this scenario:
\begin{center}
\includegraphics[width=0.45\linewidth]{dR_2}
\includegraphics[width=0.45\linewidth]{dR_2_small}
\end{center}
According to this preference for small values of $\Delta R$, we
parametrise the $\Delta R$ integrals as
\begin{equation}
\label{eq:dR_sampling}
\frac{\Delta R}{R} =
\begin{cases}
0.25\,x_R & x_R < 0.4 \\
1.5\,x_R - 0.5 & x_R \geq 0.4
\end{cases}\,.
\end{equation}
Next, we generate $\Theta_1 \equiv \Theta$ and use the constraint $\Theta_2 = \Theta
\pm \pi$. The transverse momentum of the first parton is then given by
\begin{equation}
\label{eq:delta_constraints}
p_{1\perp} =
\frac{p_{\mathcal{J} y} - \tan(\phi_2) p_{\mathcal{J} x}}{\sin(\phi_1)
- \tan(\phi_2)\cos(\phi_1)}\,.
\end{equation}
We get $p_{2\perp}$ by exchanging $1 \leftrightarrow 2$ in the
indices. To obtain the Jacobian of the transformation, we start from the
single jet phase space eq.~(\ref{eq:ps_jetparton}) with the rapidity
delta function already rewritten to be linear in the rapidity of the
last parton, i.e.
\begin{equation}
\label{eq:jet_2p}
\prod_{i=1,2} \bigg(\int
\mathrm{d}\mathbf{p}_{i\perp}\ \int \mathrm{d} y_i
\bigg)\delta^{(2)}\Big(\mathbf{p}_{1\perp} + \mathbf{p}_{2\perp} -
\mathbf{j}_{\perp}^B\Big)\delta(y_2- \dots)\,.
\end{equation}
The integral over the second parton momentum is now trivial; we can just replace
the integral over $y_2$ with the equivalent constraint
\begin{equation}
\label{eq:R2}
\int \mathrm{d}R_2 \ \delta\bigg(R_2 - \bigg[\phi_{\cal J} - \arctan
\bigg(\frac{p_{{\cal J}y} - p_{1y}}{p_{{\cal J}x} -
p_{1x}}\bigg)\bigg]/\cos \Theta\bigg) \,.
\end{equation}
In order to fix the integral over $p_{1\perp}$ instead, we rewrite this
$\delta$ function. This introduces the Jacobian
\begin{equation}
\label{eq:jac_pt1}
\bigg|\frac{\partial p_{1\perp}}{\partial R_2} \bigg| =
\frac{\cos(\Theta)\mathbf{p}_{2\perp}^2}{p_{{\cal J}\perp}\sin(\phi_{\cal J}-\phi_1)}\,.
\end{equation}
The final form of the integral over the two parton momenta is then
\begin{equation}
\label{eq:ps_jet_2p}
\int \mathrm{d}R_1\ R_1 \int \mathrm{d}R_2 \int \mathrm{d}x_\Theta\ 2\pi \int
\mathrm{d}p_{1\perp}\ p_{1\perp} \int \mathrm{d}p_{2\perp}
\ \bigg|\frac{\partial p_{1\perp}}{\partial R_2} \bigg|\delta(p_{1\perp}
-\dots) \delta(p_{2\perp} - \dots)\,.
\end{equation}
As is evident from section~\ref{sec:psp_ng_jet}, jets with three or more
constituents are rare and an efficient phase-space sampling is less
important. For such jets, we exploit the observation that partons with a
distance larger than $R_{\text{max}} = \tfrac{5}{3} R$ to
the jet centre are never clustered into the jet. Assuming $N$
constituents, we generate all components
for the first $N-1$ partons and fix the remaining parton with the
$\delta$-functional. In order to end up inside the jet, we use the
parametrisation
\begin{align}
\label{eq:ps_jet_param}
\phi_i ={}& \phi_{\cal J} + \Delta \phi_i\,, & \Delta \phi_i ={}& \Delta
R_i
\cos(\Theta_i)\,, \\
y_i ={}& y_{\cal J} + \Delta y_i\,, & \Delta y_i ={}& \Delta
R_i
\sin(\Theta_i)\,,
\end{align}
and generate $\Theta_i$ and $\Delta R_i$ randomly with $\Delta R_i \leq
R_{\text{max}}$ and the empiric value $R_{\text{max}} = 5\*R/3$. We can
then write the phase space integral for a single parton as $(p_\perp = |\mathbf{p}_\perp|)$
\begin{equation}
\label{eq:ps_jetparton_x}
\int \mathrm{d}\mathbf{p}_{\perp}\ \int
\mathrm{d} y \approx \int_{\Box} \mathrm{d}x_{\perp}
\mathrm{d}x_{ R}
\mathrm{d}x_{\theta}\
2\*\pi\,\*R_{\text{max}}^2\,\*x_{R}\,\*p_{\perp}\,\*(p_{\perp,\text{max}}
- p_{\perp,\text{min}})
\end{equation}
with
\begin{align}
\label{eq:ps_jetparton_parameters}
\Delta \phi ={}& R_{\text{max}}\*x_{R}\*\cos(2\*\pi\*x_\theta)\,,&
\Delta y ={}& R_{\text{max}}\*x_{R}\*\sin(2\*\pi\*x_\theta)\,, \\
p_{\perp} ={}& (p_{\perp,\text{max}} - p_{\perp,\text{min}})\*x_\perp +
p_{\perp,\text{min}}\,.
\end{align}
$p_{\perp,\text{max}}$ is determined from the requirement that the total
contribution from the first $n-1$ partons --- i.e. the projection onto the
jet $p_{\perp}$ axis --- must never exceed the jet $p_\perp$. This gives
\todo{This bound is too high}
\begin{equation}
\label{eq:pt_max}
p_{i\perp,\text{max}} = \frac{p_{{\cal J}\perp} - \sum_{j<i} p_{j\perp}
\cos \Delta
\phi_j}{\cos \Delta
\phi_i}\,.
\end{equation}
The $x$ and $y$ components of the last parton follow immediately from
the first $\delta$ function. The last rapidity is fixed by the condition that
the jet rapidity is kept fixed by the reshuffling, i.e.
\begin{equation}
\label{eq:yJ_delta}
y^B_{\cal J} = y_{\cal J} = \frac 1 2 \ln \frac{\sum_{i=1}^n E_i+ p_{iz}}{\sum_{i=1}^n E_i - p_{iz}}\,.
\end{equation}
With $E_n \pm p_{nz} = p_{n\perp}\exp(\pm y_n)$ this can be rewritten to
\begin{equation}
\label{eq:yn_quad_eq}
\exp(2y_{\cal J}) = \frac{\sum_{i=1}^{n-1} E_i+ p_{iz}+p_{n\perp} \exp(y_n)}{\sum_{i=1}^{n-1} E_i - p_{iz}+p_{n\perp} \exp(-y_n)}\,,
\end{equation}
which is a quadratic equation in $\exp(y_n)$. The physical solution is
\begin{align}
\label{eq:yn}
y_n ={}& \log\Big(-b + \sqrt{b^2 + \exp(2y_{\cal J})}\,\Big)\,,\\
b ={}& \bigg(\sum_{i=1}^{n-1} E_i + p_{iz} - \exp(2y_{\cal J})
\sum_{i=1}^{n-1} E_i - p_{iz}\bigg)/(2 p_{n\perp})\,.
\end{align}
\todo{what's wrong with the following?} To eliminate the remaining rapidity
integral, we transform the $\delta$ function to be linear in the
rapidity $y$ of the last parton. The corresponding Jacobian is
\begin{equation}
\label{eq:jacobian_y}
\bigg|\frac{\partial y_{\cal J}}{\partial y_n}\bigg|^{-1} = 2 \bigg( \frac{E_n +
p_{nz}}{E_{\cal J} + p_{{\cal J}z}} + \frac{E_n - p_{nz}}{E_{\cal J} -
p_{{\cal J}z}}\bigg)^{-1}\,.
\end{equation}
Finally, we check that all designated constituents are actually
clustered into the considered jet.
\subsubsection{Final steps}
\label{sec:final}
Knowing the rapidity span covered by the extremal partons, we can now
generate the rapdities for the partons outside jets. We perform jet
clustering on all partons and check in
\lstinline!PhaseSpacePoint::jets_ok! that all the following criteria are
fulfilled:
\begin{itemize}
\item The number of resummation jets must match the number of
fixed-order jets.
\item No partons designated to be outside jets may end up inside jets.
\item All other outgoing partons \emph{must} end up inside jets.
\item The extremal (in rapidity) partons must be inside the extremal
jets. If there is, for example, an unordered forward emission, the
most forward parton must end up inside the most forward jet and the
next parton must end up inside second jet.
\item The rapidities of fixed-order and resummation jets must match.
\end{itemize}
After this, we adjust the phase-space normalisation according to the
third line of eq.~(\ref{eq:resumdijetFKLmatched2}), determine the
flavours of the outgoing partons, and adopt any additional colourless
bosons from the fixed-order input event. Finally, we use momentum
conservation to reconstruct the momenta of the incoming partons.
\subsection{The matrix element }
\label{sec:ME}
The derivation of the \HEJ matrix element is explained in some detail
in~\cite{Andersen:2017kfc}, where also results for leading and
subleading matrix elements for pure multijet production and production
of a Higgs boson with at least two associated jets are listed. Matrix
elements for $W$ and $Z/\gamma^*$ production together with jets are
given in~\cite{Andersen:2012gk,Andersen:2016vkp}, but not yet included.
The matrix elements are implemented in the \lstinline!MatrixElement!
class. To discuss the structure, let us consider the matrix element for
FKL multijet production:
\begin{align}
\label{eq:ME}
\begin{split}
\overline{\left|\mathcal{M}_\text{HEJ}^{f_1 f_2 \to f_1
g\cdots g f_2}\right|}^2 = \ &\frac 1 {4\
(N_c^2-1)}\ \left\|S_{f_1 f_2\to f_1 f_2}\right\|^2\\
&\cdot\ \left(g_s^2 K_{f_1}\ \frac 1 {t_1}\right) \cdot\ \left(g_s^2\ K_{f_2}\ \frac 1
{t_{n-1}}\right)\\
& \cdot \prod_{i=1}^{n-2} \left( \frac{-g_s^2 C_A}{t_it_{i+1}}\
V^\mu(q_i,q_{i+1})V_\mu(q_i,q_{i+1}) \right)\\
& \cdot \prod_{j=1}^{n-1} \exp\left[\omega^0(q_{j\perp})(y_{j+1}-y_j)\right].
\end{split}
\end{align}
The square of the complete matrix element as given in eq.~(\ref{eq:ME})
-is calculated by \lstinline!MatrixElement::operator()!. \todo{There is
-most likely some overall factor}\todo{I haven't found any. Double-check - MH}
-The last line of eq.~\eqref{eq:ME}
-constitutes the all-order virtual correction, implemented in
-\lstinline!MatrixElement::virtual_corrections!. The remaining parts,
-which correspond to the square of the leading-order HEJ matrix element
-$\overline{\left|\mathcal{M}_\text{LO, HEJ}^{f_1f_2\to f_1g\cdots
-gf_2}\big(\big\{p^B_j\big\}\big)\right|}^{2}$, are computed in
-\lstinline!MatrixElement::tree!. We can further factor off the
+is calculated by \lstinline!MatrixElement::operator()!. The last line
+of eq.~\eqref{eq:ME} constitutes the all-order virtual correction,
+implemented in \lstinline!MatrixElement::virtual_corrections!. The
+remaining parts, which correspond to the square of the leading-order HEJ
+matrix element $\overline{\left|\mathcal{M}_\text{LO, HEJ}^{f_1f_2\to
+f_1g\cdots gf_2}\big(\big\{p^B_j\big\}\big)\right|}^{2}$, are computed
+in \lstinline!MatrixElement::tree!. We can further factor off the
scale-dependent ``parametric'' part
\lstinline!MatrixElement::tree_param! containing all factors of the
strong coupling $g_s$. Using this function saves some CPU time when
adjusting the renormalisation scale, see section~\ref{sec:resum}.
The remaining ``kinematic'' factors are calculated in
\lstinline!MatrixElement::kin!. The currents $\frac{K_{f_1}K_{f_2}}{t_1
t_{n-1}}\left\|S_{f_1 f_2\to f_1 f_2}\right\|^2$ are implemented in the
\lstinline!jM2!$\dots$ functions of \texttt{src/Currents.cc}.
\section{The fixed-order generator}
\label{sec:HEJFOG}
Even at leading order, standard fixed-order generators can only generate
events with a limited number of final-state particles within reasonable
CPU time. The purpose of the fixed-order generator is to supplement this
with high-multiplicity input events according to the first two lines of
eq.~\eqref{eq:resumdijetFKLmatched2} with the \HEJ approximation
$\mathcal{M}_\text{LO, HEJ}^{f_1f_2\to f_1g\cdots gf_2}$ instead of the
full fixed-order matrix element $\mathcal{M}_\text{LO}^{f_1f_2\to
f_1g\cdots gf_2}$. Its usage is described in the user
documentation\todo{link}.
\subsection{File structure}
\label{sec:HEJFOG_structure}
The code for the fixed-order generator is in the \texttt{FixedOrderGen}
directory, which contains the following:
\begin{description}
\item[include:] Contains the C++ header files.
\item[src:] Contains the C++ source files.
\item[t:] Contains the source code for the automated tests.
\item[CMakeLists.txt:] Configuration file for the \cmake build system.
\item[configFO.yml:] Sample configuration file for the fixed-order generator.
\end{description}
The code is generally in the \lstinline!HEJFOG! namespace. Functions and
classes \lstinline!MyClass! are usually declared in
\texttt{include/MyClass.hh} and implemented in \texttt{src/MyClass.cc}.
\subsection{Program flow}
\label{sec:prog_flow}
A single run of the fixed-order generator consists of three or four
stages.
First, we perform initialisation similar to reversed HEJ, see
section~\ref{sec:init}. Since there is a lot of overlap we frequently
reuse classes and functions from reversed HEJ, i.e. from the
\lstinline!RHEJ! namespace. The code for parsing the configuration file
is in \texttt{include/config.hh} and implemented in
\texttt{src/config.cc}.
If partial unweighting is requested in the user settings \todo{link},
the initialisation is followed by a calibration phase. We use a
\lstinline!EventGenerator! to produce a number of trial
events. We use these to calibrate the \lstinline!Unweighter! in
its constructor and produce a first batch of partially unweighted
events. This also allows us to estimate our unweighting efficiency.
In the next step, we continue to generate events and potentially
unweight them. Once the user-defined target number of events is reached,
we adjust their weights according to the number of required trials. As
in reversed HEJ (see section~\ref{sec:processing}), we pass the final
events to a \lstinline!RHEJ::Analysis! and a
\lstinline!RHEJ::CombinedEventWriter!.
\subsection{Event generation}
\label{sec:evgen}
Event generation is performed by the
\lstinline!EventGenerator::gen_event! member function. We begin by generating a
\lstinline!PhaseSpacePoint!. This is not to be confused with
the resummation phase space points represented by
\lstinline!RHEJ::PhaseSpacePoint!! After jet clustering, we compute the
leading-order matrix element (see section~\ref{sec:ME}) and pdf factors.
The phase space point generation is performed in the
\lstinline!PhaseSpacePoint! constructor. We first construct the
user-defined number of $n_p$ partons (by default gluons) in
\lstinline!PhaseSpacePoint::gen_LO_partons!. We use flat sampling in
rapidity and azimuthal angle. For the scalar transverse momenta, we
distinguish between two cases. By default, they are generated based on a
random variable $x_{p_\perp}$ according to
\begin{equation}
\label{eq:pt_sampling}
p_\perp = p_{\perp,\text{min}} +
\begin{cases}
p_{\perp,\text{par}}
\tan\left(
x_{p_\perp}
\arctan\left(
\frac{p_{\perp,\text{max}} - p_{\perp,\text{min}}}{p_{\perp,\text{par}}}
\right)
\right)
& y < y_\text{cut}
\\
- \tilde{p}_{\perp,\text{par}}\log\left(1 - x_{p_\perp}\left[1 -
\exp\left(\frac{p_{\perp,\text{min}} -
p_{\perp,\text{max}}}{\tilde{p}_{\perp,\text{par}}}\right)\right]\right)
& y \geq y_\text{cut}
\end{cases}\,,
\end{equation}
where $p_{\perp,\text{min}}$ is the minimum jet transverse momentum,
$p_{\perp,\text{max}}$ is the maximum transverse parton momentum,
tentatively set to the beam energy, and $y_\text{cut}$, $p_{\perp,\text{par}}$
and $\tilde{p}_{\perp,\text{par}}$ are generation parameters set to
heuristically determined values of
\begin{align}
y_\text{cut}&=3,\\
p_{\perp,\text{par}}&=p_{\perp,\min}+\frac{n_p}{5}, \\
\tilde{p}_{\perp,\text{par}}&=\frac{p_{\perp,\text{par}}}{1 +
5(y-y_\text{cut})}.
\end{align}
The problem with this generation is that the transverse momenta peak at
the minimum transverse momentum required for fixed-order jets. However,
if we use the generated events as input for \HEJ resummation, events
with such soft transverse momenta hardly contribute, see
section~\ref{sec:ptj_res}. To generate efficient input for resummation,
there is the user option \texttt{peak pt}, which specifies the
dominant transverse momentum for resummation jets. If this option is
set, most jets will be generated as above, but with
$p_{\perp,\text{min}}$ set to the peak transverse momentum $p_{\perp,
\text{peak}}$. In addition, there is a small chance of around $2\%$ to
generate softer jets. The heuristic ansatz for the transverse momentum
distribution in the ``soft'' region is
\begin{equation}
\label{FO_pt_soft}
\frac{\partial \sigma}{\partial p_\perp} \propto e^{n_p\frac{p_\perp- p_{\perp,
\text{peak}}}{\bar{p}_\perp}}\,,
\end{equation}
where $n_p$ is the number of partons and $\bar{p}_\perp \approx
4\,$GeV. To achieve this distribution, we use
\begin{equation}
\label{eq:FO_pt_soft_sampling}
p_\perp = p_{\perp, \text{peak}} + \bar{p}_\perp \frac{\log x_{p_\perp}}{n_p}
\end{equation}
and discard the phase space point if the parton is too soft, i.e. below the threshold for
fixed-order jets.
After ensuring that all partons form separate jets, we generate any
potential colourless emissions. We then determine the incoming momenta
and flavours in \lstinline!PhaseSpacePoint::reconstruct_incoming! and
adjust the outgoing flavours to ensure an FKL configuration. Finally, we
may reassign outgoing flavours to generate suppressed (for example
unordered) configurations.
\subsection{Unweighting}
\label{sec:unweight}
Straightforward event generation tends to produce many events with small
weights. Those events have a negligible contribution to the final
observables, but can take up considerable storage space and CPU time in
later processing stages. This problem can be addressed by unweighting.
For naive unweighting, one would determine the maximum weight
$w_\text{max}$ of all events, discard each event with weight $w$ with a
probability $p=w/w_\text{max}$, and set the weights of all remaining
events to $w_\text{max}$. The downside to this procedure is that it also
eliminates a sizeable fraction of events with moderate weight, so that
the statistical convergence deteriorates.
To ameliorate this problem, we perform unweighting only for events with
sufficiently small weights. This is done by the
\lstinline!Unweighter! class. In the constructor we estimate the
mean and width of the weight-weight distribution from a sample of
events. We use these estimates to determine the maximum weight below
which unweighting is performed. The actual unweighting is the done in
the \lstinline!Unweighter::unweight! function.
\bibliographystyle{JHEP}
\bibliography{biblio}
\end{document}
diff --git a/src/HepMCWriter.cc b/src/HepMCWriter.cc
index dc53e9d..300a464 100644
--- a/src/HepMCWriter.cc
+++ b/src/HepMCWriter.cc
@@ -1,145 +1,139 @@
#include "RHEJ/HepMCWriter.hh"
#include <cassert>
#ifdef RHEJ_BUILD_WITH_HepMC_VERSION
#if RHEJ_BUILD_WITH_HepMC_VERSION >= 3
#include "HepMC/WriterAscii.h"
#include "HepMC/LHEFAttributes.h"
#include "RHEJ/Version.hh"
#else
#include "HepMC/IO_GenEvent.h"
#endif
#include "HepMC/GenVertex.h"
#include "HepMC/GenParticle.h"
#include "RHEJ/utility.hh"
#include "RHEJ/Event.hh"
#include "RHEJ/exceptions.hh"
#include "RHEJ/HepMCInterface.hh"
#if RHEJ_BUILD_WITH_HepMC_VERSION >= 3
namespace {
- void add_generator_tag(LHEF::HEPRUP & heprup){
- heprup.generators.emplace_back(LHEF::XMLTag{});
- heprup.generators.back().name = RHEJ::Version::package_name();
- heprup.generators.back().version = RHEJ::Version::String();
- }
void reset_weight_info(LHEF::HEPRUP & heprup){
heprup.IDWTUP = 2;
// use placeholders for unknown init block values
// we can overwrite them after processing all events
heprup.XSECUP = {0.};
heprup.XERRUP = {0.};
heprup.XMAXUP = {0.};
}
HepMC::shared_ptr<HepMC::GenRunInfo> init_runinfo(LHEF::HEPRUP && heprup){
- add_generator_tag(heprup);
reset_weight_info(heprup);
HepMC::GenRunInfo runinfo;
auto hepr = HepMC::make_shared<HepMC::HEPRUPAttribute>();
hepr->heprup = heprup;
runinfo.add_attribute(std::string("HEPRUP"), hepr);
for (int i = 0, N = hepr->heprup.generators.size(); i < N; ++i ){
HepMC::GenRunInfo::ToolInfo tool;
tool.name = hepr->heprup.generators[i].name;
tool.version = hepr->heprup.generators[i].version;
tool.description = hepr->heprup.generators[i].contents;
runinfo.tools().push_back(tool);
}
return HepMC::make_shared<HepMC::GenRunInfo>(runinfo);
}
} // namespace anonymous
#endif // HepMC 3
namespace RHEJ{
struct HepMCWriter::HepMCWriterImpl{
HepMCInterface hepmc_;
HepMCWriterImpl & operator=(HepMCWriterImpl const & other) = delete;
HepMCWriterImpl(HepMCWriterImpl const & other) = delete;
HepMCWriterImpl & operator=(HepMCWriterImpl && other) = delete;
HepMCWriterImpl(HepMCWriterImpl && other) = delete;
#if RHEJ_BUILD_WITH_HepMC_VERSION >= 3
HepMC::WriterAscii writer_;
HepMCWriterImpl(
std::string const & file, LHEF::HEPRUP && heprup
):
hepmc_(),
writer_{file, init_runinfo(std::move(heprup))}
{}
~HepMCWriterImpl(){
writer_.close();
}
#else // HepMC 2
HepMC::IO_GenEvent writer_;
HepMCWriterImpl(
std::string const & file, LHEF::HEPRUP &&
):
hepmc_(),
writer_{file}
{}
#endif
void write(Event const & ev){
auto out_ev = hepmc_(ev);
#if RHEJ_BUILD_WITH_HepMC_VERSION >= 3
writer_.write_event(out_ev);
#else // HepMC 2
writer_.write_event(&out_ev);
#endif
}
};
void HepMCWriter::HepMCWriterImplDeleter::operator()(HepMCWriterImpl* p) {
delete p;
}
HepMCWriter::HepMCWriter(std::string const & file, LHEF::HEPRUP heprup):
impl_{std::unique_ptr<HepMCWriterImpl, HepMCWriterImplDeleter>{
new HepMCWriterImpl(file, std::move(heprup))
}}
{}
void HepMCWriter::write(Event const & ev){
impl_->write(ev);
}
} // namespace RHEJ
#else // no HepMC
namespace RHEJ{
class HepMCWriter::HepMCWriterImpl{};
HepMCWriter::HepMCWriter(std::string const &, LHEF::HEPRUP){
throw std::invalid_argument(
"Failed to create HepMC writer: "
"Reversed HEJ was built without HepMC support"
);
}
void HepMCWriter::write(Event const &){
assert(false);
}
void HepMCWriter::HepMCWriterImplDeleter::operator()(HepMCWriterImpl* p) {
delete p;
}
}
#endif
diff --git a/src/LesHouchesWriter.cc b/src/LesHouchesWriter.cc
index 4482b49..57bd7bd 100644
--- a/src/LesHouchesWriter.cc
+++ b/src/LesHouchesWriter.cc
@@ -1,93 +1,90 @@
#include <stdexcept>
#include <memory>
#include <cassert>
#include "RHEJ/LesHouchesWriter.hh"
#include "RHEJ/event_types.hh"
#include "RHEJ/Version.hh"
#include "RHEJ/Event.hh"
namespace RHEJ{
namespace{
template<class T, class... Args>
std::unique_ptr<T> make_unique(Args&&... a){
return std::unique_ptr<T>{new T{std::forward<Args>(a)...}};
}
}
LesHouchesWriter::LesHouchesWriter(
std::string const & file, LHEF::HEPRUP heprup
):
out_{file, std::fstream::in | std::fstream::out | std::fstream::trunc},
writer_{RHEJ::make_unique<LHEF::Writer>(out_)}
{
if(! out_.is_open()){
throw std::ios_base::failure("Failed to open " + file);
};
writer_->heprup = std::move(heprup);
// lhe Stardard: IDWTUP (negative => weights = +/-)
// 1: weighted events, xs = mean(weight), XMAXUP given
// 2: weighted events, xs = XSECUP, XMAXUP given
// 3: unweighted events, no additional information given
// 4: unweighted events, xs = mean(weight), no additional information given
writer_->heprup.IDWTUP = 2;
- writer_->heprup.generators.emplace_back(LHEF::XMLTag{});
- writer_->heprup.generators.back().name = Version::package_name();
- writer_->heprup.generators.back().version = Version::String();
// use placeholders for unknown init block values
// we can overwrite them after processing all events
writer_->heprup.XSECUP = std::vector<double>(event_type::last_type+1, 0.);
writer_->heprup.XERRUP = std::vector<double>(event_type::last_type+1, 0.);
writer_->heprup.XMAXUP = std::vector<double>(event_type::last_type+1, 0.);
write_init();
}
void LesHouchesWriter::write(Event const & ev){
assert(writer_ && out_.is_open());
const double wt = ev.central().weight;
writer_->hepeup = RHEJ::to_HEPEUP(std::move(ev), &heprup());
writer_->writeEvent();
heprup().XSECUP[ev.type()] += wt;
heprup().XERRUP[ev.type()] += wt*wt;
if(wt > heprup().XMAXUP[ev.type()]){
heprup().XMAXUP[ev.type()] = wt;
}
}
// this function is called after overwritting the Les Houches init block
// assert that we have overwritten *exactly* the init block,
// i.e. an intact event block is next
void assert_next_event_intact(std::istream & out){
(void) out; // suppress compiler warnings if not in debug mode
#ifndef NDEBUG
std::string line;
getline(out, line);
assert(line == "<event>");
#endif
}
void LesHouchesWriter::rewrite_init(){
assert(writer_ && out_.is_open());
// replace placeholder entries
const auto pos = out_.tellp();
out_.seekp(0);
writer_->init();
assert_next_event_intact(out_);
out_.seekp(pos);
}
LesHouchesWriter::~LesHouchesWriter(){
assert(writer_ && out_.is_open());
for(auto & xs_err: heprup().XERRUP)
{
xs_err = sqrt(xs_err);
}
rewrite_init();
}
}
diff --git a/src/bin/rHEJ.cc b/src/bin/rHEJ.cc
index e535ae3..ec61201 100644
--- a/src/bin/rHEJ.cc
+++ b/src/bin/rHEJ.cc
@@ -1,174 +1,178 @@
/**
* Name: main.cc
* Authors: Tuomas Hapola, Andreas Maier <andreas.maier@durham.ac.uk>
*
*/
#include <fstream>
#include <numeric>
#include <memory>
#include <chrono>
#include <iostream>
#include "yaml-cpp/yaml.h"
#include "RHEJ/CombinedEventWriter.hh"
#include "RHEJ/config.hh"
#include "RHEJ/EventReweighter.hh"
#include "RHEJ/get_analysis.hh"
#include "LHEF/LHEF.h"
#include "RHEJ/utility.hh"
#include "RHEJ/Version.hh"
#include "RHEJ/stream.hh"
#include "RHEJ/YAMLreader.hh"
#include "RHEJ/ProgressBar.hh"
#include "RHEJ/make_RNG.hh"
int event_number(std::string const & record){
size_t start = record.rfind("Number of Events");
start = record.find_first_of("123456789", start);
if(start == std::string::npos) {
throw std::invalid_argument("no event number record found");
}
const size_t end = record.find_first_not_of("0123456789", start);
return std::stoi(record.substr(start, end - start));
}
RHEJ::Config load_config(char const * filename){
try{
return RHEJ::load_config(filename);
}
catch(std::exception const & exc){
std::cerr << "Error: " << exc.what() << '\n';
std::exit(EXIT_FAILURE);
}
}
std::unique_ptr<RHEJ::Analysis> get_analysis(
YAML::Node const & parameters
){
try{
return RHEJ::get_analysis(parameters);
}
catch(std::exception const & exc){
std::cerr << "Failed to load analysis: " << exc.what() << '\n';
std::exit(EXIT_FAILURE);
}
}
// unique_ptr is a workaround:
// RHEJ::optional is a better fit, but gives spurious errors with g++ 7.3.0
std::unique_ptr<RHEJ::ProgressBar<double>> make_progress_bar(
std::vector<double> const & xs
) {
if(xs.empty()) return {};
const double Born_xs = std::accumulate(begin(xs), end(xs), 0.);
return std::make_unique<RHEJ::ProgressBar<double>>(std::cout, Born_xs);
}
std::string time_to_string(const time_t time){
char s[30];
struct tm * p = localtime(&time);
strftime(s, 30, "%a %b %d %Y %H:%M:%S", p);
return s;
}
int main(int argn, char** argv) {
using clock = std::chrono::system_clock;
if (argn < 3) {
std::cerr << "\n# Usage:\n."<< argv[0] <<" config_file input_file\n\n";
return EXIT_FAILURE;
}
const auto start_time = clock::now();
{
std::cout << "Starting " << RHEJ::Version::package_name_full()
<< ", revision " << RHEJ::Version::revision() << " ("
<< time_to_string(clock::to_time_t(start_time)) << ")" << std::endl;
}
fastjet::ClusterSequence::print_banner();
// read configuration
const RHEJ::Config config = load_config(argv[1]);
RHEJ::istream in{argv[2]};
std::unique_ptr<RHEJ::Analysis> analysis = get_analysis(
config.analysis_parameters
);
assert(analysis != nullptr);
LHEF::Reader reader{in};
- RHEJ::CombinedEventWriter writer{config.output, reader.heprup};
+ auto heprup = reader.heprup;
+ heprup.generators.emplace_back(LHEF::XMLTag{});
+ heprup.generators.back().name = RHEJ::Version::package_name();
+ heprup.generators.back().version = RHEJ::Version::String();
+ RHEJ::CombinedEventWriter writer{config.output, std::move(heprup)};
double global_reweight = 1.;
int max_events = std::numeric_limits<int>::max();
if(argn > 3){
max_events = std::stoi(argv[3]);
const int input_events = event_number(reader.headerBlock);
global_reweight = input_events/static_cast<double>(max_events);
std::cout << "Processing " << max_events
<< " out of " << input_events << " events\n";
}
RHEJ::ScaleGenerator scale_gen{
config.scales.base,
config.scales.factors,
config.scales.max_ratio
};
auto ran = RHEJ::make_RNG(config.rng.name, config.rng.seed);
assert(ran != nullptr);
RHEJ::EventReweighter rhej{
reader.heprup,
std::move(scale_gen),
to_EventReweighterConfig(config),
*ran
};
int nevent = 0;
std::array<int, RHEJ::event_type::last_type + 1>
nevent_type{0}, nfailed_type{0};
auto progress = make_progress_bar(reader.heprup.XSECUP);
// Loop over the events in the inputfile
while(reader.readEvent()){
// reweight events so that the total cross section is conserved
reader.hepeup.setWeight(0, global_reweight * reader.hepeup.weight());
if(nevent == max_events) break;
++nevent;
// calculate rHEJ weight
RHEJ::Event FO_event{
RHEJ::UnclusteredEvent{reader.hepeup},
config.fixed_order_jets.def, config.fixed_order_jets.min_pt,
};
auto resummed_events = rhej.reweight(FO_event, config.trials);
++nevent_type[FO_event.type()];
if(resummed_events.empty()) ++nfailed_type[FO_event.type()];
for(auto const & ev: resummed_events){
//TODO: move pass_cuts to after phase space point generation
if(analysis->pass_cuts(ev, FO_event)){
analysis->fill(ev, FO_event);
writer.write(ev);
}
}
if(progress) progress->increment(FO_event.central().weight);
} // main event loop
std::cout << '\n';
analysis->finalise();
using namespace RHEJ::event_type;
std::cout<< "Events processed: " << nevent << '\n';
for(size_t ev_type = first_type; ev_type <= last_type; ++ev_type){
std::cout << '\t' << names[ev_type] << ": " << nevent_type[ev_type]
<< ", failed to reconstruct " << nfailed_type[ev_type]
<< '\n';
}
std::chrono::duration<double> run_time = (clock::now() - start_time);
std::cout << "Finished " << RHEJ::Version::package_name() << " at "
<< time_to_string(clock::to_time_t(clock::now()))
<< "\n=> Runtime: " << run_time.count() << " sec ("
<< nevent/run_time.count() << " Events/sec).\n";
}
diff --git a/t/ME_data/ME_h_mt.dat b/t/ME_data/ME_h_mt.dat
new file mode 100644
index 0000000..ba23656
--- /dev/null
+++ b/t/ME_data/ME_h_mt.dat
@@ -0,0 +1,807 @@
+3.999227056e-08
+4.193300057e-05
+3.89321526e-05
+0.0168082963
+4.057962883e-06
+4.063245819e-07
+6.411897603e-07
+7.822040516e-07
+6.064514954e-07
+2.202805699e-05
+4.738208766e-05
+1.624045121e-06
+0.0005217822994
+1.536443326e-05
+9.5908372e-06
+0.0339113482
+8.152325116e-07
+608.615459
+5.962078524e-06
+0.01340321821
+1.931630897
+0.0005720381031
+0.0008001259943
+3.300880891e-05
+0.0003564506816
+0.0007037298905
+0.03792377651
+6.513991407e-07
+2.011798913
+8.842231663e-08
+0.0463578742
+0.1064067277
+9.357584012e-07
+0.006517735511
+0.02272391876
+3.606711294e-05
+6.725111513e-07
+0.005897687516
+7.581255699e-05
+0.005354681787
+0.0006088947447
+0.0083356406
+0.003352568463
+3.23438533e-05
+0.0001515363281
+0.005531947882
+0.03875973775
+2.312735565e-07
+5.479259483e-05
+8.480458133
+0.0001854378187
+4.957725793e-06
+5.091269317e-06
+0.0001308303481
+0.1918392997
+0.0238989834
+4.446354691
+0.0002597751966
+1.267033388
+0.005050109129
+2.130648865e-06
+0.05854510061
+0.00280707581
+6.289945688
+7.216599144e-08
+4.779195357e-06
+1.655604389e-05
+1.522656761e-05
+1.116772322e-06
+0.0009539322219
+9.279967349e-07
+2.28564472
+0.001331277201
+0.001585734932
+26972.25526
+8.555976389e-07
+0.0006894623725
+1.065245235e-07
+0.01968534835
+1.951627129e-06
+4.401131811e-05
+0.02561966521
+49705.99132
+5.383396351e-06
+1.07104069e-06
+2.808652999e-06
+0.0001066898908
+0.0002239236027
+7.851273644e-08
+0.002704063747
+0.0005663491795
+0.07550986067
+0.00128549481
+3.590772606e-05
+0.1149401044
+0.183403668
+0.000113873958
+1.730096598e-05
+0.02430322741
+0.001076022562
+5.066109523e-06
+0.1068111238
+9.253312815e-05
+91.21436585
+0.01945125628
+8.958261209e-05
+1.678018826e-06
+0.0002943803381
+0.0001299537598
+1.170328977e-06
+0.0001328201344
+3034.923139
+1.378838481
+0.08550576557
+36.1778481
+1.785612005e-06
+9.054761711e-06
+0.001615788007
+1.23208212
+0.0003642605012
+2.708222463e-06
+4.719890307e-06
+0.001710230342
+0.3182780512
+0.0008830648263
+0.08778053115
+0.2629019547
+0.02748771871
+1.833687248e-06
+0.002690141449
+1.89343446e-06
+1.77338727e-06
+2.291261955e-05
+1.31623411e-05
+0.02865008699
+0.0001316279098
+5.992647136
+427898745.7
+0.0002908564713
+180.1583563
+1.149638817
+4.549308837e-06
+0.02938178265
+12.41064482
+0.0005867842208
+1.235864829e-06
+4.615116294e-05
+0.7228122449
+0.0003119308893
+1.100546805e-07
+0.0007569561799
+3.456980613e-06
+0.001902875586
+1.766220045
+0.0006636538355
+435.7550934
+1.468265667e-05
+0.003868037348
+0.005375424637
+9.895862027e-06
+0.002634610055
+0.04246482577
+1.811370936e-05
+0.0008533026641
+1.935841576
+1.965419689e-06
+0.0001940609469
+6.432993714e-05
+436.2325667
+8.585463564e-05
+0.003171908449
+3.1453383e-07
+0.10708771
+1.134414955e-07
+5.434578249e-06
+6.681002905e-06
+0.00246782063
+0.00431428184
+2.421882039e-05
+3.783451257e-08
+5.144794528e-06
+8.550628039e-05
+13660.98467
+117.1131418
+387.5973894
+0.2157832996
+2.335478847e-06
+0.0001022288321
+0.004686932429
+0.4278874185
+0.000106415896
+1.985873409e-06
+2.25291051e-07
+2.945623146e-05
+0.3581923626
+6.785775177e-06
+0.004893941209
+0.0005325322641
+1.060152245e-07
+1.118231806e-05
+2.203633348e-06
+3.969211813e-06
+0.007400958058
+5.769269726e-05
+1.899358848e-06
+0.002291707993
+0.007504512648
+0.0946266303
+0.07641638751
+0.001194744559
+0.000131740617
+3.126413693e-07
+6.676941465e-07
+2.247981499e-07
+1.394670123e-07
+0.0003117083998
+0.003124792631
+3.196614452e-07
+0.0005145950558
+1.196095567
+1.715539058e-05
+2.527890627e-05
+2.665407449e-06
+4.483788358e-07
+3.597134253e-05
+6.514469651e-07
+0.002675227421
+0.1700338633
+0.0151159864
+1.877463443e-05
+0.001018686854
+0.0001769414038
+3.068873861e-06
+1176.630703
+9.378378651e-06
+7.10431401e-07
+4.541139018e-05
+3.954117424e-06
+0.003271203855
+0.01288628186
+0.01758893926
+0.0001040458517
+0.08587993013
+7.493297002e-07
+0.0303353991
+2.525780749e-07
+1.075360445e-06
+0.0005909996117
+3.980399911e-06
+67.71892274
+8.802422373e-06
+3.997018327e-05
+2.690317449e-06
+0.1177580565
+1.362864165e-07
+0.00719011227
+3.379854611e-06
+1.335802906e-08
+0.0008096413047
+8.017608649e-08
+5.475646202e-07
+9.89836826e-07
+0.0001097524935
+0.005006547178
+0.01328382633
+3.708077203e-10
+1.7738466e-09
+6.195338978e-10
+1.375292134e-08
+7.712749689e-06
+1.087672399e-10
+7.612820481e-06
+1.631362221e-06
+7.195255533e-08
+8.307139358e-09
+2.261143452e-08
+2.267767047e-07
+7.680504754e-08
+1.175875499e-07
+0.0006018013011
+1.467361257e-05
+1.350450721e-07
+2.569869939e-07
+7.718323457e-07
+0.0170653268
+0.0001792376906
+3.072166197e-06
+22266.55115
+0.06729772649
+1.065337157
+7.103072447e-07
+0.0008701526857
+0.3457109502
+0.001639388444
+6.829709377e-07
+103.0220399
+1.724918219e-06
+0.006442680294
+0.0002924073612
+0.0061407479
+3.014208435e-07
+1.065501509e-05
+5.819752501e-07
+0.002933183567
+1.390851143e-08
+1.359009867e-06
+0.02819012119
+1.30006471e-09
+78.00480652
+3.27001628e-06
+2.84175006e-07
+3.993660531e-07
+0.0003605232041
+1.871558115e-06
+2.854501301e-07
+0.02600693272
+0.0004600409654
+0.001348216483
+2.488693144e-08
+2.395438283e-11
+4.248528123e-07
+1.389769895e-05
+8.049135616e-07
+2.12824677e-07
+8.405933433e-07
+2.36621812e-06
+0.06259159217
+1.831343487e-07
+2.1663617e-08
+6.323481085e-10
+1.339773875e-07
+4.883049792e-06
+1.705369167e-05
+0.0001013417062
+0.01536206656
+1.67093181
+5.039245627e-10
+1.141488157e-05
+3.188719492e-07
+3.24135522e-10
+0.0002434524475
+1.383868594e-06
+0.498790572
+9.03612161e-05
+0.02764416817
+0.01341805272
+7.529041731e-07
+0.0002257161473
+0.0003717671499
+1.858022341e-08
+5.715357136e-07
+2.477427879e-09
+0.0001517676839
+2.326424939e-11
+1.729565889e-05
+7.354412281e-09
+6.877321994e-06
+0.0005289104709
+0.02479967619
+6.212490017e-06
+1.714682348e-08
+7.916835703e-09
+0.00934088881
+4.083237353e-07
+0.004696390685
+6.761095529e-10
+0.003004586166
+6.030026635e-05
+8.803060787e-08
+6.957661086e-08
+9.612925872e-08
+1.628346787e-05
+5.143053469e-07
+8.784807581e-08
+1.024493405e-08
+3.137825282e-07
+42.00569276
+0.004640201019
+1.146513701e-05
+4.057040434e-08
+3.381428808e-06
+2.603443408e-07
+0.002272672712
+1.276808422e-05
+0.06226820245
+1.778585408e-08
+2.378758757e-07
+9.7714181e-05
+0.387678215
+4.398556239e-07
+3.0060972e-08
+3.874778496e-05
+9.119354143e-06
+1.402325829e-06
+4.695875643e-06
+2.383017224e-07
+2.417306793e-06
+0.0002708928034
+0.01384277314
+1.710909166e-07
+1.217157554
+177.0639344
+0.03861638077
+7.433640647e-05
+0.0002263704383
+0.002948183736
+0.008773163623
+5.353914571e-08
+360.0073872
+6.478760385e-07
+2.750359008e-07
+8.589493477e-08
+4.263303371e-07
+2.005861581e-05
+5.503740271e-08
+0.00180575553
+16.03887797
+7.922076143
+5.036869078e-12
+1.277093114e-06
+0.003373959327
+1.743675895e-05
+1.778007998e-06
+0.0002212911594
+9.259789313e-05
+8.539690283e-05
+2.423945599e-06
+0.0001233083611
+6.913916933e-09
+2.570708518e-07
+0.0008554774811
+0.2986161102
+0.0006441353055
+1.203219038e-07
+4.75874322e-07
+258.1003097
+2.768523645e-09
+1.270156521e-11
+1.46097374e-09
+34.87523283
+8.506779165e-05
+4.671813569e-06
+1.907838682e-05
+3.716521365e-07
+1.888947852e-08
+0.02669197841
+15.26083816
+2.873162527e-10
+1.360352114e-06
+9.527755241e-07
+0.0003084354836
+0.0009838761157
+1.302357534e-05
+0.004623681617
+0.01694286548
+9.642283586e-07
+8.499293031e-06
+8.770386824e-06
+0.001276821989
+6.668543971e-07
+5.652013261e-11
+0.238779506
+1.624068992e-07
+0.0002438443448
+0.002313588467
+7.53704251e-09
+7.248545824e-12
+0.02160343338
+6.465648118e-09
+7.430463219e-08
+0.000611761207
+2.262200983e-07
+8.810901655e-07
+1.688360234e-07
+2.91786064e-06
+5.480913194e-08
+0.0004601409972
+3.364979337e-10
+8.967694023e-09
+0.1345541103
+2.492529354e-06
+9.220762636e-06
+7.386130159e-06
+3.598519539
+8.045779179e-06
+4.730716431e-07
+2.192459188e-07
+2.194928052e-06
+3.130712514e-07
+1.43402212e-05
+0.0002608840076
+23.2295584
+1.440575153e-08
+4.161369951e-06
+0.001762466634
+4.255417524e-09
+3320585.541
+0.0009563172933
+0.03635196251
+0.5161165209
+0.004796521016
+2.31018301e-06
+0.06003700358
+0.0004342970655
+2.009638611e-08
+0.0001808310494
+0.0003251077547
+4.320349927e-07
+2.882854303e-08
+0.0422146142
+8.524423251e-07
+19.01825412
+0.000339795268
+1.575704657e-07
+0.0004837774983
+4.132343892e-07
+5.931203648e-09
+0.01135109493
+5.877226405e-05
+1.480445352e-08
+0.3379228367
+1.786471685e-06
+1.021098483e-07
+1.079160999e-08
+9.622764911e-09
+0.002978083967
+4.588139504e-06
+6.6721226e-11
+3.261090243e-06
+0.06059954782
+4.847476637e-09
+1.104977758e-06
+0.0003141612317
+0.0006548050265
+0.000225642662
+4.294099642e-06
+1.185934783e-07
+0.001042912926
+8.03811263e-06
+1.755333168e-10
+6.77273477e-08
+598.9038787
+1.572919365e-07
+0.04266693618
+2.739698989e-08
+0.0005135423693
+1.264426626e-10
+2.608778512e-07
+1.857091679e-06
+1.075074651e-10
+1.188966384e-07
+0.03379476153
+0.0006473640268
+3.624390296e-09
+1.237883305e-05
+1.016117525e-05
+9.687850834e-06
+0.0002738964397
+9.038889191e-09
+7.786947927e-07
+3.198680693e-05
+1.608278345e-11
+0.003353487144
+6.107229663e-08
+0.0003768704931
+0.0002005945404
+3.670401818e-11
+9.030155872e-09
+2.042723216e-08
+2.386955262e-08
+4.835084799e-09
+3.685817471e-11
+1.066443786e-07
+2.652667148e-07
+3.649486136e-11
+0.02477462094
+3.517378902e-05
+3.080027841e-09
+1.018709237e-09
+2.870228082e-10
+8.744977958e-06
+5.273214292e-09
+3.168385461e-11
+0.001556612382
+3.116051041e-07
+2.492191426e-09
+7.615392241e-10
+0.01004217426
+1.101368692e-06
+3.918013244e-06
+3.713827343e-09
+6.857767402e-09
+1.290966146e-09
+6.224822208e-12
+3.376319521e-06
+3.94050205e-09
+7.633357106e-05
+6.164135969e-05
+1.630673977e-06
+7.08423096e-09
+4.080955534e-10
+232770.6483
+0.008068531153
+8.666549381e-10
+2.246457597e-10
+3.222301462e-08
+1.395629818e-06
+2.196466632e-08
+1.592802071e-08
+4.326999161e-05
+0.0005042839931
+0.00197308444
+3.032370137e-05
+0.0005857579152
+7.171987925e-06
+6.864328965e-06
+4.00109422e-06
+0.005161817275
+2.375668536e-09
+9.778482285e-06
+1.181506681e-08
+1.503749258e-08
+3.203838424e-11
+0.3676119133
+1.281881299e-09
+0.0004415836065
+1.190288394e-07
+8.77428226e-07
+7.210177682e-11
+1.45886715e-06
+0.02121994342
+0.01168544154
+3.982811877e-09
+8.600471911e-09
+1.434258247e-08
+8.098070419e-10
+0.0007724068717
+2.490831815e-12
+2.21358682e-06
+0.2873305685
+5.652944325e-06
+2.722304477e-09
+3.122833335e-06
+1.198571403e-06
+4.489385847e-09
+2.67381748e-08
+3.471377557e-09
+5.756644016e-06
+6.839839041e-05
+2.209846587e-07
+8.291924237e-08
+5.242886271e-07
+6.844201737e-11
+1.726048304e-06
+9.225948144e-07
+0.000185717301
+0.002602315394
+6.138200606e-11
+0.0001818846094
+7.733889705e-14
+1.358383639e-08
+1.07000594e-06
+4.579302999e-07
+7.849169727e-07
+2.448333585
+7.535064872e-11
+3.779457663e-06
+1.495153092e-07
+3.954406596e-09
+8.122956595e-09
+0.0001228391098
+2.269361403e-11
+0.001844497418
+0.01115645921
+2.08005184e-06
+4.660852875e-11
+4.199474743e-10
+4.628821128e-05
+0.0004980806998
+7.598116953e-09
+4.143456828e-13
+0.002644658829
+2.564668699e-08
+6.50580785e-08
+3.820040563e-11
+0.2398554094
+9.400695396e-08
+0.002358652784
+1.171877918e-11
+1.637975849e-09
+8.942959658e-08
+3.408065038e-10
+6.371021504e-07
+8.862071381e-07
+8.07800753e-13
+1.275141633e-07
+5.065134208e-08
+3.723768129e-05
+1.586955855e-06
+4.630822765e-08
+3.869149299e-05
+0.0002373420357
+8.091381421e-08
+2.031825403e-10
+1.493296083e-05
+2.814058365e-06
+4.293942142e-10
+1.227720708e-09
+9.668701181e-10
+5.082776865e-06
+0.001170366525
+2.950213679e-07
+2.929772297e-06
+3.908673483e-06
+1.641984659e-08
+9.22544372e-06
+0.004601199051
+1.245696104e-11
+1.623356072e-06
+0.0004075079271
+0.06223654452
+2.229639676e-12
+1.548827536e-10
+0.2098995653
+0.02739866019
+1.922659779e-12
+0.001586041269
+1.845383825
+0.8609993903
+0.09603688052
+5.904130818e-11
+9.958673683e-12
+0.05389921105
+0.0001869127049
+6.618563858e-08
+3.845467031e-07
+9.065990802e-07
+0.0003895671306
+2.782863325e-11
+0.0002689153495
+5.924260979e-06
+0.0005747447785
+4.117458599e-09
+5.117012888e-08
+3.396356612e-09
+8.109913215e-12
+4.329382871e-08
+4.881990436e-07
+1.276280038e-11
+9.120675139e-12
+7.192434752e-07
+3.913876447e-11
+1.9517939e-08
+0.11670582
+1.114874302e-06
+8.08838007e-07
+2.596368516e-07
+0.0001277437573
+0.0004507490086
+5.934514686e-10
+9.045316154e-06
+1.1680623e-08
+4.405302098e-07
+1.103299007e-15
+7.105364025e-08
+6.778113877e-05
+2.587602347e-06
+2.630707294e-07
+2.256922483e-10
+1.585065727e-05
+9.078985718e-05
+1.439932327e-09
+0.003315113132
+4.417287152e-05
+3.208255435e-09
+7.222413113e-10
+0.000300699488
+4.624336745e-05
+8.037199117e-10
+6.718759329e-12
+3.83547078e-05
+7.94104255e-05
+4.496517789e-05
+1.351102043e-06
+9.177810983e-10
+8.693619098e-07
+1.929511565e-06
+4.773961763e-09
+1.896551164e-08
+1.552346596e-08
+9.245502051e-06
+0.0001037407773
+6.368845512e-07
+4.041904863e-06
+1.506851808e-08
+9.717844454e-12
+2.297581007e-05
+0.0234227485
+0.00306449787
+5.437659639e-08
+0.001404494601
+1.582456182e-12
+10.74537802
+0.0006876957837
+3.080472252e-07
+4.979458387e-09
+2.191738006e-08
diff --git a/t/ME_data/ME_h_mtinf.dat b/t/ME_data/ME_h_mtinf.dat
new file mode 100644
index 0000000..73445f3
--- /dev/null
+++ b/t/ME_data/ME_h_mtinf.dat
@@ -0,0 +1,807 @@
+5.582760808e-08
+2.825811483e-05
+3.895336201e-05
+0.02046456494
+4.688432033e-06
+1.529291749e-07
+6.420743602e-07
+5.905800101e-07
+6.479004963e-07
+1.48600349e-05
+4.880154255e-05
+1.945209343e-06
+0.0004438901791
+1.600557037e-05
+1.231846376e-05
+0.03428118184
+8.553476044e-07
+579.5994998
+6.327145731e-06
+0.01311111831
+2.253068503
+0.001161627534
+0.0007930178104
+3.019642381e-05
+0.0001604361161
+0.0009782409472
+0.03878378325
+3.213295457e-07
+1.942127118
+9.002578443e-08
+0.04615078581
+0.1667432966
+1.060974436e-06
+0.004405656545
+0.02198234308
+3.75309762e-05
+6.894134543e-07
+0.006337938566
+0.000232790789
+0.004628411362
+0.0005874171572
+0.0114321818
+0.002787500992
+2.287103816e-05
+0.0001818288916
+0.009239536849
+0.03250478346
+3.095188489e-07
+5.383217293e-05
+12.03590233
+0.0001869627067
+4.288822951e-06
+9.903075819e-05
+0.0001133438964
+0.1413834117
+0.02374806527
+5.439959182
+0.0002200036282
+1.257502016
+0.00554436291
+1.492794515e-06
+0.0410145385
+0.003110926104
+4.895681474
+9.566466236e-08
+4.29754523e-06
+1.26422416e-05
+1.50686362e-05
+4.031356768e-06
+0.0009303622768
+3.581771412e-06
+2.209149529
+0.001347543465
+0.001948045992
+25948.6239
+4.798288837e-07
+0.001267876802
+2.198489834e-07
+0.02568760469
+1.738239331e-06
+5.383639367e-05
+0.03615040778
+48368.71222
+6.221468888e-06
+8.272280776e-06
+2.436735404e-06
+0.00082069939
+0.0001956879791
+8.464294968e-08
+0.002499166028
+0.0005268877995
+0.08401690492
+0.001109681865
+1.676530038e-05
+0.1210659795
+0.2021753863
+0.0001468722356
+1.49843416e-05
+0.01797157389
+0.001470636482
+6.108534705e-06
+0.1147500209
+6.044324523e-05
+83.23969297
+0.01952725493
+0.0001089245711
+5.957305171e-07
+0.00020821416
+0.0001540635133
+1.358665617e-06
+0.0001977851353
+2969.178567
+1.434474792
+0.07897582528
+33.78910043
+1.944172541e-06
+8.511420265e-06
+0.001671034953
+1.106372417
+0.000475762259
+2.254611547e-06
+5.024971782e-06
+0.002685315998
+0.3204715262
+0.0006830156865
+0.06161172452
+0.3792093477
+0.03062030789
+1.284034864e-06
+0.002336243826
+1.185574165e-05
+1.710024271e-06
+2.253363854e-05
+1.300040157e-05
+0.02939378205
+0.0001573213338
+6.400482717
+484632185.3
+0.0002916173761
+171.9267844
+1.633949616
+3.749772544e-06
+0.03071686489
+13.06467703
+0.0006848589426
+2.117832499e-06
+2.785336442e-05
+1.197689446
+0.0002238650692
+1.488235642e-07
+0.0007409854778
+6.610729458e-06
+0.002123335104
+1.564470647
+0.0006669917367
+473.2723045
+1.432945457e-05
+0.003779071212
+0.005378672443
+4.446772576e-06
+0.002341093315
+0.1693157671
+1.832833621e-05
+0.0006895620169
+1.807274869
+1.739710164e-06
+0.0002259708712
+6.536256314e-05
+438.0279281
+8.402208741e-05
+0.003349479787
+1.893640387e-07
+0.1112826341
+2.103910568e-07
+6.061483157e-06
+6.743144434e-06
+0.002108199475
+0.003016819088
+2.679236437e-05
+4.152061832e-08
+6.910237845e-06
+4.414649231e-05
+13180.47273
+112.3954293
+407.9774509
+0.1947783364
+3.633741266e-06
+8.949320507e-05
+0.004538125269
+0.437437773
+9.781170001e-05
+2.778563584e-06
+1.70122334e-06
+3.067502446e-05
+0.3572469398
+4.865503412e-06
+0.005054297163
+0.0005309402179
+2.786914424e-08
+1.278066173e-05
+9.682183899e-06
+2.538469656e-06
+0.006448672409
+6.607293322e-05
+1.895984721e-06
+0.001842602862
+0.07109867357
+0.09292074641
+0.08069062114
+0.0007841374424
+0.0001302148142
+6.96049157e-07
+4.632599032e-07
+2.355250248e-07
+1.220815109e-06
+0.0003380714588
+0.003131420086
+6.328928995e-07
+0.000531365034
+1.039688858
+1.577333165e-05
+4.032814959e-05
+2.142451659e-06
+4.632583267e-07
+3.564353255e-05
+7.645988124e-07
+0.004094577112
+0.171048394
+0.009405424559
+1.086927134e-05
+0.0009935783254
+0.0001283118369
+3.018764055e-06
+1142.596042
+6.149007778e-06
+5.946922879e-06
+4.484486221e-05
+4.875515638e-06
+0.003370278763
+0.01348136435
+0.01694870912
+8.663774537e-05
+0.07732984975
+2.265338646e-06
+0.0303766268
+6.402972472e-07
+2.703996094e-06
+0.0006768238778
+2.4208244e-06
+64.78560774
+8.344522601e-06
+3.189381378e-05
+1.712603497e-06
+0.09667429817
+1.393889614e-07
+0.007276533549
+5.008302644e-06
+1.627964527e-08
+0.0008524764764
+1.204352922e-07
+7.15763651e-07
+9.48877371e-07
+0.0001193656831
+0.005851402185
+0.01372488252
+5.72350516e-10
+3.614330908e-08
+7.320467704e-10
+1.453936566e-08
+1.727671119e-05
+1.524992372e-10
+7.316531328e-06
+2.398381457e-06
+9.455892839e-08
+1.199317127e-08
+2.857466479e-08
+2.923880608e-07
+1.635996229e-07
+5.008649952e-07
+0.004371895197
+1.533701531e-05
+5.07640554e-07
+3.104784661e-07
+2.161998882e-06
+0.0205663762
+0.0001962376019
+9.949633031e-06
+25087.75723
+0.07129826748
+1.131381297
+8.534253593e-07
+0.0008626947124
+0.3319845459
+0.001688940793
+7.75104173e-07
+135.8463341
+1.718736108e-06
+0.01768528155
+0.0003050867091
+0.009055152615
+4.871654822e-07
+9.928551431e-07
+6.761126468e-07
+0.002764873145
+4.20752883e-07
+1.839465257e-06
+0.02911407063
+1.429095869e-09
+123.8464804
+3.878907185e-06
+2.823435977e-07
+2.181998951e-07
+0.0003628671133
+1.956611944e-06
+2.925524697e-07
+0.02740505501
+0.0004835484089
+0.001294968821
+4.961467604e-08
+9.379585371e-11
+4.794611466e-07
+1.724735523e-05
+1.053115346e-06
+1.004941951e-06
+1.456188771e-06
+6.012953164e-06
+0.07435516808
+4.986554318e-07
+2.620286035e-08
+6.730186632e-10
+1.462561843e-07
+7.843083295e-06
+2.23455571e-05
+0.0001134351386
+0.0154680137
+1.735552558
+1.156605839e-08
+1.294157613e-05
+3.348717433e-07
+3.491517956e-09
+0.0002711312012
+5.584761145e-06
+0.9006208236
+0.0001003115609
+0.02237513671
+0.01338382443
+7.481229115e-07
+0.0002234987644
+0.0003705759886
+1.898351129e-08
+1.657239164e-06
+4.093878871e-09
+0.0001758241056
+8.431436166e-11
+0.0005833802259
+7.70161748e-09
+9.983626918e-06
+0.0007035775066
+0.02474478055
+6.270525808e-06
+3.468720779e-08
+1.724400474e-08
+0.009117250451
+6.078887573e-07
+0.005082678025
+8.119315225e-10
+0.004497344072
+7.071952691e-05
+1.118931716e-07
+9.463920409e-08
+1.144548901e-07
+2.499918959e-05
+1.104695967e-06
+9.186735585e-08
+1.118167162e-08
+3.166442107e-07
+47.27735999
+0.005619787542
+2.645615442e-05
+3.217410189e-08
+3.940942183e-06
+1.708751497e-06
+0.00269242283
+1.355472425e-05
+0.06048498767
+4.185594237e-08
+2.759300723e-07
+0.0001286610383
+0.4621788119
+8.258814737e-07
+1.033330814e-06
+4.526252424e-05
+2.303720626e-05
+1.472556126e-06
+4.863423273e-06
+2.234988348e-07
+5.189969176e-06
+0.0002798831217
+0.01368736644
+1.619630911e-06
+10.78852095
+178.0821043
+0.04172457533
+7.340946131e-05
+0.0002317124782
+0.00328386642
+0.01860979559
+4.402941352e-07
+839.2570616
+7.234469729e-07
+7.40724986e-07
+9.050963032e-08
+6.950885896e-07
+1.966481291e-05
+7.845545986e-08
+0.002422229196
+16.84187681
+8.373737241
+2.510688577e-11
+1.354006565e-06
+0.003916208072
+1.87834382e-05
+1.061893299e-06
+0.0002522663944
+0.0002802397367
+8.284625139e-05
+2.74659168e-06
+0.0001365424915
+4.110023632e-09
+3.185742762e-07
+0.001078329246
+0.2900316086
+0.0007813567568
+1.247378241e-07
+6.657616321e-07
+262.7194422
+2.172891383e-07
+2.392873607e-10
+1.822331974e-09
+34.82068803
+0.0001058890218
+3.541148799e-06
+1.877530071e-05
+4.850000997e-07
+2.306721765e-07
+0.02931008827
+65.25317796
+4.062413402e-10
+1.440603973e-06
+1.066449458e-06
+0.0003088258184
+0.001075411908
+1.303005239e-05
+0.006471141599
+0.01899374932
+1.062724655e-06
+8.253512942e-06
+9.334440409e-06
+0.001280639873
+2.283680721e-06
+6.260116796e-10
+0.2534558629
+1.630677018e-07
+0.0003298919389
+0.002533846757
+1.914586499e-08
+8.606769006e-12
+0.03263233412
+7.203104901e-09
+8.088413546e-08
+0.0006269777645
+2.473222665e-07
+9.968830416e-07
+1.746698633e-07
+4.753145056e-06
+1.328682579e-07
+0.000508160052
+3.992578846e-10
+1.106276429e-08
+0.1695339751
+3.770100276e-06
+1.040479918e-05
+8.063177206e-06
+3.629259261
+2.105289275e-05
+5.04901351e-07
+2.339407131e-07
+2.502888392e-06
+4.209803322e-07
+4.838866183e-05
+0.0002779278814
+23.3462823
+2.131347013e-08
+8.739796471e-06
+0.00178497071
+4.510412368e-09
+3191738.416
+0.0009705156251
+0.04676370472
+0.5174086123
+0.005058629843
+2.355889917e-06
+0.07473378013
+0.000615424774
+3.193865064e-08
+0.0002793243959
+0.0003800264825
+8.196122567e-07
+2.938519398e-08
+0.01644793231
+8.533246724e-07
+18.46549494
+0.0004370897888
+1.651976144e-07
+0.0005031457335
+5.155226872e-07
+6.138259697e-09
+0.01196182175
+7.309559373e-05
+2.265688286e-08
+0.3556622364
+4.67413735e-06
+1.135831503e-07
+1.079470448e-08
+4.325481297e-09
+0.003030054728
+4.908949019e-06
+8.440216766e-10
+3.26659775e-06
+0.07242812553
+9.046937232e-09
+2.805134548e-06
+0.0005862513052
+0.0003870837421
+0.0002306455526
+7.65295072e-06
+1.202363418e-07
+0.001047970198
+8.175645768e-06
+4.081373723e-10
+7.487964052e-08
+815.4072543
+2.865432621e-07
+0.05412369048
+1.096693013e-07
+0.0007231232968
+1.365769452e-10
+2.902344226e-07
+5.476182992e-06
+1.048398369e-09
+1.216243516e-07
+0.03523573874
+0.0006888643475
+4.389532248e-09
+1.427610434e-05
+1.094889803e-05
+1.015424181e-05
+0.0004168879878
+3.936307461e-08
+1.864671086e-06
+3.356134604e-05
+1.945892381e-11
+0.003928393406
+6.103483135e-08
+0.0004013798122
+0.0002197116329
+4.597021604e-11
+2.2516721e-08
+4.243452504e-07
+2.704134795e-08
+7.131477445e-09
+4.669986887e-11
+1.285703296e-07
+3.257431397e-07
+2.285919944e-10
+0.02508641314
+5.454810101e-05
+3.296367932e-09
+1.11729952e-09
+3.490965824e-10
+1.193576046e-05
+9.507580167e-09
+3.859691143e-11
+0.00167134939
+1.356646543e-06
+3.553523289e-09
+1.130199408e-09
+0.01208123117
+2.979263331e-06
+7.174415532e-06
+9.589015627e-09
+1.479759217e-08
+1.36245597e-09
+6.548210064e-12
+3.763473932e-06
+3.525318889e-09
+7.719401841e-05
+6.539626343e-05
+1.904220779e-06
+9.62129602e-09
+4.841186966e-10
+246593.2851
+0.008332653867
+1.290851572e-09
+2.990238898e-10
+3.448903773e-08
+1.448602664e-06
+2.559216695e-08
+1.719816211e-08
+4.225852213e-05
+0.0006274710758
+0.00204896325
+4.007808047e-05
+0.0005957886887
+1.602607821e-05
+1.074899093e-05
+4.545871079e-06
+0.005330085164
+3.842500322e-09
+1.156201498e-05
+1.349677574e-08
+1.757902852e-08
+3.523827584e-09
+0.3610832767
+1.900099274e-09
+0.0004719709839
+1.5697833e-07
+1.027727033e-06
+1.610359355e-10
+2.078150083e-06
+0.02721374189
+0.01309332711
+7.75903148e-09
+1.009860589e-08
+5.335001462e-08
+8.906171393e-10
+0.0008327997757
+4.736917232e-12
+2.6621679e-06
+0.3172551012
+1.837448918e-05
+3.339404398e-09
+4.301594676e-06
+1.494195327e-06
+1.412793449e-08
+2.895282688e-08
+4.767199796e-09
+6.659974899e-06
+7.56413762e-05
+7.076010906e-07
+2.80725435e-07
+5.262757869e-07
+7.802537307e-11
+1.930864852e-06
+9.117230152e-07
+0.0006021904412
+0.003598279106
+7.615276157e-11
+0.0003443063619
+2.836510651e-13
+1.703152726e-08
+1.122068602e-06
+5.640183119e-07
+8.033879766e-07
+2.374947276
+9.963323166e-10
+4.406676642e-05
+2.663286666e-07
+4.078038377e-09
+1.351374911e-08
+0.0001348505694
+2.402113845e-11
+0.002220938033
+0.01191439056
+2.454623495e-06
+8.750664283e-11
+5.138183276e-10
+5.493614804e-05
+0.0005301365869
+1.681039143e-08
+5.096493457e-13
+0.003476563315
+2.623036447e-08
+6.687561662e-08
+3.979738218e-11
+0.2399691631
+1.020773719e-07
+0.002623003
+1.49405648e-11
+2.26728107e-09
+1.094107299e-07
+4.247893449e-10
+8.192485481e-07
+9.679265525e-07
+1.206214751e-12
+1.740428847e-07
+9.278859806e-08
+4.42554943e-05
+2.313138874e-06
+1.035020982e-07
+5.313570353e-05
+0.0005596193144
+1.27100336e-07
+4.233962927e-10
+1.499268987e-05
+4.136080511e-06
+1.071383645e-08
+1.290378644e-09
+1.980075516e-09
+7.221820165e-06
+0.001200948065
+5.019340409e-07
+2.898262838e-06
+3.853847068e-06
+1.921129238e-08
+9.138522545e-06
+0.005026331226
+1.470543708e-11
+1.715756843e-06
+0.00157952304
+0.06229588336
+2.214867548e-11
+1.836663943e-10
+0.2079901302
+0.02670973542
+3.179425989e-12
+0.001740026058
+1.27258903
+1.251228844
+0.1029350173
+1.305960327e-10
+9.734742078e-12
+0.05480019337
+0.000234840227
+8.6608831e-08
+4.554624791e-07
+9.838980984e-07
+0.0004551488404
+2.739000093e-11
+0.0002836906549
+1.399175767e-05
+0.0006552138517
+6.975369626e-09
+5.512408759e-08
+3.376107196e-09
+1.45257794e-11
+4.720383471e-08
+4.73791695e-07
+1.933240743e-11
+8.53870322e-11
+7.401153388e-07
+1.081969852e-10
+2.285185619e-08
+0.1168762274
+1.412131993e-06
+1.149236718e-06
+2.677084682e-07
+0.000125731329
+0.0006240968511
+6.43106978e-10
+1.0667532e-05
+1.15844648e-08
+4.861485542e-07
+6.876119387e-14
+7.277574791e-08
+8.299245484e-05
+3.094145725e-05
+6.029188243e-07
+2.433731632e-10
+1.746514749e-05
+7.013149931e-05
+1.480501008e-09
+0.00343239303
+4.858960081e-05
+3.700071189e-09
+8.203470258e-10
+0.0003139786485
+4.972464081e-05
+1.461897125e-09
+2.470863429e-11
+5.239353747e-05
+0.0001213943975
+5.447430408e-05
+1.644260989e-06
+1.044601766e-09
+9.323205197e-07
+2.237012525e-06
+7.007232174e-09
+5.112958465e-08
+2.184186312e-08
+1.055926617e-05
+0.0001340518598
+7.51634369e-07
+5.178962613e-06
+1.784863902e-08
+1.076455468e-10
+2.735429644e-05
+0.0255466258
+0.003360195454
+9.962948316e-08
+0.003891911814
+1.805164492e-12
+10.92753296
+0.000689613333
+5.300043147e-07
+7.535926242e-09
+3.286303302e-08
diff --git a/t/ME_data/ME_h_mtmb.dat b/t/ME_data/ME_h_mtmb.dat
new file mode 100644
index 0000000..463c576
--- /dev/null
+++ b/t/ME_data/ME_h_mtmb.dat
@@ -0,0 +1,807 @@
+4.043634525e-08
+3.947735809e-05
+3.978333671e-05
+0.01707962353
+4.133534858e-06
+4.12758231e-07
+6.527379321e-07
+7.807924441e-07
+6.213247187e-07
+2.158745738e-05
+4.929604661e-05
+1.648537835e-06
+0.0005294944592
+1.551221705e-05
+9.732109994e-06
+0.03508523471
+8.385697667e-07
+608.0814044
+6.0749939e-06
+0.01370931025
+1.96358809
+0.0005777778828
+0.000827073758
+3.396653466e-05
+0.0003663185905
+0.000712177849
+0.03945983226
+6.694529237e-07
+2.002990871
+9.107975513e-08
+0.04804489153
+0.1078671772
+9.503786206e-07
+0.00665553606
+0.02341176253
+3.684250185e-05
+6.861318278e-07
+0.005400516499
+7.647108683e-05
+0.005354950581
+0.0006200616798
+0.008448060377
+0.003448403835
+3.215759836e-05
+0.0001546133101
+0.005591549046
+0.03949683081
+2.350530066e-07
+5.66847762e-05
+8.602097332
+0.0001911650059
+5.061383546e-06
+5.119802912e-06
+0.0001348170624
+0.1914012121
+0.02428150541
+4.539435975
+0.0002614830785
+1.314349378
+0.005127584311
+2.106730988e-06
+0.06042114066
+0.002877526755
+6.46292188
+7.320533172e-08
+4.686272323e-06
+1.674918088e-05
+1.561186199e-05
+1.127468541e-06
+0.0009956223126
+9.322344484e-07
+2.355779109
+0.001354263217
+0.001607437479
+27300.17114
+8.763842677e-07
+0.0006972261753
+1.071457504e-07
+0.02005686587
+1.961498186e-06
+4.460773766e-05
+0.02601668696
+50724.57307
+5.480882962e-06
+1.076584976e-06
+2.869362531e-06
+0.0001072518549
+0.0002293795185
+8.022114184e-08
+0.002698012227
+0.0005830909647
+0.07725070316
+0.001325583365
+3.695641615e-05
+0.1179283377
+0.187466976
+0.0001155356229
+1.778799047e-05
+0.02462347248
+0.001087874101
+5.154078743e-06
+0.1031621088
+9.484855944e-05
+93.85073623
+0.02033391802
+9.092704303e-05
+1.705551286e-06
+0.0002974887579
+0.0001322670062
+1.193514887e-06
+0.0001344144727
+3095.568362
+1.406347883
+0.08812981858
+35.91258076
+1.82576512e-06
+9.27160445e-06
+0.001666287395
+1.269161514
+0.0003682229129
+2.790042514e-06
+4.841321328e-06
+0.001729869633
+0.3293946648
+0.0008921922127
+0.08111101834
+0.2664463909
+0.02792004132
+1.869728836e-06
+0.002768079641
+1.908915912e-06
+1.7932868e-06
+2.342012655e-05
+1.3539045e-05
+0.02918304927
+0.0001341037055
+6.099926061
+438766540.3
+0.0003006556126
+182.0743425
+1.166556788
+4.685195193e-06
+0.02991092642
+12.79082303
+0.0005990054968
+1.245221356e-06
+4.607261468e-05
+0.729519575
+0.0003077751881
+1.113634813e-07
+0.0007837766949
+3.49282179e-06
+0.001938345149
+1.741954359
+0.000679566605
+447.8606334
+1.502685811e-05
+0.003972927121
+0.005572294955
+1.00407234e-05
+0.002697149585
+0.04268640856
+1.838639564e-05
+0.0008683709605
+1.984375256
+2.023641764e-06
+0.000198987372
+6.590604109e-05
+451.0223058
+8.870879716e-05
+0.003243511862
+3.187920803e-07
+0.1107987679
+1.145759171e-07
+5.532565231e-06
+6.903234805e-06
+0.002276208464
+0.004410664221
+2.467146722e-05
+3.863236208e-08
+5.225855864e-06
+8.615548536e-05
+13702.18793
+118.8790931
+398.6816605
+0.2170834277
+2.355356756e-06
+0.0001050726086
+0.004854583584
+0.4416609149
+0.0001091431342
+2.015137689e-06
+2.270813172e-07
+3.065037179e-05
+0.3695610208
+6.626431375e-06
+0.005034520677
+0.0005496323042
+1.077013565e-07
+1.134883954e-05
+2.217217266e-06
+4.053576036e-06
+0.007630107217
+5.870346284e-05
+1.937770986e-06
+0.002308560802
+0.00756714829
+0.09713373161
+0.07865638607
+0.001166988788
+0.0001360267898
+3.146930172e-07
+6.721239549e-07
+2.308081382e-07
+1.402445764e-07
+0.000317796034
+0.003250721897
+3.220983676e-07
+0.0005361385069
+1.1674911
+1.701754998e-05
+2.555659626e-05
+2.72704331e-06
+4.575225674e-07
+3.723503521e-05
+6.625793946e-07
+0.002713331741
+0.1757705927
+0.01528287892
+1.914346383e-05
+0.001027883901
+0.0001819228709
+3.169381656e-06
+1197.183436
+9.459688271e-06
+7.151860571e-07
+4.709453403e-05
+4.024117529e-06
+0.003399032872
+0.01326938543
+0.0181234889
+0.0001070400499
+0.08816817298
+7.557106621e-07
+0.03150507468
+2.546159018e-07
+1.083218854e-06
+0.0006029839411
+4.109536026e-06
+69.62539127
+9.053551382e-06
+3.976334138e-05
+2.74610264e-06
+0.1217200425
+1.39072939e-07
+0.007475099982
+3.40898863e-06
+1.355970335e-08
+0.0008211141398
+8.085185969e-08
+5.540708628e-07
+9.793920273e-07
+0.0001116678867
+0.005087064408
+0.01371958748
+3.747341913e-10
+1.784590611e-09
+6.285401932e-10
+1.403711889e-08
+7.942905268e-06
+1.099705834e-10
+7.694717507e-06
+1.645454731e-06
+7.281682073e-08
+8.437713425e-09
+2.293298229e-08
+2.337785459e-07
+7.748072516e-08
+1.184163566e-07
+0.0006048814428
+1.496340286e-05
+1.358047852e-07
+2.641795688e-07
+7.894816929e-07
+0.01733476065
+0.0001847336643
+3.09120607e-06
+22646.44954
+0.06825326182
+1.097908477
+7.194716543e-07
+0.0008831460117
+0.3510995791
+0.001703720233
+6.923341702e-07
+104.5021475
+1.817413974e-06
+0.006274458892
+0.000298133205
+0.006212019474
+2.971106209e-07
+1.073536231e-05
+6.001267547e-07
+0.002926676713
+1.394485679e-08
+1.378850246e-06
+0.02873430354
+1.324231029e-09
+78.68865757
+3.323528263e-06
+2.901477846e-07
+4.025907172e-07
+0.0003762267637
+1.888509349e-06
+2.984779601e-07
+0.02649146415
+0.0004735978501
+0.001390723358
+2.511800946e-08
+2.406174629e-11
+4.343797764e-07
+1.413682848e-05
+8.187236684e-07
+2.137756399e-07
+8.487602643e-07
+2.382894999e-06
+0.06385257472
+1.843038529e-07
+2.200340736e-08
+6.500296312e-10
+1.369543749e-07
+5.034654856e-06
+1.729239287e-05
+0.0001029880478
+0.01584013323
+1.717937137
+5.051793703e-10
+1.169719308e-05
+3.274892339e-07
+3.257575264e-10
+0.0002494213231
+1.395452311e-06
+0.503263349
+9.183628731e-05
+0.02652180806
+0.01351083505
+7.697372696e-07
+0.0002345916253
+0.0003785854477
+1.911331681e-08
+5.741981368e-07
+2.496303102e-09
+0.0001538309849
+2.342705363e-11
+1.733704579e-05
+7.440708121e-09
+6.97019919e-06
+0.0005366276141
+0.02581853732
+6.351131637e-06
+1.725302693e-08
+7.963798782e-09
+0.009484858796
+4.13192111e-07
+0.004760868685
+6.868174526e-10
+0.003035550995
+6.122424548e-05
+8.912863951e-08
+7.049197773e-08
+9.784536423e-08
+1.644442743e-05
+5.17689779e-07
+8.945719428e-08
+1.043648257e-08
+3.220171418e-07
+42.8640364
+0.004721640337
+1.17318728e-05
+4.163036826e-08
+3.442022366e-06
+2.613636766e-07
+0.002337930201
+1.311028403e-05
+0.06454239294
+1.788371853e-08
+2.431129133e-07
+9.887586539e-05
+0.393826783
+4.435492073e-07
+3.017594495e-08
+3.944030937e-05
+8.95155309e-06
+1.480534305e-06
+4.80208438e-06
+2.278546025e-07
+2.431285271e-06
+0.0002795248679
+0.01414217835
+1.721740596e-07
+1.225137873
+184.4942235
+0.03969271858
+7.5889168e-05
+0.0002307099487
+0.002991500607
+0.008823843341
+5.377454625e-08
+363.5972586
+6.577363411e-07
+2.768216121e-07
+8.808374953e-08
+4.320144869e-07
+2.087813482e-05
+5.563234756e-08
+0.001828101352
+16.31099206
+8.047143214
+5.054980175e-12
+1.31165661e-06
+0.003453525346
+1.775158542e-05
+1.735750599e-06
+0.0002258034827
+9.302350911e-05
+8.722704625e-05
+2.463340497e-06
+0.0001266405389
+7.186989683e-09
+2.616936561e-07
+0.0008671500774
+0.3034312162
+0.0006552874753
+1.250994614e-07
+4.836707652e-07
+267.0859855
+2.784590138e-09
+1.277933942e-11
+1.625022487e-09
+35.42273613
+8.653074513e-05
+4.787553737e-06
+1.988624796e-05
+3.776413822e-07
+1.896902599e-08
+0.02757075064
+15.7257863
+2.903221186e-10
+1.383656135e-06
+9.796122701e-07
+0.000313466982
+0.001008582166
+1.326553759e-05
+0.004678660451
+0.01717544418
+9.763673978e-07
+8.769786833e-06
+8.887711201e-06
+0.001304503338
+6.699929944e-07
+5.668340919e-11
+0.2432408562
+1.65579253e-07
+0.000246842309
+0.002380529019
+7.590644517e-09
+7.382339234e-12
+0.02191318318
+6.595093824e-09
+7.59955463e-08
+0.0006241836255
+2.306006909e-07
+8.949234518e-07
+1.745865901e-07
+2.984765317e-06
+5.515085903e-08
+0.0004672383281
+3.41289611e-10
+9.072028681e-09
+0.1360437848
+2.516296619e-06
+9.347242718e-06
+7.560996071e-06
+3.650323724
+8.084973509e-06
+4.821049366e-07
+2.262423437e-07
+2.237305224e-06
+3.165244113e-07
+1.440989373e-05
+0.0002661269589
+24.06039199
+1.471908477e-08
+4.189979743e-06
+0.001821016076
+4.36609224e-09
+3372180.49
+0.0009765702499
+0.03671795044
+0.5266849257
+0.004945279364
+2.349850855e-06
+0.06096774779
+0.000440582517
+2.030747703e-08
+0.0001828061992
+0.0003290502891
+4.360219646e-07
+2.945465959e-08
+0.04314824277
+8.695492983e-07
+19.3255897
+0.0003433917438
+1.623604604e-07
+0.0004985403516
+4.206613628e-07
+6.048976857e-09
+0.01167821222
+5.984173029e-05
+1.495919742e-08
+0.349858651
+1.798434603e-06
+1.039820674e-07
+1.118772128e-08
+9.889699357e-09
+0.003069395217
+4.6681898e-06
+6.697088362e-11
+3.350850437e-06
+0.06185962328
+4.883448988e-09
+1.110833085e-06
+0.0003200423438
+0.0006747162613
+0.0002346905421
+4.330939823e-06
+1.219218909e-07
+0.00107853716
+8.281691391e-06
+1.767444691e-10
+6.900019008e-08
+605.6626124
+1.58440332e-07
+0.04329139184
+2.752138714e-08
+0.0005185712214
+1.286238525e-10
+2.639855477e-07
+1.869433683e-06
+1.079551089e-10
+1.229225358e-07
+0.03468655392
+0.0006584274735
+3.689112319e-09
+1.260350348e-05
+1.040877611e-05
+9.878538403e-06
+0.00027695052
+9.100037836e-09
+8.033944341e-07
+3.315013452e-05
+1.633327055e-11
+0.003414227642
+6.245172231e-08
+0.0003805246583
+0.0002060494734
+3.729836801e-11
+9.13272934e-09
+2.051216421e-08
+2.436679394e-08
+4.8855578e-09
+3.735846979e-11
+1.088675628e-07
+2.683031236e-07
+3.706588258e-11
+0.02523141233
+3.549809347e-05
+3.149827076e-09
+1.040176277e-09
+2.910870047e-10
+8.86122235e-06
+5.308664537e-09
+3.218598972e-11
+0.001577139281
+3.19884756e-07
+2.55834443e-09
+7.691641781e-10
+0.01021181751
+1.105572103e-06
+3.947538433e-06
+3.734235002e-09
+6.899637353e-09
+1.304337436e-09
+6.538762872e-12
+3.449997548e-06
+4.050730086e-09
+7.962318534e-05
+6.295652823e-05
+1.652753715e-06
+7.186541052e-09
+4.161402504e-10
+240986.4156
+0.00837085449
+8.760965452e-10
+2.273546864e-10
+3.315418794e-08
+1.416471539e-06
+2.236464809e-08
+1.633966009e-08
+4.518163364e-05
+0.0005095243003
+0.001992318371
+3.067811242e-05
+0.000600451339
+7.236523643e-06
+6.932191963e-06
+4.091606026e-06
+0.00532740987
+2.395091742e-09
+9.944229055e-06
+1.205282983e-08
+1.590006487e-08
+3.215700583e-11
+0.3737229127
+1.296635624e-09
+0.0004534080133
+1.217859786e-07
+8.898892246e-07
+7.255073266e-11
+1.472337758e-06
+0.02150215784
+0.01192727543
+4.01006022e-09
+8.731812041e-09
+1.44194837e-08
+8.229024075e-10
+0.0007855372767
+2.506773438e-12
+2.241537936e-06
+0.2916665128
+5.828120261e-06
+2.752322662e-09
+3.158084387e-06
+1.215204067e-06
+4.56919432e-09
+2.790752035e-08
+3.511291959e-09
+5.863328277e-06
+6.94327717e-05
+2.22257262e-07
+8.350918192e-08
+5.331043956e-07
+6.969438883e-11
+1.760059719e-06
+9.217915003e-07
+0.000186591209
+0.002628593815
+6.211758257e-11
+0.0001834032575
+7.768284292e-14
+1.377546165e-08
+1.091734431e-06
+4.677996963e-07
+7.943499964e-07
+2.545817678
+7.562150253e-11
+3.799329956e-06
+1.510534633e-07
+3.983823634e-09
+8.184966815e-09
+0.000124713148
+2.309165194e-11
+0.001872087514
+0.01153490791
+2.14659923e-06
+4.693787525e-11
+4.26518767e-10
+4.684595348e-05
+0.0005076579094
+7.644830233e-09
+4.193877203e-13
+0.002676923954
+2.654050699e-08
+6.800195443e-08
+3.888064789e-11
+0.2446884347
+9.628489292e-08
+0.002537052711
+1.189694289e-11
+1.656137485e-09
+9.040798991e-08
+3.488787289e-10
+6.557255554e-07
+9.083676565e-07
+8.150458329e-13
+1.315321346e-07
+5.098465051e-08
+3.786448639e-05
+1.600711068e-06
+4.658140964e-08
+3.906381479e-05
+0.0002389891749
+8.287668941e-08
+2.051670765e-10
+1.60134009e-05
+2.847282071e-06
+4.308493985e-10
+1.238922307e-09
+9.741854048e-10
+5.146899686e-06
+0.001209019717
+3.025450861e-07
+2.974261102e-06
+3.899610282e-06
+1.695222355e-08
+9.332118489e-06
+0.004694117038
+1.266652132e-11
+1.65446818e-06
+0.0004100868952
+0.06302457519
+2.236229887e-12
+1.584413907e-10
+0.2179964362
+0.0278380659
+1.936978481e-12
+0.001609711018
+1.892891439
+0.8712003892
+0.09866862844
+5.934210083e-11
+1.002311696e-11
+0.05494208606
+0.000190347461
+6.722151007e-08
+4.030694324e-07
+9.278979899e-07
+0.0003978128465
+2.803383506e-11
+0.0002738201279
+5.960007001e-06
+0.0005814501443
+4.157124101e-09
+5.198894043e-08
+3.465884394e-09
+8.167365756e-12
+4.450019241e-08
+4.893968894e-07
+1.288321375e-11
+9.163290302e-12
+7.451831479e-07
+3.940383302e-11
+1.996508632e-08
+0.1236118645
+1.182043812e-06
+8.323870133e-07
+2.618934064e-07
+0.0001307650885
+0.0004571929478
+6.07226713e-10
+9.141057675e-06
+1.18356058e-08
+4.465375632e-07
+1.107633947e-15
+7.416788486e-08
+6.870175273e-05
+2.600618361e-06
+2.647425108e-07
+2.279494937e-10
+1.611105573e-05
+9.315734048e-05
+1.474187502e-09
+0.003419123567
+4.527389745e-05
+3.271142207e-09
+7.434473088e-10
+0.0003102710133
+4.710559791e-05
+8.139484815e-10
+6.749221433e-12
+3.882836934e-05
+8.016924727e-05
+4.578651435e-05
+1.374899461e-06
+9.31258077e-10
+9.107579154e-07
+2.005597861e-06
+4.853545363e-09
+1.906194275e-08
+1.569062034e-08
+9.37420375e-06
+0.0001051157618
+6.464238034e-07
+4.092696924e-06
+1.538436694e-08
+9.760057524e-12
+2.626691189e-05
+0.02379715614
+0.00314703493
+5.484543738e-08
+0.00141690332
+1.612112513e-12
+11.09088062
+0.0007046654234
+3.111374555e-07
+5.038033967e-09
+2.21156295e-08
diff --git a/t/ME_data/ME_jets.dat b/t/ME_data/ME_jets.dat
new file mode 100644
index 0000000..fc2be65
--- /dev/null
+++ b/t/ME_data/ME_jets.dat
@@ -0,0 +1,672 @@
+1089423.187
+16.03009374
+6634.608547
+75414606.57
+77554.50738
+197575.0454
+76983.13177
+38.35809694
+14.34964854
+2.299306716e+10
+769524899.8
+20.89671172
+76.68018946
+8308.664153
+2450736.574
+52.24359459
+497.7665557
+5051.865522
+94571.50832
+1362099715
+3690.443315
+4732694.851
+33060865.81
+38851015.2
+30289.18618
+283948.732
+14123.1659
+40.56212546
+6.779515122e+12
+25724.24084
+4.464903667
+26.15970999
+28358522.83
+481820.3073
+33.14640689
+91.72649968
+322690.8843
+39.04279288
+65541087.6
+2256.832837
+51.04918081
+21658.09425
+675.8508089
+0.003214397574
+413076.4025
+1194042.703
+15261824.18
+1.190018997e+13
+25852091.94
+765101.1953
+579.2342295
+23845182.04
+92.16369985
+705.5684735
+828764.9832
+81.971147
+1436684.808
+70.01623522
+230.5658106
+8794.725516
+890.1697598
+97617236.05
+31101.11974
+10351668.97
+187539232.1
+93875.71518
+43031260.49
+4.067125353
+6035.796288
+16569782.71
+46.2510317
+52.88862266
+0.5556956436
+606771384.7
+293585.0788
+8415963.902
+29.5069149
+18.96420915
+6.321441102e+14
+3.592998225
+129.9175932
+19.85723932
+20438567.28
+2.633206648e+10
+2051545.14
+13.98077027
+342.1399233
+17.65245465
+79250.23744
+1030.69768
+73407.42532
+237883.1282
+79.44231257
+81668.43497
+93415340.81
+430.286772
+2209.749603
+568258.2296
+173145.2499
+1620285.613
+337.6920671
+153.562771
+416.0616338
+2187.267876
+22.41576042
+49.58992905
+559801.7855
+45.11194021
+584508.9206
+35.3041943
+202.7925831
+15.00277101
+46.15144698
+4.588470323e+13
+597.3210575
+39.62999028
+5228016720
+822065.5972
+1.826581842e+10
+6.479966482
+486208349.5
+62.22613221
+21.46207828
+792.2244925
+98.46049107
+192.195726
+539903191.6
+7888.570415
+559946441.6
+99.19582481
+967.8349215
+1186287.112
+296.9970515
+393.8848181
+46.84865591
+27.79939609
+70102.51858
+68.89626885
+2702.415147
+257.708207
+155288432.1
+1225359.249
+90576.73064
+9434.676956
+88.90057162
+4378317.363
+15680420.87
+6567.795332
+94.80131393
+44.21436446
+61.79707787
+44.41095665
+7385.16094
+227.1571585
+94986.49609
+67239897.71
+96.37463492
+2911.906879
+49.889402
+6585.792276
+336579403.2
+564.3980943
+1554303.479
+33573.8975
+2.437369584e+10
+226746.4443
+75.83290452
+27.75735686
+94778193.37
+59.97102623
+64.15185493
+7.774417686e+13
+7438.779369
+22.08782146
+434613.7211
+34293520.02
+68775.95004
+307606184.4
+150615.387
+899713.4452
+9.223973371e+11
+59817.02462
+348.4932635
+2609243018
+24597079.8
+98.70527236
+83.33745442
+12560815.61
+1358729.011
+6976819.18
+7425.681281
+488.6665953
+13.67398052
+401.1059245
+1888.653087
+4009.955386
+175.1044182
+864180.0454
+37.7515354
+47187493.74
+109.4430381
+359.7148155
+2.340917207e+11
+179.4280073
+9827.615509
+51.43706602
+159.8742767
+2066998.073
+21197190.52
+154.5419162
+20.22828093
+96.58078786
+54.45061659
+65.30848921
+27.46863884
+399638.3233
+48.71138495
+1579072.132
+67717150.9
+1034311422
+2300.166812
+542.347043
+225.023696
+15.39245502
+5.045806172
+408439.9118
+2850069.206
+412.6441358
+461.2402375
+22.40894016
+37199025.96
+16799601.72
+1283263.103
+19.12993911
+4119931.525
+138990.0255
+122.7744597
+18733.45006
+56.21459033
+343399063.9
+40.69967214
+1257.945403
+2.549605794e+11
+100938.8483
+197079455.9
+18898.20284
+65974.22935
+2007851.613
+6824.687033
+4076901842
+7253481.757
+3.002058086
+1662154.497
+3.070682383e+12
+163.8653423
+560506230.7
+1510.132538
+26693.22493
+0.1285518109
+494554.148
+24594487.35
+271745.3087
+155.6774454
+0.01468054527
+3.34264067
+1588023.004
+11.02440854
+2.276304783
+2060679.506
+1.060994634
+7899476.433
+11392931.55
+0.002283073032
+223.8249674
+64.14820176
+222013139.2
+17285137.17
+4904.948169
+4.325556105e+10
+2.025060646
+64.72904144
+0.7119503103
+0.3108282312
+127.8083553
+231480500.9
+3004499.276
+35310522.85
+9501.295688
+34.48365857
+0.1237148274
+137.304491
+22906.44193
+2845.690168
+1425.623527
+8.059398488
+2836.653576
+545.9404471
+6676238.836
+453869217.1
+3046455.67
+3824.741211
+1733.364341
+163997517.5
+46856740.63
+18247.63995
+55504734.26
+1281.154525
+0.01317003733
+26796232.71
+614.7809373
+124.9917797
+13.32414443
+2.061438665
+1437324478
+7.653466138
+215267780.4
+304925.7977
+12340855.35
+1.292474046
+4012100939
+918.1699992
+0.001829689464
+266.1514907
+19586.50145
+9843.077092
+34563.10338
+2.665695419e+14
+194168.4335
+182.9208313
+1.789617285
+115011741.8
+1.600475143
+0.9039184064
+54.16404691
+3836249.251
+0.1453799136
+83576.30659
+75918.72599
+255875.4377
+0.0940918961
+8139.618634
+497813992.4
+9158737.332
+81.5174511
+85205.83347
+1827.041926
+5434428.481
+34601.43212
+645859.4651
+8.409921219
+102240.6791
+0.3455555689
+2.273623498e+11
+2139040773
+82682.82602
+164684299.3
+6840.694229
+6292.264217
+1.266028837
+23696.0817
+9776772.75
+146.8084842
+40186.198
+578656.2198
+0.05760079494
+1308.62247
+2846240637
+2269.666603
+0.05728404415
+97.64433608
+816.0835471
+0.003537766648
+715128749.1
+581519.0926
+28878.33931
+606.7268882
+8.540999968
+1.843057112
+14728.95688
+0.09287199376
+4.952792168
+39703045.24
+0.2253239938
+0.4007086942
+0.6679356978
+8985.245646
+1459.149121
+222.574255
+1410.034894
+49.66459517
+280057.7217
+11990414.03
+22.39502507
+6652.821533
+412.6309921
+3976669.212
+376.5293878
+0.04993631096
+325.2221445
+0.001056663136
+0.6094009495
+0.01755460768
+4169089.697
+120922.3418
+19.94386406
+1.034905758
+899064.7273
+13370.52361
+12.9511581
+7.43997218
+0.2195708041
+22.69922572
+1364456.758
+37.31335252
+193721994.5
+5387171012
+83382.68593
+4356.219769
+4.959390234
+280348.136
+4068.591503
+497.1128967
+1220.499616
+37070808.46
+5.029130002
+13411.495
+21606019.92
+13.50033011
+10356.03697
+0.853804882
+0.004144159842
+66.5696697
+215.1472766
+1.80856617
+47652781.08
+187.4720271
+0.02309225892
+221825.3003
+68502271.85
+89524.33999
+63.06068401
+35193.88681
+0.1324344626
+0.07415048828
+1.316536312
+1914484.13
+239.8833788
+404.8578233
+311149.7413
+31.91203802
+6765.191464
+3683.120991
+111.8561994
+622181.6379
+3.495906967
+1.324392325
+2932853.167
+1221.022853
+0.3964612919
+772.5117573
+1587.443569
+16407.82538
+1153370.819
+697976.9971
+5293.102236
+0.6273574623
+433333.6561
+667.0010151
+279.5253754
+13.71499773
+27581.98765
+0.1265026607
+0.1265655615
+0.217545965
+119235.7842
+1207.825866
+1.168584528
+442927.4532
+11025.74448
+12920.20556
+6.440028505
+0.01867492246
+0.009209754218
+3.078810838
+1.895501574
+4354.355041
+1376.801012
+281.0956164
+12.62981914
+0.01896307384
+1.469952339
+1070.389041
+0.3140297342
+17576.52171
+19.46940005
+72639.28176
+1.00064982
+10.58896856
+146.5709549
+240963.5215
+21441702.88
+9436.867787
+1807651761
+253.6931177
+33906.26088
+73290.29485
+58.88693751
+0.1874080618
+128856.1698
+37.9170344
+999.730091
+79.29102322
+4297075.439
+2815203002
+1.065757865
+41.63356752
+3132937.645
+0.006570811459
+114610718.1
+19766422.36
+104.7817608
+978.6047854
+1114.98705
+2962094.183
+4303175.509
+40918.11534
+6.458644406
+0.01513394186
+10888.55623
+1429.371661
+72470.14267
+1239031.521
+0.01533517898
+9486.889743
+9312.872953
+0.001369262882
+263.7126988
+27397.40987
+6.460887231e+11
+20.47552795
+11886.05681
+359511.4015
+16322071.14
+5.495229
+94.06819575
+11159.01276
+497280.7644
+152.8232363
+2819.864026
+1840.922093
+105.3416519
+1.51196546
+893.0699171
+274.2733387
+213.9840929
+29.84435445
+0.06575134529
+68473.33981
+0.04814411438
+95307.18205
+0.2103887883
+87533.03201
+170.1269759
+12266.74779
+4.906472106
+1289.949588
+206.2206679
+1339736.494
+3265180.524
+817632.5321
+777148.2976
+45580.38902
+112.3471097
+40.46375548
+113.0563114
+0.5022425417
+125.6469752
+1250753.214
+119.6687695
+7890.066678
+255592.6823
+0.09544112407
+80.43028795
+168.8638851
+231.3586602
+9.003352253
+3865310.035
+0.0327080534
+0.2259241849
+0.3309640842
+70628514.93
+1.342293964
+1558415395
+1.535156072
+1105643165
+11.12340492
+19553.68065
+0.03262635444
+43.60758633
+149.96132
+16.71798701
+862303.2755
+3779.707836
+227.1354158
+0.04181934127
+0.007441956106
+0.06547652512
+562.1841603
+0.1962099745
+47427.77928
+278196156.8
+3272270.994
+3.274695664
+3.671659367
+0.8773587215
+3475.635402
+150.1691185
+79249376.57
+100282855
+1034000.854
+34631.08123
+0.2021471708
+11338.25051
+365.9254869
+3614.640649
+38892.33006
+3206.818315
+495.3831508
+3131.024396
+1216.930767
+3.68554273
+0.0007587392874
+27267.82366
+66215.88611
+12270.6996
+2181.519341
+9.799873017
+84.23360373
+350624.9294
+2324.233559
+0.210264797
+1438761.868
+160.6723109
+39.13994322
+2317.668587
+32196.71663
+0.02293846799
+177.7024259
+294490.6299
+963.9483044
+683.7845476
+12.22284439
+4350.875864
+132.6662523
+800.0540763
+0.6752102469
+6384.094581
+667.7145567
+88.5979278
+65957.42139
+0.1009745291
+7074.842515
+512.2889019
+11.48214352
+2.877973951
+1.31777398
+7.553015302
+0.01771209643
+0.01420982393
+5.350111784
+5171856.193
+0.6877419322
diff --git a/t/ME_data/ME_mt.dat b/t/ME_data/ME_mt.dat
deleted file mode 100644
index 899ceb0..0000000
--- a/t/ME_data/ME_mt.dat
+++ /dev/null
@@ -1,832 +0,0 @@
-3.999227056e-08
-4.193300057e-05
-3.89321526e-05
-0.0168082963
-4.057962883e-06
-4.063245819e-07
-6.411897603e-07
-7.822040516e-07
-6.064514954e-07
-2.202805699e-05
-4.738208766e-05
-1.624045121e-06
-0.0005217822994
-1.536443326e-05
-9.5908372e-06
-0.0339113482
-8.152325116e-07
-608.615459
-5.962078524e-06
-0.01340321821
-1.931630897
-0.0005720381031
-0.0008001259943
-3.300880891e-05
-0.0003564506816
-0.0007037298905
-0.03792377651
-6.513991407e-07
-2.011798913
-8.842231663e-08
-0.0463578742
-0.1064067277
-9.357584012e-07
-0.006517735511
-0.02272391876
-3.606711294e-05
-6.725111513e-07
-0.005897687516
-7.581255699e-05
-0.005354681787
-0.0006088947447
-0.0083356406
-0.003352568463
-3.23438533e-05
-0.0001515363281
-0.005531947882
-0.03875973775
-2.312735565e-07
-5.479259483e-05
-8.480458133
-0.0001854378187
-4.957725793e-06
-5.091269317e-06
-0.0001308303481
-0.1918392997
-0.0238989834
-4.446354691
-0.0002597751966
-1.267033388
-0.005050109129
-2.130648865e-06
-0.05854510061
-0.00280707581
-6.289945688
-7.216599144e-08
-4.779195357e-06
-1.655604389e-05
-1.522656761e-05
-1.116772322e-06
-0.0009539322219
-9.279967349e-07
-2.28564472
-0.001331277201
-0.001585734932
-26972.25526
-8.555976389e-07
-0.0006894623725
-1.065245235e-07
-0.01968534835
-1.951627129e-06
-4.401131811e-05
-0.02561966521
-49705.99132
-5.383396351e-06
-1.07104069e-06
-2.808652999e-06
-0.0001066898908
-0.0002239236027
-7.851273644e-08
-0.002704063747
-0.0005663491795
-0.07550986067
-0.00128549481
-3.590772606e-05
-0.1149401044
-0.183403668
-0.000113873958
-1.730096598e-05
-0.02430322741
-0.001076022562
-5.066109523e-06
-0.1068111238
-9.253312815e-05
-91.21436585
-0.01945125628
-8.958261209e-05
-1.678018826e-06
-0.0002943803381
-0.0001299537598
-1.170328977e-06
-0.0001328201344
-3034.923139
-1.378838481
-0.08550576557
-36.1778481
-1.785612005e-06
-9.054761711e-06
-0.001615788007
-1.23208212
-0.0003642605012
-2.708222463e-06
-4.719890307e-06
-0.001710230342
-0.3182780512
-0.0008830648263
-0.08778053115
-0.2629019547
-0.02748771871
-1.833687248e-06
-0.002690141449
-1.89343446e-06
-1.77338727e-06
-2.291261955e-05
-1.31623411e-05
-0.02865008699
-0.0001316279098
-5.992647136
-427898745.7
-0.0002908564713
-180.1583563
-1.149638817
-4.549308837e-06
-0.02938178265
-12.41064482
-0.0005867842208
-1.235864829e-06
-4.615116294e-05
-0.7228122449
-0.0003119308893
-1.100546805e-07
-0.0007569561799
-3.456980613e-06
-0.001902875586
-1.766220045
-0.0006636538355
-435.7550934
-1.468265667e-05
-0.003868037348
-0.005375424637
-9.895862027e-06
-0.002634610055
-0.04246482577
-1.811370936e-05
-0.0008533026641
-1.935841576
-1.965419689e-06
-0.0001940609469
-6.432993714e-05
-436.2325667
-8.585463564e-05
-0.003171908449
-3.1453383e-07
-0.10708771
-1.134414955e-07
-5.434578249e-06
-6.681002905e-06
-0.00246782063
-0.00431428184
-2.421882039e-05
-3.783451257e-08
-5.144794528e-06
-8.550628039e-05
-13660.98467
-117.1131418
-387.5973894
-0.2157832996
-2.335478847e-06
-0.0001022288321
-0.004686932429
-0.4278874185
-0.000106415896
-1.985873409e-06
-2.25291051e-07
-2.945623146e-05
-0.3581923626
-6.785775177e-06
-0.004893941209
-0.0005325322641
-1.060152245e-07
-1.118231806e-05
-2.203633348e-06
-3.969211813e-06
-0.007400958058
-5.769269726e-05
-1.899358848e-06
-0.002291707993
-0.007504512648
-0.0946266303
-0.07641638751
-0.001194744559
-0.000131740617
-3.126413693e-07
-6.676941465e-07
-2.247981499e-07
-1.394670123e-07
-0.0003117083998
-0.003124792631
-3.196614452e-07
-0.0005145950558
-1.196095567
-1.715539058e-05
-2.527890627e-05
-2.665407449e-06
-4.483788358e-07
-3.597134253e-05
-6.514469651e-07
-0.002675227421
-0.1700338633
-0.0151159864
-1.877463443e-05
-0.001018686854
-0.0001769414038
-3.068873861e-06
-1176.630703
-9.378378651e-06
-7.10431401e-07
-4.541139018e-05
-3.954117424e-06
-0.003271203855
-0.01288628186
-0.01758893926
-0.0001040458517
-0.08587993013
-7.493297002e-07
-0.0303353991
-2.525780749e-07
-1.075360445e-06
-0.0005909996117
-3.980399911e-06
-67.71892274
-8.802422373e-06
-3.997018327e-05
-2.690317449e-06
-0.1177580565
-26.26873777
-4.682737548e-05
-0.01239567595
-16.44749478
-0.3150242187
-9.24097764e-06
-1.0502207e-05
-0.3654346078
-0.0002920077007
-0.005014120485
-1.37979468e-05
-2.610550465e-05
-0.0001273284709
-3.179634583e-05
-6.102729506e-07
-0.03043360439
-0.1941067281
-0.02626152766
-2.330385275e-07
-2.089775002
-0.005137045761
-2.078451882e-07
-0.07636633505
-6.943791403e-05
-0.08182004768
-5.307882315e-06
-0.01505065689
-4.601570229e-07
-0.001449579024
-0.0004062483064
-0.1311512449
-6.049424921e-05
-3.109242872
-0.001155891664
-842.4193918
-0.3892954656
-3.580521333e-08
-0.002489034361
-0.02284899485
-0.01965474774
-1.593013311e-06
-0.001587285435
-1.761113439e-06
-1.275562919e-05
-5.438407279e-06
-2.914573789e-06
-17.98383935
-0.001201241561
-0.001049345941
-0.0003002612256
-5.925134807e-05
-1.260239262e-06
-0.0002079300158
-3.029492204e-07
-0.0003910476504
-7917.215229
-2.968281206e-06
-0.4613426334
-9.035490342e-06
-3.515994567e-07
-0.006484242371
-1.667500961e-07
-0.1290091628
-3.723726081e-06
-0.0002047200546
-2.476739768e-07
-0.03234586295
-0.0002393935129
-0.000107753832
-0.5928508576
-0.005892399009
-0.002371984736
-5.106319636e-05
-0.005422903728
-0.005600697388
-0.0001165320035
-8.492553037e-08
-89.84623636
-3714.453168
-0.0008356034504
-2.715363407e-05
-6.770161279e-06
-0.0004350984106
-1.689150634
-0.1094566555
-0.01101851345
-0.1834614184
-3.27091145e-07
-0.004412216962
-8.597553483e-09
-0.3238843541
-1.251350652e-06
-4.686624512e-07
-0.0004352075129
-0.1273052185
-0.0004152526108
-0.05115780463
-0.002149661023
-40050.96602
-0.004341500916
-0.001584949738
-1.272734303e-07
-0.0006057175903
-7.219836807e-05
-0.01372195327
-2.613205505e-07
-0.0872581892
-0.005874236639
-2.280113068
-0.0004066694872
-4.332771697e-05
-3.053780593e-06
-4.008766439e-06
-1.901443143e-06
-0.001449495929
-5.970191013e-07
-7.557409002e-06
-0.09688734499
-2.043706434e-05
-1.169897058e-05
-3.355651365e-07
-4.838231677e-06
-40.18369063
-0.008249553704
-0.0001898868568
-3.316049325e-05
-0.003042105677
-0.0002633075178
-0.001384705507
-1.30613353e-05
-6.96022133e-05
-978473.5934
-1.932222352e-07
-0.0003933467542
-3.324339683e-06
-5.940328548
-0.01413957687
-0.04631344109
-0.001347286968
-0.0006301600371
-8.89982715e-07
-0.0001097966529
-0.02764333143
-0.002043612003
-0.2543376191
-0.02541046238
-0.000423184004
-3.153654508e-05
-4.081867364e-05
-7.869651663e-05
-0.07465825267
-8.687471319e-08
-3.622632916e-06
-8.558042998e-07
-0.2819014239
-1.090279565e-07
-3.478860107e-06
-0.0004982082598
-4.165564726e-05
-2.189001365
-0.01050915398
-0.00595918136
-0.06914156553
-1.6850014e-05
-0.0003787756001
-1.148987389e-07
-1037.221815
-5.22815566e-07
-1.816317404e-05
-0.02271436569
-0.0006177921286
-1.418450826e-05
-39.76431886
-1.133576912e-06
-2.60701124
-0.006872044385
-7.410621558e-07
-5.279974244
-1.77944165
-0.0002133587221
-2.396505542e-07
-0.02482296172
-2.072743803e-06
-2.493066005
-0.000369674698
-1.282580975e-06
-1.4494405e-08
-1.944487433e-06
-0.01613148481
-5.069732503e-09
-3.960253416e-05
-1.943795243e-09
-0.0002019146449
-0.2881271289
-2418.017717
-2.274366855e-06
-1.584923585e-05
-0.000359645391
-6.468635773e-06
-6.305981737e-06
-0.1373436532
-3.568394772e-08
-10.65229514
-0.0005703438949
-0.01088066221
-0.028168428
-5.143407419e-06
-0.002268752908
-0.0006795837317
-2.541900263e-05
-0.0008556471336
-1.07313036e-05
-9.323243019e-08
-2.267392161e-06
-3.058902445e-05
-7.266964483e-05
-3.060187251
-2.487187997e-07
-1424.492308
-0.0002402505927
-1.569564169e-06
-7.104680635e-05
-4.54353651e-07
-0.0009770970965
-0.01776977281
-7.694650607e-06
-0.01549991418
-97.59921393
-0.01273940891
-0.0002475880377
-1.067828133e-05
-0.01247198191
-2.607039514
-28.77996768
-2.075186467e-08
-1.767992115e-06
-0.000206016132
-2.363970384
-6.799545204e-05
-7.976187621e-07
-99.60625075
-2.156457939
-0.01061805752
-1.441250173e-07
-0.0001392548391
-2.601275478e-09
-2.416815175e-07
-0.7244761926
-2.811745822e-07
-2.383801797
-0.0001047482062
-2.477293909e-08
-9.76524148
-6.371548506e-05
-0.003653377205
-3401732.002
-0.01882155215
-8.251532341e-06
-1.708669822e-06
-0.09118209136
-1.334505516e-06
-0.080682921
-0.01496148473
-5.14047694e-06
-1.880492441e-07
-5.634245899e-06
-6.635265567
-2.416245647e-06
-0.08302238032
-28.73381197
-1.376907197e-05
-283.2884318
-0.02336664049
-0.4486810192
-3.220372193e-05
-0.05238590733
-5.323885071
-0.001157033333
-7.270065551e-08
-6.843106933e-07
-7.257694949e-05
-1.009824955
-0.0005874142779
-3.245241682e-08
-3.651272553e-05
-0.0003775347107
-0.1768738775
-0.002440176383
-4.218733473e-05
-2.72934499e-07
-5.253885353e-05
-7.259149485e-05
-0.0008990783601
-0.0001369735152
-2.897508431e-06
-0.0001313481627
-5.826420418e-08
-2.512574343e-06
-0.0003159295598
-0.3219634536
-2.279636105e-05
-0.0003446741397
-10.84610164
-3.150494588e-07
-6.942778454e-06
-0.0308024565
-9.161160741e-05
-3.931909141e-08
-0.03265515469
-7.470806563e-08
-0.04061292556
-0.922813306
-1.205602681e-05
-1.011128217e-05
-2.734796552e-07
-4.666601916e-07
-9.26341834e-06
-0.0001189830883
-4.79504653e-05
-9.40210997e-08
-0.003192196284
-3.878619838e-07
-0.3531284073
-2.254951774e-06
-1.165181138e-07
-731.8701372
-7.075023338e-07
-6.634794626e-06
-5.713134541e-05
-73.4710477
-0.0004575991874
-0.0007895179304
-9.52363626e-08
-1.049954932e-05
-0.02191077518
-4.662049748e-07
-0.04849234967
-20.10109824
-0.0004418549374
-0.0003945670361
-8.658307581e-07
-0.2818049741
-5.962723021e-08
-0.08473995074
-0.003109482617
-0.001143974946
-0.0002579867334
-6.16018968e-06
-0.0002824325622
-6.159790946e-06
-6.353896895e-06
-0.005333350602
-1.875537525e-05
-0.002663668472
-8.951186235e-06
-0.006645111068
-4.90250758e-07
-6.646679078e-06
-1.545898756
-0.3455288649
-0.0002526296983
-0.004192890649
-5.922090089e-05
-3.839797981e-07
-8.422502661e-05
-3.892382254e-08
-0.1593072782
-0.008688694827
-3.460498504
-0.1104609367
-1.279203683e-05
-7.244836207e-07
-1.279268032e-08
-1.587447528e-06
-0.002794843066
-0.156821856
-0.0004833674348
-0.1816771289
-3.352410551e-05
-0.0003992548523
-4.898360006e-07
-0.0003693030689
-0.0002003251286
-7.493081266e-05
-5.559559801e-08
-7.882137483e-08
-1.612910594e-07
-0.0001561218116
-8.945102825e-06
-1.360597005e-06
-0.0003763191539
-4.106279609e-06
-1.855281191e-06
-0.0004436577095
-5.777004064e-06
-0.0112368811
-0.0004092745855
-69.94157962
-0.2202152397
-1.958722773e-05
-19.3543816
-4.538214989e-07
-8.578209079e-06
-0.000116403028
-99.22741868
-0.01694011099
-62.49659164
-3.935367012e-07
-9.323404416e-05
-5.315549677e-05
-4.582641654e-05
-0.0009473078802
-6.118158747e-08
-4.77332888e-05
-2.819319402e-05
-9.534828885e-08
-3.493578518
-0.00120458802
-0.0002129189738
-0.02051978074
-1.185108825e-05
-0.0001769522074
-0.000242558443
-2.998165697e-07
-0.1445117016
-19213.05972
-2.614917799e-05
-2.154984047e-05
-1.065818621e-06
-0.006673934531
-9.598378178e-07
-1.810236331
-3.14499853e-09
-0.0003378887016
-0.007323689147
-0.003375793146
-1.509667047e-06
-0.0001369839097
-0.2097160223
-0.05477519178
-1.538107302e-07
-9.341781708e-06
-1.491658573
-50.18027213
-0.07865855488
-0.4389407401
-0.02394620579
-4.971505948e-06
-5.8158536
-4.798556085e-07
-0.05721134509
-9.535503358
-4.320083719e-05
-6.429836839e-07
-2.390184904e-07
-8.24958057e-06
-0.0002416033838
-8.042993697e-06
-0.01227745354
-2.280489052
-0.03007200558
-0.003131455203
-2.018782507e-06
-3.546187991e-05
-0.003575737125
-2.360447419
-0.0003965339519
-2.026909807e-06
-7.037425242
-0.001148872632
-5.724417221
-0.009887142971
-0.01158110536
-1.639813216e-06
-8.700043459e-06
-5.945561287e-07
-0.02355290158
-1.7933074e-06
-0.007890391396
-0.0001497141614
-2.089856231e-06
-0.4222896481
-8.237456974e-05
-2.149958432e-05
-0.0543901024
-0.001239424138
-7.844565659e-07
-4.36928814e-05
-0.0001011764803
-6.040757484e-08
-1.923331356
-13633.17886
-1.843758011e-07
-6.90380739e-05
-0.02700035281
-4.484299074e-07
-5.369259619e-05
-1.32435878e-05
-2.659066878e-05
-0.009446234067
-1.631778535e-06
-0.0001031459386
-0.009488728284
-2.531652988e-05
-0.0008410272087
-0.3799584599
-0.009542377363
-0.0001225704413
-2.099520291e-06
-126710.9971
-1.834562958
-0.003235750745
-0.811195511
-8.933276453e-06
-0.005973686144
-125.3751052
-0.0004398877841
-1.165775399e-06
-0.1949030136
-3.51445378
-0.1532698797
-0.005948526347
-8.543983996e-05
-0.0004389176051
-4.121262717e-07
-4.922714767e-06
-1.147581634e-06
-2.499966664
-2.250827253e-06
-0.00162441738
-1.574918543e-07
-0.07409277548
-0.01172590186
-5.010592812
-4.289106171e-06
-1.337447725e-07
-8.664771519
-0.0001766657626
-152.8240353
-3.573040925e-05
-5.91245704e-05
-6.295836167
-0.0006334920936
-5.724064003e-08
-7.462284813e-05
-2.188411387e-07
-0.001045839473
-3.750710804
-1.712819677e-05
-0.0008771061809
-0.0001653861352
-0.006676871113
-0.0257858681
-1.995077357e-05
-7.88511206e-06
-4.406534707e-06
-0.05273839985
-8.611156959e-06
-0.0001695003465
-2.179920575e-08
-0.5114110589
-2.461607769e-05
-7.702844117e-06
-0.2277828338
-1.58282654e-05
-9.443942351e-08
-0.1575967804
-3.414467442e-05
-6.21556477e-07
-7.899823458e-06
-7.442861636e-05
-3.687846001e-05
-0.2146650436
-53565.90315
-9.981028388e-06
-0.003946786772
-0.005806122559
-0.0001603468635
diff --git a/t/ME_data/ME_mtinf.dat b/t/ME_data/ME_mtinf.dat
deleted file mode 100644
index 6b585b5..0000000
--- a/t/ME_data/ME_mtinf.dat
+++ /dev/null
@@ -1,832 +0,0 @@
-5.582760808e-08
-2.825811483e-05
-3.895336201e-05
-0.02046456494
-4.688432033e-06
-1.529291749e-07
-6.420743602e-07
-5.905800101e-07
-6.479004963e-07
-1.48600349e-05
-4.880154255e-05
-1.945209343e-06
-0.0004438901791
-1.600557037e-05
-1.231846376e-05
-0.03428118184
-8.553476044e-07
-579.5994998
-6.327145731e-06
-0.01311111831
-2.253068503
-0.001161627534
-0.0007930178104
-3.019642381e-05
-0.0001604361161
-0.0009782409472
-0.03878378325
-3.213295457e-07
-1.942127118
-9.002578443e-08
-0.04615078581
-0.1667432966
-1.060974436e-06
-0.004405656545
-0.02198234308
-3.75309762e-05
-6.894134543e-07
-0.006337938566
-0.000232790789
-0.004628411362
-0.0005874171572
-0.0114321818
-0.002787500992
-2.287103816e-05
-0.0001818288916
-0.009239536849
-0.03250478346
-3.095188489e-07
-5.383217293e-05
-12.03590233
-0.0001869627067
-4.288822951e-06
-9.903075819e-05
-0.0001133438964
-0.1413834117
-0.02374806527
-5.439959182
-0.0002200036282
-1.257502016
-0.00554436291
-1.492794515e-06
-0.0410145385
-0.003110926104
-4.895681474
-9.566466236e-08
-4.29754523e-06
-1.26422416e-05
-1.50686362e-05
-4.031356768e-06
-0.0009303622768
-3.581771412e-06
-2.209149529
-0.001347543465
-0.001948045992
-25948.6239
-4.798288837e-07
-0.001267876802
-2.198489834e-07
-0.02568760469
-1.738239331e-06
-5.383639367e-05
-0.03615040778
-48368.71222
-6.221468888e-06
-8.272280776e-06
-2.436735404e-06
-0.00082069939
-0.0001956879791
-8.464294968e-08
-0.002499166028
-0.0005268877995
-0.08401690492
-0.001109681865
-1.676530038e-05
-0.1210659795
-0.2021753863
-0.0001468722356
-1.49843416e-05
-0.01797157389
-0.001470636482
-6.108534705e-06
-0.1147500209
-6.044324523e-05
-83.23969297
-0.01952725493
-0.0001089245711
-5.957305171e-07
-0.00020821416
-0.0001540635133
-1.358665617e-06
-0.0001977851353
-2969.178567
-1.434474792
-0.07897582528
-33.78910043
-1.944172541e-06
-8.511420265e-06
-0.001671034953
-1.106372417
-0.000475762259
-2.254611547e-06
-5.024971782e-06
-0.002685315998
-0.3204715262
-0.0006830156865
-0.06161172452
-0.3792093477
-0.03062030789
-1.284034864e-06
-0.002336243826
-1.185574165e-05
-1.710024271e-06
-2.253363854e-05
-1.300040157e-05
-0.02939378205
-0.0001573213338
-6.400482717
-484632185.3
-0.0002916173761
-171.9267844
-1.633949616
-3.749772544e-06
-0.03071686489
-13.06467703
-0.0006848589426
-2.117832499e-06
-2.785336442e-05
-1.197689446
-0.0002238650692
-1.488235642e-07
-0.0007409854778
-6.610729458e-06
-0.002123335104
-1.564470647
-0.0006669917367
-473.2723045
-1.432945457e-05
-0.003779071212
-0.005378672443
-4.446772576e-06
-0.002341093315
-0.1693157671
-1.832833621e-05
-0.0006895620169
-1.807274869
-1.739710164e-06
-0.0002259708712
-6.536256314e-05
-438.0279281
-8.402208741e-05
-0.003349479787
-1.893640387e-07
-0.1112826341
-2.103910568e-07
-6.061483157e-06
-6.743144434e-06
-0.002108199475
-0.003016819088
-2.679236437e-05
-4.152061832e-08
-6.910237845e-06
-4.414649231e-05
-13180.47273
-112.3954293
-407.9774509
-0.1947783364
-3.633741266e-06
-8.949320507e-05
-0.004538125269
-0.437437773
-9.781170001e-05
-2.778563584e-06
-1.70122334e-06
-3.067502446e-05
-0.3572469398
-4.865503412e-06
-0.005054297163
-0.0005309402179
-2.786914424e-08
-1.278066173e-05
-9.682183899e-06
-2.538469656e-06
-0.006448672409
-6.607293322e-05
-1.895984721e-06
-0.001842602862
-0.07109867357
-0.09292074641
-0.08069062114
-0.0007841374424
-0.0001302148142
-6.96049157e-07
-4.632599032e-07
-2.355250248e-07
-1.220815109e-06
-0.0003380714588
-0.003131420086
-6.328928995e-07
-0.000531365034
-1.039688858
-1.577333165e-05
-4.032814959e-05
-2.142451659e-06
-4.632583267e-07
-3.564353255e-05
-7.645988124e-07
-0.004094577112
-0.171048394
-0.009405424559
-1.086927134e-05
-0.0009935783254
-0.0001283118369
-3.018764055e-06
-1142.596042
-6.149007778e-06
-5.946922879e-06
-4.484486221e-05
-4.875515638e-06
-0.003370278763
-0.01348136435
-0.01694870912
-8.663774537e-05
-0.07732984975
-2.265338646e-06
-0.0303766268
-6.402972472e-07
-2.703996094e-06
-0.0006768238778
-2.4208244e-06
-64.78560774
-8.344522601e-06
-3.189381378e-05
-1.712603497e-06
-0.09667429817
-25.39833184
-2.743944024e-05
-0.01296043652
-20.1569871
-0.2322952925
-8.092460701e-06
-7.683409336e-06
-4.593017291
-0.0002449299413
-0.00523879726
-1.093000481e-05
-3.456966245e-05
-0.002971509085
-3.393728092e-05
-7.66629053e-07
-0.0296687921
-0.2319316424
-0.02436745947
-2.36332202e-07
-3.206954494
-0.001525352593
-2.188583307e-07
-0.06694976109
-7.241407041e-05
-0.08335446798
-3.264976398e-06
-0.01196235312
-2.498236803e-07
-0.001644171751
-0.0004035076896
-0.1306966818
-4.132360995e-05
-3.046243795
-0.003018382027
-881.3490404
-0.4546667796
-9.064466579e-08
-0.002163542438
-0.02494252294
-0.2044042461
-1.599883465e-06
-0.007429973372
-1.09914769e-05
-1.169335913e-05
-4.821816233e-06
-3.011624738e-06
-12.90229638
-0.001487606977
-0.0007102204656
-0.0003222391159
-6.637367056e-05
-9.602293783e-07
-0.00030032421
-5.892234555e-07
-0.0007800663346
-12693.45926
-1.246239219e-06
-0.4608136937
-4.686105226e-06
-3.674558139e-07
-0.007204718245
-2.497379172e-07
-0.1285504712
-6.664917332e-06
-0.0002313036825
-2.723172326e-07
-0.02540409803
-0.0002585621672
-0.0001530539048
-0.5030883756
-0.006200984495
-0.002372899322
-5.068583054e-05
-0.004070797031
-0.005964523973
-0.0003127185963
-9.497442021e-08
-96.42268432
-4083.471729
-0.0006980174959
-3.015487844e-05
-8.754876963e-06
-0.000471310047
-1.61799055
-0.09080692266
-0.01178844473
-0.2041351557
-4.108894361e-07
-0.004272505699
-8.935607949e-09
-0.2820889003
-4.921322971e-07
-5.026967052e-07
-0.0004235166152
-0.1180618535
-0.0004099511171
-0.05308356847
-0.002175398288
-39909.86411
-0.00495724079
-0.001561254132
-4.274998406e-07
-0.0005440402015
-6.670698143e-05
-0.007258668602
-2.60943183e-07
-0.08498614435
-0.003815925744
-2.222628302
-0.0003541676329
-0.0002242870094
-3.562100507e-06
-4.429739157e-06
-1.054004097e-06
-0.00149447791
-2.446511801e-07
-8.312474957e-06
-0.09692624138
-1.962072264e-05
-1.2633533e-05
-7.957449465e-07
-5.360998441e-06
-48.84222045
-0.004935712395
-0.0001615305681
-3.485868016e-05
-0.002967534836
-0.0001747606262
-0.001368673963
-1.342339291e-05
-5.663732737e-05
-951296.5058
-1.357017004e-06
-0.0002430876264
-3.339872253e-06
-6.34972199
-0.01196251522
-0.04158958011
-0.001124850802
-0.0005723482545
-3.76044586e-06
-0.0001107580785
-0.02979295371
-0.001963441807
-0.2509977533
-0.02148405854
-0.000477635843
-3.64858986e-05
-4.22822887e-05
-7.374986907e-05
-0.07628127654
-1.225401468e-07
-2.210039033e-06
-2.296422651e-06
-0.3334783046
-3.635226811e-06
-1.934035792e-06
-0.000495060149
-4.36935663e-05
-2.300727955
-0.01283761114
-0.00557883427
-0.0719699281
-1.82181339e-05
-0.0002216021515
-1.270486172e-07
-1940.25235
-5.666656484e-07
-1.860611735e-05
-0.02277334026
-0.0004230285398
-9.859060738e-06
-38.8284244
-8.550126322e-07
-3.194573173
-0.006724981284
-8.390467411e-07
-5.394729567
-1.782355254
-0.0003270406706
-6.012608249e-07
-0.02626165732
-3.670279208e-05
-2.520894603
-0.0002698573318
-1.489017532e-06
-2.274548381e-07
-8.737249392e-07
-0.01456085266
-5.457807092e-09
-3.939935595e-05
-5.512277202e-08
-0.0001532601043
-0.3195672069
-2373.003409
-2.764694852e-06
-1.305719793e-05
-0.0003978915682
-0.0004948555655
-3.280521304e-06
-0.1355744285
-9.539967905e-09
-13.00984733
-0.0006801756974
-0.009644637547
-0.01845330108
-6.988596187e-06
-0.002416821323
-0.00139251527
-2.210462606e-05
-0.0007144938597
-3.098029761e-06
-1.062795442e-07
-1.540921293e-06
-2.32085136e-05
-5.906898983e-05
-2.815067914
-2.626278681e-07
-1274.911143
-0.0002414189765
-1.41029294e-06
-4.085495754e-05
-1.469023465e-07
-0.0009099324663
-0.01571147903
-5.271086554e-06
-0.01643989004
-95.99695683
-0.01147762039
-0.0004091694655
-1.129916931e-05
-0.01252602377
-4.580344108
-28.23428139
-2.368518601e-08
-2.513102323e-06
-0.0002047039974
-4.37385486
-5.204147828e-05
-1.008501642e-06
-95.82595938
-2.175220964
-0.0115792253
-1.661729845e-07
-0.0001629816332
-3.481252916e-09
-3.932760717e-07
-0.6931415301
-2.649359888e-06
-19.89107003
-0.0001015626597
-6.352730089e-08
-9.11889889
-0.0001280302506
-0.003632084954
-3124257.466
-0.01997968917
-9.062763488e-06
-1.806685549e-06
-0.1029017089
-2.388024181e-05
-0.0775425829
-0.01376323377
-5.039646218e-06
-4.796601753e-07
-5.928185403e-06
-6.609428695
-2.071273115e-06
-0.06145615574
-34.02694087
-0.0001099850522
-309.0075882
-0.0252581086
-0.6382068759
-0.0001946000044
-0.0516754215
-7.972945986
-0.001197095771
-8.848572412e-08
-7.101584754e-07
-8.849498802e-05
-1.023229262
-0.0006030655775
-4.908547066e-08
-2.778161765e-05
-0.0003856084955
-0.1802215093
-0.00251917308
-4.852765701e-05
-9.377007406e-06
-5.223044658e-05
-7.844738702e-05
-0.000823109541
-0.0001514552108
-2.072745569e-06
-0.0002134409917
-7.112786559e-08
-2.164548157e-06
-0.0003701532283
-0.312453918
-2.434168445e-05
-0.0003196660825
-12.97570663
-1.904837416e-07
-7.687264107e-06
-0.03050279549
-4.836388196e-05
-5.182325816e-08
-0.03523445202
-7.50629725e-08
-0.06397706877
-2.980733906
-3.332054567e-05
-9.304399694e-06
-3.605909581e-07
-3.495748221e-07
-2.909153473e-05
-0.0001259318878
-3.342472966e-05
-1.190805029e-07
-0.003194334166
-3.802662978e-07
-0.3086555319
-2.265399049e-06
-1.240985629e-07
-731.5204
-3.878565182e-07
-5.265372774e-06
-5.881803437e-05
-72.22818259
-0.0005582677709
-0.0007784672343
-3.288490726e-07
-5.842140857e-06
-0.02113918753
-1.269654727e-06
-0.04675464579
-25.32179403
-0.0004668340387
-0.0002363968869
-9.693807373e-07
-0.2691575098
-7.591428395e-08
-0.1238104135
-0.002057529127
-0.001297319821
-0.0002578808711
-7.16917206e-06
-0.0003134969133
-1.89886682e-05
-6.816646149e-06
-0.005208076232
-1.553721275e-05
-0.002163306945
-2.08087712e-05
-0.005912329881
-7.100353372e-07
-6.064675674e-06
-0.9278861068
-0.5398096181
-0.0002482834202
-0.003251431474
-0.0001351409005
-2.255074238e-07
-0.0001966273358
-6.081228169e-08
-0.1621182535
-0.01356394984
-3.352599887
-0.1093060945
-2.033105897e-05
-7.790291723e-07
-2.741226766e-08
-8.197142111e-07
-0.002938152183
-0.1688720412
-0.000878388148
-0.1894727785
-3.540031572e-05
-0.004719476749
-2.505006373e-06
-0.000370334811
-0.0001955296755
-3.081872404e-05
-6.16568903e-08
-8.125050517e-08
-2.841463125e-07
-0.0001447605718
-9.723473085e-06
-3.842440928e-07
-0.0002002004501
-2.876662005e-06
-3.206621379e-06
-0.0004757973033
-4.080935553e-06
-0.01811675801
-0.0007422047088
-68.55774269
-0.1990430793
-2.759390865e-05
-19.36133318
-4.056803417e-07
-7.078017532e-06
-0.0001003816537
-95.07498601
-0.0166477342
-58.88190306
-5.649108811e-06
-7.793584968e-05
-3.33470689e-05
-4.499180952e-05
-0.001054772305
-3.624979571e-07
-3.383705855e-05
-2.096652126e-05
-1.473256752e-07
-3.532179606
-0.001202187739
-0.0002243301291
-0.03530559137
-1.42504564e-05
-0.0001771562336
-0.0005429903338
-8.539507323e-07
-0.1459649943
-20128.72412
-2.201327763e-05
-1.696658284e-05
-3.790831048e-07
-0.04339432699
-9.266881901e-07
-1.887969395
-2.176375742e-08
-0.00041568578
-0.005946060998
-0.003630083583
-1.843851889e-06
-0.0001066579874
-0.1832324847
-0.0549419187
-1.690213495e-07
-2.903264073e-05
-0.9908674856
-40.39497028
-0.06466477233
-0.4552476095
-0.02983905871
-3.007284175e-06
-7.663231527
-8.632942557e-07
-0.05647465611
-9.475595107
-3.173116103e-05
-6.500398209e-07
-4.795277195e-07
-3.084669237e-05
-0.0001714345968
-9.284790614e-06
-0.01835802147
-2.417931971
-0.03035911898
-0.002916783982
-2.055801117e-06
-5.608563674e-05
-0.003529644789
-2.402227217
-0.0004235577219
-1.616006288e-06
-8.793458292
-0.001511400213
-5.513047874
-0.0299723301
-0.02530694857
-1.717601897e-06
-7.331930421e-06
-3.358969417e-07
-0.02404885266
-4.304661703e-05
-0.01008941554
-0.0007910273119
-2.106198066e-06
-0.4584734408
-8.699683257e-05
-6.48207849e-05
-0.05740194674
-0.003273869341
-7.745530251e-07
-4.450468335e-05
-7.060712812e-05
-3.43962665e-07
-1.885858086
-12596.28093
-7.072848576e-08
-6.227410044e-05
-0.03238383706
-1.250340711e-06
-4.945891575e-05
-1.370420025e-05
-2.608359212e-05
-0.008367079389
-1.618717241e-06
-8.263294097e-05
-0.009865263551
-2.46915343e-05
-0.0006910733202
-0.3051615944
-0.006723371218
-0.0001049903836
-2.090318377e-06
-122200.0192
-6.4260476
-0.00311260923
-0.804791076
-9.274397022e-06
-0.00597987096
-120.084753
-0.0004896982366
-5.653748148e-07
-0.1951903766
-3.914175763
-0.2811999263
-0.006352658402
-0.0001248073811
-0.0004464717139
-8.169285959e-06
-5.314216957e-06
-3.487632483e-06
-2.491846834
-2.244331191e-06
-0.001698490011
-5.486417742e-07
-0.07593626881
-0.01172688569
-5.245932064
-2.769113043e-06
-1.578286837e-07
-8.819666854
-0.0001747929634
-315.309908
-3.248098689e-05
-0.0001311790738
-6.108978378
-0.000623564243
-1.994199188e-07
-0.0001065229989
-2.805236353e-07
-0.001004843236
-4.091307134
-1.921939229e-05
-0.0009479594902
-9.252608359e-05
-0.007417908515
-0.02541094922
-2.087328265e-05
-1.096806923e-05
-4.220356991e-06
-0.1765502926
-6.331905678e-06
-0.0001671692355
-3.555632432e-08
-0.5800063575
-2.038149359e-05
-5.536258117e-06
-0.2750958952
-7.26265839e-05
-3.711138323e-07
-0.153305225
-2.41543512e-05
-9.396344975e-07
-3.729430265e-06
-7.422782227e-05
-3.769490057e-05
-0.2081251952
-54520.66861
-8.707297925e-06
-0.004163637881
-0.0056787473
-0.0001293319359
diff --git a/t/ME_data/ME_mtmb.dat b/t/ME_data/ME_mtmb.dat
deleted file mode 100644
index 9c9cedf..0000000
--- a/t/ME_data/ME_mtmb.dat
+++ /dev/null
@@ -1,832 +0,0 @@
-4.043634525e-08
-3.947735809e-05
-3.978333671e-05
-0.01707962353
-4.133534858e-06
-4.12758231e-07
-6.527379321e-07
-7.807924441e-07
-6.213247187e-07
-2.158745738e-05
-4.929604661e-05
-1.648537835e-06
-0.0005294944592
-1.551221705e-05
-9.732109994e-06
-0.03508523471
-8.385697667e-07
-608.0814044
-6.0749939e-06
-0.01370931025
-1.96358809
-0.0005777778828
-0.000827073758
-3.396653466e-05
-0.0003663185905
-0.000712177849
-0.03945983226
-6.694529237e-07
-2.002990871
-9.107975513e-08
-0.04804489153
-0.1078671772
-9.503786206e-07
-0.00665553606
-0.02341176253
-3.684250185e-05
-6.861318278e-07
-0.0054005165
-7.647108683e-05
-0.005354950581
-0.0006200616798
-0.008448060377
-0.003448403835
-3.215759836e-05
-0.0001546133101
-0.005591549046
-0.03949683081
-2.350530066e-07
-5.66847762e-05
-8.602097332
-0.0001911650059
-5.061383546e-06
-5.119802912e-06
-0.0001348170624
-0.1914012121
-0.02428150541
-4.539435975
-0.0002614830785
-1.314349378
-0.005127584311
-2.106730988e-06
-0.06042114066
-0.002877526755
-6.46292188
-7.320533172e-08
-4.686272323e-06
-1.674918088e-05
-1.561186199e-05
-1.127468541e-06
-0.0009956223126
-9.322344484e-07
-2.355779109
-0.001354263217
-0.001607437479
-27300.17114
-8.763842677e-07
-0.0006972261753
-1.071457504e-07
-0.02005686587
-1.961498186e-06
-4.460773766e-05
-0.02601668696
-50724.57307
-5.480882962e-06
-1.076584976e-06
-2.869362531e-06
-0.0001072518549
-0.0002293795185
-8.022114184e-08
-0.002698012227
-0.0005830909647
-0.07725070316
-0.001325583365
-3.695641615e-05
-0.1179283377
-0.187466976
-0.0001155356229
-1.778799047e-05
-0.02462347248
-0.001087874101
-5.154078743e-06
-0.1031621088
-9.484855944e-05
-93.85073623
-0.02033391802
-9.092704303e-05
-1.705551286e-06
-0.0002974887579
-0.0001322670062
-1.193514887e-06
-0.0001344144727
-3095.568362
-1.406347883
-0.08812981858
-35.91258076
-1.82576512e-06
-9.27160445e-06
-0.001666287395
-1.269161514
-0.0003682229129
-2.790042514e-06
-4.841321328e-06
-0.001729869633
-0.3293946648
-0.0008921922127
-0.08111101834
-0.2664463909
-0.02792004132
-1.869728836e-06
-0.002768079641
-1.908915912e-06
-1.7932868e-06
-2.342012655e-05
-1.3539045e-05
-0.02918304927
-0.0001341037055
-6.099926061
-438766540.3
-0.0003006556126
-182.0743425
-1.166556788
-4.685195193e-06
-0.02991092642
-12.79082303
-0.0005990054968
-1.245221356e-06
-4.607261468e-05
-0.729519575
-0.0003077751881
-1.113634813e-07
-0.0007837766949
-3.49282179e-06
-0.001938345149
-1.741954359
-0.000679566605
-447.8606334
-1.502685811e-05
-0.003972927121
-0.005572294955
-1.00407234e-05
-0.002697149585
-0.04268640856
-1.838639564e-05
-0.0008683709605
-1.984375256
-2.023641764e-06
-0.000198987372
-6.590604109e-05
-451.0223058
-8.870879716e-05
-0.003243511862
-3.187920803e-07
-0.1107987679
-1.145759171e-07
-5.532565231e-06
-6.903234805e-06
-0.002276208464
-0.004410664221
-2.467146722e-05
-3.863236208e-08
-5.225855864e-06
-8.615548536e-05
-13702.18793
-118.8790931
-398.6816605
-0.2170834277
-2.355356756e-06
-0.0001050726086
-0.004854583584
-0.4416609149
-0.0001091431342
-2.015137689e-06
-2.270813172e-07
-3.065037179e-05
-0.3695610208
-6.626431375e-06
-0.005034520677
-0.0005496323042
-1.077013565e-07
-1.134883954e-05
-2.217217266e-06
-4.053576036e-06
-0.007630107217
-5.870346284e-05
-1.937770986e-06
-0.002308560802
-0.00756714829
-0.09713373161
-0.07865638607
-0.001166988788
-0.0001360267898
-3.146930172e-07
-6.721239549e-07
-2.308081382e-07
-1.402445764e-07
-0.000317796034
-0.003250721897
-3.220983676e-07
-0.0005361385069
-1.1674911
-1.701754998e-05
-2.555659626e-05
-2.72704331e-06
-4.575225674e-07
-3.723503521e-05
-6.625793946e-07
-0.002713331741
-0.1757705927
-0.01528287892
-1.914346383e-05
-0.001027883901
-0.0001819228709
-3.169381656e-06
-1197.183436
-9.459688271e-06
-7.151860571e-07
-4.709453403e-05
-4.024117529e-06
-0.003399032872
-0.01326938543
-0.0181234889
-0.0001070400499
-0.08816817298
-7.557106621e-07
-0.03150507468
-2.546159018e-07
-1.083218854e-06
-0.0006029839411
-4.109536026e-06
-69.62539127
-9.053551382e-06
-3.976334138e-05
-2.74610264e-06
-0.1217200425
-26.75585209
-4.684309171e-05
-0.0125877428
-16.76545876
-0.3169610574
-9.353760931e-06
-1.07865377e-05
-0.3676199238
-0.0002880060428
-0.005152194769
-1.41099509e-05
-2.650294015e-05
-0.0001281502944
-3.240562789e-05
-6.20473578e-07
-0.0308225291
-0.1976686504
-0.02685240634
-2.401876547e-07
-2.123598943
-0.004691551094
-2.124344304e-07
-0.07570215778
-7.036967951e-05
-0.0844338852
-5.219464468e-06
-0.01538384477
-4.721421083e-07
-0.001472983279
-0.0004188110921
-0.1327823736
-6.152491259e-05
-3.239586149
-0.001175629336
-871.2691721
-0.3969457222
-3.611093443e-08
-0.002547168341
-0.02355053075
-0.01978826536
-1.629432987e-06
-0.001596181487
-1.774300225e-06
-1.307034784e-05
-5.363905162e-06
-2.975238583e-06
-18.34285454
-0.001221513765
-0.0009527644543
-0.0003056723918
-6.023182126e-05
-1.298019757e-06
-0.0002110795487
-3.062568333e-07
-0.0003944099269
-8017.105462
-3.048710777e-06
-0.4712988875
-9.179254412e-06
-3.590591627e-07
-0.006691166206
-1.686803436e-07
-0.1318215783
-3.772842175e-06
-0.0002081270913
-2.51985643e-07
-0.03315604586
-0.0002442306805
-0.0001036935455
-0.6114605634
-0.005979836691
-0.002454121593
-5.221039976e-05
-0.005566052602
-0.005695987836
-0.0001177971925
-8.660478127e-08
-92.35458829
-3774.456141
-0.0008295364887
-2.790062671e-05
-6.887441946e-06
-0.0004559717764
-1.726900007
-0.1119717456
-0.01128883757
-0.1878506264
-3.315999686e-07
-0.004492437144
-8.865150661e-09
-0.3300895513
-1.275749227e-06
-4.798867243e-07
-0.000430605922
-0.1309265622
-0.0004243515512
-0.05223983499
-0.002194423992
-41439.51884
-0.004465008791
-0.001613828195
-1.282980178e-07
-0.0006190141029
-7.194498513e-05
-0.0127208867
-2.657540728e-07
-0.08995942708
-0.00587596777
-2.264314747
-0.0004128827152
-4.362055433e-05
-3.094154924e-06
-4.090167641e-06
-1.882247111e-06
-0.001494985492
-6.072358792e-07
-7.78861065e-06
-0.0996671474
-2.087743173e-05
-1.195758488e-05
-3.380560937e-07
-4.91156394e-06
-40.88650322
-0.008316202508
-0.0001878582209
-3.404003278e-05
-0.003154799666
-0.0002683422743
-0.001433303192
-1.333987217e-05
-6.541815356e-05
-1013853.193
-1.940142974e-07
-0.0004055213885
-3.394978594e-06
-6.063971952
-0.01415112473
-0.04748322673
-0.001347493205
-0.0006384625135
-8.930551094e-07
-0.0001131648923
-0.02696899175
-0.00210540585
-0.2604784401
-0.02517939777
-0.0004319754673
-3.063538811e-05
-4.193167229e-05
-8.029758676e-05
-0.07617825103
-8.947767582e-08
-3.564349846e-06
-8.615322031e-07
-0.2862634339
-1.095651653e-07
-3.556638869e-06
-0.0005160716003
-4.254724595e-05
-2.251498738
-0.01064347946
-0.00595090805
-0.07130582332
-1.715996605e-05
-0.0003787498521
-1.170286199e-07
-1048.330067
-5.346123785e-07
-1.85189906e-05
-0.02346029109
-0.0006336000299
-1.463283193e-05
-40.56959009
-1.155413601e-06
-2.65750517
-0.006991122396
-7.550294362e-07
-5.39504168
-1.819080071
-0.0002156785189
-2.41982949e-07
-0.02549539581
-2.082507467e-06
-2.550531382
-0.0003820138242
-1.299986774e-06
-1.455908814e-08
-2.003136276e-06
-0.01637378355
-5.186594387e-09
-4.050969997e-05
-1.950260502e-09
-0.0002073579324
-0.2954012611
-2507.99095
-2.303197695e-06
-1.626886935e-05
-0.000358110058
-6.506359464e-06
-6.349514789e-06
-0.141951511
-3.653167133e-08
-10.83944388
-0.000581830268
-0.01113402828
-0.02901228806
-5.200506604e-06
-0.002314138958
-0.0006876278156
-2.598198607e-05
-0.0008673190386
-1.104294394e-05
-9.472389228e-08
-2.273856546e-06
-3.102576679e-05
-7.204946629e-05
-3.149762116
-2.553239064e-07
-1447.112074
-0.0002454065232
-1.585493248e-06
-6.936515813e-05
-4.61710293e-07
-0.0009759829851
-0.01797761937
-7.919050081e-06
-0.015785171
-101.1723454
-0.01276989221
-0.0002497757192
-1.085505893e-05
-0.01278096264
-2.62760061
-29.78160572
-2.114860505e-08
-1.796977691e-06
-0.000212496441
-2.39626398
-6.650574936e-05
-8.088806448e-07
-101.0676422
-2.229055132
-0.01080386956
-1.463589916e-07
-0.0001417681527
-2.631212208e-09
-2.437682518e-07
-0.7418377718
-2.827578969e-07
-2.398079218
-0.0001065771066
-2.489806849e-08
-9.908827254
-6.309778751e-05
-0.00378003465
-3448698.134
-0.0191559099
-8.431146522e-06
-1.751755013e-06
-0.09328656347
-1.344398036e-06
-0.0832260718
-0.01527786949
-5.257717365e-06
-1.896902995e-07
-5.795109393e-06
-6.921522693
-2.430515824e-06
-0.08522895716
-29.14927122
-1.391614037e-05
-291.2904419
-0.02368583104
-0.455782889
-3.245758864e-05
-0.05414701889
-5.397475156
-0.00117384466
-7.373271822e-08
-6.998778683e-07
-7.360475302e-05
-1.031451727
-0.0006051798457
-3.280229684e-08
-3.658394412e-05
-0.0003902233216
-0.1847558715
-0.002502672198
-4.330167288e-05
-2.746526695e-07
-5.382406862e-05
-7.387092921e-05
-0.0009192813545
-0.0001391713976
-2.955690693e-06
-0.000132686335
-5.921191753e-08
-2.502885009e-06
-0.0003262384663
-0.3265859428
-2.334781115e-05
-0.0003350764479
-11.0179341
-3.140510375e-07
-7.10747644e-06
-0.03028450502
-9.387163118e-05
-3.980818682e-08
-0.03318852235
-7.700654025e-08
-0.04100966415
-0.9306026389
-1.215482723e-05
-1.034771565e-05
-2.770367845e-07
-4.819430633e-07
-9.354060954e-06
-0.0001213408411
-4.808376106e-05
-9.523603178e-08
-0.003303910055
-3.999500869e-07
-0.3387252545
-2.327975106e-06
-1.193353844e-07
-757.1807265
-7.115390345e-07
-6.560832522e-06
-5.770150753e-05
-75.80554368
-0.0004649766183
-0.0008181187127
-9.611138e-08
-1.050468883e-05
-0.02231539522
-4.714535231e-07
-0.04960630889
-20.45328123
-0.0004507112914
-0.0003943023904
-8.83939796e-07
-0.2888949087
-6.028232019e-08
-0.08199621721
-0.003189896256
-0.001178339913
-0.0002676374351
-6.300649506e-06
-0.0002877679041
-6.219862368e-06
-6.536771847e-06
-0.005420161908
-1.929584044e-05
-0.002684252623
-9.007563628e-06
-0.006849588416
-4.982293287e-07
-6.550466828e-06
-1.525211127
-0.3494894924
-0.0002577621296
-0.004105595304
-5.973345851e-05
-3.857080672e-07
-8.503467113e-05
-3.923603006e-08
-0.1589617537
-0.00876907058
-3.523322109
-0.1145258675
-1.314045265e-05
-7.418942915e-07
-1.324490352e-08
-1.636964114e-06
-0.002880065675
-0.1592007386
-0.000488013888
-0.1854220636
-3.428390405e-05
-0.000400235162
-4.927032357e-07
-0.0003805175174
-0.000206049728
-7.750331497e-05
-5.676053662e-08
-8.091087725e-08
-1.629883973e-07
-0.0001578475837
-9.210829305e-06
-1.399868064e-06
-0.0003876430699
-4.2060549e-06
-1.883530473e-06
-0.0004551528347
-5.95661238e-06
-0.01133751836
-0.0004124881431
-71.4449825
-0.2266258562
-1.982517645e-05
-20.09198674
-4.459887837e-07
-8.668133286e-06
-0.0001140358258
-102.1596809
-0.01756860925
-62.01103257
-3.965683742e-07
-9.023902295e-05
-5.460174999e-05
-4.68974667e-05
-0.0009707252228
-6.16128967e-08
-4.820773721e-05
-2.909887948e-05
-9.643723985e-08
-3.587151784
-0.001228769589
-0.0002189407164
-0.02067709918
-1.209973132e-05
-0.0001786793751
-0.0002447224315
-3.025119469e-07
-0.1497591661
-19584.06048
-2.500843214e-05
-2.04819022e-05
-1.095592696e-06
-0.006707927898
-9.546423416e-07
-1.866364403
-3.157061971e-09
-0.0003427210638
-0.007403790719
-0.003447595712
-1.541603529e-06
-0.0001414653533
-0.2149924951
-0.05682353296
-1.567122577e-07
-9.425600019e-06
-1.516658444
-50.26258345
-0.08123197798
-0.4596831029
-0.02440378929
-5.099462111e-06
-5.906994171
-4.856007707e-07
-0.05934270293
-9.872551101
-4.453291031e-05
-6.635660769e-07
-2.414882955e-07
-8.323783478e-06
-0.0002494377442
-8.203056301e-06
-0.01241361914
-2.321986401
-0.03082695358
-0.003121736152
-2.087698212e-06
-3.593358952e-05
-0.003653855273
-2.446855355
-0.0004048717468
-2.08378296e-06
-7.165260738
-0.001166207815
-5.887536838
-0.01013970123
-0.01167811614
-1.669471654e-06
-8.203051977e-06
-6.120151852e-07
-0.02407438684
-1.799870317e-06
-0.008006468116
-0.0001507123435
-2.105686938e-06
-0.43233826
-8.455083902e-05
-2.171052457e-05
-0.05586438426
-0.001252930378
-8.127056447e-07
-4.513680206e-05
-0.000101171823
-6.088181113e-08
-1.986796246
-13982.8461
-1.890254682e-07
-6.95701156e-05
-0.02773032584
-4.52781786e-07
-5.545292211e-05
-1.350514515e-05
-2.770516139e-05
-0.009515518928
-1.670611457e-06
-0.0001058171739
-0.009687328162
-2.61605968e-05
-0.0008428900179
-0.3772131567
-0.009451023908
-0.0001231850963
-2.166013419e-06
-129053.5102
-1.855117513
-0.003319355848
-0.8451507447
-9.204040714e-06
-0.006100179513
-127.1446333
-0.0004502251479
-1.200723497e-06
-0.1990753321
-3.594711114
-0.1555226613
-0.00608608373
-8.62682525e-05
-0.0004527910284
-4.13348428e-07
-5.012452903e-06
-1.156677713e-06
-2.539881224
-2.304629442e-06
-0.001603902323
-1.585646341e-07
-0.07566755294
-0.01213301266
-5.167295408
-4.417103555e-06
-1.361518308e-07
-8.948305052
-0.0001806738857
-154.3418669
-3.623518482e-05
-5.974885223e-05
-6.401585137
-0.0006471766428
-5.756911442e-08
-7.587748259e-05
-2.226631044e-07
-0.001059477349
-3.833339206
-1.745860723e-05
-0.0008961167773
-0.0001663597116
-0.006854517077
-0.02656008377
-2.034231653e-05
-7.965188898e-06
-4.2693379e-06
-0.05310976055
-8.657283413e-06
-0.0001752091525
-2.196673634e-08
-0.5195738207
-2.507876817e-05
-7.96780687e-06
-0.2311732213
-1.593154072e-05
-9.508397921e-08
-0.160706417
-3.362657365e-05
-6.278716013e-07
-8.130699775e-06
-7.350692748e-05
-3.801653842e-05
-0.2136266576
-55290.08117
-1.024416571e-05
-0.004063025436
-0.006012059135
-0.0001637062434
diff --git a/t/ME_data/PSP.lhe.gz b/t/ME_data/PSP.lhe.gz
deleted file mode 100644
index 98df72d..0000000
Binary files a/t/ME_data/PSP.lhe.gz and /dev/null differ
diff --git a/t/ME_data/PSP_h.lhe.gz b/t/ME_data/PSP_h.lhe.gz
new file mode 100644
index 0000000..02d1a23
Binary files /dev/null and b/t/ME_data/PSP_h.lhe.gz differ
diff --git a/t/ME_data/PSP_jets.lhe.gz b/t/ME_data/PSP_jets.lhe.gz
new file mode 100644
index 0000000..bd8379b
Binary files /dev/null and b/t/ME_data/PSP_jets.lhe.gz differ
File Metadata
Details
Attached
Mime Type
text/x-diff
Expires
Mon, Jan 20, 10:27 PM (1 d, 10 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4242662
Default Alt Text
(199 KB)
Attached To
rHEJ HEJ
Event Timeline
Log In to Comment