Page MenuHomeHEPForge
Files F7878943

No OneTemporary
Actions

View Options

diff --git a/src/bin/HEJ.cc b/src/bin/HEJ.cc
index 55b4071..a048f3d 100644
--- a/src/bin/HEJ.cc
+++ b/src/bin/HEJ.cc
@@ -1,363 +1,372 @@
/**
* \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/BufferedEventReader.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;
}
HEJ::Event to_event(
LHEF::HEPEUP const & hepeup,
HEJ::JetParameters const & fixed_order_jets
) {
HEJ::Event::EventData event_data{hepeup};
event_data.reconstruct_intermediate();
return HEJ::Event{
std::move(event_data).cluster(
fixed_order_jets.def, fixed_order_jets.min_pt
)
};
}
void unweight(
HEJ::Unweighter & unweighter,
HEJ::WeightType weight_type,
std::vector<HEJ::Event> & events,
HEJ::RNG & ran
) {
if(weight_type == HEJ::WeightType::unweighted_resum){
unweighter.set_cut_to_maxwt(events);
}
events.erase(
unweighter.unweight(begin(events), end(events), ran),
end(events)
);
}
// peek up to nevents events from reader
std::vector<LHEF::HEPEUP> peek_events(
HEJ::BufferedEventReader & reader,
const int nevents
) {
std::vector<LHEF::HEPEUP> events;
while(
static_cast<int>(events.size()) < nevents
&& reader.read_event()
) {
events.emplace_back(reader.hepeup());
}
// put everything back into the reader
for(auto it = rbegin(events); it != rend(events); ++it) {
reader.emplace(*it);
}
return events;
}
void append_resummed_events(
std::vector<HEJ::Event> & resummation_events,
HEJ::EventReweighter & reweighter,
LHEF::HEPEUP const & hepeup,
const int trials,
HEJ::JetParameters const & fixed_order_jets
) {
const HEJ::Event FO_event = to_event(hepeup, fixed_order_jets);
if(reweighter.treatment(FO_event.type()) != HEJ::EventTreatment::reweight) {
return;
}
const auto resummed = reweighter.reweight(FO_event, trials);
resummation_events.insert(
end(resummation_events),
begin(resummed), end(resummed)
);
}
void train(
HEJ::Unweighter & unweighter,
HEJ::BufferedEventReader & reader,
HEJ::EventReweighter & reweighter,
const int total_trials,
const double max_dev,
+ double reweight_factor,
HEJ::JetParameters const & fixed_order_jets
) {
- const auto FO_events = peek_events(reader, total_trials);
+ std::cout << "Reading up to " << total_trials << " training events...\n";
+ auto FO_events = peek_events(reader, total_trials);
const int trials = total_trials/FO_events.size();
+ // adjust reweight factor so that the overall normalisation
+ // is the same as in the full run
+ reweight_factor *= trials;
+ for(auto & hepeup: FO_events) {
+ hepeup.setWeight(0, reweight_factor * hepeup.weight());
+ }
std::cout << "Training unweighter with "
<< trials << '*' << FO_events.size() << " events\n";
auto progress = HEJ::ProgressBar<int>{
std::cout, static_cast<int>(FO_events.size())
};
std::vector<HEJ::Event> resummation_events;
for(auto const & hepeup: FO_events) {
append_resummed_events(
resummation_events,
reweighter, hepeup, trials, fixed_order_jets
);
++progress;
}
unweighter.set_cut_to_peakwt(resummation_events, max_dev);
std::cout << "\nUnweighting events with weight up to "
<< unweighter.get_cut() << '\n';
}
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::optional<HEJ::Unweighter> unweighter{};
if(config.weight_type != HEJ::WeightType::weighted) {
unweighter = HEJ::Unweighter{};
}
if(config.weight_type == HEJ::WeightType::partially_unweighted) {
// TODO: turn these into configuration options
constexpr unsigned long training_trials = 100000;
constexpr double max_dev = 0.;
HEJ::BufferedEventReader buffered_reader{std::move(reader)};
train(
*unweighter,
buffered_reader,
hej,
training_trials,
max_dev,
+ global_reweight/config.trials,
config.fixed_order_jets
);
reader = std::make_unique<HEJ::BufferedEventReader>(
std::move(buffered_reader)
);
}
// 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());
const auto FO_event = to_event(hepeup, config.fixed_order_jets);
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()];
if(unweighter) {
unweight(*unweighter, config.weight_type, resummed_events, *ran);
}
// analysis
for(auto & 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);
} else {
ev.parameters()*=0; // do not use discarded events afterwards
}
}
xs.fill_correlated(resummed_events);
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;
}

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 7:02 PM (1 d, 9 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3801153
Default Alt Text
(11 KB)

Event Timeline