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; }