diff --git a/include/HEJ/Event.hh b/include/HEJ/Event.hh index a15dbe4..e272141 100644 --- a/include/HEJ/Event.hh +++ b/include/HEJ/Event.hh @@ -1,369 +1,370 @@ /** \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 "boost/iterator/filter_iterator.hpp" #include "fastjet/ClusterSequence.hh" +#include "HEJ/EWConstants.hh" #include "HEJ/event_types.hh" #include "HEJ/Parameters.hh" #include "HEJ/Particle.hh" #include "HEJ/RNG.hh" namespace LHEF { class HEPEUP; class HEPRUP; } namespace fastjet { class JetDefinition; } namespace HEJ { struct UnclusteredEvent; /** @brief An event with clustered jets * * This is the main HEJ 2 event class. * It contains kinematic information including jet clustering, * parameter (e.g. scale) settings and the event weight. */ class Event { public: class EventData; //! Iterator over partons using ConstPartonIterator = boost::filter_iterator< bool (*)(Particle const &), std::vector::const_iterator >; //! Reverse Iterator over partons using ConstReversePartonIterator = std::reverse_iterator< ConstPartonIterator>; //! No default Constructor Event() = delete; //! Event Constructor adding jet clustering to an unclustered event //! @deprecated UnclusteredEvent will be replaced by EventData in HEJ 2.2.0 [[deprecated("UnclusteredEvent will be replaced by EventData")]] Event( UnclusteredEvent const & ev, fastjet::JetDefinition const & jet_def, double min_jet_pt ); //! @name Particle Access //! @{ //! Incoming particles std::array const & incoming() const{ return incoming_; } //! Outgoing particles std::vector const & outgoing() const{ return outgoing_; } //! Iterator to the first outgoing parton ConstPartonIterator begin_partons() const; //! Iterator to the first outgoing parton ConstPartonIterator cbegin_partons() const; //! Iterator to the end of the outgoing partons ConstPartonIterator end_partons() const; //! Iterator to the end of the outgoing partons ConstPartonIterator cend_partons() const; //! Reverse Iterator to the first outgoing parton ConstReversePartonIterator rbegin_partons() const; //! Reverse Iterator to the first outgoing parton ConstReversePartonIterator crbegin_partons() const; //! Reverse Iterator to the first outgoing parton ConstReversePartonIterator rend_partons() const; //! Reverse Iterator to the first outgoing parton ConstReversePartonIterator crend_partons() const; //! Particle decays /** * The key in the returned map corresponds to the index in the * vector returned by outgoing() */ std::unordered_map> const & decays() const{ return decays_; } //! The jets formed by the outgoing partons, sorted in rapidity std::vector const & jets() const{ return jets_; } //! @} //! @name Weight variations //! @{ //! All chosen parameter, i.e. scale choices (const version) Parameters const & parameters() const{ return parameters_; } //! All chosen parameter, i.e. scale choices Parameters & parameters(){ return parameters_; } //! Central parameter choice (const version) EventParameters const & central() const{ return parameters_.central; } //! Central parameter choice EventParameters & central(){ return parameters_.central; } //! Parameter (scale) variations (const version) std::vector const & variations() const{ return parameters_.variations; } //! Parameter (scale) variations std::vector & variations(){ return parameters_.variations; } //! Parameter (scale) variation (const version) /** * @param i Index of the requested variation */ EventParameters const & variations(size_t i) const{ return parameters_.variations.at(i); } //! Parameter (scale) variation /** * @param i Index of the requested variation */ EventParameters & variations(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(HEJ::RNG &); //! Check that current colours are leading in the high energy limit /** * @details Checks that the colour configuration can be split up in * multiple, rapidity ordered, non-overlapping ladders. Such * configurations are leading in the MRK limit, see * \cite Andersen:2011zd * * @note This is _not_ to be confused with \ref is_resummable(), however * for all resummable states it is possible to create a leading colour * configuration, see generate_colours() */ bool is_leading_colour() const; /** * @brief Check if given event could have been produced by HEJ * @details A HEJ state has to fulfil: * 1. type() has to be \ref is_resummable() "resummable" * 2. Soft radiation in the tagging jets contributes at most to * `max_ext_soft_pt_fraction` of the total jet \f$ p_\perp \f$ * * @note This is true for any resummed stated produced by the * EventReweighter or any \ref is_resummable() "resummable" Leading * Order state. * * @param max_ext_soft_pt_fraction Maximum transverse momentum fraction from * soft radiation in extremal jets * @param min_extparton_pt Absolute minimal \f$ p_\perp \f$, * \b deprecated use max_ext_soft_pt_fraction * instead * @return True if this state could have been produced by HEJ */ bool valid_hej_state( double max_ext_soft_pt_fraction, double min_extparton_pt = 0.) const; private: //! \internal //! @brief Construct Event explicitly from input. /** This is only intended to be called from EventData. * * \warning The input is taken _as is_, sorting and classification has to be * done externally, i.e. by EventData */ Event( std::array && incoming, std::vector && outgoing, std::unordered_map> && decays, Parameters && parameters, fastjet::JetDefinition const & jet_def, double const min_jet_pt ); //! Iterator over partons (non-const) using PartonIterator = boost::filter_iterator< bool (*)(Particle const &), std::vector::iterator >; //! Reverse Iterator over partons (non-const) using ReversePartonIterator = std::reverse_iterator; //! Iterator to the first outgoing parton (non-const) PartonIterator begin_partons(); //! Iterator to the end of the outgoing partons (non-const) PartonIterator end_partons(); //! Reverse Iterator to the first outgoing parton (non-const) ReversePartonIterator rbegin_partons(); //! Reverse Iterator to the first outgoing parton (non-const) ReversePartonIterator rend_partons(); std::array incoming_; std::vector outgoing_; std::unordered_map> decays_; std::vector jets_; Parameters parameters_; fastjet::ClusterSequence cs_; double min_jet_pt_; event_type::EventType type_; }; // end class Event //! Class to store general Event setup, i.e. Phase space and weights class Event::EventData { public: //! Default Constructor EventData() = default; //! Constructor from LesHouches event information EventData(LHEF::HEPEUP const & hepeup); //! Constructor with all values given EventData( std::array incoming, std::vector outgoing, std::unordered_map> decays, Parameters parameters ): incoming(std::move(incoming)), outgoing(std::move(outgoing)), decays(std::move(decays)), parameters(std::move(parameters)) {} //! Generate an Event from the stored EventData. /** * @details Do jet clustering and classification. * Use this to generate an Event. * * @note Calling this function destroys EventData * * @param jet_def Jet definition * @param min_jet_pt minimal \f$p_T\f$ for each jet * * @returns Full clustered and classified event. */ Event cluster( fastjet::JetDefinition const & jet_def, double const 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(); + void reconstruct_intermediate(EWConstants const & ew_parameters); //! Incoming particles std::array incoming; //! Outcoing particles std::vector outgoing; //! Particle decays in the format {outgoing index, decay products} std::unordered_map> decays; //! Parameters, e.g. scale or inital weight Parameters parameters; }; // end class EventData //! Print Event std::ostream& operator<<(std::ostream & os, Event const & ev); //! Square of the partonic centre-of-mass energy \f$\hat{s}\f$ double shat(Event const & ev); //! Convert an event to a LHEF::HEPEUP LHEF::HEPEUP to_HEPEUP(Event const & event, LHEF::HEPRUP *); // put deprecated warning at the end, so don't get the warning inside Event.hh, // additionally doxygen can not identify [[deprecated]] correctly struct [[deprecated("UnclusteredEvent will be replaced by EventData")]] UnclusteredEvent; //! An event before jet clustering //! @deprecated UnclusteredEvent will be replaced by EventData in HEJ 2.2.0 struct UnclusteredEvent{ //! Default Constructor UnclusteredEvent() = default; //! Constructor from LesHouches event information UnclusteredEvent(LHEF::HEPEUP const & hepeup); std::array incoming; /**< Incoming Particles */ std::vector outgoing; /**< Outgoing Particles */ //! Particle decays in the format {outgoing index, decay products} std::unordered_map> decays; //! Central parameter (e.g. scale) choice EventParameters central; std::vector variations; /**< For parameter variation */ }; } // namespace HEJ diff --git a/src/Event.cc b/src/Event.cc index 7f2d860..2ec49e9 100644 --- a/src/Event.cc +++ b/src/Event.cc @@ -1,1128 +1,1128 @@ /** * \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 "LHEF/LHEF.h" #include "fastjet/JetDefinition.hh" #include "HEJ/Constants.hh" #include "HEJ/exceptions.hh" #include "HEJ/PDG_codes.hh" namespace HEJ{ namespace { constexpr int status_in = -1; constexpr int status_decayed = 2; constexpr int status_out = 1; //! 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]); } /// @name helper functions to determine event type //@{ /** * \brief check if final state valid for HEJ * * check final state has the expected number of valid decays for bosons * and all the rest are quarks or gluons */ bool final_state_ok(Event const & ev){ size_t invalid_decays = ev.decays().size(); std::vector const & outgoing = ev.outgoing(); for( size_t i=0; i0?+1:-1, decay->second) ){ --invalid_decays; } else return false; } // Higgs decays (optional) else if(out.type == ParticleID::h){ if(decay != ev.decays().cend()) --invalid_decays; } } else if(! is_parton(out.type)) return false; } // any invalid decays? if(invalid_decays != 0) return false; return true; } /** * returns all EventTypes implemented in HEJ */ size_t implemented_types(std::vector const & bosons){ using namespace event_type; if(bosons.empty()) return FKL | unob | unof | qqxexb | qqxexf | qqxmid; else if(bosons.size() == 1) { switch (bosons[0].type) { case ParticleID::Wp: case ParticleID::Wm: return FKL | unob | unof | qqxexb | qqxexf | qqxmid; case ParticleID::h: return FKL | unob | unof; default: return non_resummable; } } else if(bosons.size() == 2) { if(bosons[0].type == ParticleID::Wp && bosons[1].type == ParticleID::Wp) { return FKL; } return non_resummable; } else { return non_resummable; } } /** * \brief function which determines if type change is consistent with Wp emission. * @param in incoming Particle id * @param out outgoing Particle id * @param qqx Current both incoming/both outgoing? * * \see is_Wm_Change */ bool is_Wp_Change(ParticleID in, ParticleID out, bool qqx){ if(!qqx && (in==-1 || in== 2 || in==-3 || in== 4)) return out== (in-1); if( qqx && (in== 1 || in==-2 || in== 3 || in==-4)) return out==-(in+1); return false; } /** * \brief function which determines if type change is consistent with Wm emission. * @param in incoming Particle id * @param out outgoing Particle id * @param qqx Current both incoming/both outgoing? * * Ensures that change type of quark line is possible by a flavour changing * Wm emission. Allows checking of qqx currents also. */ bool is_Wm_Change(ParticleID in, ParticleID out, bool qqx){ if(!qqx && (in== 1 || in==-2 || in== 3 || in==-4)) return out== (in+1); if( qqx && (in==-1 || in== 2 || in==-3 || in== 4)) return out==-(in-1); return false; } /** * \brief checks if particle type remains same from incoming to outgoing * @param in incoming Particle * @param out outgoing Particle * @param qqx Current both incoming/outgoing? */ bool no_flavour_change(ParticleID in, ParticleID out, bool qqx){ const int qqxCurrent = qqx?-1:1; if(abs(in)<=6 || in==pid::gluon) return (in==out*qqxCurrent); else return false; } bool has_2_jets(Event const & event){ return event.jets().size() >= 2; } /** * \brief check if we have a valid Impact factor * @param in incoming Particle * @param out outgoing Particle * @param qqx Current both incoming/outgoing? * @param W_change returns +1 if Wp, -1 if Wm, else 0 */ bool is_valid_impact_factor( ParticleID in, ParticleID out, bool qqx, int & W_change ){ if( no_flavour_change(in, out, qqx) ){ return true; } if( is_Wp_Change(in, out, qqx) ) { W_change+=1; return true; } if( is_Wm_Change(in, out, qqx) ) { W_change-=1; return true; } return false; } //! Returns all possible classifications from the impact factors // the beginning points are changed s.t. after the the classification they // point to the beginning of the (potential) FKL chain // sets W_change: + if Wp change // 0 if no change // - if Wm change // This function can be used with forward & backwards iterators template size_t possible_impact_factors( ParticleID incoming_id, // incoming OutIterator & begin_out, OutIterator const & end_out, // outgoing int & W_change, std::vector const & boson, bool const backward // backward? ){ using namespace event_type; // keep track of all states that we don't test size_t not_tested = qqxmid; if(backward) not_tested |= unof | qqxexf; else not_tested |= unob | qqxexb; // Is this LL current? if( is_valid_impact_factor(incoming_id, begin_out->type, false, W_change) ){ ++begin_out; return not_tested | FKL; } // or NLL current? // -> needs two partons in two different jets if( std::distance(begin_out, end_out)>=2 ){ auto next = std::next(begin_out); // Is this unordered emisson? if( incoming_id!=pid::gluon && begin_out->type==pid::gluon ){ if( is_valid_impact_factor( incoming_id, next->type, false, W_change ) ){ // veto Higgs inside uno assert(next!=end_out); if( !boson.empty() && boson.front().type == ParticleID::h ){ if( (backward && boson.front().rapidity() < next->rapidity()) ||(!backward && boson.front().rapidity() > next->rapidity())) return non_resummable; } begin_out = std::next(next); return not_tested | (backward?unob:unof); } } // Is this QQbar? else if( incoming_id==pid::gluon ){ if( is_valid_impact_factor( begin_out->type, next->type, true, W_change ) ){ // veto Higgs inside qqx assert(next!=end_out); if( !boson.empty() && boson.front().type == ParticleID::h ){ if( (backward && boson.front().rapidity() < next->rapidity()) ||(!backward && boson.front().rapidity() > next->rapidity())) return non_resummable; } begin_out = std::next(next); return not_tested | (backward?qqxexb:qqxexf); } } } return non_resummable; } //! Returns all possible classifications from central emissions // the beginning points are changed s.t. after the the classification they // point to the end of the emission chain // sets W_change: + if Wp change // 0 if no change // - if Wm change template size_t possible_central( OutIterator & begin_out, OutIterator const & end_out, int & W_change, std::vector const & boson ){ using namespace event_type; // if we already passed the central chain, // then it is not a valid all-order state if(std::distance(begin_out, end_out) < 0) return non_resummable; // keep track of all states that we don't test size_t possible = unob | unof | qqxexb | qqxexf; // Find the first non-gluon/non-FKL while( (begin_out->type==pid::gluon) && (begin_out!=end_out) ){ ++begin_out; } // end of chain -> FKL if( begin_out==end_out ){ return possible | FKL; } // is this a qqbar-pair? // needs two partons in two separate jets auto next = std::next(begin_out); if( is_valid_impact_factor( begin_out->type, next->type, true, W_change ) ){ // veto Higgs inside qqx if( !boson.empty() && boson.front().type == ParticleID::h && boson.front().rapidity() > begin_out->rapidity() && boson.front().rapidity() < next->rapidity() ){ return non_resummable; } begin_out = std::next(next); // remaining chain should be pure gluon/FKL for(; begin_out!=end_out; ++begin_out){ if(begin_out->type != pid::gluon) return non_resummable; } return possible | qqxmid; } return non_resummable; } /** * \brief Checks for all event types * @param ev Event * @returns Event Type * */ event_type::EventType classify(Event const & ev){ using namespace event_type; if(! has_2_jets(ev)) return no_2_jets; // currently we can't handle multiple boson states in the ME. So they are // considered "bad_final_state" even though the "classify" could work with // them. if(! final_state_ok(ev)) return bad_final_state; // initialise variables auto const & in = ev.incoming(); // range for current checks auto begin_out{ev.cbegin_partons()}; auto end_out{ev.crbegin_partons()}; assert(std::distance(begin(in), end(in)) == 2); assert(std::distance(begin_out, end_out.base()) >= 2); assert(std::is_sorted(begin_out, end_out.base(), rapidity_less{})); auto const bosons{ filter_AWZH_bosons(ev.outgoing()) }; // keep track of potential W couplings, at the end the sum should be 0 int remaining_Wp = 0; int remaining_Wm = 0; for(auto const & boson : bosons){ if(boson.type == ParticleID::Wp) ++remaining_Wp; else if(boson.type == ParticleID::Wm) ++remaining_Wm; } size_t final_type = ~(no_2_jets | bad_final_state); // check forward impact factor int W_change = 0; final_type &= possible_impact_factors( in.front().type, begin_out, end_out.base(), W_change, bosons, true ); if( final_type == non_resummable ) return non_resummable; if(W_change>0) remaining_Wp-=W_change; else if(W_change<0) remaining_Wm+=W_change; // check backward impact factor W_change = 0; final_type &= possible_impact_factors( in.back().type, end_out, std::make_reverse_iterator(begin_out), W_change, bosons, false ); if( final_type == non_resummable ) return non_resummable; if(W_change>0) remaining_Wp-=W_change; else if(W_change<0) remaining_Wm+=W_change; // check central emissions W_change = 0; final_type &= possible_central( begin_out, end_out.base(), W_change, bosons ); if( final_type == non_resummable ) return non_resummable; if(W_change>0) remaining_Wp-=W_change; else if(W_change<0) remaining_Wm+=W_change; // Check whether the right number of Ws are present if( remaining_Wp != 0 || remaining_Wm != 0 ) return non_resummable; // result has to be unique if( (final_type & (final_type-1)) != 0) return non_resummable; // check that each sub processes is implemented // (has to be done at the end) if( (final_type & ~implemented_types(bosons)) != 0 ) return non_resummable; return static_cast(final_type); } //@} Particle extract_particle(LHEF::HEPEUP const & hepeup, size_t i){ const ParticleID id = static_cast(hepeup.IDUP[i]); const fastjet::PseudoJet momentum{ hepeup.PUP[i][0], hepeup.PUP[i][1], hepeup.PUP[i][2], hepeup.PUP[i][3] }; if(is_parton(id)) return Particle{ id, std::move(momentum), hepeup.ICOLUP[i] }; return Particle{ id, std::move(momentum), {} }; } bool is_decay_product(std::pair const & mothers){ if(mothers.first == 0) return false; return mothers.second == 0 || mothers.first == mothers.second; } } // namespace anonymous 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] == 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 o1, Particle 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 & decays) { assert(std::is_sorted(begin(decays), end(decays), type_less{})); if(valid_W_decay(+1, decays)) { return ParticleID::Wp; } if(valid_W_decay(-1, decays)) { return ParticleID::Wm; } throw not_implemented{ "final state with decay X -> " + name(decays[0].type) + " + " + name(decays[1].type) }; } Particle reconstruct_boson(std::vector const & decays) { assert(std::is_sorted(begin(decays), end(decays), type_less{})); Particle decayed_boson; decayed_boson.p = decays[0].p + decays[1].p; decayed_boson.type = reconstruct_type(decays); assert(is_AWZH_boson(decayed_boson)); return decayed_boson; } } - void Event::EventData::reconstruct_intermediate() { + void Event::EventData::reconstruct_intermediate(EWConstants const & /* ew_parameters */) { const auto begin_leptons = std::partition( begin(outgoing), end(outgoing), [](Particle const & p) {return !is_anylepton(p);} ); if(begin_leptons == end(outgoing)) return; assert(is_anylepton(*begin_leptons)); std::vector leptons(begin_leptons, end(outgoing)); outgoing.erase(begin_leptons, end(outgoing)); if(leptons.size() != 2) { throw not_implemented{"Final states with one or more than two leptons"}; } std::sort(begin(leptons), end(leptons), type_less{}); outgoing.emplace_back(reconstruct_boson(leptons)); decays.emplace(outgoing.size()-1, std::move(leptons)); } Event Event::EventData::cluster( fastjet::JetDefinition const & jet_def, double const min_jet_pt ){ sort(); Event ev{ std::move(incoming), std::move(outgoing), std::move(decays), std::move(parameters), jet_def, min_jet_pt }; assert(std::is_sorted(begin(ev.outgoing_), end(ev.outgoing_), rapidity_less{})); ev.type_ = classify(ev); return ev; } 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_)); } 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 anonymous bool Event::is_leading_colour() const { if( !correct_colour(incoming()[0]) || !correct_colour(incoming()[1]) ) return false; Colour line_colour = *incoming()[0].colour; std::swap(line_colour.first, line_colour.second); // reasonable colour if(!std::all_of(outgoing().cbegin(), outgoing().cend(), correct_colour)) return false; for(auto it_part = cbegin_partons(); it_part!=cend_partons(); ++it_part){ switch (type()) { case event_type::FKL: if( !try_connect_t(it_part, line_colour) ) return false; break; case event_type::unob: case event_type::qqxexb: { if( !try_connect_t(it_part, line_colour) // u-channel only allowed at impact factor && (std::distance(cbegin_partons(), it_part)!=0 || !try_connect_u(it_part, line_colour))) return false; break; } case event_type::unof: case event_type::qqxexf: { if( !try_connect_t(it_part, line_colour) // u-channel only allowed at impact factor && (std::distance(it_part, cend_partons())!=2 || !try_connect_u(it_part, line_colour))) return false; break; } case event_type::qqxmid:{ auto it_next = std::next(it_part); if( !try_connect_t(it_part, line_colour) // u-channel only allowed at qqx/qxq pair && ( ( !(is_quark(*it_part) && is_antiquark(*it_next)) && !(is_antiquark(*it_part) && is_quark(*it_next))) || !try_connect_u(it_part, line_colour)) ) return false; break; } default: throw std::logic_error{"unreachable"}; } // no colour singlet exchange/disconnected diagram if(line_colour.first == line_colour.second) return false; } return (incoming()[1].colour->first == line_colour.first) && (incoming()[1].colour->second == line_colour.second); } namespace { //! connect incoming Particle to colour flow void connect_incoming(Particle & in, int & colour, int & anti_colour){ in.colour = std::make_pair(anti_colour, colour); // gluon if(in.type == pid::gluon) return; if(in.type > 0){ // quark assert(is_quark(in)); in.colour->second = 0; colour*=-1; return; } // anti-quark assert(is_antiquark(in)); in.colour->first = 0; anti_colour*=-1; return; } //! connect outgoing Particle to t-channel colour flow template void connect_tchannel( OutIterator & it_part, int & colour, int & anti_colour, RNG & ran ){ assert(colour>0 || anti_colour>0); if(it_part->type == ParticleID::gluon){ // gluon if(colour>0 && anti_colour>0){ // on g line => connect to colour OR anti-colour (random) if(ran.flat() < 0.5){ it_part->colour = std::make_pair(colour+2,colour); colour+=2; } else { it_part->colour = std::make_pair(anti_colour, anti_colour+2); anti_colour+=2; } } else if(colour > 0){ // on q line => connect to available colour it_part->colour = std::make_pair(colour+2, colour); colour+=2; } else { assert(colour<0 && anti_colour>0); // on qx line => connect to available anti-colour it_part->colour = std::make_pair(anti_colour, anti_colour+2); anti_colour+=2; } } else if(is_quark(*it_part)) { // quark assert(anti_colour>0); if(colour>0){ // on g line => connect and remove anti-colour it_part->colour = std::make_pair(anti_colour, 0); anti_colour+=2; anti_colour*=-1; } else { // on qx line => new colour colour*=-1; it_part->colour = std::make_pair(colour, 0); } } else if(is_antiquark(*it_part)) { // anti-quark assert(colour>0); if(anti_colour>0){ // on g line => connect and remove colour it_part->colour = std::make_pair(0, colour); colour+=2; colour*=-1; } else { // on q line => new anti-colour anti_colour*=-1; it_part->colour = std::make_pair(0, anti_colour); } } else { // not a parton assert(!is_parton(*it_part)); it_part->colour = {}; } } //! connect to t- or u-channel colour flow template void connect_utchannel( OutIterator & it_part, int & colour, int & anti_colour, RNG & ran ){ OutIterator it_first = it_part++; if(ran.flat()<.5) {// t-channel connect_tchannel(it_first, colour, anti_colour, ran); connect_tchannel(it_part, colour, anti_colour, ran); } else { // u-channel connect_tchannel(it_part, colour, anti_colour, ran); connect_tchannel(it_first, colour, anti_colour, ran); } } } // namespace anonymous bool Event::generate_colours(RNG & ran){ // generate only for HEJ events if(!event_type::is_resummable(type())) return false; assert(std::is_sorted( begin(outgoing()), end(outgoing()), rapidity_less{})); assert(incoming()[0].pz() < incoming()[1].pz()); // positive (anti-)colour -> can connect // negative (anti-)colour -> not available/used up by (anti-)quark int colour = COLOUR_OFFSET; int anti_colour = colour+1; // initialise first connect_incoming(incoming_[0], colour, anti_colour); // reset outgoing colours std::for_each(outgoing_.begin(), outgoing_.end(), [](Particle & part){ part.colour = {};}); for(auto it_part = begin_partons(); it_part!=end_partons(); ++it_part){ switch (type()) { // subleading can connect to t- or u-channel case event_type::unob: case event_type::qqxexb: { if( std::distance(begin_partons(), it_part)==0) connect_utchannel(it_part, colour, anti_colour, ran); else connect_tchannel(it_part, colour, anti_colour, ran); break; } case event_type::unof: case event_type::qqxexf: { if( std::distance(it_part, end_partons())==2) connect_utchannel(it_part, colour, anti_colour, ran); else connect_tchannel(it_part, colour, anti_colour, ran); break; } case event_type::qqxmid:{ auto it_next = std::next(it_part); if( std::distance(begin_partons(), it_part)>0 && std::distance(it_part, end_partons())>2 && ( (is_quark(*it_part) && is_antiquark(*it_next)) || (is_antiquark(*it_part) && is_quark(*it_next)) ) ) connect_utchannel(it_part, colour, anti_colour, ran); else connect_tchannel(it_part, colour, anti_colour, ran); break; } default: // rest has to be t-channel connect_tchannel(it_part, colour, anti_colour, ran); } } // Connect last connect_incoming(incoming_[1], anti_colour, colour); assert(is_leading_colour()); return true; } // generate_colours namespace { bool valid_parton( std::vector const & jets, Particle const & parton, int const idx, double const max_ext_soft_pt_fraction, double const min_extparton_pt ){ // TODO code overlap with PhaseSpacePoint::pass_extremal_cuts if(min_extparton_pt > parton.pt()) return false; if(idx<0) return false; assert(static_cast(jets.size())>=idx); auto const & jet{ jets[idx] }; if( (parton.p - jet).pt()/jet.pt() > max_ext_soft_pt_fraction) return false; return true; } } // this should work with multiple types bool Event::valid_hej_state(double const max_frac, double const min_pt ) const { using namespace event_type; if(!is_resummable(type())) return false; auto const & jet_idx{ particle_jet_indices() }; auto idx_begin{ jet_idx.cbegin() }; auto idx_end{ jet_idx.crbegin() }; auto part_begin{ cbegin_partons() }; auto part_end{ crbegin_partons() }; // always seperate extremal jets if( !valid_parton(jets(), *part_begin, *idx_begin, max_frac, min_pt) ) return false; ++part_begin; ++idx_begin; if( !valid_parton(jets(), *part_end, *idx_end, max_frac, min_pt) ) return false; ++part_end; ++idx_end; // unob -> second parton in own jet if( type() & (unob | qqxexb) ){ if( !valid_parton(jets(), *part_begin, *idx_begin, max_frac, min_pt) ) return false; ++part_begin; ++idx_begin; } if( type() & (unof | qqxexf) ){ if( !valid_parton(jets(), *part_end, *idx_end, max_frac, min_pt) ) return false; ++part_end; ++idx_end; } if( type() & qqxmid ){ // find qqx pair auto begin_qqx{ std::find_if( part_begin, part_end.base(), [](Particle const & part) -> bool { return part.type != ParticleID::gluon; } )}; assert(begin_qqx != part_end.base()); long int qqx_pos{ std::distance(part_begin, begin_qqx) }; assert(qqx_pos >= 0); idx_begin+=qqx_pos; if( !( valid_parton(jets(),*begin_qqx, *idx_begin, max_frac,min_pt) && valid_parton(jets(),*(++begin_qqx),*(++idx_begin),max_frac,min_pt) )) return false; } return true; } Event::ConstPartonIterator Event::begin_partons() const { return cbegin_partons(); } Event::ConstPartonIterator Event::cbegin_partons() const { return boost::make_filter_iterator( static_cast(is_parton), cbegin(outgoing()), cend(outgoing()) ); } Event::ConstPartonIterator Event::end_partons() const { return cend_partons(); } Event::ConstPartonIterator Event::cend_partons() const { return boost::make_filter_iterator( static_cast(is_parton), cend(outgoing()), cend(outgoing()) ); } Event::ConstReversePartonIterator Event::rbegin_partons() const { return crbegin_partons(); } Event::ConstReversePartonIterator Event::crbegin_partons() const { return std::reverse_iterator( cend_partons() ); } Event::ConstReversePartonIterator Event::rend_partons() const { return crend_partons(); } Event::ConstReversePartonIterator Event::crend_partons() const { return std::reverse_iterator( cbegin_partons() ); } Event::PartonIterator Event::begin_partons() { return boost::make_filter_iterator( static_cast(is_parton), begin(outgoing_), end(outgoing_) ); } Event::PartonIterator Event::end_partons() { return boost::make_filter_iterator( static_cast(is_parton), end(outgoing_), end(outgoing_) ); } Event::ReversePartonIterator Event::rbegin_partons() { return std::reverse_iterator( end_partons() ); } Event::ReversePartonIterator Event::rend_partons() { return std::reverse_iterator( begin_partons() ); } namespace { void print_momentum(std::ostream & os, fastjet::PseudoJet const & part){ const std::streamsize orig_prec = os.precision(); os < const & col){ if(!col) os << "(no color)"; // American spelling for better alignment else os << "(" <first << ", " <second << ")"; } } std::ostream& operator<<(std::ostream & os, Event const & ev){ const std::streamsize orig_prec = os.precision(); os < rapidity=" < rapidity=" < 0 ){ os << "\nDecays: " << ev.decays().size() << "\n"; for(auto const & decay: ev.decays()){ os < 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(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)?status_decayed: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(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; } } diff --git a/src/bin/HEJ.cc b/src/bin/HEJ.cc index 1ff402b..67c5c7e 100644 --- a/src/bin/HEJ.cc +++ b/src/bin/HEJ.cc @@ -1,382 +1,387 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #include #include #include #include #include #include #include "yaml-cpp/yaml.h" #include "fastjet/ClusterSequence.hh" #include "HEJ/CombinedEventWriter.hh" #include "HEJ/Config.hh" #include "HEJ/CrossSectionAccumulator.hh" #include "HEJ/Event.hh" #include "HEJ/EventReader.hh" #include "HEJ/BufferedEventReader.hh" #include "HEJ/EventReweighter.hh" #include "HEJ/get_analysis.hh" #include "HEJ/make_RNG.hh" #include "HEJ/optional.hh" #include "HEJ/ProgressBar.hh" #include "HEJ/stream.hh" #include "HEJ/Unweighter.hh" #include "HEJ/Version.hh" #include "HEJ/YAMLreader.hh" HEJ::Config load_config(char const * filename){ try{ return HEJ::load_config(filename); } catch(std::exception const & exc){ std::cerr << "Error: " << exc.what() << '\n'; std::exit(EXIT_FAILURE); } } std::vector> get_analyses( std::vector const & parameters, LHEF::HEPRUP const & heprup ){ try{ return HEJ::get_analyses(parameters, heprup); } catch(std::exception const & exc){ std::cerr << "Failed to load analysis: " << exc.what() << '\n'; std::exit(EXIT_FAILURE); } } // unique_ptr is a workaround: // HEJ::optional is a better fit, but gives spurious errors with g++ 7.3.0 std::unique_ptr> make_progress_bar( std::vector const & xs ) { if(xs.empty()) return {}; const double Born_xs = std::accumulate(begin(xs), end(xs), 0.); return std::make_unique>(std::cout, Born_xs); } std::string time_to_string(const time_t time){ char s[30]; struct tm * p = localtime(&time); strftime(s, 30, "%a %b %d %Y %H:%M:%S", p); return s; } HEJ::Event to_event( LHEF::HEPEUP const & hepeup, - HEJ::JetParameters const & fixed_order_jets + HEJ::JetParameters const & fixed_order_jets, + HEJ::EWConstants const & ew_parameters ) { HEJ::Event::EventData event_data{hepeup}; - event_data.reconstruct_intermediate(); + event_data.reconstruct_intermediate(ew_parameters); return HEJ::Event{ std::move(event_data).cluster( fixed_order_jets.def, fixed_order_jets.min_pt ) }; } void unweight( HEJ::Unweighter & unweighter, HEJ::WeightType weight_type, std::vector & events, HEJ::RNG & ran ) { if(weight_type == HEJ::WeightType::unweighted_resum){ unweighter.set_cut_to_maxwt(events); } events.erase( unweighter.unweight(begin(events), end(events), ran), end(events) ); } // peek up to nevents events from reader std::vector peek_events( HEJ::BufferedEventReader & reader, const int nevents ) { std::vector events; while( static_cast(events.size()) < nevents && reader.read_event() ) { events.emplace_back(reader.hepeup()); } // put everything back into the reader for(auto it = rbegin(events); it != rend(events); ++it) { reader.emplace(*it); } return events; } void append_resummed_events( std::vector & resummation_events, HEJ::EventReweighter & reweighter, LHEF::HEPEUP const & hepeup, const size_t trials, - HEJ::JetParameters const & fixed_order_jets + HEJ::JetParameters const & fixed_order_jets, + HEJ::EWConstants const & ew_parameters ) { - const HEJ::Event FO_event = to_event(hepeup, fixed_order_jets); + const HEJ::Event FO_event = to_event(hepeup, fixed_order_jets, ew_parameters); if(reweighter.treatment(FO_event.type()) != HEJ::EventTreatment::reweight) { return; } const auto resummed = reweighter.reweight(FO_event, trials); resummation_events.insert( end(resummation_events), begin(resummed), end(resummed) ); } void train( HEJ::Unweighter & unweighter, HEJ::BufferedEventReader & reader, HEJ::EventReweighter & reweighter, const size_t total_trials, const double max_dev, double reweight_factor, - HEJ::JetParameters const & fixed_order_jets + HEJ::JetParameters const & fixed_order_jets, + HEJ::EWConstants const & ew_parameters ) { std::cout << "Reading up to " << total_trials << " training events...\n"; auto FO_events = peek_events(reader, total_trials); if(FO_events.empty()) { throw std::runtime_error{ "No events generated to calibrate the unweighting weight!" "Please increase the number \"trials\" or deactivate the unweighting." }; } const size_t trials = total_trials/FO_events.size(); // adjust reweight factor so that the overall normalisation // is the same as in the full run reweight_factor *= trials; for(auto & hepeup: FO_events) { hepeup.XWGTUP *= reweight_factor; } std::cout << "Training unweighter with " << trials << '*' << FO_events.size() << " events\n"; auto progress = HEJ::ProgressBar{ std::cout, static_cast(FO_events.size()) }; std::vector resummation_events; for(auto const & hepeup: FO_events) { append_resummed_events( resummation_events, - reweighter, hepeup, trials, fixed_order_jets + reweighter, hepeup, trials, fixed_order_jets, + ew_parameters ); ++progress; } unweighter.set_cut_to_peakwt(resummation_events, max_dev); std::cout << "\nUnweighting events with weight up to " << unweighter.get_cut() << '\n'; } int main(int argn, char** argv) { using clock = std::chrono::system_clock; if (argn != 3) { std::cerr << "\n# Usage:\n."<< argv[0] <<" config_file input_file\n\n"; return EXIT_FAILURE; } const auto start_time = clock::now(); { std::cout << "Starting " << HEJ::Version::package_name_full() << ", revision " << HEJ::Version::revision() << " (" << time_to_string(clock::to_time_t(start_time)) << ")" << std::endl; } fastjet::ClusterSequence::print_banner(); // read configuration const HEJ::Config config = load_config(argv[1]); auto reader = HEJ::make_reader(argv[2]); assert(reader); auto heprup{ reader->heprup() }; heprup.generators.emplace_back(LHEF::XMLTag{}); heprup.generators.back().name = HEJ::Version::package_name(); heprup.generators.back().version = HEJ::Version::String(); auto analyses = get_analyses( config.analyses_parameters, heprup ); assert(analyses.empty() || analyses.front() != nullptr); HEJ::CombinedEventWriter writer{config.output, std::move(heprup)}; double global_reweight = 1.; const auto & max_events = config.max_events; // if we need the event number: if(std::abs(heprup.IDWTUP) == 4 || std::abs(heprup.IDWTUP) == 1 || max_events){ // try to read from LHE head auto input_events{reader->number_events()}; if(!input_events) { // else count manually auto t_reader = HEJ::make_reader(argv[2]); input_events = 0; while(t_reader->read_event()) ++(*input_events); } if(std::abs(heprup.IDWTUP) == 4 || std::abs(heprup.IDWTUP) == 1){ // IDWTUP 4 or 1 assume average(weight)=xs, but we need sum(weights)=xs std::cout << "Found IDWTUP " << heprup.IDWTUP << ": " << "assuming \"cross section = average weight\".\n" << "converting to \"cross section = sum of weights\" "; global_reweight /= *input_events; } if(max_events && (*input_events > *max_events)){ // maximal number of events given global_reweight *= *input_events/static_cast(*max_events); std::cout << "Processing " << *max_events << " out of " << *input_events << " events\n"; } } HEJ::ScaleGenerator scale_gen{ config.scales.base, config.scales.factors, config.scales.max_ratio }; std::shared_ptr ran{ HEJ::make_RNG(config.rng.name, config.rng.seed)}; assert(ran != nullptr); HEJ::EventReweighter hej{ reader->heprup(), std::move(scale_gen), to_EventReweighterConfig(config), ran }; HEJ::optional unweighter{}; if(config.weight_type != HEJ::WeightType::weighted) { unweighter = HEJ::Unweighter{}; } if(config.weight_type == HEJ::WeightType::partially_unweighted) { HEJ::BufferedEventReader buffered_reader{std::move(reader)}; assert(config.unweight_config); train( *unweighter, buffered_reader, hej, config.unweight_config->trials, config.unweight_config->max_dev, global_reweight/config.trials, - config.fixed_order_jets + config.fixed_order_jets, + config.ew_parameters ); reader = std::make_unique( std::move(buffered_reader) ); } // status infos & eye candy size_t nevent = 0; std::array nevent_type{0}, nfailed_type{0}; auto progress = make_progress_bar(reader->heprup().XSECUP); HEJ::CrossSectionAccumulator xs; std::map status_counter; size_t total_trials = 0; size_t total_resum = 0; // Loop over the events in the input file while(reader->read_event() && (!max_events || nevent < *max_events) ){ ++nevent; // reweight events so that the total cross section is conserved auto hepeup = reader->hepeup(); hepeup.XWGTUP *= global_reweight; - const auto FO_event = to_event(hepeup, config.fixed_order_jets); + const auto FO_event = to_event(hepeup, config.fixed_order_jets, config.ew_parameters); if(FO_event.central().weight == 0) { static const bool warned_once = [argv,nevent](){ std::cerr << "WARNING: event number " << nevent << " in " << argv[2] << " has zero weight. " "Ignoring this and all further events with vanishing weight.\n"; return true; }(); (void) warned_once; // shut up compiler warnings continue; } auto resummed_events{ hej.reweight(FO_event, config.trials) }; // some bookkeeping for(auto const & s: hej.status()) ++status_counter[s]; total_trials+=hej.status().size(); ++nevent_type[FO_event.type()]; if(resummed_events.empty()) ++nfailed_type[FO_event.type()]; if(unweighter) { unweight(*unweighter, config.weight_type, resummed_events, *ran); } // analysis for(auto & ev: resummed_events){ //TODO: move pass_cuts to after phase space point generation bool passed = analyses.empty(); for(auto const & analysis: analyses){ if(analysis->pass_cuts(ev, FO_event)){ passed = true; analysis->fill(ev, FO_event); }; } if(passed){ writer.write(ev); } else { ev.parameters()*=0; // do not use discarded events afterwards } } xs.fill_correlated(resummed_events); total_resum += resummed_events.size(); if(progress) progress->increment(FO_event.central().weight); } // main event loop std::cout << '\n'; for(auto const & analysis: analyses){ analysis->finalise(); } using namespace HEJ::event_type; std::cout<< "Events processed: " << nevent << " (" << total_resum << " resummed)"<< '\n'; std::cout << '\t' << name(EventType::first_type) << ": " << nevent_type[EventType::first_type] << ", failed to reconstruct " << nfailed_type[EventType::first_type] << '\n'; for(auto i=EventType::first_type+1; i<=EventType::last_type; i*=2){ std::cout << '\t' << name(static_cast(i)) << ": " << nevent_type[i] << ", failed to reconstruct " << nfailed_type[i] << '\n'; } std::cout << '\n' << xs << '\n'; std::cout << "Generation statistic: " << status_counter[HEJ::StatusCode::good] << "/" << total_trials << " trials successful.\n"; for(auto && entry: status_counter){ const double fraction = static_cast(entry.second)/total_trials; const int percent = std::round(100*fraction); std::cout << std::left << std::setw(17) << (to_string(entry.first) + ":") << " ["; for(int i = 0; i < percent/2; ++i) std::cout << '#'; for(int i = percent/2; i < 50; ++i) std::cout << ' '; std::cout << "] " < run_time = (clock::now() - start_time); std::cout << "\nFinished " << HEJ::Version::package_name() << " at " << time_to_string(clock::to_time_t(clock::now())) << "\n=> Runtime: " << run_time.count() << " sec (" << nevent/run_time.count() << " Events/sec).\n"; return EXIT_SUCCESS; } diff --git a/t/check_res.cc b/t/check_res.cc index 81b2c8e..b4aca0d 100644 --- a/t/check_res.cc +++ b/t/check_res.cc @@ -1,150 +1,150 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #include #include #include "LHEF/LHEF.h" #include "HEJ/CrossSectionAccumulator.hh" #include "HEJ/Event.hh" #include "HEJ/EventReweighter.hh" #include "HEJ/Mixmax.hh" #include "HEJ/stream.hh" #include "hej_test.hh" namespace{ const fastjet::JetDefinition jet_def{fastjet::kt_algorithm, 0.4}; const fastjet::JetDefinition Born_jet_def{jet_def}; constexpr double Born_jetptmin = 30; constexpr double max_ext_soft_pt_fraction = 0.1; constexpr double jetptmin = 35; constexpr bool log_corr = false; const HEJ::ParticleProperties Wprop{80.385, 2.085}; const HEJ::ParticleProperties Zprop{91.187, 2.495}; const HEJ::ParticleProperties Hprop{125, 0.004165}; constexpr double vev = 246.2196508; using EventTreatment = HEJ::EventTreatment; using namespace HEJ::event_type; HEJ::EventTreatMap treat{ {no_2_jets, EventTreatment::discard}, {bad_final_state, EventTreatment::discard}, {non_resummable, EventTreatment::discard}, {unof, EventTreatment::discard}, {unob, EventTreatment::discard}, {qqxexb, EventTreatment::discard}, {qqxexf, EventTreatment::discard}, {qqxmid, EventTreatment::discard}, {FKL, EventTreatment::reweight} }; bool correct_colour(HEJ::Event const & ev){ if(!HEJ::event_type::is_resummable(ev.type())) return true; if(!ev.is_leading_colour()) return false; return true; } } int main(int argn, char** argv) { if(argn == 5 && std::string(argv[4]) == "unof"){ --argn; treat[unof] = EventTreatment::reweight; treat[unob] = EventTreatment::discard; treat[FKL] = EventTreatment::discard; } if(argn == 5 && std::string(argv[4]) == "unob"){ --argn; treat[unof] = EventTreatment::discard; treat[unob] = EventTreatment::reweight; treat[FKL] = EventTreatment::discard; } else if(argn == 5 && std::string(argv[4]) == "splitf"){ --argn; treat[qqxexb] = EventTreatment::discard; treat[qqxexf] = EventTreatment::reweight; treat[FKL] = EventTreatment::discard; } else if(argn == 5 && std::string(argv[4]) == "splitb"){ --argn; treat[qqxexb] = EventTreatment::reweight; treat[qqxexf] = EventTreatment::discard; treat[FKL] = EventTreatment::discard; } else if(argn == 5 && std::string(argv[4]) == "qqxmid"){ --argn; treat[qqxmid] = EventTreatment::reweight; treat[FKL] = EventTreatment::discard; } if(argn != 4){ std::cerr << "Usage: check_res eventfile xsection tolerance [uno]"; return EXIT_FAILURE; } const double xsec_ref = std::stod(argv[2]); const double tolerance = std::stod(argv[3]); HEJ::istream in{argv[1]}; LHEF::Reader reader{in}; HEJ::PhaseSpacePointConfig psp_conf; psp_conf.jet_param = HEJ::JetParameters{jet_def, jetptmin}; psp_conf.max_ext_soft_pt_fraction = max_ext_soft_pt_fraction; HEJ::MatrixElementConfig ME_conf; ME_conf.log_correction = log_corr; ME_conf.Higgs_coupling = HEJ::HiggsCouplingSettings{}; ME_conf.ew_parameters.set_vevWZH(vev, Wprop, Zprop, Hprop); HEJ::EventReweighterConfig conf; conf.psp_config = std::move(psp_conf); conf.ME_config = std::move(ME_conf); conf.treat = treat; reader.readEvent(); const bool has_Higgs = std::find( begin(reader.hepeup.IDUP), end(reader.hepeup.IDUP), 25 ) != end(reader.hepeup.IDUP); const double mu = has_Higgs?125.:91.188; HEJ::ScaleGenerator scale_gen{ {{std::to_string(mu), HEJ::FixedScale{mu}}}, {}, 1. }; std::shared_ptr ran{std::make_shared()}; HEJ::EventReweighter hej{reader.heprup, std::move(scale_gen), conf, ran}; HEJ::CrossSectionAccumulator xs; do{ auto ev_data = HEJ::Event::EventData{reader.hepeup}; shuffle_particles(ev_data); - ev_data.reconstruct_intermediate(); + ev_data.reconstruct_intermediate(ME_conf.ew_parameters); HEJ::Event ev{ ev_data.cluster( Born_jet_def, Born_jetptmin ) }; auto resummed_events = hej.reweight(ev, 100); for(auto const & res_ev: resummed_events) { ASSERT(correct_colour(res_ev)); ASSERT(std::isfinite(res_ev.central().weight)); // we fill the xs uncorrelated since we only want to test the uncertainty // of the resummation xs.fill(res_ev); } } while(reader.readEvent()); const double xsec = xs.total().value; const double xsec_err = std::sqrt(xs.total().error); const double significance = std::abs(xsec - xsec_ref) / std::sqrt( xsec_err*xsec_err + tolerance*tolerance ); std::cout << xsec_ref << " +/- " << tolerance << " ~ " << xsec << " +- " << xsec_err << " => " << significance << " sigma\n"; if(significance > 3.){ std::cerr << "Cross section is off by over 3 sigma!\n"; return EXIT_FAILURE; } return EXIT_SUCCESS; } diff --git a/t/cmp_events.cc b/t/cmp_events.cc index 0a3c833..492c15d 100644 --- a/t/cmp_events.cc +++ b/t/cmp_events.cc @@ -1,53 +1,61 @@ #include "HEJ/stream.hh" #include "HEJ/Event.hh" #include "LHEF/LHEF.h" +namespace { + const HEJ::ParticleProperties Wprop{80.385, 2.085}; + const HEJ::ParticleProperties Zprop{91.187, 2.495}; + const HEJ::ParticleProperties Hprop{125, 0.004165}; + constexpr double vev = 246.2196508; + const HEJ::EWConstants ew_parameters{vev, Wprop, Zprop, Hprop}; +} + int main(int argn, char** argv) { if(argn != 3){ std::cerr << "Usage: cmp_events eventfile eventfile_no_boson "; return EXIT_FAILURE; } HEJ::istream in{argv[1]}; LHEF::Reader reader{in}; HEJ::istream in_no_boson{argv[2]}; LHEF::Reader reader_no_boson{in_no_boson}; while(reader.readEvent()) { if(! reader_no_boson.readEvent()) { std::cerr << "wrong number of events in " << argv[2] << '\n'; return EXIT_FAILURE; } const auto is_AWZH = [](HEJ::Particle const & p) { return is_AWZH_boson(p); }; const HEJ::Event::EventData event_data{reader.hepeup}; const auto boson = std::find_if( begin(event_data.outgoing), end(event_data.outgoing), is_AWZH ); if(boson == end(event_data.outgoing)) { std::cerr << "no boson in event\n"; return EXIT_FAILURE; } HEJ::Event::EventData data_no_boson{reader_no_boson.hepeup}; - data_no_boson.reconstruct_intermediate(); + data_no_boson.reconstruct_intermediate(ew_parameters); const auto new_boson = std::find_if( begin(event_data.outgoing), end(event_data.outgoing), is_AWZH ); if(new_boson == end(data_no_boson.outgoing)) { std::cerr << "failed to find reconstructed boson\n"; return EXIT_FAILURE; } if(new_boson->type != boson->type) { std::cerr << "reconstructed wrong boson type\n"; return EXIT_FAILURE; } if(!HEJ::nearby(new_boson->p, boson->p)) { std::cerr << "reconstructed boson momentum is off\n"; return EXIT_FAILURE; } } return EXIT_SUCCESS; }