diff --git a/include/HEJ/Event.hh b/include/HEJ/Event.hh index c9cd6dd..dae110a 100644 --- a/include/HEJ/Event.hh +++ b/include/HEJ/Event.hh @@ -1,403 +1,413 @@ /** \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/Constants.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 { class 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 got superseded by EventData //! and will be removed in HEJ 2.2.0 [[deprecated("UnclusteredEvent got superseded 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 * `soft_pt_regulator` 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 soft_pt_regulator Maximum transverse momentum fraction from soft * radiation in tagging jets * @param min_pt Absolute minimal \f$ p_\perp \f$, * \b deprecated use soft_pt_regulator instead * @return True if this state could have been produced by HEJ */ bool valid_hej_state( double soft_pt_regulator = DEFAULT_SOFT_PT_REGULATOR, double min_pt = 0. ) const; //! Check that the incoming momenta are valid /** * @details Checks that the incoming parton momenta are on-shell and have * vanishing transverse components. * */ bool valid_incoming() 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 //! Detect if a backward incoming gluon turns into a backward outgoing Higgs boson inline bool is_backward_g_to_h(Event const & ev) { return ev.outgoing().front().type == pid::higgs && ev.incoming().front().type == pid::gluon; } //! Detect if a forward incoming gluon turns into a forward outgoing Higgs boson inline bool is_forward_g_to_h(Event const & ev) { return ev.outgoing().back().type == pid::higgs && ev.incoming().back().type == pid::gluon; } //! 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*/); + //! Repair momenta of massless particles + /** + * This function changes the momenta of massless particles as follows: + * - Close-to-zero incoming transverse momenta are set to zero. + * - Nearly on-shell momenta are made lightlike by rescaling energy + * and spatial components. This rescaling is applied to both incoming + * and outgoing particles, including decay products. + */ + void repair_momenta(double tolerance); + //! 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); //! Tolerance parameter for validity check on incoming momenta static constexpr double TOL = 1e-6; //! 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 got superseded by EventData //! and will be removed 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/include/HEJ/PDG_codes.hh b/include/HEJ/PDG_codes.hh index 23694de..00f8af4 100644 --- a/include/HEJ/PDG_codes.hh +++ b/include/HEJ/PDG_codes.hh @@ -1,241 +1,261 @@ /** \file PDG_codes.hh * \brief Contains the Particle IDs of all relevant SM particles. * * Large enumeration included which has multiple entries for potential * alternative names of different particles. There are also functions * which can be used to determine if a particle is a parton or if * it is a non-gluon boson. * * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #pragma once #include #include namespace HEJ { //! particle ids according to PDG namespace pid { //! The possible particle identities. We use PDG IDs as standard. enum ParticleID: int{ //! Unspecified type, should never be used!, debug only unspecified = 0, d = 1, /*!< Down Quark */ down = d, /*!< Down Quark */ u = 2, /*!< Up Quark */ up = u, /*!< Up Quark */ s = 3, /*!< Strange Quark */ strange = s, /*!< Strange Quark */ c = 4, /*!< Charm Quark */ charm = c, /*!< Charm Quark */ b = 5, /*!< Bottom Quark */ bottom = b, /*!< Bottom Quark */ t = 6, /*!< Top Quark */ top = t, /*!< Top Quark */ e = 11, /*!< Electron */ electron = e, /*!< Electron */ nu_e = 12, /*!< Electron Neutrino */ electron_neutrino = nu_e, /*!< Electron neutrino */ mu = 13, /*!< Muon */ muon = mu, /*!< Muon */ nu_mu = 14, /*!< Muon Neutrino */ muon_neutrino = nu_mu, /*!< Muon Neutrino */ tau = 15, /*!< Tau */ nu_tau = 16, /*!< Tau Neutrino */ tau_neutrino = nu_tau, /*!< Tau Neutrino */ d_bar = -d, /*!< Anti-Down Quark */ antidown = d_bar, /*!< Anti-Down Quark */ u_bar = -u, /*!< Anti-Up quark */ antiup = -u, /*!< Anti-Up quark */ s_bar = -s, /*!< Anti-Strange Quark */ antistrange = -s, /*!< Anti-Strange Quark */ c_bar = -c, /*!< Anti-Charm Quark */ anticharm = -c, /*!< Anti-Charm Quark */ b_bar = -b, /*!< Anti-Bottom Quark */ antibottom = -b, /*!< Anti-Bottom Quark */ t_bar = -t, /*!< Anti-Top Quark */ antitop = -t, /*!< Anti-Top Quark */ e_bar = -e, /*!< Positron */ positron = e_bar, /*!< Positron */ antielectron = positron, /*!< Positron */ nu_e_bar = -nu_e, /*!< Electron Anti-Neutrino */ electron_antineutrino = nu_e_bar,/*!< Electron Anti-Neutrino */ mu_bar = -mu, /*!< Anti-Muon */ antimuon = -mu, /*!< Anti-Muon */ nu_mu_bar = -nu_mu, /*!< Muon Anti-Neutrino */ muon_antineutrino = nu_mu_bar, /*!< Muon Anti-Neutrino */ tau_bar = -tau, /*!< Anti-Tau */ antitau = tau_bar, /*!< Anti-Tau */ nu_tau_bar = -nu_tau, /*!< Tau Anti-Neutrino */ tau_antineutrino = nu_tau_bar, /*!< Tau Anti-Neutrino */ gluon = 21, /*!< Gluon */ g = gluon, /*!< Gluon */ photon = 22, /*!< Photon */ gamma = photon, /*!< Photon */ Z = 23, /*!< Z Boson */ Z_photon_mix = 81, /*!< Z/photon superposition */ Z_gamma_mix = Z_photon_mix, /*!< Z/photon superposition */ Wp = 24, /*!< W- Boson */ Wm = -Wp, /*!< W+ Boson */ h = 25, /*!< Higgs Boson */ Higgs = h, /*!< Higgs Boson */ higgs = h, /*!< Higgs Boson */ p = 2212, /*!< Proton */ proton = p, /*!< Proton */ p_bar = -p, /*!< Anti-Proton */ antiproton = p_bar, /*!< Anti-Proton */ }; //! Get the of the particle with the given PDG ID std::string name(ParticleID id); //! return the negative flavour of the given PDG ID ParticleID anti(ParticleID id); } // namespace pid using ParticleID = pid::ParticleID; //! Convert a particle name to the corresponding PDG particle ID ParticleID to_ParticleID(std::string const & name); /** * \brief Function to determine if particle is a quark * @param id PDG ID of particle * @returns true if the particle is a quark, false otherwise */ inline constexpr bool is_quark(ParticleID id){ return (id >= pid::down && id <= pid::top); } /** * \brief Function to determine if particle is an antiquark * @param id PDG ID of particle * @returns true if the particle is an antiquark, false otherwise */ inline constexpr bool is_antiquark(ParticleID id){ return (id <= pid::d_bar && id >= pid::t_bar); } /** * \brief Function to determine if particle is an (anti-)quark * @param id PDG ID of particle * @returns true if the particle is a quark or antiquark, false otherwise */ inline constexpr bool is_anyquark(ParticleID id){ return is_quark(id) || is_antiquark(id); } /** * \brief Function to determine if particle is a gluon * @param id PDG ID of particle * @returns true if the particle is a gluon, false otherwise */ inline constexpr bool is_gluon(ParticleID id){ return id == pid::gluon; } /** * \brief Function to determine if particle is a parton * @param id PDG ID of particle * @returns true if the particle is a parton, false otherwise */ inline constexpr bool is_parton(ParticleID id){ return is_gluon(id) || (is_anyquark(id) && std::abs(id) != pid::top); } /** * \brief function to determine if the particle is a photon, W or Z * @param id PDG ID of particle * @returns true if the partice is an A,W,Z, or H, false otherwise */ inline constexpr bool is_AWZ_boson(ParticleID id){ return id == pid::Wm || (id >= pid::photon && id <= pid::Wp) || id == pid::Z_photon_mix; } /** * \brief function to determine if the particle is a photon, W, Z, or Higgs * boson * @param id PDG ID of particle * @returns true if the partice is an A,W,Z, or H, false otherwise */ inline constexpr bool is_AWZH_boson(ParticleID id){ return is_AWZ_boson(id) || (id == pid::Higgs); } /** * \brief Function to determine if particle is a lepton * @param id PDG ID of particle * @returns true if the particle is a lepton, false otherwise */ inline constexpr bool is_lepton(ParticleID id){ return (id >= pid::electron && id <= pid::tau_neutrino); } /** * \brief Function to determine if particle is an antilepton * @param id PDG ID of particle * @returns true if the particle is an antilepton, false otherwise */ inline constexpr bool is_antilepton(ParticleID id){ return (id <= pid::positron && id >= pid::nu_tau_bar); } /** * \brief Function to determine if particle is an (anti-)lepton * @param id PDG ID of particle * @returns true if the particle is a lepton or antilepton, * false otherwise */ inline constexpr bool is_anylepton(ParticleID id){ return ( is_lepton(id) || is_antilepton(id)); } /** * \brief Function to determine if particle is a neutrino * @param id PDG ID of particle * @returns true if the particle is a neutrino, false otherwise */ inline constexpr bool is_neutrino(ParticleID id){ return (id == pid::nu_e || id == pid::tau_neutrino || id == pid::muon_neutrino); } /** * \brief Function to determine if particle is an antineutrino * @param id PDG ID of particle * @returns true if the particle is an antineutrino, false otherwise */ inline constexpr bool is_antineutrino(ParticleID id){ return (id == pid::nu_e_bar || id == pid::nu_tau_bar || id == pid::nu_mu_bar); } /** * \brief Function to determine if particle is an (anti-)neutrino * @param id PDG ID of particle * @returns true if the particle is a neutrino or antineutrino, * false otherwise */ inline constexpr bool is_anyneutrino(ParticleID id){ return ( is_neutrino(id) || is_antineutrino(id)); } + + //! Check if a particle is massless + inline + constexpr bool is_massless(ParticleID id){ + // cannot use `std::abs` because it's not `constexpr` + const int abs_id = (id >= 0)?id:-id; + switch(abs_id){ + case pid::bottom: + case pid::top: + case pid::tau: + case pid::Z: + case pid::Z_photon_mix: + case pid::Wp: + case pid::Higgs: + case pid::proton: + return false; + default: + return true; + } + } } // namespace HEJ diff --git a/include/HEJ/Particle.hh b/include/HEJ/Particle.hh index 94f0b4e..48e905f 100644 --- a/include/HEJ/Particle.hh +++ b/include/HEJ/Particle.hh @@ -1,239 +1,245 @@ /** * \file Particle.hh * \brief Contains the particle struct * * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2022 * \copyright GPLv2 or later */ #pragma once #include #include #include #include "fastjet/PseudoJet.hh" #include "HEJ/PDG_codes.hh" namespace HEJ { using Colour = std::pair; //! Class representing a particle struct Particle { //! particle type ParticleID type = pid::unspecified; //! particle momentum fastjet::PseudoJet p; //! (optional) colour & anti-colour std::optional colour; //! get rapidity double rapidity() const{ return p.rapidity(); } //! get transverse momentum double perp() const{ return p.perp(); } //! get transverse momentum double pt() const{ return perp(); } //! get momentum in x direction double px() const{ return p.px(); } //! get momentum in y direction double py() const{ return p.py(); } //! get momentum in z direction double pz() const{ return p.pz(); } //! get energy double E() const{ return p.E(); } //! get mass double m() const{ return p.m(); } }; //! Functor to compare rapidities /** * This can be used whenever a rapidity comparison function is needed, * for example in many standard library functions. * * @see pz_less */ struct rapidity_less{ template bool operator()(FourVector const & p1, FourVector const & p2){ return p1.rapidity() < p2.rapidity(); } }; //! Functor to compare momenta in z direction /** * This can be used whenever a pz comparison function is needed, * for example in many standard library functions. * * @see rapidity_less */ struct pz_less{ template bool operator()(FourVector const & p1, FourVector const & p2){ return p1.pz() < p2.pz(); } }; //! Convert a vector of Particles to a vector of particle momenta inline std::vector to_PseudoJet( std::vector const & v ){ std::vector result; result.reserve(v.size()); for(auto && sp: v) result.emplace_back(sp.p); return result; } //! Functor to compare particle type (PDG) /** * This can be used whenever a particle-type comparison is needed * for example in many standard library functions. */ struct type_less{ template bool operator()(Particle const & p1, Particle const & p2){ return p1.type < p2.type; } }; //! Check if a particle is a parton, i.e. quark, antiquark, or gluon inline constexpr bool is_parton(Particle const & p){ return is_parton(p.type); } //! Check if a particle is a quark inline constexpr bool is_quark(Particle const & p){ return is_quark(p.type); } //! Check if a particle is an anti-quark inline constexpr bool is_antiquark(Particle const & p){ return is_antiquark(p.type); } //! Check if a particle is a quark or anit-quark inline constexpr bool is_anyquark(Particle const & p){ return is_anyquark(p.type); } /** * \brief Function to determine if particle is a lepton * @param p the particle * @returns true if the particle is a lepton, false otherwise */ inline constexpr bool is_lepton(Particle const & p){ return is_lepton(p.type); } /** * \brief Function to determine if particle is an antilepton * @param p the particle * @returns true if the particle is an antilepton, false otherwise */ inline constexpr bool is_antilepton(Particle const & p){ return is_antilepton(p.type); } /** * \brief Function to determine if particle is an (anti-)lepton * @param p the particle * @returns true if the particle is a lepton or antilepton, false otherwise */ inline constexpr bool is_anylepton(Particle const & p){ return is_anylepton(p.type); } /** * \brief Function to determine if particle is a neutrino * @param p the particle * @returns true if the particle is a neutrino, false otherwise */ inline constexpr bool is_neutrino(Particle const & p){ return is_neutrino(p.type); } /** * \brief Function to determine if particle is an antineutrino * @param p the particle * @returns true if the particle is an antineutrino, false otherwise */ inline constexpr bool is_antineutrino(Particle const & p){ return is_antineutrino(p.type); } /** * \brief Function to determine if particle is an (anti-)neutrino * @param p the particle * @returns true if the particle is a neutrino or antineutrino, false otherwise */ inline constexpr bool is_anyneutrino(Particle const & p){ return is_anyneutrino(p.type); } + //! Check if a particle is massless + inline + constexpr bool is_massless(Particle const & p){ + return is_massless(p.type); + } + //! Check if a particle is a photon, W or Z boson inline constexpr bool is_AWZ_boson(Particle const & particle){ return is_AWZ_boson(particle.type); } //! Check if a particle is a photon, W, Z, or Higgs boson inline constexpr bool is_AWZH_boson(Particle const & particle){ return is_AWZH_boson(particle.type); } //! Extract all partons from a vector of particles inline std::vector filter_partons( std::vector const & v ){ std::vector result; result.reserve(v.size()); std::copy_if( begin(v), end(v), std::back_inserter(result), [](Particle const & p){ return is_parton(p); } ); return result; } //! Extract all AWZH bosons from a vector of particles inline std::vector filter_AWZH_bosons( std::vector const & v ){ std::vector result; std::copy_if( begin(v), end(v), std::back_inserter(result), [](Particle const & p){ return is_AWZH_boson(p); } ); return result; } } // namespace HEJ diff --git a/src/Event.cc b/src/Event.cc index cf62654..7e935ac 100644 --- a/src/Event.cc +++ b/src/Event.cc @@ -1,1361 +1,1394 @@ /** * \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 #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/LorentzVector.hh" #include "HEJ/utility.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; i second; fastjet::PseudoJet res = std::accumulate( progeny.cbegin(), progeny.cend(), fastjet::PseudoJet(), [](fastjet::PseudoJet & sum, Particle const & p) -> fastjet::PseudoJet { return std::move(sum) + p.p; } ); if(!nearby(out.p, res, out.E())){ return false; } } // W decays (required) if(std::abs(out.type) == ParticleID::Wp){ if( decay != ev.decays().cend() && valid_W_decay(out.type>0?+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? return invalid_decays == 0; } /** * returns all EventTypes implemented in HEJ */ size_t implemented_types(std::vector const & bosons){ using namespace event_type; // no bosons if(bosons.empty()) return FKL | unob | unof | qqbar_exb | qqbar_exf | qqbar_mid; // 1 boson if(bosons.size()== 1) { switch (bosons[0].type) { case ParticleID::Wp: case ParticleID::Wm: return FKL | unob | unof | qqbar_exb | qqbar_exf | qqbar_mid; case ParticleID::Z_photon_mix: return FKL | unob | unof; case ParticleID::h: return FKL | unob | unof; default: return non_resummable; } } // 2 bosons if(bosons.size() == 2) { // Resum only samesign W events if(bosons[0].type == ParticleID::Wp && bosons[1].type == ParticleID::Wp) { return FKL; } else if(bosons[0].type == ParticleID::Wm && bosons[1].type == ParticleID::Wm) { 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 is_qqbar Current both incoming/both outgoing? * * \see is_Wm_Change */ bool is_Wp_Change(ParticleID in, ParticleID out, bool is_qqbar){ using namespace pid; if(!is_qqbar && (in==d_bar || in==u || in==s_bar || in==c)) return out == (in-1); if( is_qqbar && (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 is_qqbar Current both incoming/both outgoing? * * Ensures that change type of quark line is possible by a flavour changing * Wm emission. Allows checking of is_qqbar currents also. */ bool is_Wm_Change(ParticleID in, ParticleID out, bool is_qqbar){ using namespace pid; if(!is_qqbar && (in==d || in==u_bar || in==s || in==c_bar)) return out == (in+1); if( is_qqbar && (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 is_qqbar Current both incoming/outgoing? */ bool no_flavour_change(ParticleID in, ParticleID out, bool is_qqbar){ const int qqbarCurrent = is_qqbar?-1:1; if(std::abs(in)<=pid::top || in==pid::gluon) return (in==out*qqbarCurrent); return false; } bool has_enough_jets(Event const & event){ if(event.jets().size() >= 2) return true; if(event.jets().empty()) return false; // check for h+jet const auto the_higgs = std::find_if( begin(event.outgoing()), end(event.outgoing()), [](const auto & particle) { return particle.type == pid::higgs; } ); return the_higgs != end(event.outgoing()); } 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 is_qqbar 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 is_qqbar, int & W_change ){ if( no_flavour_change(in, out, is_qqbar) || is_gluon_to_Higgs(in, out)) { return true; } if( is_Wp_Change(in, out, is_qqbar) ) { W_change+=1; return true; } if( is_Wm_Change(in, out, is_qqbar) ) { 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; if(begin_out == end_out) return non_resummable; // keep track of all states that we don't test size_t not_tested = qqbar_mid; if(backward) not_tested |= unof | qqbar_exf; else not_tested |= unob | qqbar_exb; // 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 const auto next = std::next(begin_out); if( // first ensure that the next particle is not part of the *other* impact factor next != end_out && is_extremal_higgs_off_quark(incoming_id, begin_out->type, next->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 qqbar 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?qqbar_exb:qqbar_exf); } } } 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 | qqbar_exb | qqbar_exf; // 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( next != end_out && is_valid_impact_factor(begin_out->type, next->type, true, W_change) ){ // veto Higgs inside qqbar 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 FKL (gluon or higgs) if(std::any_of( begin_out, end_out, [](Particle const & p) { return is_anyquark(p); } )) { return non_resummable; } return possible | qqbar_mid; } 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_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 = 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, rbegin_out.base()) >= 2); assert(std::is_sorted(begin_out, rbegin_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 = ~(not_enough_jets | bad_final_state); // check forward impact factor int W_change = 0; final_type &= possible_impact_factors( in.front().type, begin_out, rbegin_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, rbegin_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, rbegin_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]:std::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; } using GroupedParticles = std::vector >; using Decay = std::pair >; using Decays = std::vector; // 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)); if(leptons.size() != 2 || neutrinos.size() != 2) { return {}; } assert(leptons.size() == 2 && neutrinos.size() == 2); std::vector valid_groupings; GroupedParticles candidate_grouping{{leptons[0], neutrinos[0]}, {leptons[1], neutrinos[1]}}; if(valid_decay(candidate_grouping.front()) && valid_decay(candidate_grouping.back())) { valid_groupings.emplace_back(std::move(candidate_grouping)); } candidate_grouping = {{leptons[1], neutrinos[0]}, {leptons[0], neutrinos[1]}}; if(valid_decay(candidate_grouping.front()) && valid_decay(candidate_grouping.back())) { valid_groupings.emplace_back(std::move(candidate_grouping)); } return valid_groupings; } // 'best' decay ordering measure double decay_measure(const 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"}; } if (valid_decays.size() == 1) { return valid_decays[0]; } // 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);} ); std::vector leptons(begin_leptons, end(outgoing)); outgoing.erase(begin_leptons, end(outgoing)); if(leptons.empty()) { return; } // nothing to do 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) { 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" }; } } + namespace { + void repair_momentum(fastjet::PseudoJet & p, const double tolerance) { + if(p.e() > 0. && p.m2() != 0. && (p.m2() < tolerance * tolerance)) { + const double rescale = std::sqrt(p.modp() / p.e()); + const double e = p.e() * rescale; + const double px = p.px() / rescale; + const double py = p.py() / rescale; + const double pz = p.pz() / rescale; + p.reset(px, py, pz, e); + } + } + } + + void Event::EventData::repair_momenta(const double tolerance) { + for(auto & in: incoming) { + if(is_massless(in)) { + const double px = (std::abs(in.px()) < tolerance)?0.:in.px(); + const double py = (std::abs(in.py()) < tolerance)?0.:in.py(); + in.p.reset(px, py, in.p.pz(), in.p.e()); + repair_momentum(in.p, tolerance); + } + } + for(auto & out: outgoing) { + if(is_massless(out)) repair_momentum(out.p, tolerance); + } + for(auto & decay: decays) { + for(auto & out: decay.second) { + if(is_massless(out)) repair_momentum(out.p, tolerance); + } + } + } + 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::qqbar_exb: { 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::qqbar_exf: { 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::qqbar_mid:{ auto it_next = std::next(it_part); if( !try_connect_t(it_part, line_colour) // u-channel only allowed at q-qbar/qbar-q 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 qbar 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 qbar 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::qqbar_exb: { 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::qqbar_exf: { 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::qqbar_mid:{ 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(!is_backward_g_to_h(*this)) { if(! valid_parton(jets(), *part_begin, *idx_begin, soft_pt_regulator, min_pt)) { return false; } ++part_begin; ++idx_begin; // unob -> second parton in own jet if( type() & (unob | qqbar_exb) ){ if( !valid_parton(jets(), *part_begin, *idx_begin, soft_pt_regulator, min_pt) ) return false; ++part_begin; ++idx_begin; } } if(!is_forward_g_to_h(*this)) { if(!valid_parton(jets(), *part_end, *idx_end, soft_pt_regulator, min_pt)) { return false; } ++part_end; ++idx_end; if( type() & (unof | qqbar_exf) ){ 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() & qqbar_mid ){ // find qqbar pair auto begin_qqbar{ std::find_if( part_begin, part_end.base(), [](Particle const & part) -> bool { return part.type != ParticleID::gluon; } )}; assert(begin_qqbar != part_end.base()); long int qqbar_pos{ std::distance(part_begin, begin_qqbar) }; assert(qqbar_pos >= 0); idx_begin+=qqbar_pos; if( !( valid_parton(jets(), *begin_qqbar, *idx_begin, soft_pt_regulator, min_pt) && valid_parton(jets(), *std::next(begin_qqbar), *std::next(idx_begin), soft_pt_regulator, min_pt) )) return false; } return true; } bool Event::valid_incoming() const{ for(std::size_t i=0; i < incoming_.size(); ++i){ if(!(HEJ::nearby_ep(std::abs(incoming_[i].pz()), incoming_[i].E(), TOL*incoming_[i].E()) && (incoming_[i].pt()==0.))) 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