diff --git a/FixedOrderGen/src/main.cc b/FixedOrderGen/src/main.cc index 3be02f3..461202a 100644 --- a/FixedOrderGen/src/main.cc +++ b/FixedOrderGen/src/main.cc @@ -1,233 +1,233 @@ /** * \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 "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.particles_properties, config.Higgs_coupling, *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{}; 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, 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)); 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)); 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/src/bin/HEJ.cc b/src/bin/HEJ.cc index 4733eb4..eb089aa 100644 --- a/src/bin/HEJ.cc +++ b/src/bin/HEJ.cc @@ -1,239 +1,240 @@ /** * \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" 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 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{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(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() && 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 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"; + return EXIT_SUCCESS; } diff --git a/t/check_hepmc.cc b/t/check_hepmc.cc index 0f824ac..e26be89 100644 --- a/t/check_hepmc.cc +++ b/t/check_hepmc.cc @@ -1,26 +1,27 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019 * \copyright GPLv2 or later */ #include <iostream> #include <stdexcept> #include "HepMC3/ReaderAscii.h" static constexpr double ep = 1e-3; int main(int argn, char** argv) { if(argn != 2){ std::cerr << "Usage: check_hepmc hepmc_file\n"; return EXIT_FAILURE; } HepMC3::ReaderAscii input{argv[1]}; if(input.failed()) throw std::runtime_error{"failed to open HepMC file"}; while(true){ HepMC3::GenEvent ev{}; if ( !input.read_event(ev) || ev.event_number() == 0 ) break; if(input.failed()) throw std::runtime_error{"failed to read HepMC event"}; } + return EXIT_SUCCESS; } diff --git a/t/check_lhe.cc b/t/check_lhe.cc index aafc0be..b770227 100644 --- a/t/check_lhe.cc +++ b/t/check_lhe.cc @@ -1,40 +1,41 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019 * \copyright GPLv2 or later */ #include <iostream> #include <unordered_map> #include "HEJ/stream.hh" #include "LHEF/LHEF.h" static constexpr double ep = 1e-3; int main(int argn, char** argv) { if(argn != 2){ std::cerr << "Usage: check_lhe lhe_file\n"; return EXIT_FAILURE; } HEJ::istream in{argv[1]}; LHEF::Reader reader{in}; std::unordered_map<int, double> xsec_ref; for(int i=0; i < reader.heprup.NPRUP; ++i) xsec_ref[reader.heprup.LPRUP[i]] = 0.; while(reader.readEvent()){ xsec_ref[reader.hepeup.IDPRUP] += reader.hepeup.weight(); } for(size_t i = 0; i < xsec_ref.size(); ++i){ double const ref = xsec_ref[reader.heprup.LPRUP[i]]; double const calc = reader.heprup.XSECUP[i]; std::cout << ref << '\t' << calc << '\n'; if(std::abs(calc-ref) > ep*calc){ std::cerr << "Cross sections deviate substantially"; return EXIT_FAILURE; } } + return EXIT_SUCCESS; } diff --git a/t/check_res.cc b/t/check_res.cc index 0d04ea9..274427b 100644 --- a/t/check_res.cc +++ b/t/check_res.cc @@ -1,165 +1,166 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019 * \copyright GPLv2 or later */ #include <iostream> #include <math.h> #include "LHEF/LHEF.h" #include "HEJ/Event.hh" #include "HEJ/EventReweighter.hh" #include "HEJ/Mixmax.hh" #include "HEJ/stream.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 fastjet::JetDefinition Born_jet_def{jet_def}; constexpr double Born_jetptmin = 30; constexpr double extpartonptmin = 30; constexpr double max_ext_soft_pt_fraction = 0.1; constexpr double jetptmin = 35; constexpr bool log_corr = false; using EventTreatment = HEJ::EventTreatment; using namespace HEJ::event_type; HEJ::EventTreatMap treat{ {no_2_jets, EventTreatment::discard}, {bad_final_state, EventTreatment::discard}, {non_resummable, EventTreatment::discard}, {unof, EventTreatment::discard}, {unob, EventTreatment::discard}, {qqxexb, EventTreatment::discard}, {qqxexf, EventTreatment::discard}, {qqxmid, EventTreatment::discard}, {FKL, EventTreatment::reweight} }; /// true if colour is allowed for particle bool correct_colour(HEJ::Particle const & part){ if(HEJ::is_AWZH_boson(part) && !part.colour) return true; if(!part.colour) return false; int const colour = part.colour->first; int const anti_colour = part.colour->second; if(part.type == HEJ::ParticleID::gluon) return colour != anti_colour && colour > 0 && anti_colour > 0; if(HEJ::is_quark(part)) return anti_colour == 0 && colour > 0; return colour == 0 && anti_colour > 0; } bool correct_colour(HEJ::Event const & ev){ if(!HEJ::event_type::is_resummable(ev.type())) return true; for(auto const & part: ev.incoming()){ if(!correct_colour(part)) return false; } for(auto const & part: ev.outgoing()){ if(!correct_colour(part)) return false; } return true; } }; int main(int argn, char** argv) { if(argn == 5 && std::string(argv[4]) == "unof"){ --argn; treat[unof] = EventTreatment::reweight; treat[unob] = EventTreatment::discard; treat[FKL] = EventTreatment::discard; } if(argn == 5 && std::string(argv[4]) == "unob"){ --argn; treat[unof] = EventTreatment::discard; treat[unob] = EventTreatment::reweight; treat[FKL] = EventTreatment::discard; } else if(argn == 5 && std::string(argv[4]) == "splitf"){ --argn; treat[qqxexb] = EventTreatment::discard; treat[qqxexf] = EventTreatment::reweight; treat[FKL] = EventTreatment::discard; } else if(argn == 5 && std::string(argv[4]) == "splitb"){ --argn; treat[qqxexb] = EventTreatment::reweight; treat[qqxexf] = EventTreatment::discard; treat[FKL] = EventTreatment::discard; } else if(argn == 5 && std::string(argv[4]) == "qqxmid"){ --argn; treat[qqxmid] = EventTreatment::reweight; treat[FKL] = EventTreatment::discard; } if(argn != 4){ std::cerr << "Usage: check_res eventfile xsection tolerance [uno]"; return EXIT_FAILURE; } const double xsec_ref = std::stod(argv[2]); const double tolerance = std::stod(argv[3]); HEJ::istream in{argv[1]}; LHEF::Reader reader{in}; HEJ::PhaseSpacePointConfig psp_conf; psp_conf.jet_param = HEJ::JetParameters{jet_def, jetptmin}; psp_conf.min_extparton_pt = extpartonptmin; psp_conf.max_ext_soft_pt_fraction = max_ext_soft_pt_fraction; HEJ::MatrixElementConfig ME_conf; ME_conf.log_correction = log_corr; ME_conf.Higgs_coupling = HEJ::HiggsCouplingSettings{}; HEJ::EventReweighterConfig conf; conf.psp_config = std::move(psp_conf); conf.ME_config = std::move(ME_conf); conf.jet_param = psp_conf.jet_param; conf.treat = treat; reader.readEvent(); const bool has_Higgs = std::find( begin(reader.hepeup.IDUP), end(reader.hepeup.IDUP), 25 ) != end(reader.hepeup.IDUP); const double mu = has_Higgs?125.:91.188; HEJ::ScaleGenerator scale_gen{ {{std::to_string(mu), HEJ::FixedScale{mu}}}, {}, 1. }; HEJ::Mixmax ran{}; HEJ::EventReweighter hej{reader.heprup, std::move(scale_gen), conf, ran}; double xsec = 0.; double xsec_err = 0.; do{ auto ev_data = HEJ::Event::EventData{reader.hepeup}; ev_data.reconstruct_intermediate(); HEJ::Event ev{ ev_data.cluster( Born_jet_def, Born_jetptmin ) }; auto resummed_events = hej.reweight(ev, 100); for(auto const & ev: resummed_events) { ASSERT(correct_colour(ev)); ASSERT(isfinite(ev.central().weight)); xsec += ev.central().weight; xsec_err += ev.central().weight*ev.central().weight; } } while(reader.readEvent()); xsec_err = std::sqrt(xsec_err); const double significance = std::abs(xsec - xsec_ref) / std::sqrt( xsec_err*xsec_err + tolerance*tolerance ); std::cout << xsec_ref << " +/- " << tolerance << " ~ " << xsec << " +- " << xsec_err << " => " << significance << " sigma\n"; if(significance > 3.){ std::cerr << "Cross section is off by over 3 sigma!\n"; return EXIT_FAILURE; } + return EXIT_SUCCESS; } diff --git a/t/cmp_events.cc b/t/cmp_events.cc index 93f7e9f..0a3c833 100644 --- a/t/cmp_events.cc +++ b/t/cmp_events.cc @@ -1,53 +1,53 @@ #include "HEJ/stream.hh" #include "HEJ/Event.hh" #include "LHEF/LHEF.h" int main(int argn, char** argv) { if(argn != 3){ std::cerr << "Usage: cmp_events eventfile eventfile_no_boson "; return EXIT_FAILURE; } HEJ::istream in{argv[1]}; LHEF::Reader reader{in}; HEJ::istream in_no_boson{argv[2]}; LHEF::Reader reader_no_boson{in_no_boson}; while(reader.readEvent()) { if(! reader_no_boson.readEvent()) { std::cerr << "wrong number of events in " << argv[2] << '\n'; return EXIT_FAILURE; } const auto is_AWZH = [](HEJ::Particle const & p) { return is_AWZH_boson(p); }; const HEJ::Event::EventData event_data{reader.hepeup}; const auto boson = std::find_if( begin(event_data.outgoing), end(event_data.outgoing), is_AWZH ); if(boson == end(event_data.outgoing)) { std::cerr << "no boson in event\n"; return EXIT_FAILURE; } HEJ::Event::EventData data_no_boson{reader_no_boson.hepeup}; data_no_boson.reconstruct_intermediate(); const auto new_boson = std::find_if( begin(event_data.outgoing), end(event_data.outgoing), is_AWZH ); if(new_boson == end(data_no_boson.outgoing)) { std::cerr << "failed to find reconstructed boson\n"; return EXIT_FAILURE; } if(new_boson->type != boson->type) { std::cerr << "reconstructed wrong boson type\n"; return EXIT_FAILURE; } if(!HEJ::nearby(new_boson->p, boson->p)) { std::cerr << "reconstructed boson momentum is off\n"; return EXIT_FAILURE; } - } + return EXIT_SUCCESS; } diff --git a/t/test_classify_ref.cc b/t/test_classify_ref.cc index a53a0a3..7a8edb4 100644 --- a/t/test_classify_ref.cc +++ b/t/test_classify_ref.cc @@ -1,95 +1,95 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019 * \copyright GPLv2 or later */ #include <fstream> #include <random> #include <algorithm> #include "LHEF/LHEF.h" #include "HEJ/YAMLreader.hh" #include "HEJ/event_types.hh" #include "HEJ/Event.hh" #include "HEJ/EventReader.hh" namespace{ // this is deliberately chosen bigger than in the generation, // to cluster multiple partons in one jet constexpr double min_jet_pt = 40.; const fastjet::JetDefinition jet_def{fastjet::kt_algorithm, 0.4}; void shuffle_particles(HEJ::Event::EventData & ev) { static std::mt19937_64 ran{0}; // incoming std::shuffle(begin(ev.incoming), end(ev.incoming), ran); // outgoing (through index) auto old_outgoing = std::move(ev.outgoing); std::vector<size_t> idx(old_outgoing.size()); std::iota(idx.begin(), idx.end(), 0); std::shuffle(begin(idx), end(idx), ran); ev.outgoing.clear(); ev.outgoing.reserve(old_outgoing.size()); for(size_t i: idx) { ev.outgoing.emplace_back(std::move(old_outgoing[i])); } // find decays again if(!ev.decays.empty()){ auto old_decays = std::move(ev.decays); ev.decays.clear(); for(size_t i=0; i<idx.size(); ++i) { auto decay = old_decays.find(idx[i]); if(decay != old_decays.end()) ev.decays.emplace(i, std::move(decay->second)); } } } } int main(int argn, char** argv) { if(argn != 3 && argn != 4){ std::cerr << "Usage: " << argv[0] << " reference_classification input_file.lhe\n"; return EXIT_FAILURE; } bool OUTPUT_MODE = false; if(argn == 4 && std::string("OUTPUT")==std::string(argv[3])) OUTPUT_MODE = true; std::fstream ref_file; if ( OUTPUT_MODE ) { std::cout << "_______________________USING OUTPUT MODE!_______________________" << std::endl; ref_file.open(argv[1], std::fstream::out); } else { ref_file.open(argv[1], std::fstream::in); } auto reader{ HEJ::make_reader(argv[2]) }; while(reader->read_event()){ HEJ::Event::EventData data{ reader->hepeup() }; shuffle_particles(data); const HEJ::Event ev{ data.cluster( jet_def, min_jet_pt ) }; if ( OUTPUT_MODE ) { ref_file << ev.type() << std::endl; } else { std::string line; if(!std::getline(ref_file,line)) break; const auto expected{static_cast<HEJ::event_type::EventType>(std::stoi(line))}; if(ev.type() != expected){ using HEJ::event_type::name; std::cerr << "wrong classification of event:\n" << ev << "classified as " << name(ev.type()) << ", is " << name(expected) << '\n'; return EXIT_FAILURE; } } } - + return EXIT_SUCCESS; } diff --git a/t/test_descriptions.cc b/t/test_descriptions.cc index 2d3c9d2..4ac462a 100644 --- a/t/test_descriptions.cc +++ b/t/test_descriptions.cc @@ -1,67 +1,68 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019 * \copyright GPLv2 or later */ #include <iostream> #include <cstddef> #include "HEJ/Event.hh" #include "HEJ/EventReweighter.hh" #include "HEJ/ScaleFunction.hh" #define ASSERT(x) if(!(x)) { \ std::cerr << "Assertion '" #x "' failed.\n"; \ return EXIT_FAILURE; \ } int main() { constexpr double mu = 125.; HEJ::ScaleFunction fun{"125", HEJ::FixedScale{mu}}; ASSERT(fun.name() == "125"); HEJ::ScaleGenerator scale_gen{ {std::move(fun)}, {0.5, 1, 2.}, 2.1 }; HEJ::Event::EventData tmp; tmp.outgoing.push_back( {HEJ::ParticleID::gluon, fastjet::PtYPhiM(50., -1., 0.3, 0.), {}} ); tmp.outgoing.push_back( {HEJ::ParticleID::gluon, fastjet::PtYPhiM(30., 1., -0.3, 0.), {}} ); HEJ::Event ev{ tmp.cluster( fastjet::JetDefinition{fastjet::kt_algorithm, 0.4}, 20. ) }; auto rescaled = scale_gen(std::move(ev)); ASSERT(rescaled.central().description->scale_name == "125"); for(auto const & var: rescaled.variations()) { ASSERT(var.description->scale_name == "125"); } ASSERT(rescaled.central().description->mur_factor == 1.); ASSERT(rescaled.central().description->muf_factor == 1.); ASSERT(rescaled.variations(0).description->mur_factor == 1.); ASSERT(rescaled.variations(0).description->muf_factor == 1.); ASSERT(rescaled.variations(1).description->mur_factor == 0.5); ASSERT(rescaled.variations(1).description->muf_factor == 0.5); ASSERT(rescaled.variations(2).description->mur_factor == 0.5); ASSERT(rescaled.variations(2).description->muf_factor == 1.); ASSERT(rescaled.variations(3).description->mur_factor == 1.); ASSERT(rescaled.variations(3).description->muf_factor == 0.5); ASSERT(rescaled.variations(4).description->mur_factor == 1.); ASSERT(rescaled.variations(4).description->muf_factor == 2.); ASSERT(rescaled.variations(5).description->mur_factor == 2.); ASSERT(rescaled.variations(5).description->muf_factor == 1.); ASSERT(rescaled.variations(6).description->mur_factor == 2.); ASSERT(rescaled.variations(6).description->muf_factor == 2.); + return EXIT_SUCCESS; } diff --git a/t/test_hdf5.cc b/t/test_hdf5.cc index 08f7c10..3f370e3 100644 --- a/t/test_hdf5.cc +++ b/t/test_hdf5.cc @@ -1,45 +1,46 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019 * \copyright GPLv2 or later */ #include <iostream> #include "HEJ/EventReader.hh" int main(int argc, char** argv) { if(argc != 2) { std::cerr << "Usage: " << argv[0] << " file.hdf5\n"; return EXIT_FAILURE; } auto reader = HEJ::make_reader(argv[1]); if( reader->heprup().EBMUP != std::make_pair(7000., 7000.) || reader->heprup().PDFSUP != std::make_pair(13000, 13000) ) { std::cerr << "Read incorrect init parameters\n"; return EXIT_FAILURE; } int nevent = 0; while(reader->read_event()) { ++nevent; if(reader->hepeup().NUP != 13) { std::cerr << "Read wrong number of particles: " << reader->hepeup().NUP << " != 13 in event " << nevent; return EXIT_FAILURE; } for(size_t i = 0; i < 2; ++i) { for(size_t j = 0; j < 2; ++j) { if(reader->hepeup().PUP[i][j] != 0) { std::cerr << "Non-vanishing transverse momentum in incoming particle" " in event " << nevent; return EXIT_FAILURE; } } } } if(nevent != 51200) { std::cerr << "Wrong number of events " << nevent << " != 51200\n"; return EXIT_FAILURE; } + return EXIT_SUCCESS; } diff --git a/t/test_psp.cc b/t/test_psp.cc index c232813..92b6a2b 100644 --- a/t/test_psp.cc +++ b/t/test_psp.cc @@ -1,68 +1,68 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019 * \copyright GPLv2 or later */ #include "LHEF/LHEF.h" #include "HEJ/stream.hh" #include "HEJ/config.hh" #include "HEJ/event_types.hh" #include "HEJ/Event.hh" #include "HEJ/PhaseSpacePoint.hh" #include "HEJ/Ranlux64.hh" namespace{ constexpr int max_trials = 100; constexpr double extpartonptmin = 45.; constexpr double max_ext_soft_pt_fraction = std::numeric_limits<double>::infinity(); const fastjet::JetDefinition jet_def{fastjet::kt_algorithm, 0.4}; constexpr double min_jet_pt = 50; }; int main(int argn, char** argv) { if(argn != 2){ std::cerr << "Usage: " << argv[0] << " eventfile"; return EXIT_FAILURE; } HEJ::istream in{argv[1]}; LHEF::Reader reader{in}; LHEF::Writer writer{std::cerr}; writer.heprup = reader.heprup; HEJ::PhaseSpacePointConfig conf; conf.jet_param = HEJ::JetParameters{jet_def, min_jet_pt}; conf.min_extparton_pt = extpartonptmin; conf.max_ext_soft_pt_fraction = max_ext_soft_pt_fraction; HEJ::Ranlux64 ran{}; while(reader.readEvent()){ const HEJ::Event ev{ HEJ::Event::EventData{reader.hepeup}( jet_def, min_jet_pt ) }; for(int trial = 0; trial < max_trials; ++trial){ HEJ::PhaseSpacePoint psp{ev, conf, ran}; if(psp.weight() != 0){ HEJ::Event::EventData tmp_ev; tmp_ev.incoming = psp.incoming(); tmp_ev.outgoing = psp.outgoing(); tmp_ev.parameters.central = {0,0,0}; HEJ::Event out_ev{ tmp_ev(jet_def, min_jet_pt) }; if(out_ev.type() != ev.type()){ using HEJ::event_type::name; std::cerr << "Wrong class of phase space point:\n" "original event of class " << name(ev.type()) << ":\n"; writer.hepeup = reader.hepeup; writer.writeEvent(); std::cerr << "Phase space point of class " << name(out_ev.type()) << ":\n"; writer.hepeup = to_HEPEUP(out_ev, &writer.heprup); writer.writeEvent(); return EXIT_FAILURE; } } } } - + return EXIT_SUCCESS; } diff --git a/t/test_scale_arithmetics.cc b/t/test_scale_arithmetics.cc index 3a66bab..933971e 100644 --- a/t/test_scale_arithmetics.cc +++ b/t/test_scale_arithmetics.cc @@ -1,120 +1,120 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019 * \copyright GPLv2 or later */ #include <fstream> #include <random> #include <algorithm> #include "LHEF/LHEF.h" #include "HEJ/EventReweighter.hh" #include "HEJ/make_RNG.hh" #include "HEJ/Event.hh" #include "HEJ/YAMLreader.hh" #include "HEJ/stream.hh" constexpr double alpha_s = 0.118; constexpr double ep = 1e-13; void dump(HEJ::Event const & ev){ { LHEF::Writer writer{std::cout}; std::cout << std::setprecision(6); writer.hepeup = to_HEPEUP(std::move(ev), nullptr); writer.writeEvent(); } std::cout << "Rapidity ordering:\n"; for(const auto & part: ev.outgoing()){ std::cout << std::setw(2) << part.type << ": "<< std::setw(7) << part.rapidity() << std::endl; } } void shuffle_particles(HEJ::Event::EventData & ev) { static std::mt19937_64 ran{0}; // incoming std::shuffle(begin(ev.incoming), end(ev.incoming), ran); // outgoing (through index) auto old_outgoing = std::move(ev.outgoing); std::vector<size_t> idx(old_outgoing.size()); std::iota(idx.begin(), idx.end(), 0); std::shuffle(begin(idx), end(idx), ran); ev.outgoing.clear(); ev.outgoing.reserve(old_outgoing.size()); for(size_t i: idx) { ev.outgoing.emplace_back(std::move(old_outgoing[i])); } // find decays again if(!ev.decays.empty()){ auto old_decays = std::move(ev.decays); ev.decays.clear(); for(size_t i=0; i<idx.size(); ++i) { auto decay = old_decays.find(idx[i]); if(decay != old_decays.end()) ev.decays.emplace(i, std::move(decay->second)); } } } int main(int argn, char** argv){ if(argn != 3){ std::cerr << "\n# Usage:\n."<< argv[0] <<" config.yml input_file.lhe\n\n"; return EXIT_FAILURE; } HEJ::Config config = HEJ::load_config(argv[1]); config.scales = HEJ::to_ScaleConfig( YAML::Load("scales: [H_T, 1 * H_T, 2/2 * H_T, 2*H_T/2, H_T/2*2, H_T/2/2*4, H_T*H_T/H_T]") ); HEJ::istream in{argv[2]}; LHEF::Reader reader{in}; auto ran = HEJ::make_RNG(config.rng.name, config.rng.seed); HEJ::ScaleGenerator scale_gen{ config.scales.base, config.scales.factors, config.scales.max_ratio }; HEJ::EventReweighter resum{ reader.heprup, std::move(scale_gen), to_EventReweighterConfig(config), *ran }; size_t i = 0; while(reader.readEvent()){ ++i; HEJ::Event::EventData data{reader.hepeup}; shuffle_particles(data); HEJ::Event event{ data.cluster( config.resummation_jets.def, config.resummation_jets.min_pt ) }; auto resummed = resum.reweight(event, config.trials); for(auto && ev: resummed) { for(auto &&var: ev.variations()) { if(std::abs(var.muf - ev.central().muf) > ep) { std::cerr << std::setprecision(15) << "unequal scales: " << var.muf << " != " << ev.central().muf << '\n' << "in resummed event:\n"; dump(ev); std::cerr << "\noriginal event:\n"; dump(event); return EXIT_FAILURE; } } } } - + return EXIT_SUCCESS; } diff --git a/t/test_scale_import.cc b/t/test_scale_import.cc index 36b72e9..2fc7a28 100644 --- a/t/test_scale_import.cc +++ b/t/test_scale_import.cc @@ -1,34 +1,35 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019 * \copyright GPLv2 or later */ #include <stdexcept> #include <iostream> #include "HEJ/YAMLreader.hh" #include "HEJ/Event.hh" int main(int argc, char** argv) { constexpr double ep = 1e-7; if (argc != 2) { throw std::logic_error{"wrong number of args"}; } const HEJ::Config config = HEJ::load_config(argv[1]); HEJ::Event::EventData tmp; tmp.outgoing.push_back( {HEJ::ParticleID::gluon, fastjet::PtYPhiM(50., -1., 0.3, 0.), {}} ); tmp.outgoing.push_back( {HEJ::ParticleID::gluon, fastjet::PtYPhiM(30., 1., -0.3, 0.), {}} ); HEJ::Event ev{ tmp(fastjet::JetDefinition{fastjet::kt_algorithm, 0.4}, 20.) }; const double softest_pt = config.scales.base[0](ev); if(std::abs(softest_pt-30.) > ep){ throw std::logic_error{"wrong softest pt"}; } + return EXIT_SUCCESS; }