diff --git a/Changes.md b/Changes.md index aea4b72..b30bb04 100644 --- a/Changes.md +++ b/Changes.md @@ -1,33 +1,34 @@ # Version 2.0 ## 2.X.0 * Made `MatrixElement::tree_kin` and `MatrixElement::tree_param` public * Allow multiplication and division of multiple scale functions e.g. `H_T/2*m_j1j2` * Follow HepMC convention for particle Status codes: incoming = 11, decaying = 2, outgoing = 1 (unchanged) * New class `CrossSectionAccumulator` to keep track of Cross Section of the different subproccess +* Removed default constructor for `Event` ## 2.0.4 * Fixed wrong path of `HEJ_INCLUDE_DIR` in `hej-config.cmake` * Correctly include rivet headers * New `EXCLUDE_package` variable in `cmake` to not interface to specific packages * Consistent search and include for side packages in `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/FixedOrderGen/include/EventGenerator.hh b/FixedOrderGen/include/EventGenerator.hh index a039144..e2064c7 100644 --- a/FixedOrderGen/include/EventGenerator.hh +++ b/FixedOrderGen/include/EventGenerator.hh @@ -1,56 +1,57 @@ #pragma once -#include "HEJ/PDF.hh" #include "HEJ/MatrixElement.hh" +#include "HEJ/optional.hh" +#include "HEJ/PDF.hh" #include "HEJ/RNG.hh" +#include "Beam.hh" #include "JetParameters.hh" +#include "ParticleProperties.hh" #include "Process.hh" -#include "Beam.hh" #include "Status.hh" -#include "ParticleProperties.hh" namespace HEJ{ class Event; class HiggsCouplingSettings; class ScaleGenerator; } //! Namespace for HEJ Fixed Order Generator namespace HEJFOG{ class EventGenerator{ public: EventGenerator( Process process, Beam beam, HEJ::ScaleGenerator scale_gen, JetParameters jets, int pdf_id, double subl_change, unsigned int subl_channels, ParticlesPropMap particles_properties, HEJ::HiggsCouplingSettings Higgs_coupling, HEJ::RNG & ran ); - HEJ::Event gen_event(); + HEJ::optional gen_event(); Status status() const { return status_; } private: HEJ::PDF pdf_; HEJ::MatrixElement ME_; HEJ::ScaleGenerator scale_gen_; Process process_; JetParameters jets_; Beam beam_; Status status_; double subl_change_; unsigned int subl_channels_; ParticlesPropMap particles_properties_; std::reference_wrapper ran_; }; } diff --git a/FixedOrderGen/src/EventGenerator.cc b/FixedOrderGen/src/EventGenerator.cc index eb99958..8762c35 100644 --- a/FixedOrderGen/src/EventGenerator.cc +++ b/FixedOrderGen/src/EventGenerator.cc @@ -1,80 +1,80 @@ #include "EventGenerator.hh" #include "Process.hh" #include "Beam.hh" #include "JetParameters.hh" #include "PhaseSpacePoint.hh" #include "HEJ/Event.hh" #include "HEJ/config.hh" namespace HEJFOG{ EventGenerator::EventGenerator( Process process, Beam beam, HEJ::ScaleGenerator scale_gen, JetParameters jets, int pdf_id, double subl_change, unsigned int subl_channels, ParticlesPropMap particles_properties, HEJ::HiggsCouplingSettings Higgs_coupling, HEJ::RNG & ran ): pdf_{pdf_id, beam.particles[0], beam.particles[1]}, ME_{ [this](double mu){ return pdf_.Halphas(mu); }, HEJ::MatrixElementConfig{ false, std::move(Higgs_coupling) } }, scale_gen_{std::move(scale_gen)}, process_{std::move(process)}, jets_{std::move(jets)}, beam_{std::move(beam)}, subl_change_{subl_change}, subl_channels_{subl_channels}, particles_properties_{std::move(particles_properties)}, ran_{ran} { } - HEJ::Event EventGenerator::gen_event(){ + HEJ::optional EventGenerator::gen_event(){ HEJFOG::PhaseSpacePoint psp{ process_, jets_, pdf_, beam_.energy, subl_change_, subl_channels_, particles_properties_, ran_ }; status_ = psp.status(); if(status_ != good) return {}; HEJ::Event ev = scale_gen_( HEJ::Event{ to_UnclusteredEvent(std::move(psp)), jets_.def, jets_.min_pt } ); const double shat = HEJ::shat(ev); const double xa = (ev.incoming()[0].E()-ev.incoming()[0].pz())/(2.*beam_.energy); const double xb = (ev.incoming()[1].E()+ev.incoming()[1].pz())/(2.*beam_.energy); // evaluate matrix element on this point const auto ME_weights = ME_.tree(ev); ev.central().weight *= ME_weights.central/(shat*shat); ev.central().weight *= pdf_.pdfpt(0,xa,ev.central().muf, ev.incoming()[0].type); ev.central().weight *= pdf_.pdfpt(0,xb,ev.central().muf, ev.incoming()[1].type); for(size_t i = 0; i < ev.variations().size(); ++i){ auto & var = ev.variations(i); var.weight *= ME_weights.variations[i]/(shat*shat); var.weight *= pdf_.pdfpt(0,xa,var.muf, ev.incoming()[0].type); var.weight *= pdf_.pdfpt(0,xb,var.muf, ev.incoming()[1].type); } return ev; } } diff --git a/FixedOrderGen/src/main.cc b/FixedOrderGen/src/main.cc index 6d55edb..7f4b583 100644 --- a/FixedOrderGen/src/main.cc +++ b/FixedOrderGen/src/main.cc @@ -1,233 +1,235 @@ /** * Name: main.cc * Authors: Jeppe R. Andersen */ #include #include #include #include #include #include #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" " / __/ / /> 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); } } void rescale_weights(HEJ::Event & ev, double factor) { ev.central().weight *= factor; for(auto & var: ev.variations()){ var.weight *= factor; } } 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 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(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(1.); heprup.XERRUP=std::vector(1.); heprup.LPRUP=std::vector{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 unweighter{}; std::map status_counter; std::vector 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 warmup_progress{std::cout, warmup_events}; for(; events.size() < warmup_events; ++trials){ auto ev = generator.gen_event(); ++status_counter[generator.status()]; - if(generator.status() == HEJFOG::good && analysis->pass_cuts(ev, ev)) { - events.emplace_back(std::move(ev)); + 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 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(events.size())/config.events; const std::chrono::duration remaining_time = (warmup_end- warmup_start)*(1./completion - 1); const auto finish = clock::to_time_t( std::chrono::time_point_cast(warmup_end + remaining_time) ); std::cout << "Generated " << events.size() << "/" << config.events << " events (" << static_cast(std::round(100*completion)) << "%)\n" << "Estimated remaining generation time: " << remaining_time.count() << " seconds (" << std::put_time(std::localtime(&finish), "%c") << ")\n\n"; } HEJ::ProgressBar progress{std::cout, config.events}; progress.increment(events.size()); events.reserve(config.events); for(; events.size() < static_cast(config.events); ++trials){ auto ev = generator.gen_event(); ++status_counter[generator.status()]; - if(generator.status() == HEJFOG::good && analysis->pass_cuts(ev, ev)) { + 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)); if(! unweighted) continue; - ev = std::move(*unweighted); + ev = std::move(unweighted); } - events.emplace_back(std::move(ev)); + events.emplace_back(std::move(*ev)); ++progress; } } std::cout << std::endl; HEJ::CrossSectionAccumulator xs; for(auto & ev: events){ rescale_weights(ev, invGeV2_to_pb/trials); analysis->fill(ev, ev); writer.write(ev); xs.fill(ev); } analysis->finalise(); const std::chrono::duration 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(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"; } } diff --git a/FixedOrderGen/t/2j.cc b/FixedOrderGen/t/2j.cc index 80bdff3..ea68d04 100644 --- a/FixedOrderGen/t/2j.cc +++ b/FixedOrderGen/t/2j.cc @@ -1,59 +1,60 @@ #ifdef NDEBUG #undef NDEBUG #endif #include #include #include #include #include "config.hh" #include "EventGenerator.hh" #include "HEJ/Mixmax.hh" #include "HEJ/Event.hh" #include "HEJ/PDF.hh" #include "HEJ/MatrixElement.hh" using namespace HEJFOG; int main(){ constexpr double invGeV2_to_pb = 389379292.; constexpr double xs_ref = 86.42031848*1e6; //calculated with "combined" HEJ svn r3480 auto config = load_config("config_2j.yml"); HEJ::Mixmax ran{}; HEJFOG::EventGenerator generator{ config.process, config.beam, HEJ::ScaleGenerator{ config.scales.base, config.scales.factors, config.scales.max_ratio }, config.jets, config.pdf_id, config.subleading_fraction, config.subleading_channels, config.particles_properties, config.Higgs_coupling, ran }; double xs = 0., xs_err = 0.; for (int trials = 0; trials < config.events; ++trials){ auto ev = generator.gen_event(); if(generator.status() != good) continue; - ev.central().weight *= invGeV2_to_pb; - ev.central().weight /= config.events; + assert(ev); + ev->central().weight *= invGeV2_to_pb; + ev->central().weight /= config.events; - xs += ev.central().weight; - xs_err += ev.central().weight*ev.central().weight; + xs += ev->central().weight; + xs_err += ev->central().weight*ev->central().weight; } xs_err = std::sqrt(xs_err); std::cout << xs_ref << " ~ " << xs << " +- " << xs_err << '\n'; assert(std::abs(xs - xs_ref) < 3*xs_err); assert(xs_err < 0.01*xs); } diff --git a/FixedOrderGen/t/4j.cc b/FixedOrderGen/t/4j.cc index 5b4e1b4..a597242 100644 --- a/FixedOrderGen/t/4j.cc +++ b/FixedOrderGen/t/4j.cc @@ -1,60 +1,61 @@ #ifdef NDEBUG #undef NDEBUG #endif #include #include #include #include #include "config.hh" #include "EventGenerator.hh" #include "HEJ/Mixmax.hh" #include "HEJ/Event.hh" #include "HEJ/PDF.hh" #include "HEJ/MatrixElement.hh" using namespace HEJFOG; int main(){ constexpr double invGeV2_to_pb = 389379292.; constexpr double xs_ref = 0.81063619*1e6; //calculated with "combined" HEJ svn r3480 auto config = load_config("config_2j.yml"); config.process.njets = 4; HEJ::Mixmax ran{}; HEJFOG::EventGenerator generator{ config.process, config.beam, HEJ::ScaleGenerator{ config.scales.base, config.scales.factors, config.scales.max_ratio }, config.jets, config.pdf_id, config.subleading_fraction, config.subleading_channels, config.particles_properties, config.Higgs_coupling, ran }; double xs = 0., xs_err = 0.; for (int trials = 0; trials < config.events; ++trials){ auto ev = generator.gen_event(); if(generator.status() != good) continue; - ev.central().weight *= invGeV2_to_pb; - ev.central().weight /= config.events; + assert(ev); + ev->central().weight *= invGeV2_to_pb; + ev->central().weight /= config.events; - xs += ev.central().weight; - xs_err += ev.central().weight*ev.central().weight; + xs += ev->central().weight; + xs_err += ev->central().weight*ev->central().weight; } xs_err = std::sqrt(xs_err); std::cout << xs_ref << " ~ " << xs << " +- " << xs_err << '\n'; assert(std::abs(xs - xs_ref) < 3*xs_err); assert(xs_err < 0.03*xs); } diff --git a/FixedOrderGen/t/h_2j.cc b/FixedOrderGen/t/h_2j.cc index 3e3135e..4ad315d 100644 --- a/FixedOrderGen/t/h_2j.cc +++ b/FixedOrderGen/t/h_2j.cc @@ -1,67 +1,68 @@ #ifdef NDEBUG #undef NDEBUG #endif #include #include #include #include #include "config.hh" #include "EventGenerator.hh" #include "HEJ/Ranlux64.hh" #include "HEJ/Event.hh" #include "HEJ/PDF.hh" #include "HEJ/MatrixElement.hh" using namespace HEJFOG; int main(){ constexpr double invGeV2_to_pb = 389379292.; constexpr double xs_ref = 2.04928; // +- 0.00377252 //calculated with HEJ revision 9570e3809613272ac4b8bf3236279ba23cf64d20 auto config = load_config("config_h_2j.yml"); HEJ::Ranlux64 ran{}; HEJFOG::EventGenerator generator{ config.process, config.beam, HEJ::ScaleGenerator{ config.scales.base, config.scales.factors, config.scales.max_ratio }, config.jets, config.pdf_id, config.subleading_fraction, config.subleading_channels, config.particles_properties, config.Higgs_coupling, ran }; double xs = 0., xs_err = 0.; for (int trials = 0; trials < config.events; ++trials){ auto ev = generator.gen_event(); if(generator.status() != good) continue; - ev.central().weight *= invGeV2_to_pb; - ev.central().weight /= config.events; + assert(ev); + ev->central().weight *= invGeV2_to_pb; + ev->central().weight /= config.events; const auto the_Higgs = std::find_if( - begin(ev.outgoing()), end(ev.outgoing()), + begin(ev->outgoing()), end(ev->outgoing()), [](HEJ::Particle const & p){ return p.type == HEJ::ParticleID::h; } ); - assert(the_Higgs != end(ev.outgoing())); + assert(the_Higgs != end(ev->outgoing())); if(std::abs(the_Higgs->rapidity()) > 5.) continue; - xs += ev.central().weight; - xs_err += ev.central().weight*ev.central().weight; + xs += ev->central().weight; + xs_err += ev->central().weight*ev->central().weight; } xs_err = std::sqrt(xs_err); std::cout << xs_ref << " ~ " << xs << " +- " << xs_err << std::endl; assert(std::abs(xs - xs_ref) < 3*xs_err); assert(xs_err < 0.01*xs); } diff --git a/FixedOrderGen/t/h_2j_decay.cc b/FixedOrderGen/t/h_2j_decay.cc index 8454b7b..e404707 100644 --- a/FixedOrderGen/t/h_2j_decay.cc +++ b/FixedOrderGen/t/h_2j_decay.cc @@ -1,86 +1,87 @@ #ifdef NDEBUG #undef NDEBUG #endif #include #include #include #include #include "config.hh" #include "EventGenerator.hh" #include "HEJ/Event.hh" #include "HEJ/MatrixElement.hh" #include "HEJ/Particle.hh" #include "HEJ/PDF.hh" #include "HEJ/Ranlux64.hh" #include "HEJ/utility.hh" using namespace HEJFOG; bool pass_dR_cut( std::vector const & jets, std::vector const & photons ){ constexpr double delta_R_min = 0.7; for(auto const & jet: jets){ for(auto const & photon: photons){ if(jet.delta_R(photon.p) < delta_R_min) return false; } } return true; } int main(){ constexpr double invGeV2_to_pb = 389379292.; constexpr double xs_ref = 0.00429198; // +- 1.0488e-05 //calculated with HEJ revision 9570e3809613272ac4b8bf3236279ba23cf64d20 auto config = load_config("config_h_2j_decay.yml"); HEJ::Ranlux64 ran{}; HEJFOG::EventGenerator generator{ config.process, config.beam, HEJ::ScaleGenerator{ config.scales.base, config.scales.factors, config.scales.max_ratio }, config.jets, config.pdf_id, config.subleading_fraction, config.subleading_channels, config.particles_properties, config.Higgs_coupling, ran }; double xs = 0., xs_err = 0.; for (int trials = 0; trials < config.events; ++trials){ auto ev = generator.gen_event(); if(generator.status() != good) continue; - assert(ev.decays().size() == 1); - const auto decay = begin(ev.decays()); - assert(ev.outgoing().size() > decay->first); - const auto & the_Higgs = ev.outgoing()[decay->first]; + assert(ev); + assert(ev->decays().size() == 1); + const auto decay = begin(ev->decays()); + assert(ev->outgoing().size() > decay->first); + const auto & the_Higgs = ev->outgoing()[decay->first]; assert(the_Higgs.type == HEJ::pid::Higgs); assert(decay->second.size() == 2); auto const & gamma = decay->second; assert(gamma[0].type == HEJ::pid::photon); assert(gamma[1].type == HEJ::pid::photon); assert(HEJ::nearby_ep(gamma[0].p + gamma[1].p, the_Higgs.p, 1e-6)); - if(!pass_dR_cut(ev.jets(), gamma)) continue; - ev.central().weight *= invGeV2_to_pb; - ev.central().weight /= config.events; + if(!pass_dR_cut(ev->jets(), gamma)) continue; + ev->central().weight *= invGeV2_to_pb; + ev->central().weight /= config.events; - xs += ev.central().weight; - xs_err += ev.central().weight*ev.central().weight; + xs += ev->central().weight; + xs_err += ev->central().weight*ev->central().weight; } xs_err = std::sqrt(xs_err); std::cout << xs_ref << " ~ " << xs << " +- " << xs_err << std::endl; assert(std::abs(xs - xs_ref) < 3*xs_err); assert(xs_err < 0.012*xs); } diff --git a/FixedOrderGen/t/h_3j.cc b/FixedOrderGen/t/h_3j.cc index cad089a..bb2baff 100644 --- a/FixedOrderGen/t/h_3j.cc +++ b/FixedOrderGen/t/h_3j.cc @@ -1,68 +1,69 @@ #ifdef NDEBUG #undef NDEBUG #endif #include #include #include #include #include "config.hh" #include "EventGenerator.hh" #include "HEJ/Ranlux64.hh" #include "HEJ/Event.hh" #include "HEJ/MatrixElement.hh" #include "HEJ/PDF.hh" using namespace HEJFOG; int main(){ constexpr double invGeV2_to_pb = 389379292.; constexpr double xs_ref = 1.07807; // +- 0.0071 //calculated with HEJ revision 93efdc851b02a907a6fcc63956387f9f4c1111c2 +1 auto config = load_config("config_h_2j.yml"); config.process.njets = 3; HEJ::Ranlux64 ran{}; HEJFOG::EventGenerator generator{ config.process, config.beam, HEJ::ScaleGenerator{ config.scales.base, config.scales.factors, config.scales.max_ratio }, config.jets, config.pdf_id, config.subleading_fraction, config.subleading_channels, config.particles_properties, config.Higgs_coupling, ran }; double xs = 0., xs_err = 0.; for (int trials = 0; trials < config.events; ++trials){ auto ev = generator.gen_event(); if(generator.status() != good) continue; - ev.central().weight *= invGeV2_to_pb; - ev.central().weight /= config.events; + assert(ev); + ev->central().weight *= invGeV2_to_pb; + ev->central().weight /= config.events; const auto the_Higgs = std::find_if( - begin(ev.outgoing()), end(ev.outgoing()), + begin(ev->outgoing()), end(ev->outgoing()), [](HEJ::Particle const & p){ return p.type == HEJ::ParticleID::h; } ); - assert(the_Higgs != end(ev.outgoing())); + assert(the_Higgs != end(ev->outgoing())); if(std::abs(the_Higgs->rapidity()) > 5.) continue; - xs += ev.central().weight; - xs_err += ev.central().weight*ev.central().weight; + xs += ev->central().weight; + xs_err += ev->central().weight*ev->central().weight; } xs_err = std::sqrt(xs_err); std::cout << xs_ref << " ~ " << xs << " +- " << xs_err << std::endl; assert(std::abs(xs - xs_ref) < 3*xs_err); assert(xs_err < 0.02*xs); } diff --git a/FixedOrderGen/t/h_3j_uno1.cc b/FixedOrderGen/t/h_3j_uno1.cc index 4c0f6f3..c100dcb 100644 --- a/FixedOrderGen/t/h_3j_uno1.cc +++ b/FixedOrderGen/t/h_3j_uno1.cc @@ -1,72 +1,73 @@ #ifdef NDEBUG #undef NDEBUG #endif // check that adding uno emissions doesn't change the FKL cross section #include #include #include #include #include "config.hh" #include "EventGenerator.hh" #include "HEJ/Ranlux64.hh" #include "Subleading.hh" #include "HEJ/Event.hh" #include "HEJ/MatrixElement.hh" #include "HEJ/PDF.hh" using namespace HEJFOG; int main(){ constexpr double invGeV2_to_pb = 389379292.; constexpr double xs_ref = 0.0243548; // +- 0.000119862 //calculated with HEJ revision 9570e3809613272ac4b8bf3236279ba23cf64d20 auto config = load_config("config_h_2j.yml"); config.process.njets = 3; config.process.incoming = {HEJ::pid::u, HEJ::pid::u}; config.subleading_channels = HEJFOG::Subleading::uno; HEJ::Ranlux64 ran{}; HEJFOG::EventGenerator generator{ config.process, config.beam, HEJ::ScaleGenerator{ config.scales.base, config.scales.factors, config.scales.max_ratio }, config.jets, config.pdf_id, config.subleading_fraction, config.subleading_channels, config.particles_properties, config.Higgs_coupling, ran }; double xs = 0., xs_err = 0.; int uno_found = 0; for (int trials = 0; trials < config.events; ++trials){ auto ev = generator.gen_event(); if(generator.status() != good) continue; - if(ev.type() != HEJ::event_type::FKL){ + assert(ev); + if(ev->type() != HEJ::event_type::FKL){ ++uno_found; continue; } - ev.central().weight *= invGeV2_to_pb; - ev.central().weight /= config.events; + ev->central().weight *= invGeV2_to_pb; + ev->central().weight /= config.events; - xs += ev.central().weight; - xs_err += ev.central().weight*ev.central().weight; + xs += ev->central().weight; + xs_err += ev->central().weight*ev->central().weight; } xs_err = std::sqrt(xs_err); std::cout << xs_ref << " ~ " << xs << " +- " << xs_err << '\n'; std::cout << uno_found << " events with unordered emission" << std::endl; assert(uno_found > 0); assert(std::abs(xs - xs_ref) < 3*xs_err); assert(xs_err < 0.05*xs); } diff --git a/FixedOrderGen/t/h_3j_uno2.cc b/FixedOrderGen/t/h_3j_uno2.cc index 5c6c409..4dc4770 100644 --- a/FixedOrderGen/t/h_3j_uno2.cc +++ b/FixedOrderGen/t/h_3j_uno2.cc @@ -1,67 +1,68 @@ #ifdef NDEBUG #undef NDEBUG #endif // check uno cross section #include #include #include #include #include "config.hh" #include "EventGenerator.hh" #include "HEJ/Ranlux64.hh" #include "Subleading.hh" #include "HEJ/Event.hh" #include "HEJ/MatrixElement.hh" #include "HEJ/PDF.hh" using namespace HEJFOG; int main(){ constexpr double invGeV2_to_pb = 389379292.; constexpr double xs_ref = 0.00347538; // +- 3.85875e-05 //calculated with HEJ revision 9570e3809613272ac4b8bf3236279ba23cf64d20 auto config = load_config("config_h_2j.yml"); config.process.njets = 3; config.process.incoming = {HEJ::pid::u, HEJ::pid::u}; config.subleading_fraction = 1.; config.subleading_channels = HEJFOG::Subleading::uno; HEJ::Ranlux64 ran{}; HEJFOG::EventGenerator generator{ config.process, config.beam, HEJ::ScaleGenerator{ config.scales.base, config.scales.factors, config.scales.max_ratio }, config.jets, config.pdf_id, config.subleading_fraction, config.subleading_channels, config.particles_properties, config.Higgs_coupling, ran }; double xs = 0., xs_err = 0.; for (int trials = 0; trials < config.events; ++trials){ auto ev = generator.gen_event(); - if(ev.type() == HEJ::event_type::FKL) continue; if(generator.status() != good) continue; - ev.central().weight *= invGeV2_to_pb; - ev.central().weight /= config.events; + assert(ev); + if(ev->type() == HEJ::event_type::FKL) continue; + ev->central().weight *= invGeV2_to_pb; + ev->central().weight /= config.events; - xs += ev.central().weight; - xs_err += ev.central().weight*ev.central().weight; + xs += ev->central().weight; + xs_err += ev->central().weight*ev->central().weight; } xs_err = std::sqrt(xs_err); std::cout << xs_ref << " ~ " << xs << " +- " << xs_err << std::endl; assert(std::abs(xs - xs_ref) < 3*xs_err); assert(xs_err < 0.05*xs); } diff --git a/FixedOrderGen/t/h_5j.cc b/FixedOrderGen/t/h_5j.cc index 6aa3d40..a1b60e8 100644 --- a/FixedOrderGen/t/h_5j.cc +++ b/FixedOrderGen/t/h_5j.cc @@ -1,63 +1,64 @@ #ifdef NDEBUG #undef NDEBUG #endif // This is a regression test // the reference cross section has not been checked against any other program #include #include #include #include #include "config.hh" #include "EventGenerator.hh" #include "HEJ/Ranlux64.hh" #include "HEJ/Event.hh" #include "HEJ/MatrixElement.hh" #include "HEJ/PDF.hh" using namespace HEJFOG; int main(){ constexpr double invGeV2_to_pb = 389379292.; constexpr double xs_ref = 0.252273; // +- 0.00657742 //calculated with HEJ revision 9570e3809613272ac4b8bf3236279ba23cf64d20 auto config = load_config("config_h_2j.yml"); config.process.njets = 5; HEJ::Ranlux64 ran{}; HEJFOG::EventGenerator generator{ config.process, config.beam, HEJ::ScaleGenerator{ config.scales.base, config.scales.factors, config.scales.max_ratio }, config.jets, config.pdf_id, config.subleading_fraction, config.subleading_channels, config.particles_properties, config.Higgs_coupling, ran }; double xs = 0., xs_err = 0.; for (int trials = 0; trials < config.events; ++trials){ auto ev = generator.gen_event(); if(generator.status() != good) continue; - ev.central().weight *= invGeV2_to_pb; - ev.central().weight /= config.events; + assert(ev); + ev->central().weight *= invGeV2_to_pb; + ev->central().weight /= config.events; - xs += ev.central().weight; - xs_err += ev.central().weight*ev.central().weight; + xs += ev->central().weight; + xs_err += ev->central().weight*ev->central().weight; } xs_err = std::sqrt(xs_err); std::cout << xs_ref << " ~ " << xs << " +- " << xs_err << std::endl; assert(std::abs(xs - xs_ref) < 3*xs_err); assert(xs_err < 0.06*xs); } diff --git a/include/HEJ/Event.hh b/include/HEJ/Event.hh index 6b59620..e87bb10 100644 --- a/include/HEJ/Event.hh +++ b/include/HEJ/Event.hh @@ -1,199 +1,199 @@ /** \file * \brief Declares the Event class and helpers * * \authors Jeppe Andersen, Tuomas Hapola, Marian Heil, Andreas Maier, Jennifer Smillie * \date 2019 * \copyright GPLv2 or later */ #pragma once #include #include #include #include #include #include "HEJ/event_types.hh" #include "HEJ/Particle.hh" #include "fastjet/ClusterSequence.hh" namespace LHEF{ class HEPEUP; class HEPRUP; } namespace fastjet{ class JetDefinition; } namespace HEJ{ struct ParameterDescription; //! Event parameters struct EventParameters{ double mur; /**< Value of the Renormalisation Scale */ double muf; /**< Value of the Factorisation Scale */ double weight; /**< Event Weight */ //! Optional description std::shared_ptr description = nullptr; }; //! Description of event parameters struct ParameterDescription { //! Name of central scale choice (e.g. "H_T/2") std::string scale_name; //! Actual renormalisation scale divided by central scale double mur_factor; //! Actual factorisation scale divided by central scale double muf_factor; ParameterDescription() = default; ParameterDescription( std::string scale_name, double mur_factor, double muf_factor ): scale_name{scale_name}, mur_factor{mur_factor}, muf_factor{muf_factor} {}; }; //! An event before jet clustering struct UnclusteredEvent{ //! Default Constructor UnclusteredEvent() = default; //! Constructor from LesHouches event information UnclusteredEvent(LHEF::HEPEUP const & hepeup); std::array incoming; /**< Incoming Particles */ std::vector outgoing; /**< Outgoing Particles */ //! Particle decays in the format {outgoing index, decay products} std::unordered_map> decays; //! Central parameter (e.g. scale) choice EventParameters central; std::vector variations; /**< For parameter variation */ }; /** An event with clustered jets * * This is the main HEJ 2 event class. * It contains kinematic information including jet clustering, * parameter (e.g. scale) settings and the event weight. */ class Event{ public: //! Default Event Constructor - Event() = default; + Event() = delete; //! Event Constructor adding jet clustering to an unclustered event Event( UnclusteredEvent ev, fastjet::JetDefinition const & jet_def, double min_jet_pt ); //! The jets formed by the outgoing partons std::vector jets() const; //! The corresponding event before jet clustering UnclusteredEvent const & unclustered() const { return ev_; } //! Central parameter choice EventParameters const & central() const{ return ev_.central; } //! Central parameter choice EventParameters & central(){ return ev_.central; } //! Incoming particles std::array const & incoming() const{ return ev_.incoming; } //! Outgoing particles std::vector const & outgoing() const{ return ev_.outgoing; } //! Particle decays /** * The key in the returned map corresponds to the index in the * vector returned by outgoing() */ std::unordered_map> const & decays() const{ return ev_.decays; } //! Parameter (scale) variations std::vector const & variations() const{ return ev_.variations; } //! Parameter (scale) variations std::vector & variations(){ return ev_.variations; } //! Parameter (scale) variation /** * @param i Index of the requested variation */ EventParameters const & variations(size_t i) const{ return ev_.variations[i]; } //! Parameter (scale) variation /** * @param i Index of the requested variation */ EventParameters & variations(size_t i){ return ev_.variations[i]; } //! Indices of the jets the outgoing partons belong to /** * @param jets Jets to be tested * @returns A vector containing, for each outgoing parton, * the index in the vector of jets the considered parton * belongs to. If the parton is not inside any of the * passed jets, the corresponding index is set to -1. */ std::vector particle_jet_indices( std::vector const & jets ) const{ return cs_.particle_jet_indices(jets); } //! Jet definition used for clustering fastjet::JetDefinition const & jet_def() const{ return cs_.jet_def(); } //! Minimum jet transverse momentum double min_jet_pt() const{ return min_jet_pt_; } //! Event type event_type::EventType type() const{ return type_; } private: UnclusteredEvent ev_; fastjet::ClusterSequence cs_; double min_jet_pt_; event_type::EventType type_; }; //! Square of the partonic centre-of-mass energy \f$\hat{s}\f$ double shat(Event const & ev); //! Convert an event to a LHEF::HEPEUP LHEF::HEPEUP to_HEPEUP(Event const & event, LHEF::HEPRUP *); } diff --git a/t/test_ME_generic.cc b/t/test_ME_generic.cc index 00e19f0..75090d4 100644 --- a/t/test_ME_generic.cc +++ b/t/test_ME_generic.cc @@ -1,104 +1,107 @@ // Generic tester for the ME for a given set of PSP // reference weights and PSP (as LHE file) have to be given as _individual_ files #include #include "LHEF/LHEF.h" #include "HEJ/MatrixElement.hh" #include "HEJ/Event.hh" #include "HEJ/YAMLreader.hh" #include "HEJ/stream.hh" constexpr double alpha_s = 0.118; constexpr double ep = 1e-6; 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; } } int main(int argn, char** argv){ if(argn != 4 && argn != 5){ std::cerr << "\n# Usage:\n."<< argv[0] <<" config.yml ME_weights input_file.lhe\n\n"; return EXIT_FAILURE; } bool OUTPUT_MODE = false; if(argn == 5 && std::string("OUTPUT")==std::string(argv[4])) OUTPUT_MODE = true; const HEJ::Config config = HEJ::load_config(argv[1]); std::fstream wgt_file; if ( OUTPUT_MODE ) { std::cout << "_______________________USING OUTPUT MODE!_______________________" << std::endl; wgt_file.open(argv[2], std::fstream::out); wgt_file.precision(10); } else { wgt_file.open(argv[2], std::fstream::in); } HEJ::istream in{argv[3]}; LHEF::Reader reader{in}; HEJ::MatrixElement ME{ [](double){ return alpha_s; }, HEJ::to_MatrixElementConfig(config) }; double max_ratio = 0.; size_t idx_max_ratio = 0; - HEJ::Event ev_max_ratio; + HEJ::Event ev_max_ratio(HEJ::UnclusteredEvent(), + config.resummation_jets.def, + config.resummation_jets.min_pt + ); double av_ratio = 0; size_t i = 0; while(reader.readEvent()){ ++i; HEJ::Event event{ HEJ::UnclusteredEvent{reader.hepeup}, config.resummation_jets.def, config.resummation_jets.min_pt }; const double our_ME = ME.tree(event).central; if ( OUTPUT_MODE ) { wgt_file << our_ME << std::endl; } else { std::string line; if(!std::getline(wgt_file,line)) break; const double ref_ME = std::stod(line); const double diff = std::abs(our_ME/ref_ME-1.); av_ratio+=diff; if( diff > max_ratio ) { max_ratio = diff; idx_max_ratio = i; ev_max_ratio = event; } if( diff > ep ){ size_t precision(std::cout.precision()); std::cout.precision(16); std::cout<< "Large difference in PSP " << i << "\nis: "< difference: " << diff << std::endl; std::cout.precision(precision); dump(event); return EXIT_FAILURE; } } } wgt_file.close(); if ( !OUTPUT_MODE ) { size_t precision(std::cout.precision()); std::cout.precision(16); std::cout << "Avg ratio after " << i << " PSP: " << av_ratio/i << std::endl; std::cout << "maximal ratio at " << idx_max_ratio << ": " << max_ratio << std::endl; std::cout.precision(precision); } return EXIT_SUCCESS; }