diff --git a/include/HEJ/event_types.hh b/include/HEJ/event_types.hh index b53ba36..d3273b8 100644 --- a/include/HEJ/event_types.hh +++ b/include/HEJ/event_types.hh @@ -1,109 +1,109 @@ /** \file * \brief Define different types of events. * * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #pragma once #include #include "HEJ/exceptions.hh" namespace HEJ { //! Namespace for event types namespace event_type { //! Possible event types enum EventType: std::size_t { non_resummable = 0, //!< event configuration not covered by All Order resummation bad_final_state = 1, //!< event with an unsupported final state - no_2_jets = 2, //!< event with less than two jets + not_enough_jets = 2, //!< event with less than two jets FKL = 4, //!< FKL-type event unordered_backward = 8, //!< event with unordered backward emission unordered_forward = 16, //!< event with unordered forward emission extremal_qqxb = 32, //!< event with a backward extremal qqbar extremal_qqxf = 64, //!< event with a forward extremal qqbar central_qqx = 128, //!< event with a central qqbar unob = unordered_backward, //!< alias for unordered_backward unof = unordered_forward, //!< alias for unordered_forward qqxexb = extremal_qqxb, //!< alias for extremal_qqxb qqxexf = extremal_qqxf, //!< alias for extremal_qqxf qqxmid = central_qqx, //!< alias for central_qqx first_type = non_resummable, //!< alias for non_resummable last_type = central_qqx //!< alias for central_qqx }; //! Event type names /** * For example, name(FKL) is the string "FKL" */ inline std::string name(EventType type) { switch(type) { case FKL: return "FKL"; case unordered_backward: return "unordered backward"; case unordered_forward: return "unordered forward"; case extremal_qqxb: return "extremal qqbar backward"; case extremal_qqxf: return "extremal qqbar forward"; case central_qqx: return "central qqbar"; case non_resummable: return "non-resummable"; - case no_2_jets: - return "no 2 jets"; + case not_enough_jets: + return "not enough jets"; case bad_final_state: return "bad final state"; default: throw std::logic_error{"Unreachable"}; } } //! Returns True for a HEJ \ref event_type::EventType "EventType" inline constexpr bool is_resummable(EventType type) { switch(type) { case FKL: case unordered_backward: case unordered_forward: case extremal_qqxb: case extremal_qqxf: case central_qqx: return true; default: return false; } } //! Returns True for an unordered \ref event_type::EventType "EventType" inline constexpr bool is_uno(EventType type) { return type == unordered_backward || type == unordered_forward; } //! Returns True for an extremal_qqx \ref event_type::EventType "EventType" inline constexpr bool is_ex_qqx(EventType type) { return type == extremal_qqxb || type == extremal_qqxf; } //! Returns True for an central_qqx \ref event_type::EventType "EventType" inline constexpr bool is_mid_qqx(EventType type) { return type == central_qqx; } //! Returns True for any qqx event \ref event_type::EventType "EventType" inline constexpr bool is_qqx(EventType type) { return is_ex_qqx(type) || is_mid_qqx(type); } } // namespace event_type } // namespace HEJ diff --git a/src/Event.cc b/src/Event.cc index cac0fa1..b70a8a2 100644 --- a/src/Event.cc +++ b/src/Event.cc @@ -1,1175 +1,1219 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #include "HEJ/Event.hh" #include #include #include #include #include #include #include #include #include #include #include "fastjet/ClusterSequence.hh" #include "fastjet/JetDefinition.hh" #include "fastjet/PseudoJet.hh" #include "LHEF/LHEF.h" #include "HEJ/Constants.hh" #include "HEJ/PDG_codes.hh" #include "HEJ/RNG.hh" #include "HEJ/exceptions.hh" #include "HEJ/optional.hh" namespace HEJ { namespace { using std::size_t; //! LHE status codes namespace lhe_status { enum Status: int { in = -1, decay = 2, out = 1, }; } using LHE_Status = lhe_status::Status; //! true if leptonic W decay bool valid_W_decay( int const w_type, // sign of W std::vector const & decays ){ if(decays.size() != 2) // no 1->2 decay return false; const int pidsum = decays[0].type + decays[1].type; if( std::abs(pidsum) != 1 || pidsum != w_type ) // correct charge return false; // leptonic decay (only check first, second follows from pidsum) if( w_type == 1 ) // W+ return is_antilepton(decays[0]) || is_neutrino(decays[0]); // W- return is_lepton(decays[0]) || is_antineutrino(decays[0]); } //! true for Z decay to charged leptons bool valid_Z_decay(std::vector const & decays){ if(decays.size() != 2) // no 1->2 decay return false; const int pidsum = decays[0].type + decays[1].type; if( std::abs(pidsum) != 0 ) // correct charge return false; // leptonic decay (only check first, second follows from pidsum) return is_anylepton(decays[0]) && !is_anyneutrino(decays[0]); } /// @name helper functions to determine event type //@{ /** * \brief check if final state valid for HEJ * * check if there is at most one photon, W, H, Z in the final state * and all the rest are quarks or gluons */ bool final_state_ok(Event const & ev){ std::vector const & outgoing = ev.outgoing(); if(ev.decays().size() > 1) // at most one decay return false; bool has_AWZH_boson = false; for( size_t i=0; ifirst != i ) return false; if( !valid_W_decay(out.type>0?+1:-1, ev.decays().cbegin()->second) ) return false; } // valid decay for Z if(out.type == ParticleID::Z_photon_mix){ // exactly 1 decay if( ev.decays().size() != 1 || ev.decays().cbegin()->first != i ) return false; if( !valid_Z_decay(ev.decays().cbegin()->second) ) return false; } } else if(! is_parton(out.type)) return false; } return true; } /** * returns all EventTypes implemented in HEJ */ size_t implemented_types(std::vector const & bosons){ using namespace event_type; if(bosons.empty()) return FKL | unob | unof | qqxexb | qqxexf | qqxmid; if(bosons.size()>1) return non_resummable; // multi boson switch (bosons[0].type) { case ParticleID::Wp: case ParticleID::Wm: return FKL | unob | unof | qqxexb | qqxexf | qqxmid; case ParticleID::Z_photon_mix: return FKL | unob | unof; case ParticleID::h: return FKL | unob | unof; default: return non_resummable; } } /** * \brief function which determines if type change is consistent with Wp emission. * @param in incoming Particle id * @param out outgoing Particle id * @param qqx Current both incoming/both outgoing? * * \see is_Wm_Change */ bool is_Wp_Change(ParticleID in, ParticleID out, bool qqx){ using namespace pid; if(!qqx && (in==d_bar || in==u || in==s_bar || in==c)) return out == (in-1); if( qqx && (in==d || in==u_bar || in==s || in==c_bar)) return out == -(in+1); return false; } /** * \brief function which determines if type change is consistent with Wm emission. * @param in incoming Particle id * @param out outgoing Particle id * @param qqx Current both incoming/both outgoing? * * Ensures that change type of quark line is possible by a flavour changing * Wm emission. Allows checking of qqx currents also. */ bool is_Wm_Change(ParticleID in, ParticleID out, bool qqx){ using namespace pid; if(!qqx && (in==d || in==u_bar || in==s || in==c_bar)) return out == (in+1); if( qqx && (in==d_bar || in==u || in==s_bar || in==c)) return out == -(in-1); return false; } /** * \brief checks if particle type remains same from incoming to outgoing * @param in incoming Particle * @param out outgoing Particle * @param qqx Current both incoming/outgoing? */ bool no_flavour_change(ParticleID in, ParticleID out, bool qqx){ const int qqxCurrent = qqx?-1:1; if(std::abs(in)<=pid::top || in==pid::gluon) return (in==out*qqxCurrent); return false; } - bool has_2_jets(Event const & event){ - return event.jets().size() >= 2; + bool has_enough_jets(Event const & event){ + if(event.jets().size() >= 2) return true; + if(event.jets().empty()) return false; + const auto the_higgs = std::find_if( + begin(event.outgoing()), end(event.outgoing()), + [](const auto & particle) { return particle.type == pid::higgs; } + ); + if(the_higgs == end(event.outgoing())) return false; + // so we have Higgs + 1 jet + // check if we can have a g -> Higgs conversion + return event.incoming().front().type == pid::gluon + || event.incoming().back().type == pid::gluon; + } + + bool is_gluon_to_Higgs(const ParticleID in, const ParticleID out) { + return in == pid::gluon && out == pid::Higgs; } /** * \brief check if we have a valid Impact factor * @param in incoming Particle * @param out outgoing Particle * @param qqx Current both incoming/outgoing? * @param W_change returns +1 if Wp, -1 if Wm, else 0 */ bool is_valid_impact_factor( ParticleID in, ParticleID out, bool qqx, int & W_change ){ - if( no_flavour_change(in, out, qqx) ){ + if( no_flavour_change(in, out, qqx) || is_gluon_to_Higgs(in, out)) { return true; } if( is_Wp_Change(in, out, qqx) ) { W_change+=1; return true; } if( is_Wm_Change(in, out, qqx) ) { W_change-=1; return true; } return false; } + bool is_extremal_higgs_off_quark( + const ParticleID in, + const ParticleID extremal_out, + const ParticleID out + ) { + return in == out && extremal_out == pid::higgs && is_anyquark(in); + } + //! Returns all possible classifications from the impact factors // the beginning points are changed s.t. after the the classification they // point to the beginning of the (potential) FKL chain // sets W_change: + if Wp change // 0 if no change // - if Wm change // This function can be used with forward & backwards iterators template size_t possible_impact_factors( ParticleID incoming_id, // incoming OutIterator & begin_out, OutIterator const & end_out, // outgoing int & W_change, std::vector const & boson, bool const backward // backward? ){ using namespace event_type; assert(boson.size() < 2); // keep track of all states that we don't test size_t not_tested = qqxmid; if(backward) not_tested |= unof | qqxexf; else not_tested |= unob | qqxexb; // Is this LL current? if( is_valid_impact_factor(incoming_id, begin_out->type, false, W_change) ){ ++begin_out; return not_tested | FKL; } + // q -> H q and qbar -> H qbar are technically not LL, + // but we treat them as such anyway + if(is_extremal_higgs_off_quark(incoming_id, begin_out->type, std::next(begin_out)->type)) { + std::advance(begin_out, 2); + return not_tested | FKL; + } // or NLL current? // -> needs two partons in two different jets if( std::distance(begin_out, end_out)>=2 ){ auto next = std::next(begin_out); // Is this unordered emisson? if( incoming_id!=pid::gluon && begin_out->type==pid::gluon ){ if( is_valid_impact_factor( incoming_id, next->type, false, W_change ) ){ // veto Higgs inside uno assert(next!=end_out); if( !boson.empty() && boson.front().type == ParticleID::h ){ if( (backward && boson.front().rapidity() < next->rapidity()) ||(!backward && boson.front().rapidity() > next->rapidity())) return non_resummable; } begin_out = std::next(next); return not_tested | (backward?unob:unof); } } // Is this QQbar? else if( incoming_id==pid::gluon ){ if( is_valid_impact_factor( begin_out->type, next->type, true, W_change ) ){ // veto Higgs inside qqx assert(next!=end_out); if( !boson.empty() && boson.front().type == ParticleID::h ){ if( (backward && boson.front().rapidity() < next->rapidity()) ||(!backward && boson.front().rapidity() > next->rapidity())) return non_resummable; } begin_out = std::next(next); return not_tested | (backward?qqxexb:qqxexf); } } } return non_resummable; } //! Returns all possible classifications from central emissions // the beginning points are changed s.t. after the the classification they // point to the end of the emission chain // sets W_change: + if Wp change // 0 if no change // - if Wm change template size_t possible_central( OutIterator & begin_out, OutIterator const & end_out, int & W_change, std::vector const & boson ){ using namespace event_type; assert(boson.size() < 2); // if we already passed the central chain, // then it is not a valid all-order state if(std::distance(begin_out, end_out) < 0) return non_resummable; // keep track of all states that we don't test size_t possible = unob | unof | qqxexb | qqxexf; - // Find the first non-gluon/non-FKL - while( (begin_out->type==pid::gluon) && (begin_out!=end_out) ){ - ++begin_out; - } + // Find the first quark or antiquark emission + begin_out = std::find_if( + begin_out, end_out, + [](Particle const & p) { return is_anyquark(p); } + ); // end of chain -> FKL if( begin_out==end_out ){ return possible | FKL; } // is this a qqbar-pair? // needs two partons in two separate jets auto next = std::next(begin_out); if( is_valid_impact_factor( begin_out->type, next->type, true, W_change ) ){ // veto Higgs inside qqx if( !boson.empty() && boson.front().type == ParticleID::h && boson.front().rapidity() > begin_out->rapidity() && boson.front().rapidity() < next->rapidity() ){ return non_resummable; } begin_out = std::next(next); - // remaining chain should be pure gluon/FKL - for(; begin_out!=end_out; ++begin_out){ - if(begin_out->type != pid::gluon) return non_resummable; + // remaining chain should be pure FKL (gluon or higgs) + if(std::any_of( + begin_out, end_out, + [](Particle const & p) { return is_anyquark(p); } + )) { + return non_resummable; } return possible | qqxmid; } return non_resummable; } + namespace { + bool is_parton_or_higgs(Particle const & p) { + return is_parton(p) || p.type == pid::higgs; + } + } + /** * \brief Checks for all event types * @param ev Event * @returns Event Type * */ event_type::EventType classify(Event const & ev){ using namespace event_type; - if(! has_2_jets(ev)) - return no_2_jets; + if(! has_enough_jets(ev)) + return not_enough_jets; // currently we can't handle multiple boson states in the ME. So they are // considered "bad_final_state" even though the "classify" could work with // them. if(! final_state_ok(ev)) return bad_final_state; // initialise variables auto const & in = ev.incoming(); // range for current checks - auto begin_out{ev.cbegin_partons()}; - auto end_out{ev.crbegin_partons()}; + auto begin_out = boost::make_filter_iterator( + is_parton_or_higgs, cbegin(ev.outgoing()), cend(ev.outgoing()) + ); + auto rbegin_out = std::make_reverse_iterator( + boost::make_filter_iterator( + is_parton_or_higgs, cend(ev.outgoing()), cend(ev.outgoing()) + ) + ); assert(std::distance(begin(in), end(in)) == 2); - assert(std::distance(begin_out, end_out.base()) >= 2); - assert(std::is_sorted(begin_out, end_out.base(), rapidity_less{})); + assert(std::distance(begin_out, rbegin_out.base()) >= 2); + assert(std::is_sorted(begin_out, rbegin_out.base(), rapidity_less{})); auto const boson{ filter_AWZH_bosons(ev.outgoing()) }; // we only allow one boson through final_state_ok assert(boson.size()<=1); // keep track of potential W couplings, at the end the sum should be 0 int remaining_Wp = 0; int remaining_Wm = 0; if(!boson.empty() && std::abs(boson.front().type) == ParticleID::Wp ){ if(boson.front().type>0) ++remaining_Wp; else ++remaining_Wm; } int W_change = 0; - size_t final_type = ~(no_2_jets | bad_final_state); + size_t final_type = ~(not_enough_jets | bad_final_state); // check forward impact factor final_type &= possible_impact_factors( in.front().type, - begin_out, end_out.base(), + begin_out, rbegin_out.base(), W_change, boson, true ); if( final_type == non_resummable ) return non_resummable; if(W_change>0) remaining_Wp-=W_change; else if(W_change<0) remaining_Wm+=W_change; W_change = 0; // check backward impact factor final_type &= possible_impact_factors( in.back().type, - end_out, std::make_reverse_iterator(begin_out), + rbegin_out, std::make_reverse_iterator(begin_out), W_change, boson, false ); if( final_type == non_resummable ) return non_resummable; if(W_change>0) remaining_Wp-=W_change; else if(W_change<0) remaining_Wm+=W_change; W_change = 0; // check central emissions final_type &= possible_central( - begin_out, end_out.base(), W_change, boson ); + begin_out, rbegin_out.base(), W_change, boson ); if( final_type == non_resummable ) return non_resummable; if(W_change>0) remaining_Wp-=W_change; else if(W_change<0) remaining_Wm+=W_change; // Check whether the right number of Ws are present if( remaining_Wp != 0 || remaining_Wm != 0 ) return non_resummable; // result has to be unique if( (final_type & (final_type-1)) != 0) return non_resummable; // check that each sub processes is implemented // (has to be done at the end) if( (final_type & ~implemented_types(boson)) != 0 ) return non_resummable; return static_cast(final_type); } //@} Particle extract_particle(LHEF::HEPEUP const & hepeup, size_t i){ auto id = static_cast(hepeup.IDUP[i]); auto colour = is_parton(id)?hepeup.ICOLUP[i]:optional(); return { id, { hepeup.PUP[i][0], hepeup.PUP[i][1], hepeup.PUP[i][2], hepeup.PUP[i][3] }, colour }; } bool is_decay_product(std::pair const & mothers){ if(mothers.first == 0) return false; return mothers.second == 0 || mothers.first == mothers.second; } } // namespace Event::EventData::EventData(LHEF::HEPEUP const & hepeup){ parameters.central = EventParameters{ hepeup.scales.mur, hepeup.scales.muf, hepeup.XWGTUP }; size_t in_idx = 0; for (int i = 0; i < hepeup.NUP; ++i) { // skip decay products // we will add them later on, but we have to ensure that // the decayed particle is added before if(is_decay_product(hepeup.MOTHUP[i])) continue; auto particle = extract_particle(hepeup, i); // needed to identify mother particles for decay products particle.p.set_user_index(i+1); if(hepeup.ISTUP[i] == LHE_Status::in){ if(in_idx > incoming.size()) { throw std::invalid_argument{ "Event has too many incoming particles" }; } incoming[in_idx++] = std::move(particle); } else outgoing.emplace_back(std::move(particle)); } // add decay products for (int i = 0; i < hepeup.NUP; ++i) { if(!is_decay_product(hepeup.MOTHUP[i])) continue; const int mother_id = hepeup.MOTHUP[i].first; const auto mother = std::find_if( begin(outgoing), end(outgoing), [mother_id](Particle const & particle){ return particle.p.user_index() == mother_id; } ); if(mother == end(outgoing)){ throw std::invalid_argument{"invalid decay product parent"}; } const int mother_idx = std::distance(begin(outgoing), mother); assert(mother_idx >= 0); decays[mother_idx].emplace_back(extract_particle(hepeup, i)); } } Event::Event( UnclusteredEvent const & ev, fastjet::JetDefinition const & jet_def, double const min_jet_pt ): Event( Event::EventData{ ev.incoming, ev.outgoing, ev.decays, Parameters{ev.central, ev.variations} }.cluster(jet_def, min_jet_pt) ) {} //! @TODO remove in HEJ 2.2.0 UnclusteredEvent::UnclusteredEvent(LHEF::HEPEUP const & hepeup){ Event::EventData const evData{hepeup}; incoming = evData.incoming; outgoing = evData.outgoing; decays = evData.decays; central = evData.parameters.central; variations = evData.parameters.variations; } void Event::EventData::sort(){ // sort particles std::sort( begin(incoming), end(incoming), [](Particle const & o1, Particle const & o2){return o1.p.pz() idx(old_outgoing.size()); std::iota(idx.begin(), idx.end(), 0); std::sort(idx.begin(), idx.end(), [&old_outgoing](size_t i, size_t j){ return old_outgoing[i].rapidity() < old_outgoing[j].rapidity(); }); outgoing.clear(); outgoing.reserve(old_outgoing.size()); for(size_t i: idx) { outgoing.emplace_back(std::move(old_outgoing[i])); } // find decays again if(!decays.empty()){ auto old_decays = std::move(decays); decays.clear(); for(size_t i=0; isecond)); } assert(old_decays.size() == decays.size()); } } namespace { Particle reconstruct_boson(std::vector const & leptons) { Particle decayed_boson; decayed_boson.p = leptons[0].p + leptons[1].p; const int pidsum = leptons[0].type + leptons[1].type; if(pidsum == +1) { assert(is_antilepton(leptons[0])); if(is_antineutrino(leptons[0])) { throw not_implemented{"lepton-flavour violating final state"}; } assert(is_neutrino(leptons[1])); // charged antilepton + neutrino means we had a W+ decayed_boson.type = pid::Wp; } else if(pidsum == -1) { assert(is_antilepton(leptons[0])); if(is_neutrino(leptons[1])) { throw not_implemented{"lepton-flavour violating final state"}; } assert(is_antineutrino(leptons[0])); // charged lepton + antineutrino means we had a W- decayed_boson.type = pid::Wm; } else if(pidsum == 0) { assert(is_anylepton(leptons[0])); if(is_anyneutrino(leptons[0])) { throw not_implemented{"final state with two neutrinos"}; } // charged lepton-antilepton pair means we had a Z/photon decayed_boson.type = pid::Z_photon_mix; } else { throw not_implemented{ "final state with leptons " + name(leptons[0].type) + " and " + name(leptons[1].type) }; } return decayed_boson; } } // namespace void Event::EventData::reconstruct_intermediate() { const auto begin_leptons = std::partition( begin(outgoing), end(outgoing), [](Particle const & p) {return !is_anylepton(p);} ); // We can only reconstruct FS with 2 leptons if(std::distance(begin_leptons, end(outgoing)) != 2) return; std::vector leptons(begin_leptons, end(outgoing)); std::sort( begin(leptons), end(leptons), [](Particle const & p0, Particle const & p1) { assert(is_anylepton(p0) && is_anylepton(p1)); return p0.type < p1.type; } ); // `reconstruct_boson` can throw, it should therefore be called before // changing `outgoing` to allow the user to recover the original EventData auto boson = reconstruct_boson(leptons); outgoing.erase(begin_leptons, end(outgoing)); outgoing.emplace_back(std::move(boson)); decays.emplace(outgoing.size()-1, std::move(leptons)); } Event Event::EventData::cluster( fastjet::JetDefinition const & jet_def, double const min_jet_pt ){ sort(); return Event{ std::move(incoming), std::move(outgoing), std::move(decays), std::move(parameters), jet_def, min_jet_pt }; } Event::Event( std::array && incoming, std::vector && outgoing, std::unordered_map> && decays, Parameters && parameters, fastjet::JetDefinition const & jet_def, double const min_jet_pt ): incoming_{std::move(incoming)}, outgoing_{std::move(outgoing)}, decays_{std::move(decays)}, parameters_{std::move(parameters)}, cs_{ to_PseudoJet( filter_partons(outgoing_) ), jet_def }, min_jet_pt_{min_jet_pt} { jets_ = sorted_by_rapidity(cs_.inclusive_jets(min_jet_pt_)); assert(std::is_sorted(begin(outgoing_), end(outgoing_), rapidity_less{})); type_ = classify(*this); } namespace { //! check that Particles have a reasonable colour bool correct_colour(Particle const & part){ ParticleID id{ part.type }; if(!is_parton(id)) return !part.colour; if(!part.colour) return false; Colour const & col{ *part.colour }; if(is_quark(id)) return col.first != 0 && col.second == 0; if(is_antiquark(id)) return col.first == 0 && col.second != 0; assert(id==ParticleID::gluon); return col.first != 0 && col.second != 0 && col.first != col.second; } //! Connect parton to t-channel colour line & update the line //! returns false if connection not possible template bool try_connect_t(OutIterator const & it_part, Colour & line_colour){ if( line_colour.first == it_part->colour->second ){ line_colour.first = it_part->colour->first; return true; } if( line_colour.second == it_part->colour->first ){ line_colour.second = it_part->colour->second; return true; } return false; } //! Connect parton to u-channel colour line & update the line //! returns false if connection not possible template bool try_connect_u(OutIterator & it_part, Colour & line_colour){ auto it_next = std::next(it_part); if( try_connect_t(it_next, line_colour) && try_connect_t(it_part, line_colour) ){ it_part=it_next; return true; } return false; } } // namespace bool Event::is_leading_colour() const { if( !correct_colour(incoming()[0]) || !correct_colour(incoming()[1]) ) return false; Colour line_colour = *incoming()[0].colour; std::swap(line_colour.first, line_colour.second); // reasonable colour if(!std::all_of(outgoing().cbegin(), outgoing().cend(), correct_colour)) return false; for(auto it_part = cbegin_partons(); it_part!=cend_partons(); ++it_part){ switch (type()) { case event_type::FKL: if( !try_connect_t(it_part, line_colour) ) return false; break; case event_type::unob: case event_type::qqxexb: { if( !try_connect_t(it_part, line_colour) // u-channel only allowed at impact factor && (std::distance(cbegin_partons(), it_part)!=0 || !try_connect_u(it_part, line_colour))) return false; break; } case event_type::unof: case event_type::qqxexf: { if( !try_connect_t(it_part, line_colour) // u-channel only allowed at impact factor && (std::distance(it_part, cend_partons())!=2 || !try_connect_u(it_part, line_colour))) return false; break; } case event_type::qqxmid:{ auto it_next = std::next(it_part); if( !try_connect_t(it_part, line_colour) // u-channel only allowed at qqx/qxq pair && ( ( !(is_quark(*it_part) && is_antiquark(*it_next)) && !(is_antiquark(*it_part) && is_quark(*it_next))) || !try_connect_u(it_part, line_colour)) ) return false; break; } default: throw std::logic_error{"unreachable"}; } // no colour singlet exchange/disconnected diagram if(line_colour.first == line_colour.second) return false; } return (incoming()[1].colour->first == line_colour.first) && (incoming()[1].colour->second == line_colour.second); } namespace { //! connect incoming Particle to colour flow void connect_incoming(Particle & in, int & colour, int & anti_colour){ in.colour = std::make_pair(anti_colour, colour); // gluon if(in.type == pid::gluon) return; if(in.type > 0){ // quark assert(is_quark(in)); in.colour->second = 0; colour*=-1; return; } // anti-quark assert(is_antiquark(in)); in.colour->first = 0; anti_colour*=-1; } //! connect outgoing Particle to t-channel colour flow template void connect_tchannel( OutIterator & it_part, int & colour, int & anti_colour, RNG & ran ){ assert(colour>0 || anti_colour>0); if(it_part->type == ParticleID::gluon){ // gluon if(colour>0 && anti_colour>0){ // on g line => connect to colour OR anti-colour (random) if(ran.flat() < 0.5){ it_part->colour = std::make_pair(colour+2,colour); colour+=2; } else { it_part->colour = std::make_pair(anti_colour, anti_colour+2); anti_colour+=2; } } else if(colour > 0){ // on q line => connect to available colour it_part->colour = std::make_pair(colour+2, colour); colour+=2; } else { assert(colour<0 && anti_colour>0); // on qx line => connect to available anti-colour it_part->colour = std::make_pair(anti_colour, anti_colour+2); anti_colour+=2; } } else if(is_quark(*it_part)) { // quark assert(anti_colour>0); if(colour>0){ // on g line => connect and remove anti-colour it_part->colour = std::make_pair(anti_colour, 0); anti_colour+=2; anti_colour*=-1; } else { // on qx line => new colour colour*=-1; it_part->colour = std::make_pair(colour, 0); } } else if(is_antiquark(*it_part)) { // anti-quark assert(colour>0); if(anti_colour>0){ // on g line => connect and remove colour it_part->colour = std::make_pair(0, colour); colour+=2; colour*=-1; } else { // on q line => new anti-colour anti_colour*=-1; it_part->colour = std::make_pair(0, anti_colour); } } else { // not a parton assert(!is_parton(*it_part)); it_part->colour = {}; } } //! connect to t- or u-channel colour flow template void connect_utchannel( OutIterator & it_part, int & colour, int & anti_colour, RNG & ran ){ OutIterator it_first = it_part++; if(ran.flat()<.5) {// t-channel connect_tchannel(it_first, colour, anti_colour, ran); connect_tchannel(it_part, colour, anti_colour, ran); } else { // u-channel connect_tchannel(it_part, colour, anti_colour, ran); connect_tchannel(it_first, colour, anti_colour, ran); } } } // namespace bool Event::generate_colours(RNG & ran){ // generate only for HEJ events if(!event_type::is_resummable(type())) return false; assert(std::is_sorted( begin(outgoing()), end(outgoing()), rapidity_less{})); assert(incoming()[0].pz() < incoming()[1].pz()); // positive (anti-)colour -> can connect // negative (anti-)colour -> not available/used up by (anti-)quark int colour = COLOUR_OFFSET; int anti_colour = colour+1; // initialise first connect_incoming(incoming_[0], colour, anti_colour); // reset outgoing colours std::for_each(outgoing_.begin(), outgoing_.end(), [](Particle & part){ part.colour = {};}); for(auto it_part = begin_partons(); it_part!=end_partons(); ++it_part){ switch (type()) { // subleading can connect to t- or u-channel case event_type::unob: case event_type::qqxexb: { if( std::distance(begin_partons(), it_part)==0) connect_utchannel(it_part, colour, anti_colour, ran); else connect_tchannel(it_part, colour, anti_colour, ran); break; } case event_type::unof: case event_type::qqxexf: { if( std::distance(it_part, end_partons())==2) connect_utchannel(it_part, colour, anti_colour, ran); else connect_tchannel(it_part, colour, anti_colour, ran); break; } case event_type::qqxmid:{ auto it_next = std::next(it_part); if( std::distance(begin_partons(), it_part)>0 && std::distance(it_part, end_partons())>2 && ( (is_quark(*it_part) && is_antiquark(*it_next)) || (is_antiquark(*it_part) && is_quark(*it_next)) ) ) connect_utchannel(it_part, colour, anti_colour, ran); else connect_tchannel(it_part, colour, anti_colour, ran); break; } default: // rest has to be t-channel connect_tchannel(it_part, colour, anti_colour, ran); } } // Connect last connect_incoming(incoming_[1], anti_colour, colour); assert(is_leading_colour()); return true; } // generate_colours namespace { bool valid_parton( std::vector const & jets, Particle const & parton, int const idx, double const soft_pt_regulator, double const min_extparton_pt ){ // TODO code overlap with PhaseSpacePoint::pass_extremal_cuts if(min_extparton_pt > parton.pt()) return false; if(idx<0) return false; assert(static_cast(jets.size())>=idx); auto const & jet{ jets[idx] }; return (parton.p - jet).pt()/jet.pt() <= soft_pt_regulator; } } // namespace // this should work with multiple types bool Event::valid_hej_state(double const soft_pt_regulator, double const min_pt ) const { using namespace event_type; if(!is_resummable(type())) return false; auto const & jet_idx{ particle_jet_indices() }; auto idx_begin{ jet_idx.cbegin() }; auto idx_end{ jet_idx.crbegin() }; auto part_begin{ cbegin_partons() }; auto part_end{ crbegin_partons() }; // always seperate extremal jets if( !valid_parton(jets(), *part_begin, *idx_begin, soft_pt_regulator, min_pt) ) return false; ++part_begin; ++idx_begin; if( !valid_parton(jets(), *part_end, *idx_end, soft_pt_regulator, min_pt) ) return false; ++part_end; ++idx_end; // unob -> second parton in own jet if( type() & (unob | qqxexb) ){ if( !valid_parton(jets(), *part_begin, *idx_begin, soft_pt_regulator, min_pt) ) return false; ++part_begin; ++idx_begin; } if( type() & (unof | qqxexf) ){ if( !valid_parton(jets(), *part_end, *idx_end, soft_pt_regulator, min_pt) ) return false; ++part_end; // ++idx_end; // last check, we don't need idx_end afterwards } if( type() & qqxmid ){ // find qqx pair auto begin_qqx{ std::find_if( part_begin, part_end.base(), [](Particle const & part) -> bool { return part.type != ParticleID::gluon; } )}; assert(begin_qqx != part_end.base()); long int qqx_pos{ std::distance(part_begin, begin_qqx) }; assert(qqx_pos >= 0); idx_begin+=qqx_pos; if( !( valid_parton(jets(), *begin_qqx, *idx_begin, soft_pt_regulator, min_pt) && valid_parton(jets(), *std::next(begin_qqx), *std::next(idx_begin), soft_pt_regulator, min_pt) )) return false; } return true; } Event::ConstPartonIterator Event::begin_partons() const { return cbegin_partons(); } Event::ConstPartonIterator Event::cbegin_partons() const { return {HEJ::is_parton, cbegin(outgoing()), cend(outgoing())}; } Event::ConstPartonIterator Event::end_partons() const { return cend_partons(); } Event::ConstPartonIterator Event::cend_partons() const { return {HEJ::is_parton, cend(outgoing()), cend(outgoing())}; } Event::ConstReversePartonIterator Event::rbegin_partons() const { return crbegin_partons(); } Event::ConstReversePartonIterator Event::crbegin_partons() const { return std::reverse_iterator( cend_partons() ); } Event::ConstReversePartonIterator Event::rend_partons() const { return crend_partons(); } Event::ConstReversePartonIterator Event::crend_partons() const { return std::reverse_iterator( cbegin_partons() ); } Event::PartonIterator Event::begin_partons() { return {HEJ::is_parton, begin(outgoing_), end(outgoing_)}; } Event::PartonIterator Event::end_partons() { return {HEJ::is_parton, end(outgoing_), end(outgoing_)}; } Event::ReversePartonIterator Event::rbegin_partons() { return std::reverse_iterator( end_partons() ); } Event::ReversePartonIterator Event::rend_partons() { return std::reverse_iterator( begin_partons() ); } namespace { void print_momentum(std::ostream & os, fastjet::PseudoJet const & part){ constexpr int prec = 6; const std::streamsize orig_prec = os.precision(); os < const & col){ constexpr int width = 3; if(!col) os << "(no color)"; // American spelling for better alignment else os << "(" <first << ", " <second << ")"; } } // namespace std::ostream& operator<<(std::ostream & os, Event const & ev){ constexpr int prec = 4; constexpr int wtype = 3; // width for types const std::streamsize orig_prec = os.precision(); os < rapidity=" < rapidity=" < rapidity=" < incoming{ event.incoming() }; // First incoming should be positive pz according to LHE standard // (or at least most (everyone?) do it this way, and Pythia assumes it) if(incoming[0].pz() < incoming[1].pz()) std::swap(incoming[0], incoming[1]); for(Particle const & in: incoming){ result.IDUP.emplace_back(in.type); result.ISTUP.emplace_back(LHE_Status::in); result.PUP.push_back({in.p[0], in.p[1], in.p[2], in.p[3], in.p.m()}); result.MOTHUP.emplace_back(0, 0); assert(in.colour); result.ICOLUP.emplace_back(*in.colour); } // outgoing for(size_t i = 0; i < event.outgoing().size(); ++i){ Particle const & out = event.outgoing()[i]; result.IDUP.emplace_back(out.type); const int status = event.decays().count(i) != 0u ?LHE_Status::decay :LHE_Status::out; result.ISTUP.emplace_back(status); result.PUP.push_back({out.p[0], out.p[1], out.p[2], out.p[3], out.p.m()}); result.MOTHUP.emplace_back(1, 2); if(out.colour) result.ICOLUP.emplace_back(*out.colour); else{ result.ICOLUP.emplace_back(std::make_pair(0,0)); } } // decays for(auto const & decay: event.decays()){ for(auto const & out: decay.second){ result.IDUP.emplace_back(out.type); result.ISTUP.emplace_back(LHE_Status::out); result.PUP.push_back({out.p[0], out.p[1], out.p[2], out.p[3], out.p.m()}); const size_t mother_idx = 1 + event.incoming().size() + decay.first; result.MOTHUP.emplace_back(mother_idx, mother_idx); result.ICOLUP.emplace_back(0,0); } } assert(result.ICOLUP.size() == num_particles); static constexpr double unknown_spin = 9.; //per Les Houches accord result.VTIMUP = std::vector(num_particles, unknown_spin); result.SPINUP = result.VTIMUP; return result; } } // namespace HEJ diff --git a/src/YAMLreader.cc b/src/YAMLreader.cc index 2896389..cd1d08c 100644 --- a/src/YAMLreader.cc +++ b/src/YAMLreader.cc @@ -1,574 +1,574 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #include "HEJ/YAMLreader.hh" #include #include #include #include #include #include #include #include #include "HEJ/ConfigFlags.hh" #include "HEJ/Constants.hh" #include "HEJ/ScaleFunction.hh" #include "HEJ/event_types.hh" #include "HEJ/output_formats.hh" namespace HEJ { class Event; namespace { //! Get YAML tree of supported options /** * The configuration file is checked against this tree of options * in assert_all_options_known. */ YAML::Node const & get_supported_options(){ const static YAML::Node supported = [](){ YAML::Node supported; static const auto opts = { "trials", "min extparton pt", "max ext soft pt fraction", "soft pt regulator", "scales", "scale factors", "max scale ratio", "import scales", "log correction", "event output", "analysis", "analyses", "vev", "regulator parameter", "max events" }; // add subnodes to "supported" - the assigned value is irrelevant for(auto && opt: opts) supported[opt] = ""; for(auto && jet_opt: {"min pt", "algorithm", "R"}){ supported["resummation jets"][jet_opt] = ""; supported["fixed order jets"][jet_opt] = ""; } for(auto && opt: {"mt", "use impact factors", "include bottom", "mb"}){ supported["Higgs coupling"][opt] = ""; } for(auto && opt: {"name", "seed"}){ supported["random generator"][opt] = ""; } for(auto && opt: {"FKL", "unordered", "extremal qqx", "central qqx", "non-resummable"}){ supported["event treatment"][opt] = ""; } for(auto && particle_type: {"Higgs", "W", "Z"}){ for(auto && particle_opt: {"mass", "width"}){ supported["particle properties"][particle_type][particle_opt] = ""; } } for(auto && opt: {"type", "trials", "max deviation"}){ supported["unweight"][opt] = ""; } return supported; }(); return supported; } fastjet::JetAlgorithm to_JetAlgorithm(std::string const & algo){ using namespace fastjet; static const std::map known = { {"kt", kt_algorithm}, {"cambridge", cambridge_algorithm}, {"antikt", antikt_algorithm}, {"cambridge for passive", cambridge_for_passive_algorithm}, {"plugin", plugin_algorithm} }; const auto res = known.find(algo); if(res == known.end()){ throw std::invalid_argument("Unknown jet algorithm \"" + algo + "\""); } return res->second; } EventTreatment to_EventTreatment(std::string const & name){ static const std::map known = { {"reweight", EventTreatment::reweight}, {"keep", EventTreatment::keep}, {"discard", EventTreatment::discard} }; const auto res = known.find(name); if(res == known.end()){ throw std::invalid_argument("Unknown event treatment \"" + name + "\""); } return res->second; } WeightType to_weight_type(std::string const & setting){ if(setting == "weighted") return WeightType::weighted; if(setting =="resummation") return WeightType::unweighted_resum; if(setting =="partial") return WeightType::partially_unweighted; throw std::invalid_argument{"Unknown weight type \"" + setting + "\""}; } } // namespace namespace detail{ void set_from_yaml(fastjet::JetAlgorithm & setting, YAML::Node const & yaml){ setting = to_JetAlgorithm(yaml.as()); } void set_from_yaml(EventTreatment & setting, YAML::Node const & yaml){ setting = to_EventTreatment(yaml.as()); } void set_from_yaml(ParticleID & setting, YAML::Node const & yaml){ setting = to_ParticleID(yaml.as()); } void set_from_yaml(WeightType & setting, YAML::Node const & yaml){ setting = to_weight_type(yaml.as()); } } // namespace detail JetParameters get_jet_parameters( YAML::Node const & node, std::string const & entry ){ assert(node); JetParameters result; fastjet::JetAlgorithm jet_algo = fastjet::antikt_algorithm; double R = NAN; set_from_yaml_if_defined(jet_algo, node, entry, "algorithm"); set_from_yaml(R, node, entry, "R"); result.def = fastjet::JetDefinition{jet_algo, R}; set_from_yaml(result.min_pt, node, entry, "min pt"); return result; } RNGConfig to_RNGConfig( YAML::Node const & node, std::string const & entry ){ assert(node); RNGConfig result; set_from_yaml(result.name, node, entry, "name"); set_from_yaml_if_defined(result.seed, node, entry, "seed"); return result; } ParticleProperties get_particle_properties( YAML::Node const & node, std::string const & entry, std::string const & boson ){ ParticleProperties result{}; set_from_yaml(result.mass, node, entry, boson, "mass"); set_from_yaml(result.width, node, entry, boson, "width"); return result; } EWConstants get_ew_parameters(YAML::Node const & node){ EWConstants result; double vev = NAN; set_from_yaml(vev, node, "vev"); result.set_vevWZH(vev, get_particle_properties(node, "particle properties", "W"), get_particle_properties(node, "particle properties", "Z"), get_particle_properties(node, "particle properties", "Higgs") ); return result; } HiggsCouplingSettings get_Higgs_coupling( YAML::Node const & node, std::string const & entry ){ assert(node); static constexpr double mt_max = 2e4; #ifndef HEJ_BUILD_WITH_QCDLOOP if(node[entry].IsDefined()){ throw std::invalid_argument{ "Higgs coupling settings require building HEJ 2 " "with QCDloop support" }; } #endif HiggsCouplingSettings settings; set_from_yaml_if_defined(settings.mt, node, entry, "mt"); set_from_yaml_if_defined(settings.mb, node, entry, "mb"); set_from_yaml_if_defined(settings.include_bottom, node, entry, "include bottom"); set_from_yaml_if_defined(settings.use_impact_factors, node, entry, "use impact factors"); if(settings.use_impact_factors){ if(settings.mt != std::numeric_limits::infinity()){ throw std::invalid_argument{ "Conflicting settings: " "impact factors may only be used in the infinite top mass limit" }; } } else{ // huge values of the top mass are numerically unstable settings.mt = std::min(settings.mt, mt_max); } return settings; } FileFormat to_FileFormat(std::string const & name){ static const std::map known = { {"Les Houches", FileFormat::Les_Houches}, {"HepMC", FileFormat::HepMC}, {"HepMC2", FileFormat::HepMC2}, {"HepMC3", FileFormat::HepMC3}, {"HDF5", FileFormat::HDF5} }; const auto res = known.find(name); if(res == known.end()){ throw std::invalid_argument("Unknown file format \"" + name + "\""); } return res->second; } std::string extract_suffix(std::string const & filename){ size_t separator = filename.rfind('.'); if(separator == std::string::npos) return {}; return filename.substr(separator + 1); } FileFormat format_from_suffix(std::string const & filename){ const std::string suffix = extract_suffix(filename); if(suffix == "lhe") return FileFormat::Les_Houches; if(suffix == "hepmc") return FileFormat::HepMC; if(suffix == "hepmc3") return FileFormat::HepMC3; if(suffix == "hepmc2") return FileFormat::HepMC2; if(suffix == "hdf5") return FileFormat::HDF5; throw std::invalid_argument{ "Can't determine format for output file \"" + filename + "\"" }; } void assert_all_options_known( YAML::Node const & conf, YAML::Node const & supported ){ if(!conf.IsMap()) return; if(!supported.IsMap()) throw invalid_type{"must not have sub-entries"}; for(auto const & entry: conf){ const auto name = entry.first.as(); if(! supported[name]) throw unknown_option{name}; /* check sub-options, e.g. 'resummation jets: min pt' * we don't check analyses sub-options * those depend on the analysis being used and should be checked there * similar for "import scales" */ if(name != "analyses" && name != "analysis" && name != "import scales"){ try{ assert_all_options_known(conf[name], supported[name]); } catch(unknown_option const & ex){ throw unknown_option{name + ": " + ex.what()}; } catch(invalid_type const & ex){ throw invalid_type{name + ": " + ex.what()}; } } } } } // namespace HEJ namespace YAML { Node convert::encode(HEJ::OutputFile const & outfile) { Node node; node[to_string(outfile.format)] = outfile.name; return node; } bool convert::decode(Node const & node, HEJ::OutputFile & out) { switch(node.Type()){ case NodeType::Map: { YAML::const_iterator it = node.begin(); out.format = HEJ::to_FileFormat(it->first.as()); out.name = it->second.as(); return true; } case NodeType::Scalar: out.name = node.as(); out.format = HEJ::format_from_suffix(out.name); return true; default: return false; } } } // namespace YAML namespace HEJ { namespace detail{ void set_from_yaml(OutputFile & setting, YAML::Node const & yaml){ setting = yaml.as(); } } namespace { void update_fixed_order_jet_parameters( JetParameters & fixed_order_jets, YAML::Node const & yaml ){ if(!yaml["fixed order jets"]) return; set_from_yaml_if_defined( fixed_order_jets.min_pt, yaml, "fixed order jets", "min pt" ); fastjet::JetAlgorithm algo = fixed_order_jets.def.jet_algorithm(); set_from_yaml_if_defined(algo, yaml, "fixed order jets", "algorithm"); double R = fixed_order_jets.def.R(); set_from_yaml_if_defined(R, yaml, "fixed order jets", "R"); fixed_order_jets.def = fastjet::JetDefinition{algo, R}; } // like std::stod, but throw if not the whole string can be converted double to_double(std::string const & str){ std::size_t pos = 0; const double result = std::stod(str, &pos); if(pos < str.size()){ throw std::invalid_argument(str + " is not a valid double value"); } return result; } using EventScale = double (*)(Event const &); void import_scale_functions( std::string const & file, std::vector const & scale_names, std::unordered_map & known ) { void * handle = dlopen(file.c_str(), RTLD_NOW); char * error = dlerror(); if(error != nullptr) throw std::runtime_error{error}; for(auto const & scale: scale_names) { void * sym = dlsym(handle, scale.c_str()); error = dlerror(); if(error != nullptr) throw std::runtime_error{error}; known.emplace(scale, reinterpret_cast(sym)); // NOLINT } } auto get_scale_map( YAML::Node const & yaml ) { std::unordered_map scale_map; scale_map.emplace("H_T", H_T); scale_map.emplace("max jet pperp", max_jet_pt); scale_map.emplace("jet invariant mass", jet_invariant_mass); scale_map.emplace("m_j1j2", m_j1j2); if(yaml["import scales"].IsDefined()) { if(! yaml["import scales"].IsMap()) { throw invalid_type{"Entry 'import scales' is not a map"}; } for(auto const & import: yaml["import scales"]) { const auto file = import.first.as(); const auto scale_names = import.second.IsSequence() ?import.second.as>() :std::vector{import.second.as()}; import_scale_functions(file, scale_names, scale_map); } } return scale_map; } // simple (as in non-composite) scale functions /** * An example for a simple scale function would be H_T, * H_T/2 is then composite (take H_T and then divide by 2) */ ScaleFunction parse_simple_ScaleFunction( std::string const & scale_fun, std::unordered_map const & known ) { assert( scale_fun.empty() || (!std::isspace(scale_fun.front()) && !std::isspace(scale_fun.back())) ); const auto it = known.find(scale_fun); if(it != end(known)) return {it->first, it->second}; try{ const double scale = to_double(scale_fun); return {scale_fun, FixedScale{scale}}; } catch(std::invalid_argument const &){} throw std::invalid_argument{"Unknown scale choice: \"" + scale_fun + "\""}; } std::string trim_front(std::string const & str){ const auto new_begin = std::find_if( begin(str), end(str), [](char c){ return std::isspace(c) == 0; } ); return std::string(new_begin, end(str)); } std::string trim_back(std::string str){ size_t pos = str.size() - 1; // use guaranteed wrap-around behaviour to check whether we have // traversed the whole string for(; pos < str.size() && std::isspace(str[pos]); --pos) {} str.resize(pos + 1); // note that pos + 1 can be 0 return str; } ScaleFunction parse_ScaleFunction( std::string const & scale_fun, std::unordered_map const & known ){ assert( scale_fun.empty() || (!std::isspace(scale_fun.front()) && !std::isspace(scale_fun.back())) ); // parse from right to left => a/b/c gives (a/b)/c const size_t delim = scale_fun.find_last_of("*/"); if(delim == std::string::npos){ return parse_simple_ScaleFunction(scale_fun, known); } const std::string first = trim_back(std::string{scale_fun, 0, delim}); const std::string second = trim_front(std::string{scale_fun, delim+1}); if(scale_fun[delim] == '/'){ return parse_ScaleFunction(first, known) / parse_ScaleFunction(second, known); } assert(scale_fun[delim] == '*'); return parse_ScaleFunction(first, known) * parse_ScaleFunction(second, known); } EventTreatMap get_event_treatment( YAML::Node const & node, std::string const & entry ){ using namespace event_type; EventTreatMap treat { - {no_2_jets, EventTreatment::discard}, + {not_enough_jets, EventTreatment::discard}, {bad_final_state, EventTreatment::discard}, {FKL, EventTreatment::discard}, {unob, EventTreatment::discard}, {unof, EventTreatment::discard}, {qqxexb, EventTreatment::discard}, {qqxexf, EventTreatment::discard}, {qqxmid, EventTreatment::discard}, {non_resummable, EventTreatment::discard} }; set_from_yaml(treat.at(FKL), node, entry, "FKL"); set_from_yaml(treat.at(unob), node, entry, "unordered"); treat.at(unof) = treat.at(unob); set_from_yaml(treat.at(qqxexb), node, entry, "extremal qqx"); treat.at(qqxexf) = treat.at(qqxexb); set_from_yaml(treat.at(qqxmid), node, entry, "central qqx"); set_from_yaml(treat.at(non_resummable), node, entry, "non-resummable"); if(treat[non_resummable] == EventTreatment::reweight){ throw std::invalid_argument{"Cannot reweight non-resummable events"}; } return treat; } Config to_Config(YAML::Node const & yaml){ try{ assert_all_options_known(yaml, get_supported_options()); } catch(unknown_option const & ex){ throw unknown_option{std::string{"Unknown option '"} + ex.what() + "'"}; } Config config; config.resummation_jets = get_jet_parameters(yaml, "resummation jets"); config.fixed_order_jets = config.resummation_jets; update_fixed_order_jet_parameters(config.fixed_order_jets, yaml); set_from_yaml_if_defined(config.min_extparton_pt, yaml, "min extparton pt"); if(config.min_extparton_pt!=0) std::cerr << "WARNING: \"min extparton pt\" is deprecated." << " Please remove this entry or set \"soft pt regulator\" instead.\n"; set_from_yaml_if_defined( config.max_ext_soft_pt_fraction, yaml, "max ext soft pt fraction" ); if(config.max_ext_soft_pt_fraction){ std::cerr << "WARNING: \"max ext soft pt fraction\" is deprecated." << " Please remove this entry or set \"soft pt regulator\" instead.\n"; config.soft_pt_regulator = *config.max_ext_soft_pt_fraction; } else { set_from_yaml_if_defined( config.soft_pt_regulator, yaml, "soft pt regulator" ); } // Sets the standard value, then changes this if defined config.regulator_lambda=CLAMBDA; set_from_yaml_if_defined(config.regulator_lambda, yaml, "regulator parameter"); set_from_yaml_if_defined(config.max_events, yaml, "max events"); set_from_yaml(config.trials, yaml, "trials"); config.weight_type = WeightType::weighted; set_from_yaml_if_defined(config.weight_type, yaml, "unweight", "type"); if(config.weight_type == WeightType::partially_unweighted) { config.unweight_config = PartialUnweightConfig{}; set_from_yaml( config.unweight_config->trials, yaml, "unweight", "trials" ); set_from_yaml( config.unweight_config->max_dev, yaml, "unweight", "max deviation" ); } else if(yaml["unweight"].IsDefined()) { for(auto && opt: {"trials", "max deviation"}) { if(yaml["unweight"][opt].IsDefined()) { throw std::invalid_argument{ "'unweight: " + std::string{opt} + "' " "is only supported if 'unweight: type' is set to 'partial'" }; } } } set_from_yaml(config.log_correction, yaml, "log correction"); config.treat = get_event_treatment(yaml, "event treatment"); set_from_yaml_if_defined(config.output, yaml, "event output"); config.rng = to_RNGConfig(yaml, "random generator"); set_from_yaml_if_defined(config.analyses_parameters, yaml, "analyses"); if(yaml["analysis"].IsDefined()){ std::cerr << "WARNING: Configuration entry 'analysis' is deprecated. " " Use 'analyses' instead.\n"; set_from_yaml(config.analysis_parameters, yaml, "analysis"); if(!config.analysis_parameters.IsNull()){ config.analyses_parameters.push_back(config.analysis_parameters); } } config.scales = to_ScaleConfig(yaml); config.ew_parameters = get_ew_parameters(yaml); config.Higgs_coupling = get_Higgs_coupling(yaml, "Higgs coupling"); return config; } } // namespace ScaleConfig to_ScaleConfig(YAML::Node const & yaml){ ScaleConfig config; auto scale_funs = get_scale_map(yaml); std::vector scales; set_from_yaml(scales, yaml, "scales"); config.base.reserve(scales.size()); std::transform( begin(scales), end(scales), std::back_inserter(config.base), [scale_funs](auto const & entry){ return parse_ScaleFunction(entry, scale_funs); } ); set_from_yaml_if_defined(config.factors, yaml, "scale factors"); config.max_ratio = std::numeric_limits::infinity(); set_from_yaml_if_defined(config.max_ratio, yaml, "max scale ratio"); return config; } Config load_config(std::string const & config_file){ try{ return to_Config(YAML::LoadFile(config_file)); } catch(...){ std::cerr << "Error reading " << config_file << ":\n "; throw; } } } // namespace HEJ diff --git a/t/check_res.cc b/t/check_res.cc index 33b0103..bfa8d00 100644 --- a/t/check_res.cc +++ b/t/check_res.cc @@ -1,164 +1,164 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #include "hej_test.hh" #include #include #include #include #include #include #include #include #include "HEJ/Config.hh" #include "HEJ/CrossSectionAccumulator.hh" #include "HEJ/Event.hh" #include "HEJ/event_types.hh" #include "HEJ/EventReweighter.hh" #include "HEJ/EWConstants.hh" #include "HEJ/Fraction.hh" #include "HEJ/HiggsCouplingSettings.hh" #include "HEJ/Mixmax.hh" #include "HEJ/Parameters.hh" #include "HEJ/ScaleFunction.hh" #include "HEJ/stream.hh" #include "fastjet/JetDefinition.hh" #include "LHEF/LHEF.h" namespace HEJ { struct RNG; } 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 JETPTMIN = 35; constexpr bool LOG_CORR = false; constexpr std::size_t NUM_TRIES = 100; constexpr HEJ::ParticleProperties WPROP{80.385, 2.085}; constexpr HEJ::ParticleProperties ZPROP{91.187, 2.495}; constexpr HEJ::ParticleProperties HPROP{125, 0.004165}; constexpr double VEV = 246.2196508; using EventTreatment = HEJ::EventTreatment; using namespace HEJ::event_type; HEJ::EventTreatMap TREAT{ - {no_2_jets, EventTreatment::discard}, + {not_enough_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} }; bool correct_colour(HEJ::Event const & ev){ if(!HEJ::event_type::is_resummable(ev.type())) return true; return ev.is_leading_colour(); } } // namespace 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}; HEJ::MatrixElementConfig ME_conf; ME_conf.log_correction = LOG_CORR; ME_conf.Higgs_coupling = HEJ::HiggsCouplingSettings{}; ME_conf.ew_parameters.set_vevWZH(VEV, WPROP, ZPROP, HPROP); HEJ::EventReweighterConfig conf; conf.psp_config = std::move(psp_conf); conf.ME_config = std::move(ME_conf); 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. }; std::shared_ptr ran{std::make_shared()}; HEJ::EventReweighter hej{reader.heprup, std::move(scale_gen), conf, ran}; HEJ::CrossSectionAccumulator xs; do{ auto ev_data = HEJ::Event::EventData{reader.hepeup}; shuffle_particles(ev_data); ev_data.reconstruct_intermediate(); HEJ::Event ev{ ev_data.cluster( BORN_JET_DEF, BORN_JETPTMIN ) }; auto resummed_events = hej.reweight(ev, NUM_TRIES); for(auto const & res_ev: resummed_events) { ASSERT(correct_colour(res_ev)); ASSERT(std::isfinite(res_ev.central().weight)); // we fill the xs uncorrelated since we only want to test the uncertainty // of the resummation xs.fill(res_ev); } } while(reader.readEvent()); const double xsec = xs.total().value; const double xsec_err = std::sqrt(xs.total().error); 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/test_classify.cc b/t/test_classify.cc index 14aad7b..5c47fc8 100644 --- a/t/test_classify.cc +++ b/t/test_classify.cc @@ -1,525 +1,525 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #include "hej_test.hh" #include #include #include #include #include #include #include "fastjet/JetDefinition.hh" #include "HEJ/Event.hh" #include "HEJ/event_types.hh" #include "HEJ/exceptions.hh" #include "HEJ/PDG_codes.hh" namespace { const fastjet::JetDefinition JET_DEF{fastjet::JetAlgorithm::antikt_algorithm, 0.4}; const double MIN_JET_PT{30.}; const std::vector ALL_QUARKS{"-4","-1","1","2","3","4"}; const std::vector ALL_PARTONS{"g","-2","-1","1","2","3","4"}; const std::vector ALL_BOSONS{"h", "Wp", "Wm", "Z_photon_mix"}; const std::vector ALL_G_Z{"photon", "Z"}; const std::vector ALL_W{"W+", "W-"}; const std::size_t MAX_MULTI = 6; std::mt19937_64 RAN{0}; bool couple_quark(std::string const & boson, std::string & quark){ if(std::abs(HEJ::to_ParticleID(boson)) == HEJ::ParticleID::Wp){ auto qflav{ HEJ::to_ParticleID(quark) }; if(!HEJ::is_anyquark(qflav)) return false; const int W_charge = HEJ::to_ParticleID(boson)>0?1:-1; if(W_charge*qflav < 0 && !(std::abs(qflav)%2)) return false; // not anti-down if(W_charge*qflav > 0 && (std::abs(qflav)%2)) return false; // not up quark=std::to_string(qflav-W_charge); } if(HEJ::to_ParticleID(boson) == HEJ::ParticleID::Z_photon_mix){ auto qflav{ HEJ::to_ParticleID(quark) }; if(!HEJ::is_anyquark(qflav)) return false; } return true; } bool match_expectation( HEJ::event_type::EventType expected, std::array const & in, std::vector const & out, int const overwrite_boson = 0 ){ HEJ::Event ev{ parse_configuration( in,out,overwrite_boson ).cluster(JET_DEF, MIN_JET_PT)}; if(ev.type() != expected){ std::cerr << "Expected type " << name(expected) << " but found " << name(ev.type()) << "\n" << ev; auto jet_idx{ ev.particle_jet_indices() }; std::cout << "Particle Jet indices: "; for(int const i: jet_idx) std::cout << i << " "; std::cout << std::endl; return false; } return true; } //! test FKL configurations //! if implemented==false : check processes that are not in HEJ yet bool check_fkl( bool const implemented=true ){ using namespace HEJ; auto const type{ implemented?event_type::FKL:event_type::non_resummable }; std::vector bosons; if(implemented) bosons = ALL_BOSONS; else { bosons = ALL_G_Z; } for(std::string const & first: ALL_PARTONS) // all quark flavours for(std::string const & last: ALL_PARTONS){ for(std::size_t njet=2; njet<=MAX_MULTI; ++njet){ // all multiplicities if(njet==5) continue; std::array base_in{first,last}; std::vector base_out(njet, "g"); base_out.front() = first; base_out.back() = last; if(implemented && !match_expectation(type, base_in, base_out)) return false; for(auto const & boson: bosons) // any boson for(std::size_t pos=0; pos<=njet; ++pos){ // at any position auto const & in{base_in}; auto out{base_out}; // change quark flavours for W const bool couple_idx = std::uniform_int_distribution{0,1}(RAN) != 0; if(!couple_quark(boson, couple_idx?out.back():out.front())) continue; out.insert(out.begin()+pos, boson); if(!match_expectation(type, in, out)) return false; } } } return true; } //! test unordered configurations //! if implemented==false : check processes that are not in HEJ yet bool check_uno( bool const implemented=true ){ using namespace HEJ; auto const b{ implemented?event_type::unob:event_type::non_resummable }; auto const f{ implemented?event_type::unof:event_type::non_resummable }; std::vector bosons; if(implemented) { bosons = ALL_BOSONS; } else { bosons = ALL_G_Z; } for(std::string const & uno: ALL_QUARKS) // all quark flavours for(std::string const & fkl: ALL_PARTONS){ for(std::size_t njet=3; njet<=MAX_MULTI; ++njet){ // all multiplicities >2 if(njet==5) continue; for(std::size_t i=0; i<2; ++i){ // forward & backwards std::array base_in; std::vector base_out(njet, "g"); const std::size_t uno_pos = i?1:(njet-2); const std::size_t fkl_pos = i?(njet-1):0; base_in[i?0:1] = uno; base_in[i?1:0] = fkl; base_out[uno_pos] = uno; base_out[fkl_pos] = fkl; auto expectation{ i?b:f }; if( implemented && !match_expectation(expectation, base_in, base_out) ) return false; for(auto const & boson: bosons){ // any boson // at any position (higgs only inside FKL chain) std::size_t start = 0; std::size_t end = njet; if(to_ParticleID(boson) == pid::higgs){ start = i?(uno_pos+1):fkl_pos; end = i?(fkl_pos+1):uno_pos; } for(std::size_t pos=start; pos<=end; ++pos){ auto const & in{base_in}; auto out{base_out}; // change quark flavours for W const bool couple_idx = std::uniform_int_distribution{0,1}(RAN) != 0; if(!couple_quark(boson, couple_idx?out[fkl_pos]:out[uno_pos])) continue; out.insert(out.begin()+pos, boson); if(!match_expectation(expectation, in, out)) return false; } } } } } return true; } //! test extremal qqx configurations //! if implemented==false : check processes that are not in HEJ yet bool check_extremal_qqx( bool const implemented=true ){ using namespace HEJ; auto const b{ implemented?event_type::qqxexb:event_type::non_resummable }; auto const f{ implemented?event_type::qqxexf:event_type::non_resummable }; std::vector bosons; if(implemented) bosons = ALL_W; else { bosons = ALL_G_Z; bosons.emplace_back("h"); bosons.emplace_back("Z_photon_mix"); } for(std::string const & qqx: ALL_QUARKS) // all quark flavours for(std::string const & fkl: ALL_PARTONS){ std::string const qqx2{ std::to_string(HEJ::to_ParticleID(qqx)*-1) }; for(std::size_t njet=3; njet<=MAX_MULTI; ++njet){ // all multiplicities >2 if(njet==5) continue; for(std::size_t i=0; i<2; ++i){ // forward & backwards std::array base_in; std::vector base_out(njet, "g"); const std::size_t qqx_pos = i?0:(njet-2); const std::size_t fkl_pos = i?(njet-1):0; base_in[i?0:1] = "g"; base_in[i?1:0] = fkl; base_out[fkl_pos] = fkl; base_out[qqx_pos] = qqx; base_out[qqx_pos+1] = qqx2; auto expectation{ i?b:f }; if( implemented && !match_expectation(expectation, base_in, base_out) ) return false; for(auto const & boson: bosons){ // all bosons // at any position (higgs only inside FKL chain) std::size_t start = 0; std::size_t end = njet; if(to_ParticleID(boson) == pid::higgs){ start = i?(qqx_pos+2):fkl_pos; end = i?(fkl_pos+1):qqx_pos; } for(std::size_t pos=start; pos<=end; ++pos){ auto const & in{base_in}; auto out{base_out}; // change quark flavours for W const bool couple_idx = std::uniform_int_distribution{0,1}(RAN) != 0; if(couple_idx || !couple_quark(boson, out[fkl_pos]) ){ // (randomly) try couple to FKL, else fall-back to qqx if(!couple_quark(boson, out[qqx_pos])) couple_quark(boson, out[qqx_pos+1]); } out.insert(out.begin()+pos, boson); if(!match_expectation(expectation, in, out)) return false; } } } } // test allowed jet configurations if( implemented){ if( !( match_expectation(f,{fkl,"g"},{fkl,"g","g","g","g",qqx,qqx2}, -3) && match_expectation(b,{"g",fkl},{qqx,qqx2,"g","g","g","g",fkl}, -4) && match_expectation(f,{fkl,"g"},{fkl,"g","g","g","g",qqx,qqx2}, -5) && match_expectation(b,{"g",fkl},{qqx,qqx2,"g","g","g","g",fkl}, -5) && match_expectation(f,{fkl,"g"},{fkl,"g","g","g","g",qqx,qqx2}, -6) && match_expectation(f,{fkl,"g"},{fkl,"g","g","g","g",qqx,qqx2}, -7) && match_expectation(b,{"g",fkl},{qqx,qqx2,"g","g","g","g",fkl}, -7) && match_expectation(f,{fkl,"g"},{fkl,"g","g","g","g",qqx,qqx2}, -8) && match_expectation(b,{"g",fkl},{qqx,qqx2,"g","g","g","g",fkl}, -8) && match_expectation(b,{"g",fkl},{qqx,qqx2,"g","g","g","g",fkl}, -9) && match_expectation(f,{fkl,"g"},{fkl,"g","g","g","g",qqx,qqx2}, -10) && match_expectation(f,{fkl,"g"},{fkl,"g","g","g","g",qqx,qqx2}, -11) && match_expectation(b,{"g",fkl},{qqx,qqx2,"g","g","g","g",fkl}, -11) && match_expectation(f,{fkl,"g"},{fkl,"g","g","g","g",qqx,qqx2}, -12) && match_expectation(b,{"g",fkl},{qqx,qqx2,"g","g","g","g",fkl}, -12) )) return false; if (fkl == "2") { if( !( match_expectation(f,{"2","g"},{"1","Wp","g","g","g",qqx,qqx2}, -3) && match_expectation(b,{"g","2"},{qqx,qqx2,"g","Wp","g","g","1"}, -4) && match_expectation(f,{"2","g"},{"1","Wp","g","g","g",qqx,qqx2}, -5) && match_expectation(b,{"g","2"},{qqx,qqx2,"g","Wp","g","g","1"}, -5) && match_expectation(f,{"2","g"},{"1","g","Wp","g","g",qqx,qqx2}, -6) && match_expectation(f,{"2","g"},{"1","g","g","g","Wp",qqx,qqx2}, -7) && match_expectation(b,{"g","2"},{qqx,qqx2,"g","g","g","Wp","1"}, -7) && match_expectation(f,{"2","g"},{"1","Wp","g","g","g",qqx,qqx2}, -8) && match_expectation(b,{"g","2"},{qqx,qqx2,"Wp","g","g","g","1"}, -8) && match_expectation(b,{"g","2"},{qqx,qqx2,"g","Wp","g","g","1"}, -9) && match_expectation(f,{"2","g"},{"1","g","g","g","Wp",qqx,qqx2}, -10) && match_expectation(f,{"2","g"},{"1","g","g","g","Wp",qqx,qqx2}, -11) && match_expectation(b,{"g","2"},{qqx,qqx2,"g","g","g","Wp","1"}, -11) && match_expectation(f,{"2","g"},{"1","g","g","g","Wp",qqx,qqx2}, -12) && match_expectation(b,{"g","2"},{qqx,qqx2,"g","Wp","g","g","1"}, -12) )) return false; } } } return true; } //! test central qqx configurations //! if implemented==false : check processes that are not in HEJ yet bool check_central_qqx(bool const implemented=true){ using namespace HEJ; auto const t{ implemented?event_type::qqxmid:event_type::non_resummable }; std::vector bosons; if(implemented) bosons = ALL_W; else { bosons = ALL_G_Z; bosons.emplace_back("h"); bosons.emplace_back("Z_photon_mix"); } for(std::string const & qqx: ALL_QUARKS) // all quark flavours for(std::string const & fkl1: ALL_PARTONS) for(std::string const & fkl2: ALL_PARTONS){ std::string const qqx2{ std::to_string(HEJ::to_ParticleID(qqx)*-1) }; for(std::size_t njet=4; njet<=MAX_MULTI; ++njet){ // all multiplicities >3 if(njet==5) continue; for(std::size_t qqx_pos=1; qqx_pos base_in; std::vector base_out(njet, "g"); base_in[0] = fkl1; base_in[1] = fkl2; base_out.front() = fkl1; base_out.back() = fkl2; base_out[qqx_pos] = qqx; base_out[qqx_pos+1] = qqx2; if( implemented && !match_expectation(t, base_in, base_out) ) return false; for(auto const & boson: bosons) // any boson for(std::size_t pos=0; pos<=njet; ++pos){ // at any position if( to_ParticleID(boson) == pid::higgs && (pos==qqx_pos || pos==qqx_pos+1) ) continue; auto const & in{base_in}; auto out{base_out}; // change quark flavours for W const int couple_idx{ std::uniform_int_distribution{0,2}(RAN) }; // (randomly) try couple to FKL, else fall-back to qqx if( couple_idx == 0 && couple_quark(boson, out.front()) ){} else if( couple_idx == 1 && couple_quark(boson, out.back()) ){} else { if(!couple_quark(boson, out[qqx_pos])) couple_quark(boson, out[qqx_pos+1]); } out.insert(out.begin()+pos, boson); if(!match_expectation(t, in, out)) return false; } } } } return true; } // this checks a (non excessive) list of non-resummable states bool check_non_resummable(){ auto type{ HEJ::event_type::non_resummable}; return // 2j - crossing lines match_expectation(type, {"g","2"}, {"2","g"}) && match_expectation(type, {"-1","g"}, {"g","-1"}) && match_expectation(type, {"1","-1"}, {"-1","1"}) && match_expectation(type, {"g","2"}, {"2","g","h"}) && match_expectation(type, {"1","2"}, {"2","h","1"}) && match_expectation(type, {"1","-1"}, {"h","-1","1"}) && match_expectation(type, {"g","2"}, {"Wp","1","g"}) && match_expectation(type, {"1","-1"}, {"-2","Wp","1"}) && match_expectation(type, {"4","g"}, {"g","3","Wp"}) && match_expectation(type, {"1","-2"}, {"-1","Wm","1"}) && match_expectation(type, {"g","3"}, {"4","g","Wm"}) && match_expectation(type, {"1","3"}, {"Wm","4","1"}) && match_expectation(type, {"g","2"}, {"Z_photon_mix","2","g"}) && match_expectation(type, {"1","-1"}, {"-1","Z_photon_mix","1"}) && match_expectation(type, {"4","g"}, {"g","4","Z_photon_mix"}) // 2j - qqx && match_expectation(type, {"g","g"}, {"1","-1"}) && match_expectation(type, {"g","g"}, {"-2","2","h"}) && match_expectation(type, {"g","g"}, {"-4","Wp","3"}) && match_expectation(type, {"g","g"}, {"Wm","-1","2"}) && match_expectation(type, {"g","g"}, {"-3","Z_photon_mix","3"}) // 3j - crossing lines && match_expectation(type, {"g","4"}, {"4","g","g"}) && match_expectation(type, {"-1","g"}, {"g","g","-1"}) && match_expectation(type, {"1","3"}, {"3","g","1"}) && match_expectation(type, {"-2","2"}, {"2","g","-2","h"}) && match_expectation(type, {"-3","g"}, {"g","g","Wp","-4"}) && match_expectation(type, {"1","-2"}, {"Wm","-1","g","1"}) && match_expectation(type, {"-1","g"}, {"1","-1","-1"}) && match_expectation(type, {"1","-4"}, {"Z_photon_mix","-4","g","1"}) // higgs inside uno && match_expectation(type, {"-1","g"}, {"g","h","-1","g"}) && match_expectation(type, {"-1","1"}, {"g","h","-1","1"}) && match_expectation(type, {"g","2"}, {"g","2","h","g"}) && match_expectation(type, {"-1","1"}, {"-1","1","h","g"}) // higgs outside uno && match_expectation(type, {"-1","g"}, {"h","g","-1","g"}) && match_expectation(type, {"-1","1"}, {"-1","1","g","h"}) // higgs inside qqx && match_expectation(type, {"g","g"}, {"-1","h","1","g","g"}) && match_expectation(type, {"g","g"}, {"g","-1","h","1","g"}) && match_expectation(type, {"g","g"}, {"g","g","2","h","-2"}) // higgs outside qqx && match_expectation(type, {"g","g"}, {"h","-1","1","g","g"}) && match_expectation(type, {"g","g"}, {"g","g","2","-2","h"}) // 4j - two uno && match_expectation(type, {"-2","2"}, {"g","-2","2","g"}) && match_expectation(type, {"1","3"}, {"g","1","h","3","g"}) && match_expectation(type, {"1","2"}, {"g","1","3","Wp","g"}) && match_expectation(type, {"1","-2"}, {"g","Wm","1","-1","g"}) && match_expectation(type, {"3","2"}, {"g","3","Z_photon_mix","2","g"}) // 4j - two gluon outside && match_expectation(type, {"g","4"}, {"g","4","g","g"}) && match_expectation(type, {"1","3"}, {"1","3","h","g","g"}) && match_expectation(type, {"1","2"}, {"1","3","g","Wp","g"}) && match_expectation(type, {"1","-2"}, {"1","Wm","-1","g","g"}) && match_expectation(type, {"-1","g"}, {"g","g","-1","g"}) && match_expectation(type, {"1","3"}, {"g","g","1","3","h"}) && match_expectation(type, {"1","2"}, {"g","g","1","Wp","3"}) && match_expectation(type, {"1","-2"}, {"Wm","g","g","1","-1"}) && match_expectation(type, {"-1","2"}, {"g","g","-1","Z_photon_mix","2"}) // 4j - ggx+uno && match_expectation(type, {"g","4"}, {"1","-1","4","g"}) && match_expectation(type, {"2","g"}, {"g","2","-3","3"}) && match_expectation(type, {"g","4"}, {"1","-1","h","4","g"}) && match_expectation(type, {"2","g"}, {"g","2","-3","3","h"}) && match_expectation(type, {"g","4"}, {"Wp","1","-1","3","g"}) && match_expectation(type, {"2","g"}, {"g","2","-4","Wp","3"}) && match_expectation(type, {"g","4"}, {"2","Wm","-1","4","g"}) && match_expectation(type, {"2","g"}, {"g","2","Wp","-3","4"}) && match_expectation(type, {"-4","g"}, {"g","-4","-3","3","Z_photon_mix"}) // 3j - crossing+uno && match_expectation(type, {"1","4"}, {"g","4","1"}) && match_expectation(type, {"1","4"}, {"4","1","g"}) && match_expectation(type, {"1","4"}, {"g","h","4","1"}) && match_expectation(type, {"-1","-3"},{"Wm","g","-4","-1"}) && match_expectation(type, {"1","4"}, {"3","1","Wp","g"}) && match_expectation(type, {"1","4"}, {"3","1","g","h"}) && match_expectation(type, {"2","3"}, {"3","2","Z_photon_mix","g"}) // 3j - crossing+qqx && match_expectation(type, {"1","g"}, {"-1","1","g","1"}) && match_expectation(type, {"1","g"}, {"-1","1","1","g"}) && match_expectation(type, {"g","1"}, {"1","g","1","-1"}) && match_expectation(type, {"g","1"}, {"g","1","1","-1"}) && match_expectation(type, {"1","g"}, {"2","-2","g","1"}) && match_expectation(type, {"1","g"}, {"2","-2","1","g"}) && match_expectation(type, {"g","1"}, {"1","g","-2","2"}) && match_expectation(type, {"g","1"}, {"g","1","-2","2"}) && match_expectation(type, {"1","g"}, {"-1","1","h","g","1"}) && match_expectation(type, {"1","g"}, {"-1","h","1","1","g"}) && match_expectation(type, {"g","1"}, {"1","g","1","h","-1"}) && match_expectation(type, {"g","1"}, {"h","g","1","1","-1"}) && match_expectation(type, {"1","g"}, {"2","-2","1","g","h"}) && match_expectation(type, {"g","1"}, {"g","h","1","-2","2"}) && match_expectation(type, {"1","g"}, {"Wp","3","-4","g","1"}) && match_expectation(type, {"3","g"}, {"-2","Wm","1","3","g"}) && match_expectation(type, {"g","1"}, {"1","g","Wm","-3","4"}) && match_expectation(type, {"g","-3"}, {"g","-3","-1","Wp","2"}) && match_expectation(type, {"g","2"}, {"2","g","Z_photon_mix","4","-4"}) // 4j- gluon in qqx && match_expectation(type, {"g","1"}, {"1","g","-1","1"}) && match_expectation(type, {"1","g"}, {"1","-1","g","1"}) && match_expectation(type, {"g","1"}, {"1","g","Wm","-2","1"}) && match_expectation(type, {"2","g"}, {"2","-2","g","Wp","1"}) && match_expectation(type, {"g","g"}, {"Wp","3","g","-4","g"}) && match_expectation(type, {"1","g"}, {"1","h","-1","g","1"}) && match_expectation(type, {"3","g"}, {"3","1","g","Z_photon_mix","-1"}) // 6j - two qqx && match_expectation(type, {"g","g"}, {"1","-1","g","g","1","-1"}) && match_expectation(type, {"g","g"}, {"1","-1","g","1","-1","g"}) && match_expectation(type, {"g","g"}, {"g","1","-1","g","1","-1"}) && match_expectation(type, {"g","g"}, {"g","1","-1","1","-1","g"}) && match_expectation(type, {"g","g"}, {"g","1","1","-1","-1","g"}) && match_expectation(type, {"g","g"}, {"h","1","-1","g","g","1","-1"}) && match_expectation(type, {"g","g"}, {"1","Wp","-2","g","1","-1","g"}) && match_expectation(type, {"g","g"}, {"g","1","Wp","-1","g","1","-2"}) && match_expectation(type, {"g","g"}, {"g","1","-1","Wm","2","-1","g"}) && match_expectation(type, {"g","g"}, {"g","1","2","-1","Wm","-1","g"}) && match_expectation(type, {"g","g"}, {"2","-2","g","-1","1","Z_photon_mix","g"}) // random stuff (can be non-physical) && match_expectation(type, {"g","g"}, {"1","-2","2","-1"}) // != 2 qqx && match_expectation(type, {"g","g"}, {"1","-2","2","g"}) // could be qqx && match_expectation(type, {"e+","e-"},{"1","-1"}) // bad initial state && match_expectation(type, {"1","e-"}, {"g","1","Wm"}) // bad initial state && match_expectation(type, {"h","g"}, {"g","g"}) // bad initial state && match_expectation(type, {"-1","g"}, {"-1","1","1"}) // bad qqx && match_expectation(type, {"-1","g"}, {"1","1","-1"}) // crossing in bad qqx && match_expectation(type, {"-1","g"}, {"-2","1","1","Wp"}) // bad qqx && match_expectation(type, {"1","2"}, {"1","-1","g","g","g","2"}) // bad qqx && match_expectation(type, {"1","2"}, {"1","-1","-2","g","g","2"}) // gluon in bad qqx && match_expectation(type, {"g","g"}, {"-1","2","g","g"}) // wrong back qqx && match_expectation(type, {"g","g"}, {"g","g","2","1"}) // wrong forward qqx && match_expectation(type, {"g","g"}, {"g","-2","1","g"}) // wrong central qqx && match_expectation(type, {"1","g"}, {"1","-2","g","g","Wp"}) // extra quark && match_expectation(type, {"g","1"}, {"g","g","-2","1","Wp"}) // extra quark && match_expectation(type, {"g","1"}, {"g","g","Wp","-2","1"}) // extra quark && match_expectation(type, {"g","1"}, {"g","-2","1","g","Wp"}) // extra quark && match_expectation(type, {"g","g"}, {"g","g","g","-2","1","-1","Wp"}) // extra quark && match_expectation(type, {"1","g"}, {"g","Wp","1","-2","g"}) // extra quark && match_expectation(type, {"g","g"}, {"1","-1","-2","g","g","g","Wp"}) // extra quark ; } // Two boson states, that are currently not implemented bool check_bad_FS(){ auto type{ HEJ::event_type::bad_final_state}; return match_expectation(type, {"g","g"}, {"g","h","h","g"}) && match_expectation(type, {"g","g"}, {"h","g","h","g"}) && match_expectation(type, {"g","-1"}, {"g","h","Wp","-2"}) && match_expectation(type, {"-3","-1"},{"-4","g","Wp","Wp","-2"}) && match_expectation(type, {"-4","-1"},{"-3","Wp","g","Wm","-2"}) && match_expectation(type, {"-4","-1"},{"g","-3","Wp","Wm","-2"}) && match_expectation(type, {"-4","-1"},{"-4","g","11","-11","-2"}) && match_expectation(type, {"-4","-1"},{"-4","g","-13","g","-2"}) && match_expectation(type, {"3","-2"}, {"Wp","3","Wm","g","g","g","-2"}, -13) && match_expectation(type, {"1","2"},{"1","g","Z_photon_mix","Z_photon_mix","2"}) && match_expectation(type, {"-4","g"},{"-4","Z_photon_mix","g","Z_photon_mix","g"}) && match_expectation(type, {"g","-2"}, {"g","Z_photon_mix","Wm","-1"}) && match_expectation(type, {"3","1"},{"3","g","h","Z_photon_mix","1"}) ; } - // not 2 jets - bool check_not_2_jets(){ - auto type{ HEJ::event_type::no_2_jets}; + // not enough jets + bool check_not_enough_jets(){ + auto type{ HEJ::event_type::not_enough_jets}; return match_expectation(type, {"g","g"}, {}) && match_expectation(type, {"1","-1"}, {}) && match_expectation(type, {"g","-1"}, {"-1"}) && match_expectation(type, {"g","g"}, {"g"}) ; } // not implemented processes bool check_not_implemented(){ return check_fkl(false) && check_uno(false) && check_extremal_qqx(false) && check_central_qqx(false); } } // namespace int main() { // tests for "no false negatives" // i.e. all HEJ-configurations get classified correctly if(!check_fkl()) return EXIT_FAILURE; if(!check_uno()) return EXIT_FAILURE; if(!check_extremal_qqx()) return EXIT_FAILURE; if(!check_central_qqx()) return EXIT_FAILURE; // test for "no false positive" // i.e. non-resummable gives non-resummable if(!check_non_resummable()) return EXIT_FAILURE; if(!check_bad_FS()) return EXIT_FAILURE; - if(!check_not_2_jets()) return EXIT_FAILURE; + if(!check_not_enough_jets()) return EXIT_FAILURE; if(!check_not_implemented()) return EXIT_FAILURE; return EXIT_SUCCESS; }