diff --git a/Changes-API.md b/Changes-API.md index 921cc3c..134b89e 100644 --- a/Changes-API.md +++ b/Changes-API.md @@ -1,134 +1,137 @@ # Changelog for HEJ API This log lists only changes on the HEJ API. These are primarily code changes relevant for calling HEJ as an API. This file should only be read as an addition to [`Changes.md`](Changes.md), where the main features are documented. ## Version 2.2 ### 2.2.0 #### New and updated functions * New functions to decide whether particles are massive, charged, charged leptons, antiparticles. * Added `to_string` function for `Event`. * Updated function `Event::EventData.reconstruct_intermediate()` - Now requires an `EWConstants` argument - Added support for WpWp and WmWm events. - In the case of WW same-flavour the reconstruction will minimise the total off-shell momentum. * Replaced `no_2_jets` and `bad_final_state` enumerators in `EventType` with `invalid` and `unknown`. * Added helper functions `is_backward_g_to_h`, `is_forward_g_to_h` to detect incoming gluon to outgoing Higgs transitions. * Updated function `implemented_types()` to now be in HEJ namespace. * `EventWriter` has a new `set_xs_scale` member function to set the ratio between the sum of event weights and the cross section. +#### Removed functions and classes +* The `UnclusteredEvent` class was removed following its deprecation in 2.1.0. + ## Version 2.1 ### 2.1.0 #### Changes to Event class * Restructured `Event` class - `Event` can now only be build from a (new) `Event::EventData` class - Removed default constructor for `Event` - `Event::EventData` replaces the old `UnclusteredEvent` struct. - `UnclusteredEvent` is now **deprecated**, and will be removed in version 2.2.0 - Removed `Event.unclustered()` function - Added new member function `Events.parameters()`, to directly access (underlying) `Parameters` - New member functions `begin_partons`, `end_partons` (`rbegin_partons`, `rend_partons`) with aliases `cbegin_partons`, `cend_partons` (`crbegin_partons`, `crend_partons`) for constant (reversed) iterators over outgoing partons. * New function `Event::EventData.reconstruct_intermediate()` to reconstruct bosons from decays, e.g. `positron + nu_e => Wp` * Added optional Colour charges to particles (`Particle.colour`) - Colour connection in the HEJ limit can be generated via `Event.generate_colours` (automatically done in the resummation) - Colour configuration of input events can be checked with `Event.is_leading_colour` * Added function `Event.valid_hej_state` to check if event _could have_ been produced by `HEJ` according to the `soft pt regulator` cut on the jets #### New and updated functions * Renamed `EventType::nonHEJ` to `EventType::non_resummable` and `is_HEJ()` to `is_resummable()` such that run card is consistent with internal workings * Made `MatrixElement.tree_kin(...)` and `MatrixElement.tree_param(...)` public * New `EventReweighter` member function `treatment` to query the treatment with respect to resummation for the various event types. * Added auxiliary functions to obtain a `std::string` from `EventDescription` (`to_string` for human readable, and `to_simple_string` for easy parsable string) * New `get_analyses` function to read in multiple `HEJ::Analysis` at once, similar to `get_analysis` * New `get_ew_parameters` function to extract electroweak parameters from YAML configuration. #### New classes * New class `Unweighter` to do unweighting * New class `CrossSectionAccumulator` to keep track of Cross Section of the different subprocess * New template struct `Parameters` similar to old `Weights` - `Weights` are now an alias for `Parameters`. Calling `Weights` did not change - `Weights.hh` was replaced by `Parameters.hh`. The old `Weights.hh` header will be removed in version 2.2.0 * Function to multiplication and division of `EventParameters.weight` by double - This can be combined with `Parameters`, e.g. `Parameters*Weights`, see also `Events.parameters()` - Moved `EventParameters` to `Parameters.hh` header * new class `EWConstants` replaces previously hard coded `vev` - `EWConstants` have to be set in the general `Config` and the `MatrixElementConfig` #### Input/Output * New abstract `EventReader` class, as base for reading events from files - Moved LHE file reader to `HEJ::LesHouchesReader` * New (optional) function `finish()` in abstract class `EventWriter`. `finish()` is called _after_ all events are written. * Support reading (`HDF5Reader`) and writing (`HDF5Writer`) `hdf5` files * New `BufferedEventReader` class that allows to put events back into the reader. * New `SherpaLHEReader` to read Sherpa LHE files with correct weights * `get_analysis` now requires `YAML::Node` and `LHEF::HEPRUP` as arguments * Replaced `HepMCInterface` and `HepMCWriter` by `HepMCInterfaceX` and `HepMCWriterX` respectively, with `X` being the major version of HepMC (2 or 3) - Renamed `HepMCInterfaceX::init_kinematics` to `HepMCInterfaceX::init_event` and protected it, use `HepMCInterfaceX::operator()` instead - Removed redundant `HepMCInterfaceX::add_variation` function #### Linking * Export cmake target `HEJ::HEJ` to link directly against `libHEJ` * Preprocessor flags (`HEJ_BUILD_WITH_XYZ`) for enabled/disabled dependencies are now written to `ConfigFlags.hh` * Provide links to version specific object files, e.g. `libHEJ.so -> libHEJ.so.2.1 (soname) -> libHEJ.so.2.1.0` * Removed `LHAPDF` from public interface #### Miscellaneous * Capitalisation of `Config.hh` now matches class `Config` (was `config.hh`) * Renamed `Config::max_ext_soft_pt_fraction` to `Config::soft_pt_regulator`. The new `Config::soft_pt_regulator` parameter is optional and the old `Config::max_ext_soft_pt_fraction` is **deprecated**. * Replaced redundant member `EventReweighterConfig::jet_param` with getter function `EventReweighter.jet_param()` (`JetParameters` are already in `PhaseSpacePointConfig`) * Avoid storing reference to the Random Number Generator inside classes - Constructors of `EventReweighter` now expect `std::shared_ptr` (was reference) - Moved reference to `HEJ::RNG` from constructor of `JetSplitter` to `JetSplitter.split` ## Version 2.0 ### 2.0.4 * Fixed wrong path of `HEJ_INCLUDE_DIR` in `hej-config.cmake` ### 2.0.0 * First release diff --git a/include/HEJ/Event.hh b/include/HEJ/Event.hh index aa45c24..1c0e848 100644 --- a/include/HEJ/Event.hh +++ b/include/HEJ/Event.hh @@ -1,421 +1,390 @@ /** \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 #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; //! Method for accessing implemented types size_t implemented_types(std::vector const & bosons); /** @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$ * 3. Partons related to subleading configurations (uno gluon, qqbar) * must be in separate jets. * * @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); std::string to_string(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/src/Event.cc b/src/Event.cc index 198e4aa..ad9d9f2 100644 --- a/src/Event.cc +++ b/src/Event.cc @@ -1,1509 +1,1489 @@ /** * \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 #include #include "HEJ/event_types.hh" #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 { /** * 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 | UNO | QQBAR; // 1 boson if(bosons.size()== 1) { switch (bosons[0].type) { case ParticleID::Wp: case ParticleID::Wm: return FKL | UNO | QQBAR; case ParticleID::Z_photon_mix: return FKL | UNO; case ParticleID::h: return FKL | UNO; 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; } 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_charge, std::vector const & decays ){ assert(std::abs(w_charge) == 1); 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_charge ) // correct charge return false; // leptonic decay (only check first, second follows from pidsum) if( w_charge == 1 ) // W+ return is_charged_antilepton(decays[0]) || is_neutrino(decays[0]); // W- return is_charged_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; if(decays[0].type != anti(decays[1].type)) { return false; } // leptonic decay (only check first, second follows from above) return is_charged_anylepton(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 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 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 = UNO | EXTREMAL_QQBAR; // 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; } bool decay_conserves_charge( Particle const &parent, std::vector const &products ) { auto charge_diff = charge(parent); for (auto const &p : products) { charge_diff -= charge(p); } return charge_diff == 0; } bool charge_conserved(Event const &ev) { boost::rational charge_diff{0}; for (auto const &in : ev.incoming()) { charge_diff += charge(in); } for (auto const &out : ev.outgoing()) { charge_diff -= charge(out); } if (charge_diff != 0) return false; return std::all_of( ev.decays().begin(), ev.decays().end(), [&ev](auto const &decay) { auto const &[parent, products] = decay; return decay_conserves_charge(ev.outgoing()[parent], products); }); } bool decay_conserves_momentum( Particle const &parent, std::vector const &products, const double tolerance ) { fastjet::PseudoJet total_p; for (auto const &p : products) total_p += p.p; return nearby_ep(parent.p, total_p, tolerance); } bool event_momentum_conserved(Event const &ev, const double tolerance) { return momentum_conserved(ev, tolerance) && std::all_of( ev.decays().begin(), ev.decays().end(), [&ev, tolerance](auto const &decay) { auto const &[parent, products] = decay; return decay_conserves_momentum( ev.outgoing()[parent], products, tolerance); }); } template bool massless_particles_onshell( Container const &c, const double tolerance ) { return std::all_of( c.begin(), c.end(), [tolerance](Particle const & p) { return is_massive(p) || p.m() < tolerance * std::max(p.E(), 1.0); } ); } bool all_massless_particles_onshell( Event const &ev, const double tolerance ) { return massless_particles_onshell(ev.incoming(), tolerance) && massless_particles_onshell(ev.outgoing(), tolerance) && std::all_of( ev.decays().begin(), ev.decays().end(), [tolerance](auto const &decay) { return massless_particles_onshell(decay.second, tolerance); }); } bool no_incoming_pt(Event const &ev, const double tolerance) { return std::all_of( ev.incoming().cbegin(), ev.incoming().end(), [tolerance](Particle const &p) { return std::abs(p.px()) < tolerance && std::abs(p.py()) < tolerance; }); } bool is_invalid(Event const &ev, const double tolerance) { return !( charge_conserved(ev) && event_momentum_conserved(ev, tolerance) && no_incoming_pt(ev, tolerance) && all_massless_particles_onshell(ev, tolerance) ); } // TODO: choose reasonable value or make configurable constexpr double TOLERANCE = 1e-3; bool incoming_are_partons(Event const &ev) { return std::all_of( ev.incoming().begin(), ev.incoming().end(), [](Particle const &p) { return is_parton(p); } ); } bool known_outgoing(Event const &ev) { return std::all_of( ev.outgoing().begin(), ev.outgoing().end(), [](Particle const &p) { return is_parton(p) || p.type == pid::Higgs || std::abs(p.type) == pid::Wp || p.type == pid::Z_photon_mix; }); } bool is_same_sign_WW(std::vector const &particles) { return particles.size() == 2 && std::abs(particles.front().type) == pid::Wp && particles.front().type == particles.back().type; } bool all_W_Zphoton_decay(Event const &ev) { auto const &out = ev.outgoing(); for (std::size_t i = 0; i < out.size(); ++i) { if ( (std::abs(out[i].type) == pid::Wp || out[i].type == pid::Z_photon_mix) && ev.decays().count(i) == 0 ) { return false; } } return true; } bool decay_known( Particle const &parent, std::vector const &products ) { if (parent.type == pid::Higgs) return true; if (parent.type == pid::Z_photon_mix) return valid_Z_decay(products); if (std::abs(parent.type) == pid::Wp) { assert(charge(parent).denominator() == 1); return valid_W_decay(charge(parent).numerator(), products); } return false; } bool all_decays_known(Event const &ev) { return std::all_of( ev.decays().begin(), ev.decays().end(), [&ev](auto const &decay) { auto const &[parent, products] = decay; return decay_known(ev.outgoing()[parent], products); }); } bool is_known_process_type(Event const &ev) { if (!incoming_are_partons(ev)) return false; if (!known_outgoing(ev)) return false; if (!all_W_Zphoton_decay(ev)) return false; if (!all_decays_known(ev)) return false; auto const bosons = filter_AWZH_bosons(ev.outgoing()); if (bosons.size() > 2) return false; if (bosons.size() == 2 && !is_same_sign_WW(bosons)) { return false; } if (bosons.size() == 1 && bosons.front().type == pid::Higgs) { return !ev.jets().empty(); } return ev.jets().size() >= 2; } } /** * \brief Checks for all event types * @param ev Event * @returns Event Type * */ event_type::EventType classify(Event const & ev){ using namespace event_type; if(is_invalid(ev, TOLERANCE)) return invalid; if(! is_known_process_type(ev)) return unknown; // 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 = VALID; // 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 bool Event::valid_hej_state(double const soft_pt_regulator, double const min_pt ) const { using namespace event_type; const auto is_valid_parton = [&](Particle const & parton, int const jet_idx) { return valid_parton(jets(), parton, jet_idx, soft_pt_regulator, min_pt); }; if(!is_resummable(type())) return false; auto const & jet_indices{ particle_jet_indices() }; auto jet_idx_begin{ jet_indices.cbegin() }; auto jet_idx_end{ jet_indices.crbegin() }; auto part_begin{ cbegin_partons() }; auto part_end{ crbegin_partons() }; if(!is_backward_g_to_h(*this)) { const int first_jet_idx = *jet_idx_begin; if(! is_valid_parton(*part_begin, first_jet_idx)) { return false; } ++part_begin; ++jet_idx_begin; // unob -> second parton in own jet if( type() & (unob | qqbar_exb) ){ if( (*jet_idx_begin == first_jet_idx) || !is_valid_parton(*part_begin, *jet_idx_begin) ) { return false; } ++part_begin; ++jet_idx_begin; } } if(!is_forward_g_to_h(*this)) { const int last_jet_idx = *jet_idx_end; if(!is_valid_parton(*part_end, last_jet_idx)) { return false; } ++part_end; ++jet_idx_end; if( type() & (unof | qqbar_exf) ){ if( (*jet_idx_end == last_jet_idx) || !is_valid_parton(*part_end, *jet_idx_end) ) { return false; } ++part_end; // ++jet_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); jet_idx_begin += qqbar_pos; const int next_jet_idx = *std::next(jet_idx_begin); if( (*jet_idx_begin == next_jet_idx) || ! is_valid_parton(*begin_qqbar, *jet_idx_begin) || ! is_valid_parton(*std::next(begin_qqbar), next_jet_idx) ) { 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