diff --git a/include/HEJ/Event.hh b/include/HEJ/Event.hh index 336f49b..1a82773 100644 --- a/include/HEJ/Event.hh +++ b/include/HEJ/Event.hh @@ -1,372 +1,372 @@ /** \file * \brief Declares the Event class and helpers * * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #pragma once #include #include #include #include #include #include #include #include "boost/iterator/filter_iterator.hpp" #include "fastjet/ClusterSequence.hh" #include "fastjet/PseudoJet.hh" -#include "HEJ/EWConstants.hh" #include "HEJ/Parameters.hh" #include "HEJ/Particle.hh" #include "HEJ/event_types.hh" namespace LHEF { class HEPEUP; class HEPRUP; } namespace fastjet { class JetDefinition; } namespace HEJ { + struct EWConstants; struct RNG; struct UnclusteredEvent; /** @brief 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: class EventData; //! Iterator over partons using ConstPartonIterator = boost::filter_iterator< bool (*)(Particle const &), std::vector::const_iterator >; //! Reverse Iterator over partons using ConstReversePartonIterator = std::reverse_iterator< ConstPartonIterator>; //! No default Constructor Event() = delete; //! Event Constructor adding jet clustering to an unclustered event //! @deprecated UnclusteredEvent will be replaced by EventData in HEJ 2.2.0 [[deprecated("UnclusteredEvent will be replaced by EventData")]] Event( UnclusteredEvent const & ev, fastjet::JetDefinition const & jet_def, double min_jet_pt ); //! @name Particle Access //! @{ //! Incoming particles std::array const & incoming() const{ return incoming_; } //! Outgoing particles std::vector const & outgoing() const{ return outgoing_; } //! Iterator to the first outgoing parton ConstPartonIterator begin_partons() const; //! Iterator to the first outgoing parton ConstPartonIterator cbegin_partons() const; //! Iterator to the end of the outgoing partons ConstPartonIterator end_partons() const; //! Iterator to the end of the outgoing partons ConstPartonIterator cend_partons() const; //! Reverse Iterator to the first outgoing parton ConstReversePartonIterator rbegin_partons() const; //! Reverse Iterator to the first outgoing parton ConstReversePartonIterator crbegin_partons() const; //! Reverse Iterator to the first outgoing parton ConstReversePartonIterator rend_partons() const; //! Reverse Iterator to the first outgoing parton ConstReversePartonIterator crend_partons() const; //! 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 decays_; } //! The jets formed by the outgoing partons, sorted in rapidity std::vector const & jets() const{ return jets_; } //! @} //! @name Weight variations //! @{ //! All chosen parameter, i.e. scale choices (const version) Parameters const & parameters() const{ return parameters_; } //! All chosen parameter, i.e. scale choices Parameters & parameters(){ return parameters_; } //! Central parameter choice (const version) EventParameters const & central() const{ return parameters_.central; } //! Central parameter choice EventParameters & central(){ return parameters_.central; } //! Parameter (scale) variations (const version) std::vector const & variations() const{ return parameters_.variations; } //! Parameter (scale) variations std::vector & variations(){ return parameters_.variations; } //! Parameter (scale) variation (const version) /** * @param i Index of the requested variation */ EventParameters const & variations(std::size_t i) const{ return parameters_.variations.at(i); } //! Parameter (scale) variation /** * @param i Index of the requested variation */ EventParameters & variations(std::size_t i){ return parameters_.variations.at(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); } //! particle_jet_indices() of the Event jets() std::vector particle_jet_indices() const { return 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_; } //! Give colours to each particle /** * @returns true if new colours are generated, i.e. same as is_resummable() * @details Colour ordering is done according to leading colour in the MRK * limit, see \cite Andersen:2011zd. This only affects \ref * is_resummable() "HEJ" configurations, all other \ref event_type * "EventTypes" will be ignored. * @note This overwrites all previously set colours. */ bool generate_colours(RNG & /*ran*/); //! Check that current colours are leading in the high energy limit /** * @details Checks that the colour configuration can be split up in * multiple, rapidity ordered, non-overlapping ladders. Such * configurations are leading in the MRK limit, see * \cite Andersen:2011zd * * @note This is _not_ to be confused with \ref is_resummable(), however * for all resummable states it is possible to create a leading colour * configuration, see generate_colours() */ bool is_leading_colour() const; /** * @brief Check if given event could have been produced by HEJ * @details A HEJ state has to fulfil: * 1. type() has to be \ref is_resummable() "resummable" * 2. Soft radiation in the tagging jets contributes at most to * `max_ext_soft_pt_fraction` of the total jet \f$ p_\perp \f$ * * @note This is true for any resummed stated produced by the * EventReweighter or any \ref is_resummable() "resummable" Leading * Order state. * * @param max_frac Maximum transverse momentum fraction from soft radiation * in extremal jets * @param min_pt Absolute minimal \f$ p_\perp \f$, * \b deprecated use max_frac instead * @return True if this state could have been produced by HEJ */ bool valid_hej_state( double max_frac, double min_pt = 0.) const; private: //! \internal //! @brief Construct Event explicitly from input. /** This is only intended to be called from EventData. * * \warning The input is taken _as is_, sorting and classification has to be * done externally, i.e. by EventData */ Event( std::array && incoming, std::vector && outgoing, std::unordered_map> && decays, Parameters && parameters, fastjet::JetDefinition const & jet_def, double min_jet_pt ); //! Iterator over partons (non-const) using PartonIterator = boost::filter_iterator< bool (*)(Particle const &), std::vector::iterator >; //! Reverse Iterator over partons (non-const) using ReversePartonIterator = std::reverse_iterator; //! Iterator to the first outgoing parton (non-const) PartonIterator begin_partons(); //! Iterator to the end of the outgoing partons (non-const) PartonIterator end_partons(); //! Reverse Iterator to the first outgoing parton (non-const) ReversePartonIterator rbegin_partons(); //! Reverse Iterator to the first outgoing parton (non-const) ReversePartonIterator rend_partons(); std::array incoming_; std::vector outgoing_; std::unordered_map> decays_; std::vector jets_; Parameters parameters_; fastjet::ClusterSequence cs_; double min_jet_pt_; event_type::EventType type_; }; // end class Event //! Class to store general Event setup, i.e. Phase space and weights class Event::EventData { public: //! Default Constructor EventData() = default; //! Constructor from LesHouches event information EventData(LHEF::HEPEUP const & hepeup); //! Constructor with all values given EventData( std::array incoming, std::vector outgoing, std::unordered_map> decays, Parameters parameters ): incoming(std::move(incoming)), outgoing(std::move(outgoing)), decays(std::move(decays)), parameters(std::move(parameters)) {} //! Generate an Event from the stored EventData. /** * @details Do jet clustering and classification. * Use this to generate an Event. * * @note Calling this function destroys EventData * * @param jet_def Jet definition * @param min_jet_pt minimal \f$p_T\f$ for each jet * * @returns Full clustered and classified event. */ Event cluster( fastjet::JetDefinition const & jet_def, double min_jet_pt); //! Alias for cluster() Event operator()( fastjet::JetDefinition const & jet_def, double const min_jet_pt){ return cluster(jet_def, min_jet_pt); } //! Sort particles in rapidity void sort(); //! Reconstruct intermediate particles from final-state leptons /** * Final-state leptons are created from virtual photons, W, or Z bosons. * This function tries to reconstruct such intermediate bosons if they * are not part of the event record. */ - void reconstruct_intermediate(EWConstants const & ew_parameters); + void reconstruct_intermediate(EWConstants const &); //! Incoming particles std::array incoming; //! Outcoing particles std::vector outgoing; //! Particle decays in the format {outgoing index, decay products} std::unordered_map> decays; //! Parameters, e.g. scale or inital weight Parameters parameters; }; // end class EventData //! Print Event std::ostream& operator<<(std::ostream & os, Event const & ev); //! 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 * /*heprup*/); // put deprecated warning at the end, so don't get the warning inside Event.hh, // additionally doxygen can not identify [[deprecated]] correctly struct [[deprecated("UnclusteredEvent will be replaced by EventData")]] UnclusteredEvent; //! An event before jet clustering //! @deprecated UnclusteredEvent will be replaced by EventData in HEJ 2.2.0 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 */ }; } // namespace HEJ diff --git a/src/Event.cc b/src/Event.cc index 44e8f79..57cd345 100644 --- a/src/Event.cc +++ b/src/Event.cc @@ -1,1303 +1,1304 @@ /** * \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/EWConstants.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]); } //! true if supported decay bool valid_decay(std::vector const & decays){ return valid_W_decay(+1, decays) || // Wp valid_W_decay(-1, decays) || // Wm valid_Z_decay( decays) // Z/gamma ; } /// @name helper functions to determine event type //@{ /** * \brief check if final state valid for HEJ * * check final state has the expected number of valid decays for bosons * and all the rest are quarks or gluons */ bool final_state_ok(Event const & ev){ size_t invalid_decays = ev.decays().size(); std::vector const & outgoing = ev.outgoing(); for( size_t i=0; i0?+1:-1, decay->second) ){ --invalid_decays; } else return false; } // Higgs decays (optional) else if(out.type == ParticleID::h){ if(decay != ev.decays().cend()) --invalid_decays; } // Z decays (required) else if(out.type == ParticleID::Z_photon_mix){ if( decay != ev.decays().cend() && valid_Z_decay(decay->second) ){ --invalid_decays; } else return false; } } else if(! is_parton(out.type)) return false; } // any invalid decays? if(invalid_decays != 0) 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) { 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; } } if(bosons.size() == 2) { if(bosons[0].type == ParticleID::Wp && bosons[1].type == ParticleID::Wp) { return FKL; } } 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){ if(!qqx && (in==-1 || in== 2 || in==-3 || in== 4)) return out== (in-1); if( qqx && (in== 1 || in==-2 || in== 3 || in==-4)) 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){ if(!qqx && (in== 1 || in==-2 || in== 3 || in==-4)) return out== (in+1); if( qqx && (in==-1 || in== 2 || in==-3 || in== 4)) 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; } /** * \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) ){ 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; } //! 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; // 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; } // 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; // 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; } // 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; } return possible | qqxmid; } return non_resummable; } /** * \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; // 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()}; 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{})); auto const bosons{ filter_AWZH_bosons(ev.outgoing()) }; // keep track of potential W couplings, at the end the sum should be 0 int remaining_Wp = 0; int remaining_Wm = 0; for(auto const & boson : bosons){ if(boson.type == ParticleID::Wp) ++remaining_Wp; else if(boson.type == ParticleID::Wm) ++remaining_Wm; } size_t final_type = ~(no_2_jets | bad_final_state); // check forward impact factor int W_change = 0; final_type &= possible_impact_factors( in.front().type, begin_out, end_out.base(), W_change, bosons, 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; // check backward impact factor W_change = 0; final_type &= possible_impact_factors( in.back().type, end_out, std::make_reverse_iterator(begin_out), W_change, bosons, 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; // check central emissions W_change = 0; final_type &= possible_central( begin_out, end_out.base(), W_change, bosons ); 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(bosons)) != 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 { // use valid_X_decay to determine boson type ParticleID reconstruct_type(std::vector const & progeny) { if(valid_W_decay(+1, progeny)) { return ParticleID::Wp; } if(valid_W_decay(-1, progeny)) { return ParticleID::Wm; } if(valid_Z_decay(progeny)) { return ParticleID::Z_photon_mix; } throw not_implemented{ "final state with decay X -> " + name(progeny[0].type) + " + " + name(progeny[1].type) }; } // reconstruct particle with explicit ParticleID Particle reconstruct_boson( std::vector const & progeny, ParticleID const & type ) { Particle progenitor; progenitor.p = progeny[0].p + progeny[1].p; progenitor.type = type; return progenitor; } // reconstruct via call to reconstruct_type Particle reconstruct_boson(std::vector const & progeny) { Particle progenitor {reconstruct_boson(progeny, reconstruct_type(progeny))}; assert(is_AWZH_boson(progenitor)); return progenitor; } typedef std::vector< std::vector > GroupedParticles; typedef std::pair > Decay; typedef std::vector< Decay > Decays; // return groups of reconstructable progeny std::vector group_progeny ( std::vector & leptons ) { /** Warning: The partition in to charged/neutral leptons is valid ONLY for WW. **/ assert(leptons.size() == 4); auto const begin_neutrino = std::partition( begin(leptons), end(leptons), [](Particle const & p) {return !is_anyneutrino(p);} ); std::vector neutrinos (begin_neutrino, end(leptons)); leptons.erase(begin_neutrino, end(leptons)); std::sort(begin(leptons), end(leptons), type_less{}); std::sort(begin(neutrinos), end(neutrinos), type_less{}); if(leptons.size() != 2 && neutrinos.size() != 2) { return {}; } assert(leptons.size() == 2 && neutrinos.size() == 2); std::vector< GroupedParticles > candidate_grouping { {{leptons[0], neutrinos[0]}, {leptons[1], neutrinos[1]}}, {{leptons[1], neutrinos[0]}, {leptons[0], neutrinos[1]}} }; // erase groupings containing invalid decays candidate_grouping.erase( std::remove_if(begin(candidate_grouping), end(candidate_grouping), [] (GroupedParticles & candidate) -> bool { return ! std::accumulate( cbegin(candidate), cend(candidate), true, [] (bool acc_valid, std::vector const & decay) -> bool { return acc_valid && valid_decay(decay); } ); } ), candidate_grouping.end() ); return candidate_grouping; } // 'best' decay ordering measure double decay_measure(Particle reconstructed, EWConstants const & params) { ParticleProperties ref = params.prop(reconstructed.type); return std::abs(reconstructed.p.m() - ref.mass); } // decay_measure accumulated over decays double decay_measure(Decays const & decays, EWConstants const & params) { return std::accumulate( cbegin(decays), cend(decays), 0., [¶ms] (double dm, Decay const & decay) -> double { return dm + decay_measure(decay.first, params); } ); } // select best combination of decays for the event Decays select_decays ( std::vector & leptons, EWConstants const & ew_parameters ) { std::vector groupings = group_progeny(leptons); std::vector valid_decays; valid_decays.reserve(groupings.size()); // Reconstruct all groupings for(GroupedParticles const & group : groupings) { Decays decays; for(auto const & progeny : group) { decays.emplace_back(make_pair(reconstruct_boson(progeny), progeny)); } valid_decays.emplace_back(decays); } if (valid_decays.empty()) { throw not_implemented{"No supported intermediate reconstruction available"}; } else if (valid_decays.size() == 1) { return valid_decays[0]; } else { // select decay with smallest decay_measure auto selected = std::min_element(cbegin(valid_decays), cend(valid_decays), [&ew_parameters] (auto const & d1, auto const & d2) -> bool { return decay_measure(d1, ew_parameters) < decay_measure(d2, ew_parameters); } ); return *selected; } } } // namespace void Event::EventData::reconstruct_intermediate(EWConstants const & ew_parameters) { auto const begin_leptons = std::partition( begin(outgoing), end(outgoing), [](Particle const & p) {return !is_anylepton(p);} ); if(begin_leptons == end(outgoing)) return; assert(is_anylepton(*begin_leptons)); std::vector leptons(begin_leptons, end(outgoing)); outgoing.erase(begin_leptons, end(outgoing)); if(leptons.empty()) { return; } // nothing to do else if(leptons.size() == 2) { outgoing.emplace_back(reconstruct_boson(leptons)); std::sort(begin(leptons), end(leptons), type_less{}); decays.emplace(outgoing.size()-1, std::move(leptons)); } else if(leptons.size() == 4) { // select_decays only supports WpWp Decays valid_decays = select_decays(leptons, ew_parameters); for(auto &decay : valid_decays) { outgoing.emplace_back(decay.first); std::sort(begin(decay.second), end(decay.second), type_less{}); decays.emplace(outgoing.size()-1, std::move(decay.second)); } } else { throw not_implemented { std::to_string(leptons.size()) + " leptons in the final state" }; } } 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 max_ext_soft_pt_fraction, 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() <= max_ext_soft_pt_fraction; } } // namespace // this should work with multiple types bool Event::valid_hej_state(double const max_frac, 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, max_frac, min_pt) ) return false; ++part_begin; ++idx_begin; if( !valid_parton(jets(), *part_end, *idx_end, max_frac, 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, max_frac, min_pt) ) return false; ++part_begin; ++idx_begin; } if( type() & (unof | qqxexf) ){ if( !valid_parton(jets(), *part_end, *idx_end, max_frac, 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, max_frac,min_pt) && valid_parton(jets(),*(++begin_qqx),*(++idx_begin),max_frac,min_pt) )) return false; } return true; } Event::ConstPartonIterator Event::begin_partons() const { return cbegin_partons(); } Event::ConstPartonIterator Event::cbegin_partons() const { return boost::make_filter_iterator( static_cast(is_parton), cbegin(outgoing()), cend(outgoing()) ); } Event::ConstPartonIterator Event::end_partons() const { return cend_partons(); } Event::ConstPartonIterator Event::cend_partons() const { return boost::make_filter_iterator( static_cast(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 boost::make_filter_iterator( static_cast(is_parton), begin(outgoing_), end(outgoing_) ); } Event::PartonIterator Event::end_partons() { return boost::make_filter_iterator( static_cast(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/t/cmp_events.cc b/t/cmp_events.cc index bd4c540..057678f 100644 --- a/t/cmp_events.cc +++ b/t/cmp_events.cc @@ -1,67 +1,68 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #include "HEJ/Event.hh" +#include "HEJ/EWConstants.hh" #include "HEJ/stream.hh" #include "HEJ/utility.hh" #include "LHEF/LHEF.h" namespace { const HEJ::ParticleProperties Wprop{80.385, 2.085}; const HEJ::ParticleProperties Zprop{91.187, 2.495}; const HEJ::ParticleProperties Hprop{125, 0.004165}; constexpr double vev = 246.2196508; const HEJ::EWConstants ew_parameters{vev, Wprop, Zprop, Hprop}; } int main(int argn, char** argv) { if(argn != 3){ std::cerr << "Usage: cmp_events eventfile eventfile_no_boson "; return EXIT_FAILURE; } HEJ::istream in{argv[1]}; LHEF::Reader reader{in}; HEJ::istream in_no_boson{argv[2]}; LHEF::Reader reader_no_boson{in_no_boson}; while(reader.readEvent()) { if(! reader_no_boson.readEvent()) { std::cerr << "wrong number of events in " << argv[2] << '\n'; return EXIT_FAILURE; } const auto is_AWZH = [](HEJ::Particle const & p) { return is_AWZH_boson(p); }; const HEJ::Event::EventData event_data{reader.hepeup}; const auto boson = std::find_if( begin(event_data.outgoing), end(event_data.outgoing), is_AWZH ); if(boson == end(event_data.outgoing)) { std::cerr << "no boson in event\n"; return EXIT_FAILURE; } HEJ::Event::EventData data_no_boson{reader_no_boson.hepeup}; data_no_boson.reconstruct_intermediate(ew_parameters); const auto new_boson = std::find_if( begin(event_data.outgoing), end(event_data.outgoing), is_AWZH ); if(new_boson == end(data_no_boson.outgoing)) { std::cerr << "failed to find reconstructed boson\n"; return EXIT_FAILURE; } if(new_boson->type != boson->type) { std::cerr << "reconstructed wrong boson type\n"; return EXIT_FAILURE; } if(!HEJ::nearby(new_boson->p, boson->p)) { std::cerr << "reconstructed boson momentum is off\n"; return EXIT_FAILURE; } } return EXIT_SUCCESS; } diff --git a/t/hej_test.cc b/t/hej_test.cc index ecf2dd2..389b7ed 100644 --- a/t/hej_test.cc +++ b/t/hej_test.cc @@ -1,624 +1,625 @@ /** * \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 +#include "HEJ/EWConstants.hh" #include "HEJ/Particle.hh" #include "HEJ/PDG_codes.hh" namespace { const HEJ::ParticleProperties Wprop{80.385, 2.085}; const HEJ::ParticleProperties Zprop{91.187, 2.495}; const HEJ::ParticleProperties Hprop{125, 0.004165}; constexpr double vev = 246.2196508; const HEJ::EWConstants ew_parameters{vev, Wprop, Zprop, Hprop}; } HEJ::Event::EventData get_process(int const njet, int const pos_boson){ using namespace HEJ::pid; HEJ::Event::EventData ev; if(njet == 0){ // jet idx: -1 -1 ev.outgoing.push_back({gluon, { -24, 12, -57, 63}, {}}); ev.outgoing.push_back({gluon, { 24, -12, 41, 49}, {}}); ev.incoming[0] = {gluon, { 0, 0, -64, 64}, {}}; ev.incoming[1] = {gluon, { 0, 0, 48, 48}, {}}; return ev; } else if(njet == 1){ // jet idx: 0 -1 -1 ev.outgoing.push_back({gluon, { 23, 28, -44, 57}, {}}); ev.outgoing.push_back({gluon, { -11, -24, -12, 29}, {}}); ev.outgoing.push_back({gluon, { -12, -4, 39, 41}, {}}); ev.incoming[0] = {gluon, { 0, 0, -72, 72}, {}}; ev.incoming[1] = {gluon, { 0, 0, 55, 55}, {}}; return ev; } else if(njet == 2){ switch(pos_boson){ case 0: ev.outgoing.push_back({higgs, { 198, 33, -170, 291}, {}}); ev.outgoing.push_back({gluon, {-154, 68, 44, 174}, {}}); ev.outgoing.push_back({gluon, { -44, -101, 88, 141}, {}}); ev.incoming[0] = {gluon, { 0, 0, -322, 322}, {}}; ev.incoming[1] = {gluon, { 0, 0, 284, 284}, {}}; return ev; case 1: ev.outgoing.push_back({gluon, { -6, 82, -159, 179}, {}}); ev.outgoing.push_back({higgs, { 195, -106, 74, 265}, {}}); ev.outgoing.push_back({gluon, {-189, 24, 108, 219}, {}}); ev.incoming[0] = {gluon, { 0, 0, -320, 320}, {}}; ev.incoming[1] = {gluon, { 0, 0, 343, 343}, {}}; return ev; case 2: ev.outgoing.push_back({gluon, { -80, -80, -140, 180}, {}}); ev.outgoing.push_back({gluon, { -60, -32, 0, 68}, {}}); ev.outgoing.push_back({higgs, { 140, 112, 177, 281}, {}}); ev.incoming[0] = {gluon, { 0, 0, -246, 246}, {}}; ev.incoming[1] = {gluon, { 0, 0, 283, 283}, {}}; return ev; default: ev.outgoing.push_back({gluon, { -72, 24, 18, 78}, {}}); ev.outgoing.push_back({gluon, { 72, -24, 74, 106}, {}}); ev.incoming[0] = {gluon, { 0, 0, -46, 46}, {}}; ev.incoming[1] = {gluon, { 0, 0, 138, 138}, {}}; return ev; } } if(njet == 3){ switch(pos_boson){ case 0: ev.outgoing.push_back({higgs, { 152, -117, -88, 245}, {}}); ev.outgoing.push_back({gluon, {-146, 62, -11, 159}, {}}); ev.outgoing.push_back({gluon, { 126, -72, 96, 174}, {}}); ev.outgoing.push_back({gluon, {-132, 127, 144, 233}, {}}); ev.incoming[0] = {gluon, { 0, 0, -335, 335}, {}}; ev.incoming[1] = {gluon, { 0, 0, 476, 476}, {}}; return ev; case 1: ev.outgoing.push_back({gluon, {-191, 188, -128, 297}, {}}); ev.outgoing.push_back({higgs, { 199, 72, -76, 257}, {}}); ev.outgoing.push_back({gluon, { 184, -172, -8, 252}, {}}); ev.outgoing.push_back({gluon, {-192, -88, 54, 218}, {}}); ev.incoming[0] = {gluon, { 0, 0, -591, 591}, {}}; ev.incoming[1] = {gluon, { 0, 0, 433, 433}, {}}; return ev; case 2: ev.outgoing.push_back({gluon, { -42, 18, -49, 67}, {}}); ev.outgoing.push_back({gluon, { -12, -54, -28, 62}, {}}); ev.outgoing.push_back({higgs, { 99, 32, -16, 163}, {}}); ev.outgoing.push_back({gluon, { -45, 4, 72, 85}, {}}); ev.incoming[0] = {gluon, { 0, 0, -199, 199}, {}}; ev.incoming[1] = {gluon, { 0, 0, 178, 178}, {}}; return ev; case 3: ev.outgoing.push_back({gluon, { -65, -32, -76, 105}, {}}); ev.outgoing.push_back({gluon, { -22, 31, -34, 51}, {}}); ev.outgoing.push_back({gluon, { -12, -67, -36, 77}, {}}); ev.outgoing.push_back({higgs, { 99, 68, -4, 173}, {}}); ev.incoming[0] = {gluon, { 0, 0, -278, 278}, {}}; ev.incoming[1] = {gluon, { 0, 0, 128, 128}, {}}; return ev; default: ev.outgoing.push_back({gluon, { -90, -135, 30, 165}, {}}); ev.outgoing.push_back({gluon, {-108, 198, 76, 238}, {}}); ev.outgoing.push_back({gluon, { 198, -63, 126, 243}, {}}); ev.incoming[0] = {gluon, { 0, 0, -207, 207}, {}}; ev.incoming[1] = {gluon, { 0, 0, 439, 439}, {}}; return ev; } } if(njet == 4){ switch(pos_boson){ case 0: ev.outgoing.push_back({higgs, { 199, 72, -76, 257}, {}}); ev.outgoing.push_back({gluon, {-200, -155, -64, 261}, {}}); ev.outgoing.push_back({gluon, { 198, 194, 57, 283}, {}}); ev.outgoing.push_back({gluon, { 1, 32, 8, 33}, {}}); ev.outgoing.push_back({gluon, {-198, -143, 186, 307}, {}}); ev.incoming[0] = {gluon, { 0, 0, -515, 515}, {}}; ev.incoming[1] = {gluon, { 0, 0, 626, 626}, {}}; return ev; case 1: ev.outgoing.push_back({gluon, { 198, 61, -162, 263}, {}}); ev.outgoing.push_back({higgs, { 199, 72, -76, 257}, {}}); ev.outgoing.push_back({gluon, {-200, 135, 144, 281}, {}}); ev.outgoing.push_back({gluon, {-198, -186, 171, 321}, {}}); ev.outgoing.push_back({gluon, { 1, -82, 122, 147}, {}}); ev.incoming[0] = {gluon, { 0, 0, -535, 535}, {}}; ev.incoming[1] = {gluon, { 0, 0, 734, 734}, {}}; return ev; case 2: ev.outgoing.push_back({gluon, {-180, -27, -164, 245}, {}}); ev.outgoing.push_back({gluon, {-108, 78, -36, 138}, {}}); ev.outgoing.push_back({higgs, { 196, -189, 68, 307}, {}}); ev.outgoing.push_back({gluon, {-107, 136, 76, 189}, {}}); ev.outgoing.push_back({gluon, { 199, 2, 178, 267}, {}}); ev.incoming[0] = {gluon, { 0, 0, -512, 512}, {}}; ev.incoming[1] = {gluon, { 0, 0, 634, 634}, {}}; return ev; case 3: ev.outgoing.push_back({gluon, { -12, -30, -84, 90}, {}}); ev.outgoing.push_back({gluon, { -72, 22, -96, 122}, {}}); ev.outgoing.push_back({gluon, { 68, 0, -51, 85}, {}}); ev.outgoing.push_back({higgs, { 64, 72, -81, 177}, {}}); ev.outgoing.push_back({gluon, { -48, -64, 84, 116}, {}}); ev.incoming[0] = {gluon, { 0, 0, -409, 409}, {}}; ev.incoming[1] = {gluon, { 0, 0, 181, 181}, {}}; return ev; case 4: ev.outgoing.push_back({gluon, { -72, -49, -72, 113}, {}}); ev.outgoing.push_back({gluon, { -48, 0, -36, 60}, {}}); ev.outgoing.push_back({gluon, { -12, 54, -36, 66}, {}}); ev.outgoing.push_back({gluon, { 68, -77, -56, 117}, {}}); ev.outgoing.push_back({higgs, { 64, 72, -81, 177}, {}}); ev.incoming[0] = {gluon, { 0, 0, -407, 407}, {}}; ev.incoming[1] = {gluon, { 0, 0, 126, 126}, {}}; return ev; default: ev.outgoing.push_back({gluon, { 248, -56, -122, 282}, {}}); ev.outgoing.push_back({gluon, { 249, 30, -10, 251}, {}}); ev.outgoing.push_back({gluon, {-249, -18, 26, 251}, {}}); ev.outgoing.push_back({gluon, {-248, 44, 199, 321}, {}}); ev.incoming[0] = {gluon, { 0, 0, -506, 506}, {}}; ev.incoming[1] = {gluon, { 0, 0, 599, 599}, {}}; return ev; } } if(njet == 6){ switch(pos_boson){ case 0: ev.outgoing.push_back({higgs, { 349, 330, -94, 505}, {}}); ev.outgoing.push_back({gluon, {-315, -300, 0, 435}, {}}); ev.outgoing.push_back({gluon, { 347, 306, 18, 463}, {}}); ev.outgoing.push_back({gluon, {-249, -342, 162, 453}, {}}); ev.outgoing.push_back({gluon, { 345, 312, 284, 545}, {}}); ev.outgoing.push_back({gluon, {-324, -126, 292, 454}, {}}); ev.outgoing.push_back({gluon, {-153, -180, 304, 385}, {}}); ev.incoming[0] = {gluon, { 0, 0, -1137, 1137}, {}}; ev.incoming[1] = {gluon, { 0, 0, 2103, 2103}, {}}; return ev; case 1: ev.outgoing.push_back({gluon, { 242, 241, -182, 387}, {}}); ev.outgoing.push_back({higgs, { 243, 238, -190, 409}, {}}); ev.outgoing.push_back({gluon, {-218, -215, -74, 315}, {}}); ev.outgoing.push_back({gluon, {-224, -224, 112, 336}, {}}); ev.outgoing.push_back({gluon, { 241, 182, 154, 339}, {}}); ev.outgoing.push_back({gluon, { -53, -234, 126, 271}, {}}); ev.outgoing.push_back({gluon, {-231, 12, 156, 279}, {}}); ev.incoming[0] = {gluon, { 0, 0, -1117, 1117}, {}}; ev.incoming[1] = {gluon, { 0, 0, 1219, 1219}, {}}; return ev; case 2: ev.outgoing.push_back({gluon, { 151, 102, -42, 187}, {}}); ev.outgoing.push_back({gluon, { -86, -46, -17, 99}, {}}); ev.outgoing.push_back({higgs, { 152, 153, 0, 249}, {}}); ev.outgoing.push_back({gluon, { -60, -135, 64, 161}, {}}); ev.outgoing.push_back({gluon, { 150, 123, 110, 223}, {}}); ev.outgoing.push_back({gluon, {-154, -49, 98, 189}, {}}); ev.outgoing.push_back({gluon, {-153, -148, 144, 257}, {}}); ev.incoming[0] = {gluon, { 0, 0, -504, 504}, {}}; ev.incoming[1] = {gluon, { 0, 0, 861, 861}, {}}; return ev; case 3: ev.outgoing.push_back({gluon, { 198, 197, -66, 287}, {}}); ev.outgoing.push_back({gluon, {-198, -189, -54, 279}, {}}); ev.outgoing.push_back({gluon, {-200, -64, 2, 210}, {}}); ev.outgoing.push_back({higgs, { 199, 158, 6, 283}, {}}); ev.outgoing.push_back({gluon, {-199, -184, 172, 321}, {}}); ev.outgoing.push_back({gluon, { 196, 168, 177, 313}, {}}); ev.outgoing.push_back({gluon, { 4, -86, 92, 126}, {}}); ev.incoming[0] = {gluon, { 0, 0, -745, 745}, {}}; ev.incoming[1] = {gluon, { 0, 0, 1074, 1074}, {}}; return ev; case 4: ev.outgoing.push_back({gluon, { 151, 102, -42, 187}, {}}); ev.outgoing.push_back({gluon, { -86, -133, -14, 159}, {}}); ev.outgoing.push_back({gluon, {-154, -104, -8, 186}, {}}); ev.outgoing.push_back({gluon, { -60, 11, 0, 61}, {}}); ev.outgoing.push_back({higgs, { 152, 153, 0, 249}, {}}); ev.outgoing.push_back({gluon, { 150, 125, 90, 215}, {}}); ev.outgoing.push_back({gluon, {-153, -154, 126, 251}, {}}); ev.incoming[0] = {gluon, { 0, 0, -578, 578}, {}}; ev.incoming[1] = {gluon, { 0, 0, 730, 730}, {}}; return ev; case 5: ev.outgoing.push_back({gluon, { -15, -90, -94, 131}, {}}); ev.outgoing.push_back({gluon, { -11, 82, -74, 111}, {}}); ev.outgoing.push_back({gluon, { 23, -80, -64, 105}, {}}); ev.outgoing.push_back({gluon, { -48, -25, -36, 65}, {}}); ev.outgoing.push_back({gluon, { -12, 99, -16, 101}, {}}); ev.outgoing.push_back({higgs, { 68, 92, -18, 170}, {}}); ev.outgoing.push_back({gluon, { -5, -78, 54, 95}, {}}); ev.incoming[0] = {gluon, { 0, 0, -513, 513}, {}}; ev.incoming[1] = {gluon, { 0, 0, 265, 265}, {}}; return ev; case 6: ev.outgoing.push_back({gluon, { 198, 197, -66, 287}, {}}); ev.outgoing.push_back({gluon, { 4, -84, -18, 86}, {}}); ev.outgoing.push_back({gluon, {-198, -60, -36, 210}, {}}); ev.outgoing.push_back({gluon, { 196, -78, -36, 214}, {}}); ev.outgoing.push_back({gluon, {-200, 45, 0, 205}, {}}); ev.outgoing.push_back({gluon, {-199, -178, 2, 267}, {}}); ev.outgoing.push_back({higgs, { 199, 158, 6, 283}, {}}); ev.incoming[0] = {gluon, { 0, 0, -850, 850}, {}}; ev.incoming[1] = {gluon, { 0, 0, 702, 702}, {}}; return ev; default: ev.outgoing.push_back({gluon, {-350, -112, -280, 462}, {}}); ev.outgoing.push_back({gluon, { 347, 266, -322, 543}, {}}); ev.outgoing.push_back({gluon, {-349, -314, -38, 471}, {}}); ev.outgoing.push_back({gluon, { 349, 348, 12, 493}, {}}); ev.outgoing.push_back({gluon, {-342, -54, 23, 347}, {}}); ev.outgoing.push_back({gluon, { 345, -134, 138, 395}, {}}); ev.incoming[0] = {gluon, { 0, 0, -1589, 1589}, {}}; ev.incoming[1] = {gluon, { 0, 0, 1122, 1122}, {}}; return ev; } } if(njet == 7){ switch(pos_boson){ case -1: // jet idx: -1 0 1 2 3 4 5 ev.outgoing.push_back({gluon, { -15, -18, -54, 59}, {}}); ev.outgoing.push_back({gluon, { -11, 98, -70, 121}, {}}); ev.outgoing.push_back({gluon, { 23, -100, -64, 121}, {}}); ev.outgoing.push_back({gluon, { 68, 93, -20, 117}, {}}); ev.outgoing.push_back({gluon, { -5, -92, -12, 93}, {}}); ev.outgoing.push_back({gluon, { -48, -76, -2, 90}, {}}); ev.outgoing.push_back({gluon, { -12, 95, 56, 111}, {}}); ev.incoming[0] = {gluon, { 0, 0, -439, 439}, {}}; ev.incoming[1] = {gluon, { 0, 0, 273, 273}, {}}; return ev; case -2: // jet idx: 0 1 2 3 4 -1 -1 ev.outgoing.push_back({gluon, { -5, -86, -82, 119}, {}}); ev.outgoing.push_back({gluon, { 68, 93, 20, 117}, {}}); ev.outgoing.push_back({gluon, { -48, -14, 20, 54}, {}}); ev.outgoing.push_back({gluon, { 23, -50, 26, 61}, {}}); ev.outgoing.push_back({gluon, { -12, 95, 56, 111}, {}}); ev.outgoing.push_back({gluon, { -15, -18, 54, 59}, {}}); ev.outgoing.push_back({gluon, { -11, -20, 88, 91}, {}}); ev.incoming[0] = {gluon, { 0, 0, -215, 215}, {}}; ev.incoming[1] = {gluon, { 0, 0, 397, 397}, {}}; return ev; case -3: // jet idx: 0 0 1 2 2 3 4 // jet pt fraction: 0.6 0.38 1 0.49 0.51 1 1 ev.outgoing.push_back({gluon, { 23, -94, -62, 1.2e+02}, {}}); ev.outgoing.push_back({gluon, { -5, -62, -34, 71}, {}}); ev.outgoing.push_back({gluon, { 68, 93, 20, 1.2e+02}, {}}); ev.outgoing.push_back({gluon, { -12, 95, 56, 1.1e+02}, {}}); ev.outgoing.push_back({gluon, { -11, 98, 70, 1.2e+02}, {}}); ev.outgoing.push_back({gluon, { -48, -1e+02, 82, 1.4e+02}, {}}); ev.outgoing.push_back({gluon, { -15, -30, 78, 85}, {}}); ev.incoming[0] = {gluon, { 0, 0, -2.7e+02, 2.7e+02}, {}}; ev.incoming[1] = {gluon, { 0, 0, 4.8e+02, 4.8e+02}, {}}; return ev; case -4: // jet idx: 0 1 1 2 3 4 4 // jet pt fraction: 1 0.51 0.49 1 1 0.25 0.75 ev.outgoing.push_back({gluon, { -5, -88, -64, 109}, {}}); ev.outgoing.push_back({gluon, { -11, 98, -70, 121}, {}}); ev.outgoing.push_back({gluon, { -12, 95, -56, 111}, {}}); ev.outgoing.push_back({gluon, { 68, 93, 20, 117}, {}}); ev.outgoing.push_back({gluon, { 23, -70, 22, 77}, {}}); ev.outgoing.push_back({gluon, { -15, -32, 16, 39}, {}}); ev.outgoing.push_back({gluon, { -48, -96, 75, 131}, {}}); ev.incoming[0] = {gluon, { 0, 0, -381, 381}, {}}; ev.incoming[1] = {gluon, { 0, 0, 324, 324}, {}}; return ev; case -5: // jet idx: 0 1 -1 -1 2 3 4 ev.outgoing.push_back({gluon, { -15, -26, -62, 69}, {}}); ev.outgoing.push_back({gluon, { -48, -60, -54, 94}, {}}); ev.outgoing.push_back({gluon, { 23, 10, -14, 29}, {}}); ev.outgoing.push_back({gluon, { -5, -20, 0, 21}, {}}); ev.outgoing.push_back({gluon, { 68, 93, 20, 117}, {}}); ev.outgoing.push_back({gluon, { -11, -92, 40, 101}, {}}); ev.outgoing.push_back({gluon, { -12, 95, 56, 111}, {}}); ev.incoming[0] = {gluon, { 0, 0, -278, 278}, {}}; ev.incoming[1] = {gluon, { 0, 0, 264, 264}, {}}; return ev; case -6: // jet idx: 0 1 1 2 -1 2 3 // jet pt fraction: 1 0.63 0.36 0.49 1 0.51 1 ev.outgoing.push_back({gluon, { 68, 93, -20, 117}, {}}); ev.outgoing.push_back({gluon, { -48, -100, 26, 114}, {}}); ev.outgoing.push_back({gluon, { -15, -62, 26, 69}, {}}); ev.outgoing.push_back({gluon, { -12, 95, 56, 111}, {}}); ev.outgoing.push_back({gluon, { -5, -28, 20, 35}, {}}); ev.outgoing.push_back({gluon, { -11, 98, 70, 121}, {}}); ev.outgoing.push_back({gluon, { 23, -96, 92, 135}, {}}); ev.incoming[0] = {gluon, { 0, 0, -216, 216}, {}}; ev.incoming[1] = {gluon, { 0, 0, 486, 486}, {}}; return ev; case -7: // jet idx: 0 1 2 2 3 3 4 // jet pt fraction: 1 1 0.51 0.49 0.18 0.82 1 ev.outgoing.push_back({gluon, { -15, -80, -100, 129}, {}}); ev.outgoing.push_back({gluon, { 23, -96, -92, 135}, {}}); ev.outgoing.push_back({gluon, { -11, 98, -70, 121}, {}}); ev.outgoing.push_back({gluon, { -12, 95, -56, 111}, {}}); ev.outgoing.push_back({gluon, { -5, -22, -10, 25}, {}}); ev.outgoing.push_back({gluon, { -48, -88, -31, 105}, {}}); ev.outgoing.push_back({gluon, { 68, 93, 20, 117}, {}}); ev.incoming[0] = {gluon, { 0, 0, -541, 541}, {}}; ev.incoming[1] = {gluon, { 0, 0, 202, 202}, {}}; return ev; case -8: // jet idx: 0 1 2 2 2 3 4 // jet pt fraction: 1 1 0.21 0.37 0.41 1 1 ev.outgoing.push_back({gluon, { -48, -44, -62, 90}, {}}); ev.outgoing.push_back({gluon, { -12, 95, -56, 111}, {}}); ev.outgoing.push_back({gluon, { -5, -50, -22, 55}, {}}); ev.outgoing.push_back({gluon, { 23, -90, -34, 99}, {}}); ev.outgoing.push_back({gluon, { -15, -100, -28, 105}, {}}); ev.outgoing.push_back({gluon, { 68, 93, -20, 117}, {}}); ev.outgoing.push_back({gluon, { -11, 96, 76, 123}, {}}); ev.incoming[0] = {gluon, { 0, 0, -423, 423}, {}}; ev.incoming[1] = {gluon, { 0, 0, 277, 277}, {}}; return ev; case -9: // jet idx: 0 1 2 1 3 0 4 // jet pt fraction: 0.72 0.51 1 0.49 1 0.28 1 ev.outgoing.push_back({gluon, { -15, -98, -62, 117}, {}}); ev.outgoing.push_back({gluon, { -12, 95, -56, 111}, {}}); ev.outgoing.push_back({gluon, { 23, -76, -40, 89}, {}}); ev.outgoing.push_back({gluon, { -11, 92, -40, 101}, {}}); ev.outgoing.push_back({gluon, { -48, -68, -34, 90}, {}}); ev.outgoing.push_back({gluon, { -5, -38, -14, 41}, {}}); ev.outgoing.push_back({gluon, { 68, 93, 20, 117}, {}}); ev.incoming[0] = {gluon, { 0, 0, -446, 446}, {}}; ev.incoming[1] = {gluon, { 0, 0, 220, 220}, {}}; return ev; case -10: // jet idx: 0 1 3 2 4 3 1 // jet pt fraction: 1 0.33 0.51 1 1 0.49 0.67 ev.outgoing.push_back({gluon, { 68, 93, 20, 117}, {}}); ev.outgoing.push_back({gluon, { -5, -48, 16, 51}, {}}); ev.outgoing.push_back({gluon, { -12, 95, 56, 111}, {}}); ev.outgoing.push_back({gluon, { 23, -76, 52, 95}, {}}); ev.outgoing.push_back({gluon, { -48, -60, 54, 94}, {}}); ev.outgoing.push_back({gluon, { -11, 92, 68, 115}, {}}); ev.outgoing.push_back({gluon, { -15, -96, 72, 121}, {}}); ev.incoming[0] = {gluon, { 0, 0, -183, 183}, {}}; ev.incoming[1] = {gluon, { 0, 0, 521, 521}, {}}; return ev; case -11: // jet idx: 0 1 2 3 4 -1 5 // jet pt fraction: 1 1 1 1 1 1 1 ev.outgoing.push_back({gluon, { -11, 98, -70, 121}, {}}); ev.outgoing.push_back({gluon, { -15, -98, -62, 117}, {}}); ev.outgoing.push_back({gluon, { 23, -90, -2, 93}, {}}); ev.outgoing.push_back({gluon, { -48, -76, 2, 90}, {}}); ev.outgoing.push_back({gluon, { 68, 93, 20, 117}, {}}); ev.outgoing.push_back({gluon, { -5, -22, 10, 25}, {}}); ev.outgoing.push_back({gluon, { -12, 95, 56, 111}, {}}); ev.incoming[0] = {gluon, { 0, 0, -360, 360}, {}}; ev.incoming[1] = {gluon, { 0, 0, 314, 314}, {}}; return ev; case -12: // jet idx: 0 1 -1 2 3 4 3 // jet pt fraction: 1 1 1 1 0.35 1 0.65 ev.outgoing.push_back({gluon, { 23, -94, -62, 115}, {}}); ev.outgoing.push_back({gluon, { -12, 95, -56, 111}, {}}); ev.outgoing.push_back({gluon, { -5, -28, 4, 29}, {}}); ev.outgoing.push_back({gluon, { 68, 93, 20, 117}, {}}); ev.outgoing.push_back({gluon, { -15, -58, 34, 69}, {}}); ev.outgoing.push_back({gluon, { -11, 92, 68, 115}, {}}); ev.outgoing.push_back({gluon, { -48, -100, 82, 138}, {}}); ev.incoming[0] = {gluon, { 0, 0, -302, 302}, {}}; ev.incoming[1] = {gluon, { 0, 0, 392, 392}, {}}; return ev; case -13: // jet idx: 0 1 2 3 3 4 2 // jet pt fraction: 1 1 0.5 0.35 0.65 1 0.5 ev.outgoing.push_back({gluon, { -15, -98, -62, 117}, {}}); ev.outgoing.push_back({gluon, { 68, 93, 20, 117}, {}}); ev.outgoing.push_back({gluon, { -12, 95, 56, 111}, {}}); ev.outgoing.push_back({gluon, { -5, -28, 20, 35}, {}}); ev.outgoing.push_back({gluon, { -48, -96, 75, 131}, {}}); ev.outgoing.push_back({gluon, { 23, -62, 50, 83}, {}}); ev.outgoing.push_back({gluon, { -11, 96, 76, 123}, {}}); ev.incoming[0] = {gluon, { 0, 0, -241, 241}, {}}; ev.incoming[1] = {gluon, { 0, 0, 476, 476}, {}}; return ev; case -14: // jet idx: 0 1 2 3 3 4 2 // jet pt fraction: 1 1 0.52 0.35 0.65 1 0.48 ev.outgoing.push_back({gluon, { 23, -94, -62, 115}, {}}); ev.outgoing.push_back({gluon, { 68, 93, 20, 117}, {}}); ev.outgoing.push_back({gluon, { -12, 95, 56, 111}, {}}); ev.outgoing.push_back({gluon, { -15, -96, 72, 121}, {}}); ev.outgoing.push_back({gluon, { -5, -42, 38, 57}, {}}); ev.outgoing.push_back({gluon, { -48, -44, 62, 90}, {}}); ev.outgoing.push_back({gluon, { -11, 88, 88, 125}, {}}); ev.incoming[0] = {gluon, { 0, 0, -231, 231}, {}}; ev.incoming[1] = {gluon, { 0, 0, 505, 505}, {}}; return ev; case -15: // jet idx: 0 -1 1 0 2 3 4 // jet pt fraction: 0.51 1 1 0.49 1 1 1 ev.outgoing.push_back({gluon, { -11, 98, -70, 121}, {}}); ev.outgoing.push_back({gluon, { -5, -16, -12, 21}, {}}); ev.outgoing.push_back({gluon, { 23, -94, -62, 115}, {}}); ev.outgoing.push_back({gluon, { -12, 95, -56, 111}, {}}); ev.outgoing.push_back({gluon, { 68, 93, 20, 117}, {}}); ev.outgoing.push_back({gluon, { -48, -76, 70, 114}, {}}); ev.outgoing.push_back({gluon, { -15, -100, 80, 129}, {}}); ev.incoming[0] = {gluon, { 0, 0, -379, 379}, {}}; ev.incoming[1] = {gluon, { 0, 0, 349, 349}, {}}; return ev; } } throw HEJ::unknown_option{"unknown process"}; } HEJ::Event::EventData parse_configuration( std::array const & in, std::vector const & out, int const overwrite_boson ){ auto boson = std::find_if(out.cbegin(), out.cend(), [](std::string id){ return !HEJ::is_parton(HEJ::to_ParticleID(id)); }); int const pos_boson = (overwrite_boson!=0)?overwrite_boson: ((boson == out.cend())?-1:std::distance(out.cbegin(), boson)); std::size_t njets = out.size(); if( (overwrite_boson == 0) && boson != out.cend()) --njets; HEJ::Event::EventData ev{get_process(njets, pos_boson)}; ASSERT((pos_boson<0) || (ev.outgoing[pos_boson].type == HEJ::ParticleID::higgs)); for(std::size_t i=0; i const & in, std::vector const & out, bool reconstruct /* = false */, std::unordered_map< size_t, std::vector > decays /* = {} */ ) { using namespace HEJ::pid; HEJ::Event::EventData evd; const size_t n_out {out.size()}; // Parameters double pT {40.}; std::array cs_phi {1, 0}; std::array cs_dphi {cos(2.*M_PI/n_out), sin(2.*M_PI/n_out)}; double y {-4.5}, dy {std::max(-2.*y/(n_out-1), 0.5)}; double sum_pE {0}, sum_pz {0}; // Outgoing evd.outgoing.reserve(n_out); for(size_t i = 0; i < n_out; ++i) { auto pid = HEJ::to_ParticleID(out[i]); // Kinematics double mT2 = pT*pT; if(HEJ::is_AWZH_boson(pid)) { auto const ref = ew_parameters.prop(pid); mT2 += ref.mass*ref.mass; } double pz {sqrt(mT2)*sinh(y)}, pE {sqrt(mT2)*cosh(y)}; evd.outgoing.push_back({pid, {pT*cs_phi[0], pT*cs_phi[1], pz, pE}, {}}); sum_pE+=pE; sum_pz+=pz; // Update cos,sin, y cs_phi = {cs_phi[0] * cs_dphi[0] - cs_phi[1] * cs_dphi[1], cs_phi[1] * cs_dphi[0] + cs_phi[0] * cs_dphi[1]}; y += dy; } // Specified decays evd.decays.reserve(decays.size()); for(auto const & decay : decays) { auto const parent = evd.outgoing.at(decay.first); std::vector progeny = decay_kinematics(parent); for(std::size_t i = 0; i < decay.second.size(); ++i) { progeny[i].type = HEJ::to_ParticleID(decay.second[i]); } evd.decays.emplace(decay.first, std::move(progeny)); } // Reconstruct decays if(reconstruct) { evd.reconstruct_intermediate(ew_parameters); } // Incoming double p0 {0.5*(sum_pE + sum_pz)}; evd.incoming[0] = {HEJ::to_ParticleID(in[0]),{0,0,-p0,p0},{}}; double p1 {0.5*(sum_pE - sum_pz)}; evd.incoming[1] = {HEJ::to_ParticleID(in[1]),{0,0,+p1,p1},{}}; return evd; } namespace { static std::mt19937_64 ran{0}; } void shuffle_particles(HEJ::Event::EventData & ev) { // incoming std::shuffle(ev.incoming.begin(), ev.incoming.end(), ran); // outgoing (through index) auto old_outgoing = std::move(ev).outgoing; std::vector idx(old_outgoing.size()); std::iota(idx.begin(), idx.end(), 0); std::shuffle(idx.begin(), idx.end(), ran); ev.outgoing.clear(); ev.outgoing.reserve(old_outgoing.size()); for(std::size_t i: idx) { ev.outgoing.emplace_back(std::move(old_outgoing[i])); } // find decays again if(!ev.decays.empty()){ auto old_decays = std::move(ev).decays; ev.decays.clear(); for(std::size_t i=0; isecond)); } for(auto & decay: ev.decays){ std::shuffle(decay.second.begin(), decay.second.end(), ran); } } } std::vector decay_kinematics( HEJ::Particle const & parent ) { std::vector decay_products(2); const double E = parent.m()/2; const double theta = 2.*M_PI*ran()/static_cast(ran.max()); const double cos_phi = 2.*ran()/static_cast(ran.max())-1.; const double sin_phi = std::sqrt(1. - cos_phi*cos_phi); // Know 0 < phi < pi const double px = E*std::cos(theta)*sin_phi; const double py = E*std::sin(theta)*sin_phi; const double pz = E*cos_phi; decay_products[0].p.reset(px, py, pz, E); decay_products[1].p.reset(-px, -py, -pz, E); for(auto & particle: decay_products) particle.p.boost(parent.p); return decay_products; } std::vector decay_W( HEJ::Particle const & parent ){ if(parent.m() == 0.) // we can't decay massless partons return {}; std::array decays; if(parent.type==HEJ::ParticleID::Wp){ // order matters: first particle, second anti decays[0] = HEJ::ParticleID::nu_e; decays[1] = HEJ::ParticleID::e_bar; } else { // this function is for testing: we don't check that parent==W boson decays[0] = HEJ::ParticleID::e; decays[1] = HEJ::ParticleID::nu_e_bar; } std::vector decay_products = decay_kinematics(parent); for(std::size_t i = 0; i < decay_products.size(); ++i){ decay_products[i].type = decays[i]; } return decay_products; } std::vector decay_Z( HEJ::Particle const & parent ){ if(parent.m() == 0.) // we can't decay massless partons return {}; std::array decays; // order matters: first particle, second anti decays[0] = HEJ::ParticleID::electron; decays[1] = HEJ::ParticleID::positron; std::vector decay_products = decay_kinematics(parent); for(std::size_t i = 0; i < decay_products.size(); ++i){ decay_products[i].type = decays[i]; } return decay_products; } diff --git a/t/test_classify_ref.cc b/t/test_classify_ref.cc index bdb6bf4..d5256a3 100644 --- a/t/test_classify_ref.cc +++ b/t/test_classify_ref.cc @@ -1,85 +1,86 @@ /** * \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 "HEJ/event_types.hh" #include "HEJ/Event.hh" #include "HEJ/EventReader.hh" +#include "HEJ/EWConstants.hh" namespace { // this is deliberately chosen bigger than in the generation, // to cluster multiple partons in one jet constexpr double min_jet_pt = 40.; const fastjet::JetDefinition jet_def{fastjet::kt_algorithm, 0.6}; const HEJ::ParticleProperties Wprop{80.385, 2.085}; const HEJ::ParticleProperties Zprop{91.187, 2.495}; const HEJ::ParticleProperties Hprop{125, 0.004165}; constexpr double vev = 246.2196508; const HEJ::EWConstants ew_parameters{vev, Wprop, Zprop, Hprop}; } int main(int argn, char** argv) { if(argn != 3 && argn != 4){ std::cerr << "Usage: " << argv[0] << " reference_classification input_file.lhe\n"; return EXIT_FAILURE; } bool OUTPUT_MODE = false; if(argn == 4 && std::string("OUTPUT")==std::string(argv[3])) OUTPUT_MODE = true; std::fstream ref_file; if ( OUTPUT_MODE ) { std::cout << "_______________________USING OUTPUT MODE!_______________________" << std::endl; ref_file.open(argv[1], std::fstream::out); } else { ref_file.open(argv[1], std::fstream::in); } auto reader{ HEJ::make_reader(argv[2]) }; std::size_t nevent{0}; while(reader->read_event()){ ++nevent; // We don't need to test forever, the first "few" are enough if(nevent>4000) break; HEJ::Event::EventData data{ reader->hepeup() }; shuffle_particles(data); data.reconstruct_intermediate(ew_parameters); const HEJ::Event ev{ data.cluster( jet_def, min_jet_pt ) }; if ( OUTPUT_MODE ) { ref_file << ev.type() << std::endl; } else { std::string line; if(!std::getline(ref_file,line)) break; const auto expected{static_cast(std::stoi(line))}; if(ev.type() != expected){ std::cerr << "wrong classification of event " << nevent << ":\n" << ev << "classified as " << name(ev.type()) << ", expected " << name(expected) << "\nJet indices: "; for(auto const & idx: ev.particle_jet_indices()) std::cerr << idx << " "; std::cerr << "\n"; return EXIT_FAILURE; } } } return EXIT_SUCCESS; }