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/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";
}