Page MenuHomeHEPForge
Files F8725676

No OneTemporary
Actions

View Options

diff --git a/FixedOrderGen/CMakeLists.txt b/FixedOrderGen/CMakeLists.txt
index f2d3ace..da601e4 100644
--- a/FixedOrderGen/CMakeLists.txt
+++ b/FixedOrderGen/CMakeLists.txt
@@ -1,72 +1,72 @@
cmake_minimum_required(VERSION 2.8 FATAL_ERROR)
set(CMAKE_LEGACY_CYGWIN_WIN32 0)
set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
project("Fixed Order Generation" C CXX)
## Flags for the compiler. No warning allowed.
if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU")
set(warnings "-Wall -Wextra -Werror -std=c++1y")
elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Clang")
set(warnings "-Wall -Wextra -Werror -std=c++11")
elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "MSVC")
set(warnings "/W4 /WX /EHsc")
endif()
if (NOT CONFIGURED_ONCE)
set(CMAKE_CXX_FLAGS "${warnings}"
CACHE STRING "Flags used by the compiler during all build types." FORCE)
set(CMAKE_C_FLAGS "${warnings}"
CACHE STRING "Flags used by the compiler during all build types." FORCE)
endif()
## Add directories and find dependences
-include_directories(${CMAKE_CURRENT_SOURCE_DIR}/../include)
+include_directories(${CMAKE_CURRENT_SOURCE_DIR}/include ${CMAKE_CURRENT_SOURCE_DIR}/../include)
set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}/../cmake/Modules/")
#find_library(rhej rhej_)
find_package(fastjet REQUIRED)
find_package(clhep REQUIRED)
find_package(lhapdf REQUIRED)
find_package(gsl REQUIRED)
find_package(Boost REQUIRED COMPONENTS iostreams)
find_package(HepMC 3)
if(${HepMC_FOUND})
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DRHEJ_BUILD_WITH_HepMC")
endif()
find_package(QCDloop 2)
if(${QCDloop_FOUND})
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DRHEJ_BUILD_WITH_QCDLOOP")
include_directories(SYSTEM ${QCDloop_INCLUDE_PATH})
endif()
include_directories(SYSTEM ${lhapdf_INCLUDE_PATH})
include_directories(SYSTEM ${fastjet_INCLUDE_PATH})
include_directories(SYSTEM ${clhep_INCLUDE_PATH})
include_directories(SYSTEM ${gsl_INCLUDE_PATH})
#add_subdirectory(src)
## define executable
file(GLOB FOgen_source ${CMAKE_CURRENT_SOURCE_DIR}/src/*.cc)
-add_executable(FOgen src/main.cc ${FOgen_source})
+add_executable(FOgen ${FOgen_source})
## link libraries
set(libraries ${CMAKE_DL_LIBS} ${LHAPDF_LIBRARIES} ${CLHEP_LIBRARIES} ${FASTJET_LIBRARIES} ${GSL_LIBRARIES} ${Boost_LIBRARIES} ${HepMC_LIBRARIES} yaml-cpp)
if(${QCDloop_FOUND})
list(APPEND libraries ${QCDloop_LIBRARIES} quadmath)
endif()
# add libraries for reversed HEJ <by hand>
list(APPEND libraries "/home/andersen/HEJ/PURE/GitReverse/run/lib/librhej.so")
target_link_libraries(FOgen ${libraries})
file(GLOB FOgen_headers ${CMAKE_CURRENT_SOURCE_DIR}/include/*.hh)
install(TARGETS FOgen DESTINATION bin)
enable_testing()
set(tst_dir "${CMAKE_CURRENT_SOURCE_DIR}/t")
set(CONFIGURED_ONCE TRUE CACHE INTERNAL
"A flag showing that CMake has configured at least once.")
diff --git a/FixedOrderGen/include/PhaseSpacePoint.hh b/FixedOrderGen/include/PhaseSpacePoint.hh
new file mode 100644
index 0000000..74257c9
--- /dev/null
+++ b/FixedOrderGen/include/PhaseSpacePoint.hh
@@ -0,0 +1,92 @@
+/** \file PhaseSpacePoint.hh
+ * \brief Contains the PhaseSpacePoint Class
+ */
+
+#pragma once
+
+#include <vector>
+
+#include <CLHEP/Random/Randomize.h>
+#include <CLHEP/Random/RanluxEngine.h>
+
+#include "fastjet/JetDefinition.hh"
+#include "RHEJ/utility.hh"
+#include "RHEJ/Event.hh"
+
+namespace HEJFOG{
+ using RHEJ::Sparticle;
+ //! A point in resummation phase space
+ class PhaseSpacePoint{
+ public:
+ //! Default PhaseSpacePoint Constructor
+ PhaseSpacePoint() = default;
+
+ //! PhaseSpacePoint Constructor
+ /**
+ * @param ev Clustered Jet Event
+ * @param jet_def The Jet Definition Used
+ * @param jetptmin Minimum Jet Transverse Momentum
+ * @param extpartonptmin Minimum Parton Transverse Momentum
+ * @param max_ext_soft_pt Maximum Parton Transverse Momentum?
+ */
+ PhaseSpacePoint(
+ int nj, fastjet::JetDefinition jet_def,double jetptmin, double rapmax
+ );
+
+ //! Get Weight Function
+ /**
+ * @returns Weight of Event
+ */
+ double weight() const{
+ return weight_;
+ }
+
+
+ //! Get Incoming Function
+ /**
+ * @returns Incoming Particles
+ */
+ std::array<Sparticle, 2> const & incoming() const{
+ return incoming_;
+ }
+
+ //! Get Outgoing Function
+ /**
+ * @returns Outgoing Particles
+ */
+ std::vector<Sparticle> const & outgoing() const{
+ return outgoing_;
+ }
+
+ static void reset_ranlux(std::string const & init_file);
+ static void reset_ranlux(char const * init_file);
+
+ private:
+
+ void rescale_rapidities(
+ std::vector<fastjet::PseudoJet> & partons,
+ double ymin, double ymax
+ );
+ bool jets_ok(
+ std::vector<fastjet::PseudoJet> const & Born_jets,
+ std::vector<fastjet::PseudoJet> const & partons
+ ) const;
+ void reconstruct_incoming(std::array<Sparticle, 2> const & Born_incoming);
+ double phase_space_normalisation(
+ int num_Born_jets
+ ) const;
+
+ bool momentum_conserved(double ep) const;
+
+ double weight_;
+
+ double jetptmin_;
+ fastjet::JetDefinition jet_def_;
+
+ std::array<Sparticle, 2> incoming_;
+ std::vector<Sparticle> outgoing_;
+
+ static CLHEP::Ranlux64Engine ran_;
+ };
+ RHEJ::UnclusteredEvent to_UnclusteredEvent(PhaseSpacePoint const & psp);
+}
diff --git a/FixedOrderGen/src/PhaseSpacePoint.cc b/FixedOrderGen/src/PhaseSpacePoint.cc
new file mode 100644
index 0000000..25bbc97
--- /dev/null
+++ b/FixedOrderGen/src/PhaseSpacePoint.cc
@@ -0,0 +1,474 @@
+#include "PhaseSpacePoint.hh"
+
+#include <random>
+
+// #include "RHEJ/resummation_jet_momenta.hh"
+// #include "RHEJ/Jacobian.hh"
+// #include "RHEJ/RNGWrapper.hh"
+// #include "RHEJ/uno.hh"
+#include "RHEJ/debug.hh"
+#include "RHEJ/kinematics.hh"
+#include "RHEJ/utility.hh"
+
+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),
+ [](Sparticle o1, Sparticle 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.central.mur = NaN;
+ result.central.muf = NaN;
+ result.central.weight = psp.weight();
+ return result;
+ }
+
+
+
+namespace{
+ //generate Ranlux64Engine with fixed, predefined state
+ /*
+ * some (all?) of the Ranlux64Engine constructors leave fields
+ * uninitialised, invoking undefined behaviour. This can be
+ * circumvented by restoring the state from a file
+ */
+ CLHEP::Ranlux64Engine gen_Ranlux64Engine(){
+ static const std::string state =
+ "9876\n"
+ "0.91280703978419097666\n"
+ "0.41606065829518357191\n"
+ "0.99156342622341142601\n"
+ "0.030922955274050423213\n"
+ "0.16206278421638486975\n"
+ "0.76151768001958330956\n"
+ "0.43765760066092695979\n"
+ "0.42904698253748563275\n"
+ "0.11476317525663759511\n"
+ "0.026620053590963976831\n"
+ "0.65953715764414511114\n"
+ "0.30136722624439826745\n"
+ "3.5527136788005009294e-15 4\n"
+ "1 202\n";
+ const std::string file = std::tmpnam(nullptr);
+ {
+ std::ofstream out{file};
+ out << state;
+ }
+ CLHEP::Ranlux64Engine result;
+ result.restoreStatus(file.c_str());
+ return result;
+ }
+ }
+
+ CLHEP::Ranlux64Engine PhaseSpacePoint::ran_{gen_Ranlux64Engine()};
+
+
+ void PhaseSpacePoint::reset_ranlux(std::string const & init_file){
+ reset_ranlux(init_file.c_str());
+ }
+
+ void PhaseSpacePoint::reset_ranlux(char const * init_file){
+ ran_.restoreStatus(init_file);
+ }
+
+ namespace {
+ // constexpr int max_jet_user_idx = PhaseSpacePoint::ng_max;
+
+ // bool is_nonjet_parton(fastjet::PseudoJet const & parton){
+ // assert(parton.user_index() != -1);
+ // return parton.user_index() > max_jet_user_idx;
+ // }
+
+ // bool is_jet_parton(fastjet::PseudoJet const & parton){
+ // assert(parton.user_index() != -1);
+ // return parton.user_index() <= max_jet_user_idx;
+ // }
+
+ // user indices for partons with extremal rapidity
+ constexpr int y_min_user_idx = -3;
+ constexpr int y_max_user_idx = -2;
+ }
+
+ // namespace{
+ // double estimate_ng_mean(std::vector<fastjet::PseudoJet> const & Born_jets){
+ // const double delta_y =
+ // Born_jets.back().rapidity() - Born_jets.front().rapidity();
+ // assert(delta_y > 0);
+ // // Formula derived from fit to 2 jet data
+ // // for 3 jets: -0.0131111 + (1.39385 + 0.050085*delta_y)*delta_y
+ // return -0.0213569 + (1.39765 + 0.0498387*delta_y)*delta_y;
+ // }
+ // }
+
+ // std::vector<fastjet::PseudoJet> PhaseSpacePoint::cluster_jets(
+ // std::vector<fastjet::PseudoJet> const & partons
+ // ) const{
+ // fastjet::ClusterSequence cs(partons, jet_def_);
+ // return cs.inclusive_jets(jetptmin_);
+ // }
+
+ // bool PhaseSpacePoint::pass_resummation_cuts(
+ // std::vector<fastjet::PseudoJet> const & jets
+ // ) const{
+ // return cluster_jets(jets).size() == jets.size();
+ // }
+
+ // int PhaseSpacePoint::sample_ng(std::vector<fastjet::PseudoJet> const & Born_jets){
+ // const double ng_mean = estimate_ng_mean(Born_jets);
+ // std::poisson_distribution<int> dist(ng_mean);
+ // RNGWrapper<CLHEP::Ranlux64Engine> rng{ran_};
+ // const int ng = dist(rng);
+ // assert(ng >= 0);
+ // assert(ng < ng_max);
+ // weight_ *= std::tgamma(ng + 1)*std::exp(ng_mean)*std::pow(ng_mean, -ng);
+ // return ng;
+ // }
+
+ // namespace {
+ // void copy_AWZH_boson(
+ // std::vector<Sparticle> const & from,
+ // std::vector<Sparticle> & to
+ // ){
+ // const auto AWZH_boson = std::find_if(
+ // begin(from), end(from),
+ // [](Sparticle const & p){ return is_AWZH_boson(p); }
+ // );
+ // if(AWZH_boson == end(from)) return;
+ // auto insertion_point = std::lower_bound(
+ // begin(to), end(to), *AWZH_boson, rapidity_less{}
+ // );
+ // to.insert(insertion_point, *AWZH_boson);
+ // assert(std::is_sorted(begin(to), end(to), rapidity_less{}));
+ // }
+ // }
+
+ PhaseSpacePoint::PhaseSpacePoint(
+ int nj, fastjet::JetDefinition jet_def,double jetptmin, double rapmax
+ )
+ {
+ weight_ = 1;
+ weight_ /= std::tgamma(nj + 1);
+ {
+
+ outgoing_.reserve(nj + 1); // one slot for possible A, W, Z, H
+
+ // use a few variables to check for compilation
+ {
+ fastjet::JetDefinition jet_def2(jet_def);
+ std::cout << " nj: " <<nj<<" jetptmin: "<<jetptmin<<" rapmax : "<<rapmax<<std::endl;
+
+ // if(&jet_def==NULL)
+
+ }
+ // for(auto & p: out_partons){
+ // outgoing_.emplace_back(Sparticle{pid::gluon, std::move(p)});
+ // }
+ // most_backward_FKL(outgoing_).type = ev.incoming().front().type;
+ // most_forward_FKL(outgoing_).type = ev.incoming().back().type;
+ // copy_AWZH_boson(ev.outgoing(), outgoing_);
+
+ // assert(!outgoing_.empty());
+ // reconstruct_incoming(ev.incoming());
+
+ }
+ }
+
+
+ void PhaseSpacePoint::rescale_rapidities(
+ std::vector<fastjet::PseudoJet> & partons,
+ double ymin, double ymax
+ ){
+ constexpr double ep = 1e-7;
+ for(auto & parton: partons){
+ assert(0 <= parton.rapidity() && parton.rapidity() <= 1);
+ const double dy = ymax - ymin - 2*ep;
+ const double y = ymin + ep + dy*parton.rapidity();
+ parton.reset_momentum_PtYPhiM(parton.pt(), y, parton.phi());
+ weight_ *= dy;
+ assert(ymin <= parton.rapidity() && parton.rapidity() <= ymax);
+ }
+ }
+
+ // int PhaseSpacePoint::sample_ng_jets(
+ // int ng, std::vector<fastjet::PseudoJet> const & Born_jets
+ // ){
+ // RNGWrapper<CLHEP::Ranlux64Engine> rng{ran_};
+ // const double p_J = probability_in_jet(Born_jets);
+ // std::binomial_distribution<> bin_dist(ng, p_J);
+ // const int ng_J = bin_dist(rng);
+ // weight_ *= std::pow(p_J, -ng_J)*std::pow(1 - p_J, ng_J - ng);
+ // return ng_J;
+ // }
+
+ // std::vector<fastjet::PseudoJet>
+ // PhaseSpacePoint::reshuffle(
+ // std::vector<fastjet::PseudoJet> const & Born_jets,
+ // fastjet::PseudoJet const & q
+ // ){
+ // if(q == fastjet::PseudoJet{0, 0, 0, 0}) return Born_jets;
+ // std::vector<fastjet::PseudoJet> jets = resummation_jet_momenta(Born_jets, q);
+ // if(jets.empty()){
+ // weight_ = 0;
+ // return {};
+ // }
+ // // transform delta functions to integration over resummation momenta
+ // weight_ /= Jacobian(jets, q);
+ // return jets;
+ // }
+
+ // std::vector<int> PhaseSpacePoint::distribute_jet_partons(
+ // int ng_jets, std::vector<fastjet::PseudoJet> const & jets
+ // ){
+ // size_t first_valid_jet = 0;
+ // size_t num_valid_jets = jets.size();
+ // const double R_eff = 5./3.*jet_def_.R();
+ // // if there is an unordered jet too far away from the FKL jets
+ // // then extra gluon constituents of the unordered jet would
+ // // violate the FKL rapidity ordering
+ // if(unob_ && jets[0].delta_R(jets[1]) > R_eff){
+ // ++first_valid_jet;
+ // --num_valid_jets;
+ // }
+ // else if(unof_ && jets[jets.size()-1].delta_R(jets[jets.size()-2]) > R_eff){
+ // --num_valid_jets;
+ // }
+ // std::vector<int> np(jets.size(), 1);
+ // for(int i = 0; i < ng_jets; ++i){
+ // ++np[first_valid_jet + ran_.flat() * num_valid_jets];
+ // }
+ // weight_ *= std::pow(num_valid_jets, ng_jets);
+ // return np;
+ // }
+
+
+ // std::vector<fastjet::PseudoJet> PhaseSpacePoint::split(
+ // std::vector<fastjet::PseudoJet> const & jets,
+ // int ng_jets
+ // ){
+ // return split(jets, distribute_jet_partons(ng_jets, jets));
+ // }
+
+ // bool PhaseSpacePoint::pass_extremal_cuts(
+ // fastjet::PseudoJet const & ext_parton,
+ // fastjet::PseudoJet const & jet
+ // ) const{
+ // if(ext_parton.pt() < extpartonptmin_) return false;
+ // return (ext_parton - jet).pt()/ext_parton.pt() < max_ext_soft_pt_fraction_;
+ // }
+
+ // std::vector<fastjet::PseudoJet> PhaseSpacePoint::split(
+ // std::vector<fastjet::PseudoJet> const & jets,
+ // std::vector<int> const & np
+ // ){
+ // assert(! jets.empty());
+ // assert(jets.size() == np.size());
+ // assert(pass_resummation_cuts(jets));
+
+ // const size_t most_backward_FKL_idx = 0 + unob_;
+ // const size_t most_forward_FKL_idx = jets.size() - 1 - unof_;
+
+ // std::vector<fastjet::PseudoJet> jet_partons;
+ // // randomly distribute jet gluons among jets
+ // for(size_t i = 0; i < jets.size(); ++i){
+ // weight_ *= splitter_.Split(jets[i], np[i]);
+ // if(weight_ == 0) return {};
+ // assert(
+ // std::all_of(
+ // begin(splitter_.get_jcons()), end(splitter_.get_jcons()),
+ // is_jet_parton
+ // )
+ // );
+ // const auto first_new_parton = jet_partons.insert(
+ // end(jet_partons),
+ // begin(splitter_.get_jcons()), end(splitter_.get_jcons())
+ // );
+ // auto extremal = end(jet_partons);
+ // if((unob_ && i == 0) || i == most_backward_FKL_idx){
+ // // unordered or FKL backward emission
+ // extremal = std::min_element(
+ // first_new_parton, end(jet_partons), rapidity_less{}
+ // );
+ // if(i == most_backward_FKL_idx) extremal->set_user_index(y_min_user_idx);
+ // }
+ // else if((unof_ && i == jets.size() - 1) || i == most_forward_FKL_idx){
+ // // unordered or FKL forward emission
+ // extremal = std::max_element(
+ // first_new_parton, end(jet_partons), rapidity_less{}
+ // );
+ // if(i == most_forward_FKL_idx) extremal->set_user_index(y_max_user_idx);
+ // }
+ // if(
+ // extremal != end(jet_partons)
+ // && !pass_extremal_cuts(*extremal, jets[i])
+ // ){
+ // weight_ = 0;
+ // return {};
+ // }
+ // }
+ // assert(tagged_FKL_extremal(jet_partons));
+ // std::sort(begin(jet_partons), end(jet_partons), rapidity_less{});
+ // if(
+ // !extremal_FKL_ok(jet_partons)
+ // || !split_preserved_jets(jets, jet_partons)
+ // ){
+ // weight_ = 0.;
+ // return {};
+ // }
+ // return jet_partons;
+ // }
+
+ // bool PhaseSpacePoint::extremal_FKL_ok(
+ // std::vector<fastjet::PseudoJet> const & partons
+ // ) const{
+ // assert(std::is_sorted(begin(partons), end(partons), rapidity_less{}));
+ // return
+ // most_backward_FKL(partons).user_index() == y_min_user_idx
+ // && most_forward_FKL(partons).user_index() == y_max_user_idx;
+ // }
+
+ // bool PhaseSpacePoint::split_preserved_jets(
+ // std::vector<fastjet::PseudoJet> const & jets,
+ // std::vector<fastjet::PseudoJet> const & jet_partons
+ // ) const{
+ // assert(std::is_sorted(begin(jets), end(jets), rapidity_less{}));
+
+ // const auto split_jets = sorted_by_rapidity(cluster_jets(jet_partons));
+ // // this can happen if two overlapping jets
+ // // are both split into more than one parton
+ // if(split_jets.size() != jets.size()) return false;
+ // for(size_t i = 0; i < split_jets.size(); ++i){
+ // // this can happen if there are two overlapping jets
+ // // and a parton is assigned to the "wrong" jet
+ // if(!nearby_ep(jets[i].rapidity(), split_jets[i].rapidity(), 1e-2)){
+ // return false;
+ // }
+ // }
+ // return true;
+ // }
+
+ // template<class Particle>
+ // Particle const & PhaseSpacePoint::most_backward_FKL(
+ // std::vector<Particle> const & partons
+ // ) const{
+ // return partons[0 + unob_];
+ // }
+
+ // template<class Particle>
+ // Particle const & PhaseSpacePoint::most_forward_FKL(
+ // std::vector<Particle> const & partons
+ // ) const{
+ // const size_t idx = partons.size() - 1 - unof_;
+ // assert(idx < partons.size());
+ // return partons[idx];
+ // }
+
+ // template<class Particle>
+ // Particle & PhaseSpacePoint::most_backward_FKL(
+ // std::vector<Particle> & partons
+ // ) const{
+ // return partons[0 + unob_];
+ // }
+
+ // template<class Particle>
+ // Particle & PhaseSpacePoint::most_forward_FKL(
+ // std::vector<Particle> & partons
+ // ) const{
+ // const size_t idx = partons.size() - 1 - unof_;
+ // assert(idx < partons.size());
+ // return partons[idx];
+ // }
+
+ // namespace{
+ // bool contains_idx(
+ // fastjet::PseudoJet const & jet, fastjet::PseudoJet const & parton
+ // ){
+ // auto const & constituents = jet.constituents();
+ // const int idx = parton.user_index();
+ // return std::find_if(
+ // begin(constituents), end(constituents),
+ // [idx](fastjet::PseudoJet const & con){return con.user_index() == idx;}
+ // ) != end(constituents);
+ // }
+ // }
+
+ // /**
+ // * final jet test:
+ // * - number of jets must match Born kinematics
+ // * - no partons designated as nonjet may end up inside jets
+ // * - all other outgoing partons *must* end up inside jets
+ // * - the extremal (in rapidity) partons must be inside the extremal jets
+ // * - rapidities must be the same (by construction)
+ // */
+ // bool PhaseSpacePoint::jets_ok(
+ // std::vector<fastjet::PseudoJet> const & Born_jets,
+ // std::vector<fastjet::PseudoJet> const & partons
+ // ) const{
+ // fastjet::ClusterSequence cs(partons, jet_def_);
+ // const auto jets = sorted_by_rapidity(cs.inclusive_jets(jetptmin_));
+ // if(jets.size() != Born_jets.size()) return false;
+ // int in_jet = 0;
+ // for(size_t i = 0; i < jets.size(); ++i){
+ // assert(jets[i].has_constituents());
+ // for(auto && parton: jets[i].constituents()){
+ // if(is_nonjet_parton(parton)) return false;
+ // }
+ // in_jet += jets[i].constituents().size();
+ // }
+ // const int expect_in_jet = std::count_if(
+ // partons.cbegin(), partons.cend(), is_jet_parton
+ // );
+ // if(in_jet != expect_in_jet) return false;
+ // // note that PseudoJet::contains does not work here
+ // if(! (
+ // contains_idx(most_backward_FKL(jets), most_backward_FKL(partons))
+ // && contains_idx(most_forward_FKL(jets), most_forward_FKL(partons))
+ // )) return false;
+ // for(size_t i = 0; i < jets.size(); ++i){
+ // assert(nearby_ep(jets[i].rapidity(), Born_jets[i].rapidity(), 1e-2));
+ // }
+ // return true;
+ // }
+
+ namespace{
+ }
+
+ void PhaseSpacePoint::reconstruct_incoming(
+ std::array<Sparticle, 2> const & Born_incoming
+ ){
+ std::tie(incoming_[0].p, incoming_[1].p) = incoming_momenta(outgoing_);
+ for(size_t i = 0; i < incoming_.size(); ++i){
+ incoming_[i].type = Born_incoming[i].type;
+ }
+ assert(momentum_conserved(1e-7));
+ }
+
+ double PhaseSpacePoint::phase_space_normalisation(
+ int num_Born_jets
+ ) const{
+ return pow(16*pow(M_PI,3), num_Born_jets);
+ }
+
+ 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);
+ }
+
+}
diff --git a/FixedOrderGen/src/main.cc b/FixedOrderGen/src/main.cc
index 3fa280d..60203c1 100644
--- a/FixedOrderGen/src/main.cc
+++ b/FixedOrderGen/src/main.cc
@@ -1,113 +1,124 @@
/**
* Name: main.cc
* Authors: Jeppe R. Andersen
*/
#include <fstream>
#include <algorithm>
#include <memory>
#include <chrono>
#include <iostream>
#include "yaml-cpp/yaml.h"
#include "LHEF/LHEF.h"
#include "RHEJ/CombinedEventWriter.hh"
#include "RHEJ/get_analysis.hh"
#include "RHEJ/utility.hh"
-#include "RHEJ/EventReweighter.hh"
+//#include "RHEJ/EventReweighter.hh"
#include "RHEJ/config.hh"
#include "RHEJ/stream.hh"
+#include "RHEJ/MatrixElement.hh"
+#include "PhaseSpacePoint.hh"
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);
}
}
int main(int argn, char** argv) {
std::cout << " __ ___ __ ______ __ __ \n";
std::cout << " / / / (_)___ _/ /_ / ____/___ ___ _________ ___ __ / /__ / /______ \n";
std::cout << " / /_/ / / __ `/ __ \\ / __/ / __ \\/ _ \\/ ___/ __ `/ / / / __ / / _ \\/ __/ ___/ \n";
std::cout << " / __ / / /_/ / / / / / /___/ / / / __/ / / /_/ / /_/ / / /_/ / __/ /_(__ ) \n";
std::cout << " /_/ /_/_/\\__, /_/ /_/ /_____/_/ /_/\\___/_/ \\__, /\\__, / \\____/\\___/\\__/____/ \n";
std::cout << " ____///__/ __ ____ ///__//____/ ______ __ \n";
std::cout << " / ____(_) _____ ____/ / / __ \\_________/ /__ _____ / ____/__ ____ ___ _________ _/ /_____ _____\n";
std::cout << " / /_ / / |/_/ _ \\/ __ / / / / / ___/ __ / _ \\/ ___/ / / __/ _ \\/ __ \\/ _ \\/ ___/ __ `/ __/ __ \\/ ___/\n";
std::cout << " / __/ / /> </ __/ /_/ / / /_/ / / / /_/ / __/ / / /_/ / __/ / / / __/ / / /_/ / /_/ /_/ / / \n";
std::cout << " /_/ /_/_/|_|\\___/\\__,_/ \\____/_/ \\__,_/\\___/_/ \\____/\\___/_/ /_/\\___/_/ \\__,_/\\__/\\____/_/ \n";
using clock = std::chrono::system_clock;
if (argn < 2) {
std::cerr << "\n# Usage:\n.FOgen config_file\n";
return EXIT_FAILURE;
}
const auto start_time = clock::now();
// 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);
+ RHEJ::MatrixElement ME(config.fixed_order_jets.def,config.fixed_order_jets.min_pt,config.log_correction,config.Higgs_coupling);
+
// Perform N trial generations
int nevent = 0;
while(nevent< config.trials){
- // // reweight events so that the total cross section is conserved
- // reader.hepeup.setWeight(0, global_reweight * reader.hepeup.weight());
+ if (int(10000.*nevent/config.trials) % 100 == 0)
+ std::cout << ".";
+ std::cout.flush();
+
+ // Generate phase space point
+ HEJFOG::PhaseSpacePoint psp{2,config.fixed_order_jets.def,config.fixed_order_jets.min_pt,5.0};
+ std::cout<<psp.weight()<<std::endl;
+ RHEJ::Event ev{to_UnclusteredEvent(std::move(psp)), config.fixed_order_jets.def, config.fixed_order_jets.min_pt};
+
+
+ // evaluate matrix element on this point
+ ME.tree(pdf_.Halphas(psp.central().mur), psp.central().mur,
+ psp.incoming(), psp.outgoing(),false);
+
+
+
- // if(nevent == max_events) break;
- // nevent++;
-
- // if (nevent % 10000 == 0)
- // std::cout << "event number " << nevent << std::endl;
-
- // // calculate rHEJ weight
// RHEJ::Event FO_event{
// RHEJ::UnclusteredEvent{reader.hepeup},
- // config.fixed_order_jets.def, config.fixed_order_jets.min_pt,
- // };
+ // config.fixed_order_jets.def, config.fixed_order_jets.min_pt,
+ // };
// auto resummed_events = rhej.reweight(
// FO_event,
// config.trials, config.scale_gen
// );
// const auto ev_class = rhej.last_event_class();
// ++nevent_class[ev_class];
// if(resummed_events.empty()) ++nfailed_class[ev_class];
// for(auto const & ev: resummed_events){
// analysis->fill(ev);
// writer.write(ev);
// }
- ++nevent;
+ nevent++;
} // main event loop
-
+ std::cout << std::endl;
std::chrono::duration<double> run_time = (clock::now() - start_time);
std::cout << "\nTask Runtime: " << run_time.count() << " seconds.\n";
return 0;
}

File Metadata

Mime Type
text/x-diff
Expires
Tue, Jan 21, 2:17 AM (1 d, 18 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4229916
Default Alt Text
(26 KB)

Event Timeline