diff --git a/doc/doxygen/biblio.bib b/doc/doxygen/biblio.bib index 5fce287..3e960be 100644 --- a/doc/doxygen/biblio.bib +++ b/doc/doxygen/biblio.bib @@ -1,56 +1,86 @@ @article{Andersen:2011hs, author = "Andersen, Jeppe R. and Smillie, Jennifer M.", title = "{Multiple Jets at the LHC with High Energy Jets}", journal = "JHEP", volume = "06", year = "2011", pages = "010", doi = "10.1007/JHEP06(2011)010", eprint = "1101.5394", archivePrefix = "arXiv", primaryClass = "hep-ph", reportNumber = "CP3-ORIGINS-2011-02, EDINBURGH-2011-03", SLACcitation = "%%CITATION = ARXIV:1101.5394;%%" } @article{James:1993np, author = "James, F.", title = "{RANLUX: A FORTRAN implementation of the high quality pseudorandom number generator of Luscher}", journal = "Comput. Phys. Commun.", volume = "79", year = "1994", pages = "111-114", doi = "10.1016/0010-4655(94)90233-X", note = "[Erratum: Comput. Phys. Commun.97,357(1996)]", reportNumber = "CERN-CN-93-13", SLACcitation = "%%CITATION = CPHCB,79,111;%%" } @article{Luscher:1993dy, author = "Luscher, Martin", title = "{A Portable high quality random number generator for lattice field theory simulations}", journal = "Comput. Phys. Commun.", volume = "79", year = "1994", pages = "100-110", doi = "10.1016/0010-4655(94)90232-1", eprint = "hep-lat/9309020", archivePrefix = "arXiv", primaryClass = "hep-lat", reportNumber = "DESY-93-133", SLACcitation = "%%CITATION = HEP-LAT/9309020;%%" } @article{Savvidy:2014ana, author = "Savvidy, Konstantin G.", title = "{The MIXMAX random number generator}", journal = "Comput. Phys. Commun.", volume = "196", year = "2015", pages = "161-165", doi = "10.1016/j.cpc.2015.06.003", eprint = "1403.5355", archivePrefix = "arXiv", primaryClass = "hep-lat", reportNumber = "NITS-PHY-2014, NITS-PHY-2014003", SLACcitation = "%%CITATION = ARXIV:1403.5355;%%" } +@inproceedings{Boos:2001cv, + author = "Boos, E. and others", + title = "{Generic user process interface for event generators}", + booktitle = "{Physics at TeV colliders. Proceedings, Euro Summer + School, Les Houches, France, May 21-June 1, 2001}", + url = "http://lss.fnal.gov/archive/preprint/fermilab-conf-01-496-t.shtml", + year = "2001", + eprint = "hep-ph/0109068", + archivePrefix = "arXiv", + primaryClass = "hep-ph", + reportNumber = "FERMILAB-CONF-01-496-T", + SLACcitation = "%%CITATION = HEP-PH/0109068;%%" +} +@article{Andersen:2011zd, + author = "Andersen, Jeppe R. and Lonnblad, Leif and Smillie, + Jennifer M.", + title = "{A Parton Shower for High Energy Jets}", + journal = "JHEP", + volume = "07", + year = "2011", + pages = "110", + doi = "10.1007/JHEP07(2011)110", + eprint = "1104.1316", + archivePrefix = "arXiv", + primaryClass = "hep-ph", + reportNumber = "CERN-PH-TH-2011-072, CP3-ORIGINS-2011-14, + EDINBURGH-2011-16, LU-TP-11-15, MCNET-11-12, LU-TP + --11-15", + SLACcitation = "%%CITATION = ARXIV:1104.1316;%%" +} diff --git a/include/HEJ/Constants.hh b/include/HEJ/Constants.hh index c06932d..bce7476 100644 --- a/include/HEJ/Constants.hh +++ b/include/HEJ/Constants.hh @@ -1,29 +1,34 @@ /** \file * \brief Header file defining all global constants used for HEJ * * \authors Jeppe Andersen, Tuomas Hapola, Marian Heil, Andreas Maier, Jennifer Smillie * \date 2019 * \copyright GPLv2 or later */ #pragma once namespace HEJ{ /// @name QCD parameters //@{ constexpr double N_C = 3.; //!< number of Colours constexpr double C_A = N_C; //!< \f$C_A\f$ constexpr double C_F = (N_C*N_C - 1.)/(2.*N_C); //!< \f$C_F\f$ constexpr double t_f = 0.5; //!< \f$t_f\f$ constexpr double n_f = 5.; //!< number light flavours constexpr double beta0 = 11./3.*C_A - 4./3.*t_f*n_f; //!< \f$\beta_0\f$ //@} /// @name QFT parameters //@{ constexpr double vev = 246.2196508; //!< Higgs vacuum expectation value in GeV //@} /// @name Generation Parameters //@{ constexpr double CLAMBDA = 0.2; //!< Scale for virtual correction, \f$\lambda\f$ cf. eq. (20) in \cite Andersen:2011hs constexpr double CMINPT = CLAMBDA; //!< minimal \f$p_t\f$ of all partons //@} +/// @name Conventional Parameters +//@{ + //! Value of first colour for colour dressing, according to LHE convention \cite Boos:2001cv + constexpr int COLOUR_OFFSET = 501; +//@} } diff --git a/include/HEJ/Event.hh b/include/HEJ/Event.hh index 6b59620..66e6c55 100644 --- a/include/HEJ/Event.hh +++ b/include/HEJ/Event.hh @@ -1,199 +1,201 @@ /** \file * \brief Declares the Event class and helpers * * \authors Jeppe Andersen, Tuomas Hapola, Marian Heil, Andreas Maier, Jennifer Smillie * \date 2019 * \copyright GPLv2 or later */ #pragma once #include #include #include #include #include #include "HEJ/event_types.hh" #include "HEJ/Particle.hh" #include "fastjet/ClusterSequence.hh" namespace LHEF{ class HEPEUP; class HEPRUP; } namespace fastjet{ class JetDefinition; } namespace HEJ{ struct ParameterDescription; //! Event parameters struct EventParameters{ double mur; /**< Value of the Renormalisation Scale */ double muf; /**< Value of the Factorisation Scale */ double weight; /**< Event Weight */ //! Optional description std::shared_ptr description = nullptr; }; //! Description of event parameters struct ParameterDescription { //! Name of central scale choice (e.g. "H_T/2") std::string scale_name; //! Actual renormalisation scale divided by central scale double mur_factor; //! Actual factorisation scale divided by central scale double muf_factor; ParameterDescription() = default; ParameterDescription( std::string scale_name, double mur_factor, double muf_factor ): scale_name{scale_name}, mur_factor{mur_factor}, muf_factor{muf_factor} {}; }; - //! An event before jet clustering + //! An event before jet clustering & classification 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 */ + //! Generate the particle colour leading in the MRK limit, see \cite Andersen:2011zd + void generate_colour_flow(); }; /** 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: //! Default Event Constructor Event() = default; //! Event Constructor adding jet clustering to an unclustered event Event( UnclusteredEvent ev, fastjet::JetDefinition const & jet_def, double min_jet_pt ); //! The jets formed by the outgoing partons std::vector jets() const; //! The corresponding event before jet clustering UnclusteredEvent const & unclustered() const { return ev_; } //! Central parameter choice EventParameters const & central() const{ return ev_.central; } //! Central parameter choice EventParameters & central(){ return ev_.central; } //! Incoming particles std::array const & incoming() const{ return ev_.incoming; } //! Outgoing particles std::vector const & outgoing() const{ return ev_.outgoing; } //! 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 ev_.decays; } //! Parameter (scale) variations std::vector const & variations() const{ return ev_.variations; } //! Parameter (scale) variations std::vector & variations(){ return ev_.variations; } //! Parameter (scale) variation /** * @param i Index of the requested variation */ EventParameters const & variations(size_t i) const{ return ev_.variations[i]; } //! Parameter (scale) variation /** * @param i Index of the requested variation */ EventParameters & variations(size_t i){ return ev_.variations[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); } //! 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_; } private: UnclusteredEvent ev_; fastjet::ClusterSequence cs_; double min_jet_pt_; event_type::EventType type_; }; //! 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 *); } diff --git a/src/Event.cc b/src/Event.cc index 8cb2a08..9f4edfd 100644 --- a/src/Event.cc +++ b/src/Event.cc @@ -1,384 +1,389 @@ /** * \authors Jeppe Andersen, Tuomas Hapola, Marian Heil, Andreas Maier, Jennifer Smillie * \date 2019 * \copyright GPLv2 or later */ #include "HEJ/Event.hh" #include #include #include #include #include "LHEF/LHEF.h" #include "fastjet/JetDefinition.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; /// @name helper functions to determine event type //@{ /** * \brief check if final state valid for HEJ * * check if there is at most one photon, W, H, Z in the final state * and all the rest are quarks or gluons */ bool final_state_ok(std::vector const & outgoing){ bool has_AWZH_boson = false; for(auto const & out: outgoing){ if(is_AWZH_boson(out.type)){ if(has_AWZH_boson) return false; has_AWZH_boson = true; } else if(! is_parton(out.type)) return false; } return true; } template Iterator remove_AWZH(Iterator begin, Iterator end){ return std::remove_if( begin, end, [](Particle const & p){return is_AWZH_boson(p);} ); } template bool valid_outgoing(Iterator begin, Iterator end){ return std::distance(begin, end) >= 2 && std::is_sorted(begin, end, rapidity_less{}) && std::count_if( begin, end, [](Particle const & s){return is_AWZH_boson(s);} ) < 2; } /// @note that this changes the outgoing range! template bool is_FKL( ConstIterator begin_incoming, ConstIterator end_incoming, Iterator begin_outgoing, Iterator end_outgoing ){ assert(std::distance(begin_incoming, end_incoming) == 2); assert(std::distance(begin_outgoing, end_outgoing) >= 2); // One photon, W, H, Z in the final state is allowed. // Remove it for remaining tests, end_outgoing = remove_AWZH(begin_outgoing, end_outgoing); // Test if this is a standard FKL configuration. return (begin_incoming->type == begin_outgoing->type) && ((end_incoming-1)->type == (end_outgoing-1)->type) && std::all_of( begin_outgoing + 1, end_outgoing - 1, [](Particle const & p){ return p.type == pid::gluon; } ); } bool is_FKL( std::array const & incoming, std::vector outgoing ){ assert(std::is_sorted(begin(incoming), end(incoming), pz_less{})); assert(valid_outgoing(begin(outgoing), end(outgoing))); return is_FKL( begin(incoming), end(incoming), begin(outgoing), end(outgoing) ); } bool has_2_jets(Event const & event){ return event.jets().size() >= 2; } /** * \brief Checks whether event is unordered backwards * @param ev Event * @returns Is Event Unordered Backwards * * - Checks there is more than 3 constuents in the final state * - Checks there is more than 3 jets * - Checks the most backwards parton is a gluon * - Checks the most forwards jet is not a gluon * - Checks the rest of the event is FKL * - Checks the second most backwards is not a different boson * - Checks the unordered gluon actually forms a jet */ bool is_unordered_backward(Event const & ev){ auto const & in = ev.incoming(); auto const & out = ev.outgoing(); assert(std::is_sorted(begin(in), end(in), pz_less{})); assert(valid_outgoing(begin(out), end(out))); if(out.size() < 3) return false; if(ev.jets().size() < 3) return false; if(in.front().type == pid::gluon) return false; if(out.front().type != pid::gluon) return false; // When skipping the unordered emission // the remainder should be a regular FKL event, // except that the (new) first outgoing particle must not be a A,W,Z,H. const auto FKL_begin = next(begin(out)); if(is_AWZH_boson(*FKL_begin)) return false; if(!is_FKL(in, {FKL_begin, end(out)})) return false; // check that the unordered gluon forms an extra jet const auto jets = sorted_by_rapidity(ev.jets()); const auto indices = ev.particle_jet_indices({jets.front()}); return indices[0] >= 0 && indices[1] == -1; } /** * \brief Checks for a forward unordered gluon emission * @param ev Event * @returns Is the event a forward unordered emission * * \see is_unordered_backward */ bool is_unordered_forward(Event const & ev){ auto const & in = ev.incoming(); auto const & out = ev.outgoing(); assert(std::is_sorted(begin(in), end(in), pz_less{})); assert(valid_outgoing(begin(out), end(out))); if(out.size() < 3) return false; if(ev.jets().size() < 3) return false; if(in.back().type == pid::gluon) return false; if(out.back().type != pid::gluon) return false; // When skipping the unordered emission // the remainder should be a regular FKL event, // except that the (new) last outgoing particle must not be a A,W,Z,H. const auto FKL_end = prev(end(out)); if(is_AWZH_boson(*prev(FKL_end))) return false; if(!is_FKL(in, {begin(out), FKL_end})) return false; // check that the unordered gluon forms an extra jet const auto jets = sorted_by_rapidity(ev.jets()); const auto indices = ev.particle_jet_indices({jets.back()}); return indices.back() >= 0 && indices[indices.size()-2] == -1; } using event_type::EventType; EventType classify(Event const & ev){ if(! final_state_ok(ev.outgoing())) return EventType::bad_final_state; if(! has_2_jets(ev)) return EventType::no_2_jets; if(is_FKL(ev.incoming(), ev.outgoing())) return EventType::FKL; if(is_unordered_backward(ev)){ return EventType::unordered_backward; } if(is_unordered_forward(ev)){ return EventType::unordered_forward; } return EventType::nonHEJ; } //@} Particle extract_particle(LHEF::HEPEUP const & hepeup, int i){ - return Particle{ - static_cast(hepeup.IDUP[i]), - fastjet::PseudoJet{ - hepeup.PUP[i][0], hepeup.PUP[i][1], - hepeup.PUP[i][2], hepeup.PUP[i][3] - } + 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 UnclusteredEvent::UnclusteredEvent(LHEF::HEPEUP const & hepeup): central(EventParameters{ hepeup.scales.mur, hepeup.scales.muf, hepeup.weight() }) { 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)); } } + void UnclusteredEvent::generate_colour_flow(){ + //TODO implement + } + Event::Event( UnclusteredEvent ev, fastjet::JetDefinition const & jet_def, double min_jet_pt ): ev_{std::move(ev)}, cs_{to_PseudoJet(filter_partons(ev_.outgoing)), jet_def}, min_jet_pt_{min_jet_pt} { // sort particles std::sort( begin(ev_.incoming), end(ev_.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(); }); ev_.outgoing.clear(); ev_.outgoing.reserve(old_outgoing.size()); for(size_t i: idx) { ev_.outgoing.emplace_back(std::move(old_outgoing[i])); } // find decays again if(!ev_.decays.empty()){ auto old_decays = std::move(ev_.decays); ev_.decays.clear(); for(size_t i=0; isecond)); } assert(old_decays.size() == ev_.decays.size()); } // classify event type_ = classify(*this); assert(std::is_sorted(begin(outgoing()), end(outgoing()), rapidity_less{})); } std::vector Event::jets() const{ return cs_.inclusive_jets(min_jet_pt_); } double shat(Event const & ev){ return (ev.incoming()[0].p + ev.incoming()[1].p).m2(); } namespace{ // colour flow according to Les Houches standard // TODO: stub std::vector> colour_flow( std::array const & incoming, std::vector const & outgoing ){ std::vector> result( incoming.size() + outgoing.size() ); for(auto & col: result){ col = std::make_pair(-1, -1); } return result; } } LHEF::HEPEUP to_HEPEUP(Event const & event, LHEF::HEPRUP * heprup){ LHEF::HEPEUP result; result.heprup = heprup; result.weights = {{event.central().weight, nullptr}}; for(auto const & var: event.variations()){ result.weights.emplace_back(var.weight, nullptr); } size_t num_particles = event.incoming().size() + event.outgoing().size(); for(auto const & decay: event.decays()) num_particles += decay.second.size(); result.NUP = num_particles; // the following entries are pretty much meaningless result.IDPRUP = event.type()+1; // event ID result.AQEDUP = 1./128.; // alpha_EW //result.AQCDUP = 0.118 // alpha_QCD // end meaningless part result.XWGTUP = event.central().weight; result.SCALUP = event.central().muf; result.scales.muf = event.central().muf; result.scales.mur = event.central().mur; result.scales.SCALUP = event.central().muf; result.pdfinfo.p1 = event.incoming().front().type; result.pdfinfo.p2 = event.incoming().back().type; result.pdfinfo.scale = event.central().muf; for(Particle const & in: event.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); } 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); } result.ICOLUP = colour_flow( event.incoming(), filter_partons(event.outgoing()) ); if(result.ICOLUP.size() < num_particles){ const size_t AWZH_boson_idx = std::find_if( begin(event.outgoing()), end(event.outgoing()), [](Particle const & s){ return is_AWZH_boson(s); } ) - begin(event.outgoing()) + event.incoming().size(); assert(AWZH_boson_idx <= result.ICOLUP.size()); result.ICOLUP.insert( begin(result.ICOLUP) + AWZH_boson_idx, std::make_pair(0,0) ); } 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; } }