No OneTemporary
Actions

Size

16 KB

Subscribers

None

View Options

	diff --git a/include/RHEJ/Version.hh b/include/RHEJ/Version.hh
	new file mode 100644
	index 0000000..95da3e7
	--- /dev/null
	+++ b/include/RHEJ/Version.hh
	@@ -0,0 +1,49 @@
	+/** \file Version.hh
	+ * \brief The file gives the current rHEJ Version
	+ */
	+
	+#pragma once
	+
	+#include <string>
	+
	+/// @brief Full name of this package.
	+#define RHEJ_PACKAGE_NAME "rHEJ"
	+
	+/// @brief rHEJ version string
	+#define RHEJ_VERSION "0.0.1"
	+
	+/// @brief Full name and version of this package.
	+#define RHEJ_PACKAGE_STRING "rHEJ 0.0.1"
	+
	+/// @brief rHEJ version as an integer, rHEJ XXYYZZ = XXYYZZ
	+///
	+/// Use like "#if RHEJ_VERSION < 010203" for < 1.2.3
	+#define RHEJ_VERSION_CODE 000001
	+
	+/// @brief Major version of this package
	+#define RHEJ_VERSION_MAJOR 0
	+
	+/// @brief Minor version of this package
	+#define RHEJ_VERSION_MINOR 0
	+
	+/// @brief Patch version of this package
	+#define RHEJ_VERSION_PATCHLEVEL 1
	+
	+
	+namespace RHEJ {
	+
	+ /// Get the rHEJ library version string
	+ inline std::string version() {
	+ return RHEJ_VERSION;
	+ }
	+
	+ /// Get the rHEJ library package name
	+ inline std::string package_name() {
	+ return RHEJ_PACKAGE_NAME;
	+ }
	+
	+ // Get the full rHEJ library name and version
	+ inline std::string package_name_full() {
	+ return RHEJ_PACKAGE_STRING;
	+ }
	+}
	diff --git a/src/HepMCWriter.cc b/src/HepMCWriter.cc
	index 82b2957..b1eb236 100644
	--- a/src/HepMCWriter.cc
	+++ b/src/HepMCWriter.cc
	@@ -1,241 +1,242 @@
	#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"
	#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/Version.hh"

	namespace RHEJ{

	#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 = "HEJ";
	- heprup.generators.back().version = "0.0.1";
	+ heprup.generators.back().name = package_name();
	+ heprup.generators.back().version = version();
	}
	void reset_weight_info(LHEF::HEPRUP & heprup){
	heprup.IDWTUP = -4;
	// 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::FourVector to_FourVector(Sparticle const & sp){
	return {sp.px(), sp.py(), sp.pz(), sp.E()};
	}

	HepMC::shared_ptr<HepMC::GenRunInfo> & add_HEPRUP(
	HepMC::shared_ptr<HepMC::GenRunInfo> & runinfo,
	LHEF::HEPRUP heprup
	){
	add_generator_tag(heprup);
	reset_weight_info(heprup);

	auto hepr = HepMC::make_shared<HepMC::HEPRUPAttribute>();
	hepr->heprup = std::move(heprup);
	runinfo->add_attribute("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 runinfo;
	}
	} // namespace anonymous

	struct HepMCWriter::HepMCWriterImpl: public EventWriter{

	HepMC::shared_ptr<HepMC::GenRunInfo> runinfo_;
	HepMC::WriterAscii writer_;

	HepMCWriterImpl(
	std::string const & file, Config const & conf, LHEF::HEPRUP && heprup
	):
	EventWriter(conf),
	runinfo_{HepMC::make_shared<HepMC::GenRunInfo>()},
	writer_{file, add_HEPRUP(runinfo_, std::move(heprup))}
	{}

	HepMCWriterImpl & operator=(HepMCWriterImpl const & other) = delete;
	HepMCWriterImpl(HepMCWriterImpl const & other) = delete;
	HepMCWriterImpl & operator=(HepMCWriterImpl && other) = delete;
	HepMCWriterImpl(HepMCWriterImpl && other) = delete;

	~HepMCWriterImpl(){
	writer_.close();
	}

	void write(Event const & ev){
	if(!ev.decays().empty()){
	throw not_implemented{"HepMC output for events with decays"};
	}
	auto vx = HepMC::make_shared<HepMC::GenVertex>();
	for(auto const & in: ev.incoming()){
	vx->add_particle_in(
	HepMC::make_shared<HepMC::GenParticle>(
	to_FourVector(in), static_cast<int>(in.type), -1
	)
	);
	}
	for(auto const & out: ev.outgoing()){
	vx->add_particle_out(
	HepMC::make_shared<HepMC::GenParticle>(
	to_FourVector(out), static_cast<int>(out.type), +1
	)
	);
	}
	HepMC::GenEvent out_ev{runinfo_, HepMC::Units::GEV, HepMC::Units::MM};
	out_ev.add_vertex(vx);
	out_ev.weights().reserve(ev.variations().size() + 1u);
	out_ev.weights().emplace_back(ev.central().weight);
	for(auto const & var: ev.variations()){
	out_ev.weights().emplace_back(var.weight);
	}
	writer_.write_event(out_ev);
	}
	};

	void HepMCWriter::HepMCWriterImplDeleter::operator()(HepMCWriterImpl* p) {
	delete p;
	}

	#else // HepMC 2

	namespace{
	HepMC::FourVector to_FourVector(Sparticle const & sp){
	return {sp.px(), sp.py(), sp.pz(), sp.E()};
	}
	}

	struct HepMCWriter::HepMCWriterImpl: public EventWriter{

	struct HepMCInfo{
	size_t event_count = 0.;
	double tot_weight = 0.;
	double tot_weight2 = 0.;

	void add_event(Event const & ev){
	double wt = ev.central().weight;
	tot_weight += wt;
	tot_weight2 += wt * wt;
	++event_count;
	}
	};

	HepMC::IO_GenEvent writer_;
	HepMCInfo genInfo_;

	HepMCWriterImpl(
	std::string const & file, Config const & conf, LHEF::HEPRUP &&
	):
	EventWriter(conf),
	writer_{file}
	{}

	HepMCWriterImpl & operator=(HepMCWriterImpl const & other) = delete;
	HepMCWriterImpl(HepMCWriterImpl const & other) = delete;
	HepMCWriterImpl & operator=(HepMCWriterImpl && other) = delete;
	HepMCWriterImpl(HepMCWriterImpl && other) = delete;

	void write(Event const & ev){
	if(!ev.decays().empty()){
	throw not_implemented{"HepMC output for events with decays"};
	}
	auto vx = new HepMC::GenVertex();
	for(auto const & in: ev.incoming()){
	vx->add_particle_in(
	new HepMC::GenParticle{
	to_FourVector(in), static_cast<int>(in.type), -1
	}
	);
	}
	for(auto const & out: ev.outgoing()){
	vx->add_particle_out(
	new HepMC::GenParticle{
	to_FourVector(out), static_cast<int>(out.type), +1
	}
	);
	}
	HepMC::GenEvent out_ev{HepMC::Units::GEV, HepMC::Units::MM};
	out_ev.add_vertex(vx);
	/// weights
	out_ev.weights().push_back(ev.central().weight);
	for(auto const & var: ev.variations()){
	out_ev.weights().push_back(var.weight);
	}
	/// @TODO add name list
	/// general informations
	genInfo_.add_event(ev);
	out_ev.set_signal_process_id(ev.type()+1); /// "+1": conistent with lhe
	out_ev.set_event_scale(ev.central().mur); /// event scale
	out_ev.set_event_number(genInfo_.event_count); /// event number
	/// cross section
	HepMC::GenCrossSection t_xs;
	t_xs.set_cross_section( genInfo_.tot_weight , sqrt(genInfo_.tot_weight2) );
	out_ev.set_cross_section( t_xs );

	/// @TODO add alphaQCD (need function) and alphaQED
	/// @TODO output pdf (currently not avaiable from event alone)
	writer_.write_event(&out_ev);
	}
	};

	void HepMCWriter::HepMCWriterImplDeleter::operator()(HepMCWriterImpl* p) {
	delete p;
	}

	#endif // HepMC 2

	HepMCWriter::HepMCWriter(std::string const & file, Config const & conf, LHEF::HEPRUP heprup):
	EventWriter(conf),
	impl_{std::unique_ptr<HepMCWriterImpl, HepMCWriterImplDeleter>{
	new HepMCWriterImpl(file, conf, 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 &, Config const & conf, LHEF::HEPRUP):
	EventWriter(conf){
	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 3c2a4b4..8532b41 100644
	--- a/src/LesHouchesWriter.cc
	+++ b/src/LesHouchesWriter.cc
	@@ -1,90 +1,91 @@
	#include <stdexcept>
	#include <memory>
	#include <cassert>

	#include "RHEJ/LesHouchesWriter.hh"
	#include "RHEJ/event_types.hh"

	#include "RHEJ/Event.hh"
	+#include "RHEJ/Version.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, Config const & conf, LHEF::HEPRUP heprup
	):
	EventWriter(conf),
	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 = +/-)
	// 3: weight=+/-, xs given in head (same as default MG)
	// 4: weight=+/-, xs = avg(weights)
	writer_->heprup.IDWTUP = -3;
	writer_->heprup.generators.emplace_back(LHEF::XMLTag{});
	- writer_->heprup.generators.back().name = "HEJ";
	- writer_->heprup.generators.back().version = "0.0.1";
	+ writer_->heprup.generators.back().name = package_name();
	+ writer_->heprup.generators.back().version = version();
	// 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/main.cc b/src/main.cc
	index 5bc57b3..050924c 100644
	--- a/src/main.cc
	+++ b/src/main.cc
	@@ -1,142 +1,162 @@
	/**
	* Name: main.cc
	* Authors: Tuomas Hapola, Andreas Maier <andreas.maier@durham.ac.uk>
	*
	*/
	#include <fstream>
	#include <algorithm>
	#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"

	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);
	}
	}

	+std::string print_time(const time_t time){
	+ char s[1000];
	+ struct tm * p = localtime(&time);
	+ strftime(s, 1000, "%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./HEJ config_file input_file\n\n";
	+ std::cerr << "\n# Usage:\n."<< argv[0] <<" config_file input_file\n\n";
	return EXIT_FAILURE;
	}

	const auto start_time = clock::now();
	+ {
	+ // char tmp_s[30];
	+ // const auto tmp_time = clock::to_time_t(start_time);
	+ // strftime(tmp_s, 30, "%a %b %d %Y %H:%M:%S", localtime(&(tmp_time)));
	+ std::cout << "Starting " << RHEJ::package_name_full() << " ("
	+ << print_time(clock::to_time_t(start_time)) << ")" << std::endl;
	+ }

	// 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, reader.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
	+ std::cout << "Processing " << max_events
	<< " out of " << input_events << " events\n";
	}
	RHEJ::EventReweighter rhej{
	reader.heprup,
	config
	};
	if(config.RanLux_init){
	RHEJ::PhaseSpacePoint::reset_ranlux(*config.RanLux_init);
	}

	int nevent = 0;
	std::array<int, RHEJ::event_type::last_type + 1>
	nevent_type{0}, nfailed_type{0};

	// 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++;
	+ ++nevent;

	- if (nevent % 10000 == 0)
	- std::cout << "event number " << nevent << std::endl;
	+ if (nevent % 10000 == 0){
	+ std::cout << "Passed " << nevent << " events ("
	+ << print_time(clock::to_time_t(clock::now())) << ")"<< std::endl;
	+ }

	// 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, config.scale_gen
	);
	++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);
	}
	}
	} // main event loop

	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 << "\nTask Runtime: " << run_time.count() << " seconds.\n";
	+ std::cout << "Finished " << RHEJ::package_name() << " at "
	+ << print_time(clock::to_time_t(clock::now()))
	+ << "\n=> Runtime: " << run_time.count() << " sec ("
	+ << nevent/run_time.count() << " Events/sec).\n";
	+
	}

File Metadata

Mime Type: text/x-diff
Expires: Tue, Jan 21, 1:46 AM (1 d, 11 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4243525
Default Alt Text: (16 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions