diff --git a/FixedOrderGen/include/Unweighter.hh b/FixedOrderGen/include/Unweighter.hh
deleted file mode 100644
index de8d9b5..0000000
--- a/FixedOrderGen/include/Unweighter.hh
+++ /dev/null
@@ -1,75 +0,0 @@
-/**
- * \authors The HEJ collaboration (see AUTHORS for details)
- * \date 2019
- * \copyright GPLv2 or later
- */
-#pragma once
-
-#include <cmath>
-
-#include "HEJ/optional.hh"
-#include "HEJ/RNG.hh"
-
-namespace HEJ {
- class Event;
-}
-
-namespace HEJFOG {
- namespace detail {
-
- bool has_jet_softer_than(HEJ::Event const & ev, double pt);
-
- template<typename Iterator>
- double calc_cut(
- Iterator begin, Iterator end, double max_dev,
- double min_pt
- ) {
- double mean = 0.;
- double err = 0.;
- double awt_sum = 0.;
- for(; begin != end; ++begin){
- if(has_jet_softer_than(*begin, min_pt)) continue;
- const double awt = std::abs(begin->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;
- return std::exp(mean + max_dev*err);
- }
- }
-
- class Unweighter {
-
- public:
- template<typename Iterator>
- Unweighter(
- Iterator begin, Iterator end, double max_dev,
- HEJ::RNG & ran,
- /* minimum pt of jets for an event to be considered for unweighting
- *
- * If the 'jets: peak pt' option is set to the *resummation* jet
- * threshold, events with softer jets 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.
- */
- double min_unweight_pt = 0.
- ):
- cut_{detail::calc_cut(begin, end, max_dev, min_unweight_pt)},
- min_unweight_pt_{min_unweight_pt},
- ran_{ran}
- {}
-
- HEJ::optional<HEJ::Event> unweight(HEJ::Event ev) const;
-
- private:
- double cut_;
- double min_unweight_pt_;
- std::reference_wrapper<HEJ::RNG> ran_;
- };
-
-}
diff --git a/FixedOrderGen/src/Unweighter.cc b/FixedOrderGen/src/Unweighter.cc
deleted file mode 100644
index 16996e5..0000000
--- a/FixedOrderGen/src/Unweighter.cc
+++ /dev/null
@@ -1,31 +0,0 @@
-/**
- * \authors The HEJ collaboration (see AUTHORS for details)
- * \date 2019
- * \copyright GPLv2 or later
- */
-#include "Unweighter.hh"
-
-#include <cassert>
-
-#include "HEJ/Event.hh"
-
-namespace HEJFOG {
-
- namespace detail {
- bool has_jet_softer_than(HEJ::Event const & ev, double pt) {
- assert(! ev.jets().empty());
- const auto softest_jet = fastjet::sorted_by_pt(ev.jets()).back();
- return softest_jet.pt() < pt;
- }
- }
-
- HEJ::optional<HEJ::Event> Unweighter::unweight(HEJ::Event ev) const {
- if(detail::has_jet_softer_than(ev, min_unweight_pt_)) return ev;
- const double awt = std::abs(ev.central().weight);
- if(ran_.get().flat() < awt/cut_) {
- if(awt < cut_) ev.parameters() *= cut_/awt;
- return ev;
- }
- return {};
- }
-}
diff --git a/FixedOrderGen/src/main.cc b/FixedOrderGen/src/main.cc
index d0b83a3..4138fec 100644
--- a/FixedOrderGen/src/main.cc
+++ b/FixedOrderGen/src/main.cc
@@ -1,234 +1,235 @@
/**
* \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 "yaml-cpp/yaml.h"
#include "LHEF/LHEF.h"
#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 "Unweighter.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);
}
}
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<HEJFOG::Unweighter> unweighter{};
+ 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 = HEJFOG::Unweighter{
- begin(events), end(events), config.unweight->max_dev, *ran,
+ unweighter = HEJ::Unweighter();
+ unweighter->set_middle(
+ events, config.unweight->max_dev,
config.jets.peak_pt?(*config.jets.peak_pt):0.
- };
+ );
std::vector<HEJ::Event> unweighted_events;
for(auto && ev: events) {
- auto unweighted = unweighter->unweight(std::move(ev));
+ 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));
+ 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 6850364..b6c3a0e 100644
--- a/include/HEJ/Unweighter.hh
+++ b/include/HEJ/Unweighter.hh
@@ -1,46 +1,59 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#pragma once
#include <functional>
#include <vector>
#include "HEJ/optional.hh"
#include "HEJ/RNG.hh"
namespace HEJ {
class Event;
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);
- //! @TODO estimate some reasonable cut for partial unweighting
- // void set_middle(std::vector<Event> const & evts)
+ //! 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,
+ double min_unweight_pt = 0
+ );
//! 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;
private:
double cut_;
};
}
diff --git a/src/Unweighter.cc b/src/Unweighter.cc
index 0cde3d0..1baec4b 100644
--- a/src/Unweighter.cc
+++ b/src/Unweighter.cc
@@ -1,50 +1,74 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "HEJ/Unweighter.hh"
#include "HEJ/Event.hh"
namespace HEJ {
namespace {
-
// returns true if element can be removed/gets discarded
// directly corrects weight if is accepted (not removed)
bool correct_weight(double cut, RNG & ran, Event & ev){
const double awt = std::abs(ev.central().weight);
if(awt >= cut ) return false;
if(ran.flat() > awt/cut) return true;
ev.parameters() *= cut/awt;
return false;
}
+
+ bool has_jet_softer_than(HEJ::Event const & ev, double pt) {
+ assert(! ev.jets().empty());
+ const auto softest_jet = fastjet::sorted_by_pt(ev.jets()).back();
+ return softest_jet.pt() < pt;
+ }
}
optional<Event> Unweighter::unweight(Event ev, RNG & ran) const {
if(correct_weight(get_cut(), ran, ev))
return {};
return ev;
}
std::vector<Event> Unweighter::unweight(
std::vector<Event> events, RNG & ran
) const {
events.erase( std::remove_if(events.begin(), events.end(),
[&](auto & ev){
return correct_weight(get_cut(), ran, ev);
}),
events.end());
return events;
}
void Unweighter::set_max(std::vector<Event> const & events){
set_cut(-1);
for(auto const & ev: events){
const double awt = std::abs(ev.central().weight);
if( awt > get_cut())
set_cut(awt);
}
}
+
+ void Unweighter::set_middle(
+ std::vector<Event> const & events, const double max_dev, const double min_pt
+ ){
+ double mean = 0.;
+ double err = 0.;
+ double awt_sum = 0.;
+ for(auto const & ev: events){
+ if(has_jet_softer_than(ev, min_pt)) continue;
+ 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));
+ }
}