Page MenuHomeHEPForge
Files F8723680

No OneTemporary
Actions

View Options

diff --git a/FixedOrderGen/CMakeLists.txt b/FixedOrderGen/CMakeLists.txt
index 29fd80d..9c8030a 100644
--- a/FixedOrderGen/CMakeLists.txt
+++ b/FixedOrderGen/CMakeLists.txt
@@ -1,68 +1,68 @@
cmake_minimum_required(VERSION 3.1 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_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 11)
+set(CMAKE_CXX_STANDARD 14)
## Add directories and find dependences
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_package(fastjet REQUIRED)
include_directories(${fastjet_INCLUDE_PATH})
find_package(clhep 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 REQUIRED)
include_directories(${YAML_CPP_INCLUDE_DIR})
find_package(HepMC 2)
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} ${GSL_LIBRARIES} ${Boost_LIBRARIES} ${YAML_CPP_LIBRARIES})
if(${QCDloop_FOUND})
list(APPEND libraries ${QCDloop_LIBRARIES} quadmath)
endif()
if(${HepMC_FOUND})
list(APPEND libraries ${HepMC_LIBRARIES})
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)
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/CrossSectionAccumulator.hh b/FixedOrderGen/include/CrossSectionAccumulator.hh
new file mode 100644
index 0000000..b7e9e86
--- /dev/null
+++ b/FixedOrderGen/include/CrossSectionAccumulator.hh
@@ -0,0 +1,43 @@
+#pragma once
+
+#include <map>
+
+#include "RHEJ/Event.hh"
+#include "RHEJ/event_types.hh"
+
+namespace HEJFOG {
+ template<typename T>
+ struct WithError {
+ T value = T{};
+ T error = T{};
+ };
+
+ class CrossSectionAccumulator {
+ public:
+ void fill(RHEJ::Event const & ev) {
+ const double wt = ev.central().weight;
+ auto & entry = xs_[ev.type()];
+ entry.value += wt;
+ entry.error += wt*wt;
+ total_.value += wt;
+ total_.error += wt*wt;
+ }
+
+ auto begin() const {
+ return std::begin(xs_);
+ }
+
+ auto end() const {
+ return std::end(xs_);
+ }
+
+ WithError<double> total() const {
+ return total_;
+ }
+
+ private:
+ std::map<RHEJ::event_type::EventType, WithError<double>> xs_;
+ WithError<double> total_;
+ };
+
+}
diff --git a/FixedOrderGen/src/main.cc b/FixedOrderGen/src/main.cc
index ed3c018..0ae7337 100644
--- a/FixedOrderGen/src/main.cc
+++ b/FixedOrderGen/src/main.cc
@@ -1,243 +1,233 @@
/**
* 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 "EventGenerator.hh"
#include "PhaseSpacePoint.hh"
#include "Ranlux64Engine.hh"
#include "ProgressBar.hh"
#include "Unweighter.hh"
+#include "CrossSectionAccumulator.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_test_trials = 100000;
}
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;
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 = HEJFOG::make_Ranlux64Engine();
if(config.RanLux_init) ran.restoreStatus(config.RanLux_init->c_str());
HEJFOG::EventGenerator generator{
config.process,
config.beam,
config.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};
std::map<HEJFOG::Status, int> status_counter;
const int test_trials = std::min(config.events, max_test_trials);
- double FKL_xs = 0., FKL_xs_err = 0.;
- double uno_xs = 0., uno_xs_err = 0.;
std::vector<RHEJ::Event> events;
// warm-up phase: check how efficient we are
std::cout << "Warm-up ...\n";
HEJFOG::ProgressBar<int> warmup_progress{std::cout, test_trials};
const auto warmup_start = clock::now();
for(int trial = 0; trial < test_trials; ++trial){
++warmup_progress;
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));
}
}
std::cout << std::endl;
RHEJ::optional<HEJFOG::Unweighter> unweighter{};
if(config.unweight) {
unweighter = HEJFOG::Unweighter{
begin(events), end(events), config.unweight->max_dev, ran
};
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);
}
const double efficiency = static_cast<double>(events.size())/test_trials;
assert(0. <= efficiency && efficiency <= 1.);
std::cout << "Efficiency: " << efficiency << '\n';
constexpr auto max_trials = std::numeric_limits<long long>::max();
if(max_trials*efficiency < test_trials) {
std::cerr << "Number of required configurations exceeds " << max_trials
<< "\nPlease reduce the number of requested events.\n";
std::cout.setstate(std::ios_base::badbit); // prevent LHAPDF output
return EXIT_FAILURE;
}
const long long total_trials = test_trials/efficiency;
+ HEJFOG::CrossSectionAccumulator xs;
for(auto & ev: events) {
rescale_weights(ev, invGeV2_to_pb/total_trials);
analysis->fill(ev, ev);
writer.write(ev);
- if(ev.type() == RHEJ::event_type::FKL){
- FKL_xs += ev.central().weight;
- FKL_xs_err += ev.central().weight*ev.central().weight;
- }
- else {
- uno_xs += ev.central().weight;
- uno_xs_err += ev.central().weight*ev.central().weight;
- }
+ xs.fill(ev);
}
const auto warmup_end = clock::now();
const std::chrono::duration<double> warmup_time = warmup_end - warmup_start;
std::cout
<< "Generating " << (total_trials - test_trials) << " additional configurations\n"
"Estimated remaining time: "
<< warmup_time.count()*(1./efficiency - 1) << " seconds.\n\n";
HEJFOG::ProgressBar<long long> progress{std::cout, total_trials};
progress.increment(test_trials);
for(int trial = test_trials - 1; trial < total_trials; ++trial){
++progress;
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);
}
rescale_weights(ev, invGeV2_to_pb/total_trials);
analysis->fill(ev, ev);
writer.write(ev);
-
- if(ev.type() == RHEJ::event_type::FKL){
- FKL_xs += ev.central().weight;
- FKL_xs_err += ev.central().weight*ev.central().weight;
- }
- else {
- uno_xs += ev.central().weight;
- uno_xs_err += ev.central().weight*ev.central().weight;
- }
+ xs.fill(ev);
}
}
std::cout << std::endl;
const std::chrono::duration<double> run_time = (clock::now() - start_time);
- std::cout << "\nTask Runtime: " << run_time.count() << " seconds.\n";
+ std::cout << "\nTask Runtime: " << run_time.count() << " seconds.\n\n";
std::cout.precision(10);
- std::cout.width(12);
std::cout.setf(std::ios::fixed);
std::cout
- << "\nCross section: " << (FKL_xs + uno_xs)
- << " +- " << std::sqrt(FKL_xs_err + uno_xs_err) << " (pb)"
- "\n FKL: " << FKL_xs << " +- " << std::sqrt(FKL_xs_err) << " (pb)"
- "\n unordered: " << uno_xs << " +- " << std::sqrt(uno_xs_err) << " (pb)"
- "\n";
+ << " " << 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)/total_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";
}
}

File Metadata

Mime Type
text/x-diff
Expires
Mon, Jan 20, 9:11 PM (1 d, 1 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4242394
Default Alt Text
(12 KB)

Event Timeline