diff --git a/Changes.md b/Changes.md index 212bb7b..176f662 100644 --- a/Changes.md +++ b/Changes.md @@ -1,53 +1,54 @@ # Changelog This is the log for changes to the HEJ program. Further changes to the HEJ API are documented in `Changes-API.md`. If you are using HEJ as a library, please also read the changes there. ## Version 2.X ### 2.X.0 * Resummation for W bosons with jets - New subleading processes `extremal qqx` & `central qqx` for a quark and anti-quark in the final state, e.g. `g g => u d_bar Wm g` (the other subleading processes also work with W's) - `HEJFOG` can generate mutliple jets together with a (off-shell) W bosons decaying into lepton & neutrino * Allow multiplication and division of multiple scale functions e.g. `H_T/2*m_j1j2` * Print cross sections at end of run * Follow HepMC convention for particle Status codes: incoming = 11, decaying = 2, outgoing = 1 (unchanged) * Partons now have a Colour charge - Colours are read from and written to LHE files - For reweighted events the colours are created according to leading colour in the FKL limit * Allow changing the regulator lambda in input (`regulator parameter`, only for advanced users) * Use `git-lfs` for raw data in test (`make test` now requires `git-lfs`) * Added support to read `hdf5` event files suggested in [arXiv:1905.05120](https://arxiv.org/abs/1905.05120) (needs [HighFive](https://github.com/BlueBrain/HighFive)) +* Support input with avarage weight equal to the cross section (`IDWTUP=1 or 4`) ## 2.0.5 * Fixed event classification for input not ordered in rapidity ### 2.0.4 * Fixed wrong path of `HEJ_INCLUDE_DIR` in `hej-config.cmake` ### 2.0.3 * Fixed parsing of (numerical factor) * (base scale) in configuration * Don't change scale names, but sanitise Rivet output file names instead ### 2.0.2 * Changed scale names to `"_over_"` and `"_times_"` for proper file names (was `"/"` and `"*"` before) ### 2.0.1 * Fixed name of fixed-order generator in error message. diff --git a/src/bin/HEJ.cc b/src/bin/HEJ.cc index b4f1a12..2472ff3 100644 --- a/src/bin/HEJ.cc +++ b/src/bin/HEJ.cc @@ -1,230 +1,250 @@ /** * \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/Version.hh" #include "HEJ/YAMLreader.hh" -int event_number(std::string const & record){ +HEJ::optional<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"); + return {}; } const size_t end = record.find_first_not_of("0123456789", start); return std::stoi(record.substr(start, end - start)); } 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.; int max_events = std::numeric_limits<int>::max(); - if(argn > 3){ - max_events = std::stoi(argv[3]); - const int input_events = event_number(reader->header()); - global_reweight = input_events/static_cast<double>(max_events); - std::cout << "Processing " << max_events - << " out of " << input_events << " events\n"; + // if we need the event number: + if(std::abs(heprup.IDWTUP) == 4 || std::abs(heprup.IDWTUP) == 1 || argn > 3){ + // try to read from LHE head + auto input_events{event_number(reader->header())}; + 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(argn > 3){ + // maximal number of events given + max_events = std::stoi(argv[3]); + 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 }; // status infos & eye candy int 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; // Loop over the events in the input file while(reader->read_event()){ // reweight events so that the total cross section is conserved auto hepeup = reader->hepeup(); hepeup.setWeight(0, global_reweight * hepeup.weight()); if(nevent == max_events) break; ++nevent; HEJ::Event::EventData event_data{hepeup}; event_data.reconstruct_intermediate(); // calculate HEJ weight 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); 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()]; 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); } } 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 << '\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"; }