No OneTemporary
Actions

Size

26 KB

Subscribers

None

View Options

	diff --git a/FixedOrderGen/src/main.cc b/FixedOrderGen/src/main.cc
	index f32d717..ae102f7 100644
	--- a/FixedOrderGen/src/main.cc
	+++ b/FixedOrderGen/src/main.cc
	@@ -1,251 +1,251 @@
	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019
	* \copyright GPLv2 or later
	*/
	#include <algorithm>
	#include <chrono>
	#include <fstream>
	#include <iostream>
	#include <map>
	#include <memory>

	#include "LHEF/LHEF.h"

	#include "yaml-cpp/yaml.h"

	#include <boost/iterator/filter_iterator.hpp>

	#include "HEJ/CombinedEventWriter.hh"
	#include "HEJ/CrossSectionAccumulator.hh"
	#include "HEJ/get_analysis.hh"
	#include "HEJ/LesHouchesWriter.hh"
	#include "HEJ/make_RNG.hh"
	#include "HEJ/ProgressBar.hh"
	#include "HEJ/stream.hh"
	#include "HEJ/Unweighter.hh"

	#include "config.hh"
	#include "EventGenerator.hh"
	#include "PhaseSpacePoint.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<HEJ::Analysis> get_analysis(
	YAML::Node const & parameters
	){
	try{
	return HEJ::get_analysis(parameters);
	}
	catch(std::exception const & exc){
	std::cerr << "Failed to load analysis: " << exc.what() << '\n';
	std::exit(EXIT_FAILURE);
	}
	}

	template<class Iterator>
	auto make_lowpt_filter(Iterator begin, Iterator end, HEJ::optional<double> peak_pt){
	return boost::make_filter_iterator(
	[peak_pt](HEJ::Event const & ev){
	assert(! ev.jets().empty());
	double min_pt = peak_pt?(*peak_pt):0.;
	const auto softest_jet = fastjet::sorted_by_pt(ev.jets()).back();
	return softest_jet.pt() > min_pt;
	},
	begin, end
	);
	}

	int main(int argn, char** argv) {
	using namespace std::string_literals;
	if (argn < 2) {
	std::cerr << "\n# Usage:\n." << argv[0] << " 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<HEJ::Analysis> analysis = get_analysis(
	config.analysis_parameters
	);
	assert(analysis != nullptr);

	auto ran = HEJ::make_RNG(config.rng.name, config.rng.seed);
	assert(ran != nullptr);

	HEJ::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.subleading_fraction,
	config.subleading_channels,
	config.particle_decays,
	config.Higgs_coupling,
	config.ew_parameters,
	*ran
	};

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

	HEJ::CombinedEventWriter writer{config.output, heprup};

	HEJ::optional<HEJ::Unweighter> unweighter{};
	std::map<HEJFOG::Status, int> status_counter;

	std::vector<HEJ::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;
	HEJ::ProgressBar<size_t> warmup_progress{std::cout, warmup_events};
	for(; events.size() < warmup_events; ++trials){
	auto ev = generator.gen_event();
	++status_counter[generator.status()];
	assert( (generator.status() == HEJFOG::good) == bool(ev) );
	if(generator.status() == HEJFOG::good && analysis->pass_cuts(ev, ev)) {
	events.emplace_back(std::move(*ev));
	++warmup_progress;
	}
	}
	std::cout << std::endl;
	unweighter = HEJ::Unweighter();
	- unweighter->set_middle(
	+ unweighter->set_cut_to_peakwt(
	make_lowpt_filter(events.cbegin(), events.cend(), config.jets.peak_pt),
	make_lowpt_filter(events.cend(), events.cend(), config.jets.peak_pt),
	config.unweight->max_dev
	);
	std::vector<HEJ::Event> unweighted_events;
	for(auto && ev: events) {
	auto unweighted = unweighter->unweight(std::move(ev), *ran);
	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";
	}

	HEJ::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()];
	assert( (generator.status() == HEJFOG::good) == bool(ev) );
	if(generator.status() == HEJFOG::good && analysis->pass_cuts(ev, ev)) {
	if(unweighter) {
	auto unweighted = unweighter->unweight(std::move(ev), ran);
	if(! unweighted) continue;
	ev = std::move(unweighted);
	}
	events.emplace_back(std::move(*ev));
	++progress;
	}
	}
	std::cout << std::endl;

	HEJ::CrossSectionAccumulator xs;
	for(auto & ev: events){
	ev.parameters() *= 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 << xs << '\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";
	}
	return EXIT_SUCCESS;
	}
	diff --git a/include/HEJ/Unweighter.hh b/include/HEJ/Unweighter.hh
	index fa86891..0ab3217 100644
	--- a/include/HEJ/Unweighter.hh
	+++ b/include/HEJ/Unweighter.hh
	@@ -1,82 +1,82 @@
	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019
	* \copyright GPLv2 or later
	*/
	#pragma once

	#include <functional>
	#include <vector>

	#include "HEJ/optional.hh"

	namespace HEJ {
	class Event;
	class RNG;

	class Unweighter {
	public:
	Unweighter(double cut = -1.):
	cut_{cut}{}

	//! explicitly set cut
	void set_cut(double max_weight){
	cut_ = max_weight;
	}
	//! set cut as max(weight) of events
	- void set_max(std::vector<Event> const & events){
	- set_max(events.cbegin(), events.cend());
	+ void set_cut_to_maxwt(std::vector<Event> const & events){
	+ set_cut_to_maxwt(events.cbegin(), events.cend());
	}
	//! iterator version of set_max()
	- template<class constIt>
	- void set_max(constIt begin, constIt end);
	+ template<class ConstIt>
	+ void set_cut_to_maxwt(ConstIt begin, ConstIt end);

	//! estimate some reasonable cut for partial unweighting
	/**
	* @param max_dev standard derivation to include above mean weight
	* @param min_unweight_pt minimum pt of jets for an event to be considered
	*
	* @note Events with softer jets than the resummation jet threshold will
	* have a spurious large weight, although they hardly contribute after
	* resummation. This destroys the unweighting efficiency. By setting
	* min_unweight_pt to the same threshold, we can exclude these events
	* from unweighting.
	*/
	- void set_middle( std::vector<Event> const & events, double max_dev){
	- set_middle(events.cbegin(), events.cend(), max_dev);
	+ void set_cut_to_peakwt( std::vector<Event> const & events, double max_dev){
	+ set_cut_to_peakwt(events.cbegin(), events.cend(), max_dev);
	}
	//! iterator version of set_middle()
	- template<class constIt>
	- void set_middle(constIt begin, constIt end, double max_dev);
	+ template<class ConstIt>
	+ void set_cut_to_peakwt(ConstIt begin, ConstIt end, double max_dev);

	//! return current value of the cut
	double get_cut() const {
	return cut_;
	}

	//! unweight one event, returns original event if weight > get_cut()
	optional<Event> unweight(Event ev, RNG & ran) const;
	//! unweight for multiple events at once
	std::vector<Event> unweight(
	std::vector<Event> events, RNG & ran
	) const;
	//! @brief iterator implementation of unweight(),
	/**
	* usage similar to std::remove(), i.e. use with erase()
	*
	* @return beginning of "discarded" range
	*/
	template<class Iterator>
	Iterator unweight(
	Iterator begin, Iterator end, RNG & ran
	) const;
	private:
	double cut_;
	//! returns true if element can be removed/gets discarded
	//! directly corrects weight if is accepted (not removed)
	bool discard(RNG & ran, Event & ev) const;
	};
	}

	// implementation of template functions
	#include "HEJ/detail/Unweighter_impl.hh"
	diff --git a/include/HEJ/detail/Unweighter_impl.hh b/include/HEJ/detail/Unweighter_impl.hh
	index d42be92..f8b49ea 100644
	--- a/include/HEJ/detail/Unweighter_impl.hh
	+++ b/include/HEJ/detail/Unweighter_impl.hh
	@@ -1,53 +1,53 @@
	/** \file
	* \brief Unweighter Class Implementation for template functions
	*
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019
	* \copyright GPLv2 or later
	*/
	#pragma once

	#include "HEJ/Unweighter.hh"

	#include "HEJ/Event.hh"
	#include "HEJ/RNG.hh"

	namespace HEJ {
	- template<class constIt>
	- void Unweighter::set_middle(constIt begin, constIt end, double max_dev){
	+ template<class ConstIt>
	+ void Unweighter::set_cut_to_peakwt(ConstIt begin, ConstIt end, double max_dev){
	double mean = 0.;
	double err = 0.;
	double awt_sum = 0.;
	for(; begin!=end; ++begin){
	auto const & ev = *begin;
	const double awt = std::abs(ev.central().weight);
	const double tmp = awt*std::log(awt);
	mean += tmp;
	err += tmp*tmp;
	awt_sum += awt;
	}
	mean /= awt_sum;
	err = std::sqrt(err)/awt_sum;
	set_cut(std::exp(mean + max_dev*err));
	}

	- template<class constIt>
	- void Unweighter::set_max(constIt begin, constIt end){
	+ template<class ConstIt>
	+ void Unweighter::set_cut_to_maxwt(ConstIt begin, ConstIt end){
	set_cut(-1);
	for(; begin!=end; ++begin){
	const double awt = std::abs(begin->central().weight);
	if( awt > get_cut())
	set_cut(awt);
	}
	}
	template<class Iterator>
	Iterator Unweighter::unweight(
	Iterator begin, Iterator end, RNG & ran
	) const {
	if(get_cut() < 0) return end;
	return std::remove_if(begin, end,
	[&](auto & ev){
	return discard(ran, ev);
	});
	}
	}
	diff --git a/src/bin/HEJ.cc b/src/bin/HEJ.cc
	index 2db7982..569767f 100644
	--- a/src/bin/HEJ.cc
	+++ b/src/bin/HEJ.cc
	@@ -1,255 +1,255 @@
	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019
	* \copyright GPLv2 or later
	*/
	#include <array>
	#include <chrono>
	#include <iostream>
	#include <limits>
	#include <memory>
	#include <numeric>

	#include "yaml-cpp/yaml.h"

	#include "fastjet/ClusterSequence.hh"

	#include "HEJ/CombinedEventWriter.hh"
	#include "HEJ/config.hh"
	#include "HEJ/CrossSectionAccumulator.hh"
	#include "HEJ/Event.hh"
	#include "HEJ/EventReader.hh"
	#include "HEJ/EventReweighter.hh"
	#include "HEJ/get_analysis.hh"
	#include "HEJ/make_RNG.hh"
	#include "HEJ/optional.hh"
	#include "HEJ/ProgressBar.hh"
	#include "HEJ/stream.hh"
	#include "HEJ/Unweighter.hh"
	#include "HEJ/Version.hh"
	#include "HEJ/YAMLreader.hh"

	HEJ::Config load_config(char const * filename){
	try{
	return HEJ::load_config(filename);
	}
	catch(std::exception const & exc){
	std::cerr << "Error: " << exc.what() << '\n';
	std::exit(EXIT_FAILURE);
	}
	}

	std::unique_ptr<HEJ::Analysis> get_analysis(
	YAML::Node const & parameters
	){
	try{
	return HEJ::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:
	// HEJ::optional is a better fit, but gives spurious errors with g++ 7.3.0
	std::unique_ptr<HEJ::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<HEJ::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 " << HEJ::Version::package_name_full()
	<< ", revision " << HEJ::Version::revision() << " ("
	<< time_to_string(clock::to_time_t(start_time)) << ")" << std::endl;
	}
	fastjet::ClusterSequence::print_banner();

	// read configuration
	const HEJ::Config config = load_config(argv[1]);
	auto reader = HEJ::make_reader(argv[2]);
	assert(reader);

	std::unique_ptr<HEJ::Analysis> analysis = get_analysis(
	config.analysis_parameters
	);
	assert(analysis != nullptr);

	auto heprup{ reader->heprup() };
	heprup.generators.emplace_back(LHEF::XMLTag{});
	heprup.generators.back().name = HEJ::Version::package_name();
	heprup.generators.back().version = HEJ::Version::String();

	analysis->initialise(heprup);
	HEJ::CombinedEventWriter writer{config.output, std::move(heprup)};

	double global_reweight = 1.;
	const auto & max_events = config.max_events;
	// if we need the event number:
	if(std::abs(heprup.IDWTUP) == 4 \|\| std::abs(heprup.IDWTUP) == 1 \|\| max_events){
	// try to read from LHE head
	auto input_events{reader->number_events()};
	if(!input_events) {
	// else count manually
	auto t_reader = HEJ::make_reader(argv[2]);
	input_events = 0;
	while(t_reader->read_event()) ++(*input_events);
	}
	if(std::abs(heprup.IDWTUP) == 4 \|\| std::abs(heprup.IDWTUP) == 1){
	// IDWTUP 4 or 1 assume average(weight)=xs, but we need sum(weights)=xs
	std::cout << "Found IDWTUP " << heprup.IDWTUP << ": "
	<< "assuming \"cross section = average weight\".\n"
	<< "converting to \"cross section = sum of weights\" ";
	global_reweight /= *input_events;
	}
	if(max_events && (input_events > max_events)){
	// maximal number of events given
	global_reweight = input_events/static_cast<double>(*max_events);
	std::cout << "Processing " << *max_events
	<< " out of " << *input_events << " events\n";
	}
	}

	HEJ::ScaleGenerator scale_gen{
	config.scales.base,
	config.scales.factors,
	config.scales.max_ratio
	};
	auto ran = HEJ::make_RNG(config.rng.name, config.rng.seed);
	assert(ran != nullptr);
	HEJ::EventReweighter hej{
	reader->heprup(),
	std::move(scale_gen),
	to_EventReweighterConfig(config),
	*ran
	};

	HEJ::Unweighter unweighter;

	// status infos & eye candy
	size_t nevent = 0;
	std::array<int, HEJ::event_type::last_type + 1>
	nevent_type{0}, nfailed_type{0};
	auto progress = make_progress_bar(reader->heprup().XSECUP);
	HEJ::CrossSectionAccumulator xs;
	std::map<HEJ::StatusCode, int> status_counter;
	size_t total_trials = 0;
	size_t total_resum = 0;

	// Loop over the events in the input file
	while(reader->read_event() && (!max_events \|\| nevent < *max_events) ){
	++nevent;

	// reweight events so that the total cross section is conserved
	auto hepeup = reader->hepeup();
	hepeup.setWeight(0, global_reweight * hepeup.weight());

	HEJ::Event::EventData event_data{hepeup};
	event_data.reconstruct_intermediate();

	// calculate HEJ weight
	const HEJ::Event FO_event{
	std::move(event_data).cluster(
	config.fixed_order_jets.def, config.fixed_order_jets.min_pt
	)
	};
	if(FO_event.central().weight == 0) {
	static const bool warned_once = [argv,nevent](){
	std::cerr
	<< "WARNING: event number " << nevent
	<< " in " << argv[2] << " has zero weight. "
	"Ignoring this and all further events with vanishing weight.\n";
	return true;
	}();
	(void) warned_once; // shut up compiler warnings
	continue;
	}

	auto resummed_events{ hej.reweight(FO_event, config.trials) };

	// some bookkeeping
	for(auto const & s: hej.status())
	++status_counter[s];
	total_trials+=hej.status().size();
	++nevent_type[FO_event.type()];

	if(resummed_events.empty()) ++nfailed_type[FO_event.type()];

	// unweight
	if(config.weight_type == HEJ::WeightType::unweighted_resum){
	- unweighter.set_max(resummed_events);
	+ unweighter.set_cut_to_maxwt(resummed_events);
	}
	resummed_events.erase( unweighter.unweight(
	resummed_events.begin(), resummed_events.end(), *ran ),
	resummed_events.end() );

	// analysis
	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);
	xs.fill(ev);
	}
	}
	total_resum += resummed_events.size();
	if(progress) progress->increment(FO_event.central().weight);
	} // main event loop
	std::cout << '\n';
	analysis->finalise();

	using namespace HEJ::event_type;
	std::cout<< "Events processed: " << nevent << " (" << total_resum << " resummed)"<< '\n';
	std::cout << '\t' << name(EventType::first_type) << ": "
	<< nevent_type[EventType::first_type]
	<< ", failed to reconstruct " << nfailed_type[EventType::first_type]
	<< '\n';
	for(auto i=EventType::first_type+1; i<=EventType::last_type; i*=2){
	std::cout << '\t' << name(static_cast<EventType>(i)) << ": "
	<< nevent_type[i]
	<< ", failed to reconstruct " << nfailed_type[i]
	<< '\n';
	}

	std::cout << '\n' << xs << '\n';

	std::cout << "Generation statistic: "
	<< status_counter[HEJ::StatusCode::good] << "/" << total_trials
	<< " trials successful.\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 << std::left << std::setw(17) << (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 << "] " <<std::setw(2)<<std::right<< percent << "%\n";
	}

	std::chrono::duration<double> run_time = (clock::now() - start_time);
	std::cout << "\nFinished " << HEJ::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";

	return EXIT_SUCCESS;
	}
	diff --git a/t/test_unweighter.cc b/t/test_unweighter.cc
	index 1685595..8fa957a 100644
	--- a/t/test_unweighter.cc
	+++ b/t/test_unweighter.cc
	@@ -1,150 +1,150 @@
	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019
	* \copyright GPLv2 or later
	*/

	#include "HEJ/CrossSectionAccumulator.hh"
	#include "HEJ/Event.hh"
	#include "HEJ/EventReader.hh"
	#include "HEJ/Mixmax.hh"
	#include "HEJ/Unweighter.hh"

	#define ASSERT(x) if(!(x)) { \
	std::cerr << "Assertion '" #x "' failed.\n"; \
	return EXIT_FAILURE; \
	}

	namespace{
	const fastjet::JetDefinition jet_def{fastjet::kt_algorithm, 0.4};
	const double min_pt{30.};
	const double sample_ratio{10.};
	const double max_dev{2.};

	bool compare_xs(
	HEJ::XSWithError<double> const & xs1, HEJ::XSWithError<double> const & xs2
	){
	std::cout << xs1.value << "+/-" << xs1.error << " vs. "
	<< xs2.value << "+/-" << xs2.error << std::endl;
	return std::abs(xs1.value/xs2.value-1.)<xs1.error;
	}
	}
	int main(int argc, char** argv) {
	if(argc != 2) {
	std::cerr << "Usage: " << argv[0] << " event_file\n";
	return EXIT_FAILURE;
	}
	std::string file{argv[1]};
	auto reader = HEJ::make_reader(file);

	// number of events
	auto nevents{reader->number_events()};
	if(!nevents) {
	auto t_reader = HEJ::make_reader(file);
	nevents = 0;
	while(t_reader->read_event()) ++(*nevents);
	}
	ASSERT(*nevents>sample_ratio);
	const size_t size_sample = floor(*nevents/sample_ratio);
	HEJ::Mixmax ran{};

	// no unweighting
	HEJ::CrossSectionAccumulator xs_base;

	std::vector<HEJ::Event> all_evts;
	// full unweighting
	HEJ::CrossSectionAccumulator xs_max;
	HEJ::Unweighter unw_max;
	size_t n_max{0};
	// midpoint on full sample
	HEJ::CrossSectionAccumulator xs_mid;
	HEJ::Unweighter unw_mid;
	size_t n_mid{0};

	// calc max from partial sample
	HEJ::CrossSectionAccumulator xs_pmax;
	HEJ::Unweighter unw_pmax;
	size_t n_pmax{0};
	// midpoint on partial sample
	HEJ::CrossSectionAccumulator xs_pmid;
	HEJ::Unweighter unw_pmid;
	size_t n_pmid{0};

	// initialise sample
	for(size_t n = 0; n < size_sample; ++n){
	if(!reader->read_event()){
	std::cerr << "Sample size bigger than event sample\n";
	return EXIT_FAILURE;
	}
	const HEJ::Event event{
	HEJ::Event::EventData{reader->hepeup()}.cluster(jet_def, min_pt)
	};

	xs_base.fill(event);
	all_evts.push_back(event);
	}

	// calculate partial settings
	- unw_pmax.set_max(all_evts);
	- unw_pmid.set_middle(all_evts, max_dev);
	+ unw_pmax.set_cut_to_maxwt(all_evts);
	+ unw_pmid.set_cut_to_peakwt(all_evts, max_dev);

	for(auto const & ev: unw_pmax.unweight(all_evts, ran)){
	xs_pmax.fill(ev);
	++n_pmax;
	}
	for(auto const & ev: unw_pmid.unweight(all_evts, ran)){
	xs_pmid.fill(ev);
	++n_pmid;
	}

	while(reader->read_event()){
	const HEJ::Event event{
	HEJ::Event::EventData{reader->hepeup()}.cluster(jet_def, min_pt)
	};
	xs_base.fill(event);
	auto ev{ unw_pmid.unweight(event, ran) };
	if(ev){
	xs_pmid.fill(*ev);
	++n_pmid;
	}
	ev = unw_pmax.unweight(event, ran);
	if(ev){
	xs_pmax.fill(*ev);
	++n_pmax;
	}
	all_evts.push_back(event);
	}

	- unw_max.set_max(all_evts);
	- unw_mid.set_middle(all_evts, max_dev);
	+ unw_max.set_cut_to_maxwt(all_evts);
	+ unw_mid.set_cut_to_peakwt(all_evts, max_dev);

	for(auto const & ev: unw_max.unweight(all_evts, ran)){
	// make sure all the events have the same weight
	ASSERT( std::abs( std::abs(unw_max.get_cut()/ev.central().weight)-1. ) < 10e-16);
	xs_max.fill(ev);
	++n_max;
	}
	for(auto const & ev: unw_mid.unweight(all_evts, ran)){
	xs_mid.fill(ev);
	++n_mid;
	}

	// sanity check number of events
	ASSERT( n_pmax < all_evts.size() );
	ASSERT( n_max < n_pmax );
	ASSERT( n_pmid < all_evts.size() );
	ASSERT( n_mid < all_evts.size() );
	ASSERT( n_max < n_mid );

	std::cout << "all_evts.size() " << all_evts.size() << " n_pmax " << n_pmax <<
	" n_max " << n_max << " n_pmid " << n_pmid << " n_mid " << n_mid << std::endl;

	// control xs (in circle)
	ASSERT(compare_xs( xs_base.total(), xs_pmax.total() ));
	ASSERT(compare_xs( xs_pmax.total(), xs_max.total() ));
	ASSERT(compare_xs( xs_max.total() , xs_pmid.total() ));
	ASSERT(compare_xs( xs_pmid.total(), xs_mid.total() ));
	ASSERT(compare_xs( xs_mid.total() , xs_base.total() ));

	return EXIT_SUCCESS;
	}

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 8:18 PM (1 d, 3 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3804465
Default Alt Text: (26 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions