From 062eb92a7054c2f7f811d81c2eb456c44e8b6f0a Mon Sep 17 00:00:00 2001 From: Andy Buckley Date: Fri, 21 May 2021 13:17:20 +0100 Subject: [PATCH 01/27] Pull in new external version of HEPUtils, matching Gambit's constness scheme, and introducing all-particles storage and multiple jet collections --- contrib/heputils/include/HEPUtils/BinnedFn.h | 2 +- contrib/heputils/include/HEPUtils/Event.h | 312 +++++++++--------- contrib/heputils/include/HEPUtils/FastJet.h | 44 +-- contrib/heputils/include/HEPUtils/Jet.h | 53 ++- contrib/heputils/include/HEPUtils/MathUtils.h | 29 +- contrib/heputils/include/HEPUtils/Particle.h | 80 ++++- .../heputils/include/HEPUtils/PhaseSpace.h | 8 +- contrib/heputils/include/HEPUtils/Utils.h | 16 +- contrib/heputils/include/HEPUtils/Vectors.h | 21 +- 9 files changed, 311 insertions(+), 254 deletions(-) diff --git a/contrib/heputils/include/HEPUtils/BinnedFn.h b/contrib/heputils/include/HEPUtils/BinnedFn.h index 6e52ee473b..45f4b513f2 100644 --- a/contrib/heputils/include/HEPUtils/BinnedFn.h +++ b/contrib/heputils/include/HEPUtils/BinnedFn.h @@ -1,7 +1,7 @@ // -*- C++ -*- // // This file is part of HEPUtils -- https://bitbucket.org/andybuckley/heputils -// Copyright (C) 2013-2018 Andy Buckley +// Copyright (C) 2013-2021 Andy Buckley // // Embedding of HEPUtils code in other projects is permitted provided this // notice is retained and the HEPUtils namespace and include path are changed. diff --git a/contrib/heputils/include/HEPUtils/Event.h b/contrib/heputils/include/HEPUtils/Event.h index 987e11832b..b4d123f1da 100644 --- a/contrib/heputils/include/HEPUtils/Event.h +++ b/contrib/heputils/include/HEPUtils/Event.h @@ -1,7 +1,7 @@ // -*- C++ -*- // // This file is part of HEPUtils -- https://bitbucket.org/andybuckley/heputils -// Copyright (C) 2013-2018 Andy Buckley +// Copyright (C) 2013-2021 Andy Buckley // // Embedding of HEPUtils code in other projects is permitted provided this // notice is retained and the HEPUtils namespace and include path are changed. @@ -11,6 +11,7 @@ #include "HEPUtils/Particle.h" #include "HEPUtils/Jet.h" #include +#include namespace HEPUtils { @@ -20,26 +21,26 @@ namespace HEPUtils { private: /// @name Internal particle / vector containers - //@{ + /// @{ /// Event weights std::vector _weights; std::vector _weight_errs; /// @name Separate particle collections - //@{ - std::vector _photons, _electrons, _muons, _taus, _invisibles; - std::vector _cphotons, _celectrons, _cmuons, _ctaus, _cinvisibles; - //@} + /// @{ + std::vector _allparticles; + mutable std::vector _visibles, _invisibles, _photons, _electrons, _muons, _taus; + /// @} - /// Jets collection (mutable to allow sorting) - mutable std::vector _jets; - mutable std::vector _cjets; + /// Jets collection(s) (mutable to allow sorting) + mutable std::map> _jets; /// Missing momentum vector P4 _pmiss; - //@} + /// @} + private: @@ -49,27 +50,22 @@ namespace HEPUtils { clear(); //< Delete current particles _weights = e._weights; _weight_errs = e._weight_errs; + // + _allparticles = e._allparticles; + _visibles = e._visibles; + _invisibles = e._invisibles; _photons = e._photons; - _cphotons = e._cphotons; _electrons = e._electrons; - _celectrons = e._celectrons; _muons = e._muons; - _cmuons = e._cmuons; _taus = e._taus; - _ctaus = e._ctaus; - _invisibles = e._invisibles; - _cinvisibles = e._cinvisibles; + // _jets = e._jets; - _cjets = e._cjets; _pmiss = e._pmiss; } public: - /// @name Constructors - //@{ - /// Default constructor Event() { clear(); } @@ -89,7 +85,8 @@ namespace HEPUtils { } - public: + /// @name Cloning (= deep copy) + /// @{ /// Clone a copy on the heap Event* clone() const { @@ -119,52 +116,50 @@ namespace HEPUtils { e._pmiss = _pmiss; } - //@} + /// @} - /// Empty the event's particle, jet and MET collections + /// Empty the event's weight, particle, jet, and MET collections void clear() { + // Weights _weights.clear(); _weight_errs.clear(); - // TODO: indexed loop -> for (Particle* p : particles()) delete p; - #define DELCLEAR(v) do { if (!v.empty()) for (size_t i = 0; i < v.size(); ++i) delete v[i]; v.clear(); } while (0) - DELCLEAR(_photons); - DELCLEAR(_electrons); - DELCLEAR(_muons); - DELCLEAR(_taus); - DELCLEAR(_invisibles); - DELCLEAR(_jets); - #undef DELCLEAR + // Particles, first the canonical collection + for (const Particle* p : _allparticles) delete p; + _allparticles.clear(); + // Then the caches + _visibles.clear(); + _invisibles.clear(); _photons.clear(); - _cphotons.clear(); _electrons.clear(); - _celectrons.clear(); _muons.clear(); - _cmuons.clear(); _taus.clear(); - _ctaus.clear(); - _invisibles.clear(); - _cinvisibles.clear(); + + // Jets + for (auto& js : _jets) { + for (const Jet* j : js.second) delete j; + } _jets.clear(); - _cjets.clear(); + // MET _pmiss.clear(); } - /////////////////////// + /// @name Weights + /// @{ /// Set the event weights (also possible directly via non-const reference) void set_weights(const std::vector& ws) { _weights = ws; } + /// Set the event weight errors (also possible directly via non-const reference) void set_weight_errs(const std::vector& werrs) { _weight_errs = werrs; } - /// Set the event weights to the single given weight void set_weight(double w) { _weights.clear(); @@ -196,7 +191,6 @@ namespace HEPUtils { std::vector& weight_errs() { return _weight_errs; } - /// Get a single event weight -- the nominal, by default double weight(size_t i=0) const { if (_weights.empty()) { @@ -215,11 +209,13 @@ namespace HEPUtils { return _weight_errs[i]; } + /// @} - ///////////////// + /// @name Particles + /// @{ - /// Add a particle to the event + /// @brief Add a particle to the event /// /// Supplied particle should be new'd, and Event will take ownership. /// @@ -228,190 +224,188 @@ namespace HEPUtils { /// immediately deleted. Accordingly, the pointer passed by user code /// must be considered potentially invalid from the moment this function is called. /// + /// @note pT-sorting can be disabled, so + /// /// @todo "Lock" at some point so that jet finding etc. only get done once - void add_particle(Particle* p) { - if (!p->is_prompt()) - delete p; - else if (p->pid() == 22) - { - _photons.push_back(p); - _cphotons.push_back(p); - } - else if (p->abspid() == 11) - { - _electrons.push_back(p); - _celectrons.push_back(p); - } - else if (p->abspid() == 13) - { - _muons.push_back(p); - _cmuons.push_back(p); + void add_particle(Particle* p, bool ptsort=true) { + + // All particles (canonical collection) + _allparticles.push_back(p); + + // Caching collections + if (!p->is_visible()) { + if (p->is_prompt()) _invisibles.push_back(p); + } else { + _visibles.push_back(p); + if (p->is_prompt()) { + if (p->pid() == 22) _photons.push_back(p); + else if (p->abspid() == 11) _electrons.push_back(p); + else if (p->abspid() == 13) _muons.push_back(p); + else if (p->abspid() == 15) _taus.push_back(p); + } } - else if (p->abspid() == 15) - { - _taus.push_back(p); - _ctaus.push_back(p); - } - else if (p->abspid() == 12 || p->abspid() == 14 || p->abspid() == 16 || - p->pid() == 1000022 || p->pid() == 1000039 || - in_range(p->abspid(), 50, 60)) //< invert definition to specify all *visibles*? - { - _invisibles.push_back(p); - _cinvisibles.push_back(p); - } - else - delete p; + + // Sort the collections + if (ptsort) _sort_particles(); } /// Add a collection of final state particles to the event /// - /// Supplied particles should be new'd, and Event will take ownership. - void add_particles(const std::vector& ps) { - for (size_t i = 0; i < ps.size(); ++i) add_particle(ps[i]); - } - - - /// Get all final state particles - /// @todo Note the return by value: it's not efficient yet! - /// @note Overlap of taus and e/mu - std::vector particles() const { - // Add together all the vectors of the different particle types - std::vector rtn; - // rtn.reserve(visible_particles().size() + _cinvisibles.size()); - rtn.reserve(_cphotons.size() + _celectrons.size() + _cmuons.size() + _ctaus.size() + _cinvisibles.size()); - #define APPEND_VEC(vec) rtn.insert(rtn.end(), vec.begin(), vec.end()) - // APPEND_VEC(visible_particles()); - APPEND_VEC(_cphotons); - APPEND_VEC(_celectrons); - APPEND_VEC(_cmuons); - APPEND_VEC(_ctaus); - APPEND_VEC(_cinvisibles); - #undef APPEND_VEC - return rtn; - /// @todo Or use Boost range join to iterate over separate containers transparently... I like this - /// @todo Cache, or otherwise think about efficiency since this gets called by the destructor - } - - - /// Get visible state particles - /// @todo Note the return by value: it's not efficient yet! - /// @note Overlap of taus and e/mu - std::vector visible_particles() const { - // Add together all the vectors of the different particle types - std::vector rtn; - rtn.reserve(_cphotons.size() + _celectrons.size() + _cmuons.size() + _ctaus.size()); - #define APPEND_VEC(vec) rtn.insert(rtn.end(), vec.begin(), vec.end() ) - APPEND_VEC(_cphotons); - APPEND_VEC(_celectrons); - APPEND_VEC(_cmuons); - APPEND_VEC(_ctaus); - #undef APPEND_VEC - return rtn; - /// @todo Add together all the vectors of the different particle types - /// @todo Or use Boost range join to iterate over separate containers transparently... I like this + /// @warning Supplied particles should be new'd, and Event will take ownership. + void add_particles(const std::vector& ps, bool ptsort=true) { + // Add each particle, without sorting each time + for (const Particle* p : ps) add_particle(new Particle(*p), false); + + // Finally sort the collections, once all new particles are in place + if (ptsort) _sort_particles(); + } + + + /// A mostly-internal function to sort the particle-vector caches + void _sort_particles() { + std::sort(_allparticles.begin(), _allparticles.end(), _cmpPtDescPtr); + std::sort(_invisibles.begin(), _invisibles.end(), _cmpPtDescPtr); + std::sort(_photons.begin(), _photons.end(), _cmpPtDescPtr); + std::sort(_electrons.begin(), _electrons.end(), _cmpPtDescPtr); + std::sort(_muons.begin(), _muons.end(), _cmpPtDescPtr); + std::sort(_taus.begin(), _taus.end(), _cmpPtDescPtr); + } + + + /// @brief Get all final state particles + /// + /// @note Small overlap of taus and e/mu? + const std::vector& particles() const { + return _allparticles; } - /// Get invisible final state particles + /// @brief Get visible state particles + /// + /// @note Small overlap of taus and e/mu? + const std::vector& visible_particles() const { + return _visibles; + } + + + /// @brief Get invisible final state particles + /// + /// @note Both prompt and non-prompt... correct? const std::vector& invisible_particles() const { - return _cinvisibles; + return _invisibles; } /// Get invisible final state particles (non-const) std::vector& invisible_particles() { - return _invisibles; + return mkunconst(_invisibles); } /// Get prompt electrons const std::vector& electrons() const { - return _celectrons; + return _electrons; } /// Get prompt electrons (non-const) std::vector& electrons() { - return _electrons; + return mkunconst(_electrons); } /// Get prompt muons const std::vector& muons() const { - return _cmuons; + return _muons; } /// Get prompt muons (non-const) std::vector& muons() { - return _muons; + return mkunconst(_muons); } /// Get prompt (hadronic) taus const std::vector& taus() const { - return _ctaus; + return _taus; } /// Get prompt (hadronic) taus (non-const) std::vector& taus() { - return _taus; + return mkunconst(_taus); } /// Get prompt photons const std::vector& photons() const { - return _cphotons; + return _photons; } /// Get prompt photons (non-const) std::vector& photons() { - return _photons; + return mkunconst(_photons); } + /// @} + /// @name Jets - /// @todo Why the new'ing? Can we use references more? - //@{ + /// + /// @{ + + // Implementation function, to avoid const/non-const duplication below + std::vector& _get_jets(const std::string& key) const { + std::vector& rtn = _jets[key]; + // std::sort(rtn.begin(), rtn.end(), _cmpPtDescPtr); //< not thread-safe! + return rtn; + } - /// @brief Get anti-kT 0.4 jets (not including charged leptons or photons) - const std::vector& jets() const { - std::sort(_cjets.begin(), _cjets.end(), _cmpPtDesc); - return _cjets; + /// @brief Get a jet collection (not including charged leptons or photons) + const std::vector& jets(const std::string& key="CANONICAL") const { + return _get_jets(key); } - /// @brief Get anti-kT 0.4 jets (not including charged leptons or photons) (non-const) - std::vector& jets() { - std::sort(_jets.begin(), _jets.end(), _cmpPtDesc); - return _jets; + /// @brief Get a jet collection (not including charged leptons or photons) (non-const) + std::vector& jets(const std::string& key="CANONICAL") { + return mkunconst(_get_jets(key)); } + /// @brief Set the jets collection /// - /// The Jets should be new'd; Event will take ownership. - void set_jets(const std::vector& jets) { - _jets = jets; - _cjets.clear(); - for (Jet* j : jets ) _cjets.push_back(j); + /// @warning The Jets should be new'd; Event will take ownership. + /// @todo "Lock" at some point so that jet finding etc. only get done once + void set_jets(const std::vector& jets, const std::string& key="CANONICAL") { + _jets[key] = jets; + std::sort(_jets[key].begin(), _jets[key].end(), _cmpPtDescPtr); } + // /// @brief Set the jets collection (non-const input) + // void set_jets(const std::vector& jets, const std::string& key="CANONICAL") { + // set_jets(mkconst(jets), key); + // } + /// @brief Add a jet to the jets collection /// - /// The Jet should be new'd; Event will take ownership. - void add_jet(Jet* j) { - _jets.push_back(j); - _cjets.push_back(j); + /// @warning The Jet should be new'd; Event will take ownership. + /// @todo "Lock" at some point so that jet finding etc. only get done once + void add_jet(const Jet* j, const std::string& key="CANONICAL") { + _jets[key].push_back(j); + std::sort(_jets[key].begin(), _jets[key].end(), _cmpPtDescPtr); } + // /// @brief Add a jet to the jets collection (non-const input) + // void add_jet(Jet* j, const std::string& key="CANONICAL") { + // add_jet(mkconst(j), key); + // } - //@} + /// @} - /// @name Missing energy - //@{ + /// @name Missing momentum + /// @{ - /// @brief Get the missing momentum vector - /// - /// Not _necessarily_ the sum over momenta of final state invisibles + /// Get the missing momentum vector + /// @note Not _necessarily_ the sum over momenta of final state invisibles const P4& missingmom() const { return _pmiss; } - /// @brief Set the missing momentum vector - /// - /// Not _necessarily_ the sum over momenta of final state invisibles + /// Set the missing momentum vector + /// @todo Not _necessarily_ the sum over momenta of final state invisibles void set_missingmom(const P4& pmiss) { _pmiss = pmiss; } @@ -421,7 +415,7 @@ namespace HEPUtils { return missingmom().pT(); } - //@} + /// @} }; diff --git a/contrib/heputils/include/HEPUtils/FastJet.h b/contrib/heputils/include/HEPUtils/FastJet.h index 8f9ec119ec..b73f180bd9 100644 --- a/contrib/heputils/include/HEPUtils/FastJet.h +++ b/contrib/heputils/include/HEPUtils/FastJet.h @@ -1,7 +1,7 @@ // -*- C++ -*- // // This file is part of HEPUtils -- https://bitbucket.org/andybuckley/heputils -// Copyright (C) 2013-2018 Andy Buckley +// Copyright (C) 2013-2021 Andy Buckley // // Embedding of HEPUtils code in other projects is permitted provided this // notice is retained and the HEPUtils namespace and include path are changed. @@ -32,18 +32,11 @@ #endif #endif -//#define JETCLUSTER_DEBUG - -/*inline bool compare_particles_by_pz(FJNS::PseudoJet jet1, FJNS::PseudoJet jet2) -{ - return (jet1.pz() > jet2.pz()); -}*/ - namespace HEPUtils { - /// @name Converters between HEPUtils and FastJet momentum types - //@{ + /// @defgroup fastjet_cnv Converters between HEPUtils and FastJet momentum types + /// @{ /// @todo Enable... conditionally on FJ version? // /// For attaching the GenParticle provenance info to a PseudoJet @@ -71,43 +64,22 @@ namespace HEPUtils { return P4::mkXYZM(p.px(), p.py(), p.pz(), (m >= 0) ? m : 0); } - //@} + /// @} - /// @name Jet builders - //@{ + /// @defgroup fastjet_mk Jet builders + /// @{ /// Construct pT-sorted jets using the @a alg measure with jet @a R parameter, and min pT @a ptmin (in MeV) inline std::vector get_jets(const std::vector& particles, double R, double ptmin, - FJNS::JetAlgorithm alg=FJNS::antikt_algorithm) - { - + FJNS::JetAlgorithm alg=FJNS::antikt_algorithm) { const FJNS::JetDefinition jet_def(alg, R); - - //std::sort(particles.begin(), particles.end(), compare_particles_by_pz); - - //int TP_TEMP_COUNTER = 0; - //for (auto particle : particles) - //{ - //std::cout << "Particle Number: " << TP_TEMP_COUNTER++ << "; Pz: " << particle.pz() << "; Px: " << particle.px() < sorting (cf. std::less) - inline bool _cmpPtDesc(const Jet* a, const Jet* b) { - return a->pT2() >= b->pT2(); + /// @defgroup jet_const Jet constness conversion + /// @{ + + /// Convenience Jet cast to const + inline const Jet* mkconst(Jet* jet) { + return const_cast(jet); + } + + /// Convenience Jet cast to non-const + inline Jet* mkunconst(const Jet* cjet) { + return const_cast(cjet); } + /// Get a reference to a vector of Jets, with each member const + inline std::vector& mkconst(const std::vector& jets) { + return * (std::vector*)(void*) (&jets); + } + + /// Get a reference to a vector of Jets, with each member non-const + inline std::vector& mkunconst(const std::vector& cjets) { + return * (std::vector*) (void*) (&cjets); + } + + /// @} + + + + // /// Function/functor for container sorting (cf. std::less) + // inline bool _cmpPtDescJet(const Jet* a, const Jet* b) { + // return a->pT2() >= b->pT2(); + // } + } diff --git a/contrib/heputils/include/HEPUtils/MathUtils.h b/contrib/heputils/include/HEPUtils/MathUtils.h index 1a10e76745..30cafab0a7 100644 --- a/contrib/heputils/include/HEPUtils/MathUtils.h +++ b/contrib/heputils/include/HEPUtils/MathUtils.h @@ -1,7 +1,7 @@ // -*- C++ -*- // // This file is part of HEPUtils -- https://bitbucket.org/andybuckley/heputils -// Copyright (C) 2013-2018 Andy Buckley +// Copyright (C) 2013-2021 Andy Buckley // // Embedding of HEPUtils code in other projects is permitted provided this // notice is retained and the HEPUtils namespace and include path are changed. @@ -25,8 +25,11 @@ namespace HEPUtils { - /// @name Numerical helpers - //@{ + /// @defgroup mathutils Mathematical utilities + /// @{ + + /// @defgroup mathutils_numhelp Numerical helpers + /// @{ /// Convenience function for getting the sign of a number template @@ -58,11 +61,11 @@ namespace HEPUtils { return sqrt(a*a + b*b + c*c); } - //@} + /// @} - /// @name Range helpers - //@{ + /// @defgroup mathutils_rangehelp Range helpers + /// @{ /// @brief Boolean function to determine if @a value is within the given range /// @@ -128,11 +131,11 @@ namespace HEPUtils { return rtn; } - //@} + /// @} - /// @name Physics maths utils - //@{ + /// @defgroup mathutils_physics Physics maths utils + /// @{ inline double deltaphi(double a, double b) { double rtn = a - b; @@ -147,17 +150,17 @@ namespace HEPUtils { return rtn; } - //@} + /// @} - /// @name Random numbers and sampling - //@{ + /// @defgroup mathutils_random Random numbers and sampling + /// @{ inline double rand01() { return rand() / (double)RAND_MAX; } - //@} + /// @} } diff --git a/contrib/heputils/include/HEPUtils/Particle.h b/contrib/heputils/include/HEPUtils/Particle.h index fd7c412a67..6bbbe45716 100644 --- a/contrib/heputils/include/HEPUtils/Particle.h +++ b/contrib/heputils/include/HEPUtils/Particle.h @@ -1,7 +1,7 @@ // -*- C++ -*- // // This file is part of HEPUtils -- https://bitbucket.org/andybuckley/heputils -// Copyright (C) 2013-2018 Andy Buckley +// Copyright (C) 2013-2021 Andy Buckley // // Embedding of HEPUtils code in other projects is permitted provided this // notice is retained and the HEPUtils namespace and include path are changed. @@ -14,14 +14,14 @@ namespace HEPUtils { - /// Simple particle class, encapsulating a momentum 4-vector and adding some extra ID info - /// @todo Derive from a PhysObj base class to centralise the momentum handling - /// @todo Provide cast operators to P4 and P4* + /// @brief Simple particle class, encapsulating a momentum 4-vector and adding some extra ID info + /// + /// @todo Derive from a PhysObj base class to centralise the momentum handling? class Particle { private: /// @name Storage - //@{ + /// @{ /// Momentum vector P4 _p4; /// PDG ID code @@ -30,13 +30,13 @@ namespace HEPUtils { bool _prompt; // /// Isolation value // double _isol4; - //@} + /// @} public: /// @name Constructors - //@{ + /// @{ /// Default constructor Particle() @@ -72,23 +72,23 @@ namespace HEPUtils { return *this; } - //@} + /// @} /// @name Implicit casts - //@{ + /// @{ operator const P4& () const { return mom(); } operator const P4* () const { return &mom(); } - //@} + /// @} /// @name Momentum /// /// Access to the P4 object, plus convenience mapping of a few popular properties - //@{ + /// @{ /// Get the 4 vector const P4& mom() const { return _p4; } @@ -126,22 +126,22 @@ namespace HEPUtils { /// Get the squared transverse momentum double pT() const { return mom().pT(); } - //@} + /// @} /// @name Promptness - //@{ + /// @{ /// Is this particle connected to the hard process or from a hadron/tau decay? bool is_prompt() const { return _prompt; } /// Set promptness void set_prompt(bool isprompt=true) { _prompt = isprompt; } - //@} + /// @} /// @name Particle ID - //@{ + /// @{ /// Get PDG ID code int pid() const { return _pdgId; } @@ -150,20 +150,64 @@ namespace HEPUtils { /// Set PDG ID code void set_pid(int pid) { _pdgId = pid; } - //@} + /// Is this particle usually visible in a detector? + bool is_visible() { + if (abspid() == 12 || abspid() == 14 || abspid() == 16) return false; + if (pid() == 1000022 || pid() == 1000039) return false; + if (in_range(abspid(), 50, 60)) return false; //< abspid zealousness since some -ve DM PIDs seen + return true; + } + + /// Is this particle usually invisible in a detector? + bool is_invisible() { + return !is_visible(); + } + /// @} // /// @name Isolation of particle - // //@{ + // /// @{ // /// Get isolation // double isol() const { return _isol4;} // void set_isol(double isol) { _isol4 = isol;} - // //@} + // /// @} }; + /// @defgroup particle_const Particle constness conversion + /// @{ + + /// Convenience Particle cast to const + inline const Particle* mkconst(Particle* particle) { + return const_cast(particle); + } + + /// Convenience Particle cast to non-const + inline Particle* mkunconst(const Particle* cparticle) { + return const_cast(cparticle); + } + + /// Get a reference to a vector of Particles, with each member const + inline std::vector& mkconst(const std::vector& particles) { + return * (std::vector*)(void*) (&particles); + } + + /// Get a reference to a vector of Particles, with each member non-const + inline std::vector& mkunconst(const std::vector& cparticles) { + return * (std::vector*) (void*) (&cparticles); + } + + /// @} + + + // /// Function/functor for container sorting (cf. std::less) + // inline bool _cmpPtDescParticle(const Particle* a, const Particle* b) { + // return a->pT2() >= b->pT2(); + // } + + } diff --git a/contrib/heputils/include/HEPUtils/PhaseSpace.h b/contrib/heputils/include/HEPUtils/PhaseSpace.h index c68888a43f..31596e7aef 100644 --- a/contrib/heputils/include/HEPUtils/PhaseSpace.h +++ b/contrib/heputils/include/HEPUtils/PhaseSpace.h @@ -1,7 +1,7 @@ // -*- C++ -*- // // This file is part of HEPUtils -- https://bitbucket.org/andybuckley/heputils -// Copyright (C) 2013-2018 Andy Buckley +// Copyright (C) 2013-2021 Andy Buckley // // Embedding of HEPUtils code in other projects is permitted provided this // notice is retained and the HEPUtils namespace and include path are changed. @@ -20,6 +20,10 @@ namespace HEPUtils { + /// @defgroup phasespace Phase-space helpers + /// @{ + + /// Integral of acos(x) constexpr double iacos(double x) { return sqrt(1 - x*x) - x*acos(x); } @@ -74,5 +78,7 @@ namespace HEPUtils { return psPhiIntegral; } + /// @} + } diff --git a/contrib/heputils/include/HEPUtils/Utils.h b/contrib/heputils/include/HEPUtils/Utils.h index b40f5497b2..fe5c54aaa8 100644 --- a/contrib/heputils/include/HEPUtils/Utils.h +++ b/contrib/heputils/include/HEPUtils/Utils.h @@ -1,28 +1,26 @@ // -*- C++ -*- // // This file is part of HEPUtils -- https://bitbucket.org/andybuckley/heputils -// Copyright (C) 2013-2018 Andy Buckley +// Copyright (C) 2013-2021 Andy Buckley // // Embedding of HEPUtils code in other projects is permitted provided this // notice is retained and the HEPUtils namespace and include path are changed. // #pragma once -#if defined(__cplusplus) && __cplusplus < 201103L -// #define XSTR(x) STR(x) -// #define STR(x) #x -// #pragma message STR(__cplusplus) -#error "This library needs at least a C++11 compliant compiler" +#if __cplusplus < 201103L +#pragma message "This library needs at least a C++11 compliant compiler" #endif + /// @file Utility functions /// @author Andy Buckley namespace HEPUtils { - /// @name Container utils - //@{ + /// @defgroup utils_container Container utils + /// @{ /// Return true if f(x) is true for any x in container c, otherwise false. template @@ -52,7 +50,7 @@ namespace HEPUtils { } - //@} + /// @} } diff --git a/contrib/heputils/include/HEPUtils/Vectors.h b/contrib/heputils/include/HEPUtils/Vectors.h index 77f32701f2..e4264f32a2 100644 --- a/contrib/heputils/include/HEPUtils/Vectors.h +++ b/contrib/heputils/include/HEPUtils/Vectors.h @@ -1,7 +1,7 @@ // -*- C++ -*- // // This file is part of HEPUtils -- https://bitbucket.org/andybuckley/heputils -// Copyright (C) 2013-2018 Andy Buckley +// Copyright (C) 2013-2021 Andy Buckley // // Embedding of HEPUtils code in other projects is permitted provided this // notice is retained and the HEPUtils namespace and include path are changed. @@ -371,9 +371,9 @@ namespace HEPUtils { double phi() const { if (rho2() == 0) return 0; else return atan2(py(),px()); } /// Get the spatial theta double theta() const { if (p2() == 0) return 0; else if (pz() == 0) return M_PI; else return atan2(rho(),pz()); } - /// Get the spatial vector pseudorapidity - double eta() const { return -log(tan( 0.5 * theta() )); } //< Optimise with a trig reln on tan(x/2) to avoid tan(atan(..)/2)? - /// Get the spatial vector absolute pseudorapidity + /// Get the spatial-vector pseudorapidity + double eta() const { return std::copysign(log((p() + fabs(pz())) / pT()), pz()); } + /// Get the spatial-vector absolute pseudorapidity double abseta() const { return fabs(eta()); } /// Get the 4-momentum rapidity double rap() const { return 0.5 * (E() + pz()) / (E() - pz()); } @@ -505,4 +505,17 @@ namespace HEPUtils { //@} + /// Function/functor for container sorting (cf. std::less) + template + inline bool _cmpPtDesc(const T& a, const T& b) { + return a.pT2() >= b.pT2(); + } + + /// Function/functor for container sorting (cf. std::less) + template + inline bool _cmpPtDescPtr(const T* a, const T* b) { + return _cmpPtDesc(*a, *b); + } + + } From 0ae0da9e777677674dad485b842fedf2dddd2009 Mon Sep 17 00:00:00 2001 From: Andy Buckley Date: Fri, 21 May 2021 13:17:57 +0100 Subject: [PATCH 02/27] Fix compiler warnings (one potentially significant) in CB analyses --- .../src/analyses/Analysis_CMS_13TeV_1LEPStop_36invfb.cpp | 4 ++-- .../src/analyses/Analysis_CMS_13TeV_2SSLEP_Stop_137invfb.cpp | 2 +- .../src/analyses/Analysis_CMS_13TeV_2SSLEP_Stop_36invfb.cpp | 2 +- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/ColliderBit/src/analyses/Analysis_CMS_13TeV_1LEPStop_36invfb.cpp b/ColliderBit/src/analyses/Analysis_CMS_13TeV_1LEPStop_36invfb.cpp index 265e290fe5..eea60e0dac 100644 --- a/ColliderBit/src/analyses/Analysis_CMS_13TeV_1LEPStop_36invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_CMS_13TeV_1LEPStop_36invfb.cpp @@ -294,7 +294,7 @@ namespace Gambit { // Mlb double deltaRlb=9999.; - double Mlb; + double Mlb = -1; ///< what if no real value is found??? for (const HEPUtils::Jet* bj :mediumbJets) { if (deltaRlb > bj->mom().deltaR_eta(Leptons.at(0)->mom())){ deltaRlb = bj->mom().deltaR_eta(Leptons.at(0)->mom()); @@ -395,7 +395,7 @@ namespace Gambit { if (baselineJets.size()>=4 and tmod>10 and Mlb<=175 and met>=450) _counters.at("aggregateSR2").add_event(event); if (baselineJets.size()>=4 and tmod<=0 and Mlb> 175 and met>=450) _counters.at("aggregateSR3").add_event(event); if (baselineJets.size()>=4 and tmod> 0 and Mlb> 175 and met>=450) _counters.at("aggregateSR4").add_event(event); - if(baselineJets.size()>=5 and leadjet_nob and deltaPhi_j12 >0.5 and Leptons.at(0)->pT() < 150 and Leptons.at(0)->mom().deltaPhi(ptot)<2. ){ + if (baselineJets.size()>=5 and leadjet_nob and deltaPhi_j12 >0.5 and Leptons.at(0)->pT() < 150 and Leptons.at(0)->mom().deltaPhi(ptot)<2. ){ if( met>=450 ) _counters.at("aggregateSR5").add_event(event); } return; diff --git a/ColliderBit/src/analyses/Analysis_CMS_13TeV_2SSLEP_Stop_137invfb.cpp b/ColliderBit/src/analyses/Analysis_CMS_13TeV_2SSLEP_Stop_137invfb.cpp index 1f11e256f2..06685b92c6 100644 --- a/ColliderBit/src/analyses/Analysis_CMS_13TeV_2SSLEP_Stop_137invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_CMS_13TeV_2SSLEP_Stop_137invfb.cpp @@ -272,8 +272,8 @@ namespace Gambit { muons.push_back(muon); } - double HT; // Jets + double HT = 0; vector candJets; for (const HEPUtils::Jet* jet : event->jets()) { if (jet->pT() > 25. && fabs(jet->eta()) < 2.4){ diff --git a/ColliderBit/src/analyses/Analysis_CMS_13TeV_2SSLEP_Stop_36invfb.cpp b/ColliderBit/src/analyses/Analysis_CMS_13TeV_2SSLEP_Stop_36invfb.cpp index 631c39b756..a1797d79a8 100644 --- a/ColliderBit/src/analyses/Analysis_CMS_13TeV_2SSLEP_Stop_36invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_CMS_13TeV_2SSLEP_Stop_36invfb.cpp @@ -290,8 +290,8 @@ namespace Gambit { muons.push_back(muon); } - double HT; // Jets + double HT = 0; vector candJets; for (const HEPUtils::Jet* jet : event->jets()) { if (jet->pT() > 25. && fabs(jet->eta()) < 2.4){ From 1e88674b44bbcaf70f51512de870c162bf5ec6e1 Mon Sep 17 00:00:00 2001 From: Andy Buckley Date: Tue, 17 Aug 2021 19:56:02 +0100 Subject: [PATCH 03/27] Massively reduce the scheduled CI frequency (to one per day at 5 am) --- .github/workflows/ci.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index abc6a9e744..d790c75839 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -2,11 +2,11 @@ name: Gambit CI on: push: - branches: [ master, gitlab-ci ] + branches: [ master ] pull_request: branches: [ master ] schedule: - - cron: '* * * * 1' + - cron: '0 5 * * *' jobs: gambit_build: From 7e37cf54336ae68222fe8ca800b53c1559dbab4e Mon Sep 17 00:00:00 2001 From: Andy Buckley Date: Tue, 17 Aug 2021 20:00:57 +0100 Subject: [PATCH 04/27] Sync CI with master --- .github/workflows/ci.yml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index d790c75839..dc284d4a2a 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -2,7 +2,7 @@ name: Gambit CI on: push: - branches: [ master ] + branches: [ master, ci-* ] pull_request: branches: [ master ] schedule: @@ -10,7 +10,7 @@ on: jobs: gambit_build: - runs-on: ubuntu-latest + runs-on: [docker, self-hosted] strategy: fail-fast: false matrix: @@ -44,7 +44,7 @@ jobs: echo "export PYTHON_INCLUDE_DIR=$PYTHON_INCLUDE_DIR" >> buildenv.sh echo "export LD_LIBRARY_PATH=/usr/lib/x86_64-linux-gnu:$LD_LIBRARY_PATH" >> buildenv.sh cat buildenv.sh - pip install pyyaml pybind11 h5py + pip install pyyaml pybind11 h5py scipy - name: Configure with cmake run: | cd BUILD/ && . buildenv.sh From bc7d34fc118cb8f738ffd6c31fe8be4761f75c50 Mon Sep 17 00:00:00 2001 From: Andy Buckley Date: Thu, 3 Nov 2022 21:24:54 +0000 Subject: [PATCH 05/27] HEPUtils update (to prototype version) to add ClusterSequence access --- contrib/heputils/include/HEPUtils/BinnedFn.h | 4 +-- contrib/heputils/include/HEPUtils/Event.h | 27 ++++++-------- contrib/heputils/include/HEPUtils/FastJet.h | 4 +-- contrib/heputils/include/HEPUtils/Jet.h | 36 +++++++++++++++++-- contrib/heputils/include/HEPUtils/MathUtils.h | 4 +-- contrib/heputils/include/HEPUtils/Particle.h | 4 +-- .../heputils/include/HEPUtils/PhaseSpace.h | 6 ++-- contrib/heputils/include/HEPUtils/Utils.h | 4 +-- contrib/heputils/include/HEPUtils/Vectors.h | 8 ++--- 9 files changed, 61 insertions(+), 36 deletions(-) diff --git a/contrib/heputils/include/HEPUtils/BinnedFn.h b/contrib/heputils/include/HEPUtils/BinnedFn.h index 45f4b513f2..72558e1051 100644 --- a/contrib/heputils/include/HEPUtils/BinnedFn.h +++ b/contrib/heputils/include/HEPUtils/BinnedFn.h @@ -1,7 +1,7 @@ // -*- C++ -*- // -// This file is part of HEPUtils -- https://bitbucket.org/andybuckley/heputils -// Copyright (C) 2013-2021 Andy Buckley +// This file is part of HEPUtils -- https://gitlab.com/hepcedar/heputils/ +// Copyright (C) 2013-2022 Andy Buckley // // Embedding of HEPUtils code in other projects is permitted provided this // notice is retained and the HEPUtils namespace and include path are changed. diff --git a/contrib/heputils/include/HEPUtils/Event.h b/contrib/heputils/include/HEPUtils/Event.h index f3bd45a978..6b556984d9 100644 --- a/contrib/heputils/include/HEPUtils/Event.h +++ b/contrib/heputils/include/HEPUtils/Event.h @@ -1,7 +1,7 @@ // -*- C++ -*- // -// This file is part of HEPUtils -- https://bitbucket.org/andybuckley/heputils -// Copyright (C) 2013-2021 Andy Buckley +// This file is part of HEPUtils -- https://gitlab.com/hepcedar/heputils/ +// Copyright (C) 2013-2022 Andy Buckley // // Embedding of HEPUtils code in other projects is permitted provided this // notice is retained and the HEPUtils namespace and include path are changed. @@ -242,26 +242,11 @@ namespace HEPUtils { else if (p->abspid() == 11) _electrons.push_back(p); else if (p->abspid() == 13) _muons.push_back(p); else if (p->abspid() == 15) _taus.push_back(p); - // } else { - // delete p; } } // Sort the collections if (ptsort) _sort_particles(); - - // else if (p->abspid() == 15) - // { - // _taus.push_back(p); - // _ctaus.push_back(p); - // } - // else if (p->abspid() == 12 || p->abspid() == 14 || p->abspid() == 16 || - // p->pid() == 1000022 || p->pid() == 1000039 || - // in_range(p->abspid(), 50, 60)) //< invert definition to specify all *visibles*? - // { - // _invisibles.push_back(p); - // _cinvisibles.push_back(p); - // } } @@ -407,6 +392,14 @@ namespace HEPUtils { // add_jet(mkconst(j), key); // } + /// @brief Access the jets' ClusterSequence object if possible (can be null) + /// + /// Optional template arg can be used to cast to a specific derived CS type if wanted. + template + const CS* clusterseq(const std::string& key="CANONICAL") { + return dynamic_cast(jets(key).front().clusterseq()); + } + /// @} diff --git a/contrib/heputils/include/HEPUtils/FastJet.h b/contrib/heputils/include/HEPUtils/FastJet.h index b73f180bd9..d902feaa73 100644 --- a/contrib/heputils/include/HEPUtils/FastJet.h +++ b/contrib/heputils/include/HEPUtils/FastJet.h @@ -1,7 +1,7 @@ // -*- C++ -*- // -// This file is part of HEPUtils -- https://bitbucket.org/andybuckley/heputils -// Copyright (C) 2013-2021 Andy Buckley +// This file is part of HEPUtils -- https://gitlab.com/hepcedar/heputils/ +// Copyright (C) 2013-2022 Andy Buckley // // Embedding of HEPUtils code in other projects is permitted provided this // notice is retained and the HEPUtils namespace and include path are changed. diff --git a/contrib/heputils/include/HEPUtils/Jet.h b/contrib/heputils/include/HEPUtils/Jet.h index 5282956dcb..40af96062c 100644 --- a/contrib/heputils/include/HEPUtils/Jet.h +++ b/contrib/heputils/include/HEPUtils/Jet.h @@ -1,13 +1,14 @@ // -*- C++ -*- // -// This file is part of HEPUtils -- https://bitbucket.org/andybuckley/heputils -// Copyright (C) 2013-2017 Andy Buckley +// This file is part of HEPUtils -- https://gitlab.com/hepcedar/heputils/ +// Copyright (C) 2013-2022 Andy Buckley // // Embedding of HEPUtils code in other projects is permitted provided this // notice is retained and the HEPUtils namespace and include path are changed. // #pragma once +#include "HEPUtils/FastJet.h" #include "HEPUtils/MathUtils.h" #include "HEPUtils/Vectors.h" @@ -15,6 +16,7 @@ namespace HEPUtils { /// Simple jet class, encapsulating a momentum 4-vector and with some extra b-tag info + /// /// @todo Derive from a PhysObj base class to centralise the momentum handling class Jet { @@ -24,6 +26,9 @@ namespace HEPUtils { P4 _p4; /// B and C tags bool _isB, _isC; + /// Optional FastJet PJ (contains link to ClusterSeq) + /// @todo Use std::optional when C++17 allowed + FJNS::PseudoJet _pj; /// @} @@ -40,6 +45,10 @@ namespace HEPUtils { Jet(double px, double py, double pz, double E, bool isB=false, bool isC=false) : _p4(px, py, pz, E), _isB(isB), _isC(isC) { } + /// "PseudoJet" constructor + Jet(const FJNS::PseudoJet& pj, bool isB=false, bool isC=false) + : _p4(mk_p4(pj)), _isB(isB), _isC(isC) { } + /// @} @@ -113,6 +122,29 @@ namespace HEPUtils { /// @} + + /// @name FastJet information + /// @{ + + /// Get the contained PseudoJet object (const) + const FJNS::PseudoJet& pseudojet() const { return _pj; } + + /// Get the contained PseudoJet object + FJNS::PseudoJet& pseudojet() { return _pj; } + + /// Set the contained PseudoJet object + void set_pseudojet(const FJNS::PseudoJet& pj) { _pj = pj; } + + /// @brief Access the ClusterSequence object if possible (can be null) + /// + /// Optional template arg can be used to cast to a specific derived CS type if wanted. + template + const CS* clusterseq() { + return dynamic_cast(_pj.associated_cs()); + } + + /// @} + }; diff --git a/contrib/heputils/include/HEPUtils/MathUtils.h b/contrib/heputils/include/HEPUtils/MathUtils.h index 30cafab0a7..850285ff46 100644 --- a/contrib/heputils/include/HEPUtils/MathUtils.h +++ b/contrib/heputils/include/HEPUtils/MathUtils.h @@ -1,7 +1,7 @@ // -*- C++ -*- // -// This file is part of HEPUtils -- https://bitbucket.org/andybuckley/heputils -// Copyright (C) 2013-2021 Andy Buckley +// This file is part of HEPUtils -- https://gitlab.com/hepcedar/heputils/ +// Copyright (C) 2013-2022 Andy Buckley // // Embedding of HEPUtils code in other projects is permitted provided this // notice is retained and the HEPUtils namespace and include path are changed. diff --git a/contrib/heputils/include/HEPUtils/Particle.h b/contrib/heputils/include/HEPUtils/Particle.h index 6bbbe45716..3e7d6d8c0a 100644 --- a/contrib/heputils/include/HEPUtils/Particle.h +++ b/contrib/heputils/include/HEPUtils/Particle.h @@ -1,7 +1,7 @@ // -*- C++ -*- // -// This file is part of HEPUtils -- https://bitbucket.org/andybuckley/heputils -// Copyright (C) 2013-2021 Andy Buckley +// This file is part of HEPUtils -- https://gitlab.com/hepcedar/heputils/ +// Copyright (C) 2013-2022 Andy Buckley // // Embedding of HEPUtils code in other projects is permitted provided this // notice is retained and the HEPUtils namespace and include path are changed. diff --git a/contrib/heputils/include/HEPUtils/PhaseSpace.h b/contrib/heputils/include/HEPUtils/PhaseSpace.h index 31596e7aef..a1501facf2 100644 --- a/contrib/heputils/include/HEPUtils/PhaseSpace.h +++ b/contrib/heputils/include/HEPUtils/PhaseSpace.h @@ -1,7 +1,7 @@ // -*- C++ -*- // -// This file is part of HEPUtils -- https://bitbucket.org/andybuckley/heputils -// Copyright (C) 2013-2021 Andy Buckley +// This file is part of HEPUtils -- https://gitlab.com/hepcedar/heputils/ +// Copyright (C) 2013-2022 Andy Buckley // // Embedding of HEPUtils code in other projects is permitted provided this // notice is retained and the HEPUtils namespace and include path are changed. @@ -25,7 +25,7 @@ namespace HEPUtils { /// Integral of acos(x) - constexpr double iacos(double x) { return sqrt(1 - x*x) - x*acos(x); } + constexpr double iacos(double x) { return x*acos(x) - sqrt(1 - x*x); } /// Inverse secant constexpr double asec(double x) { return acos(1/x); } diff --git a/contrib/heputils/include/HEPUtils/Utils.h b/contrib/heputils/include/HEPUtils/Utils.h index fe5c54aaa8..c6f762e466 100644 --- a/contrib/heputils/include/HEPUtils/Utils.h +++ b/contrib/heputils/include/HEPUtils/Utils.h @@ -1,7 +1,7 @@ // -*- C++ -*- // -// This file is part of HEPUtils -- https://bitbucket.org/andybuckley/heputils -// Copyright (C) 2013-2021 Andy Buckley +// This file is part of HEPUtils -- https://gitlab.com/hepcedar/heputils/ +// Copyright (C) 2013-2022 Andy Buckley // // Embedding of HEPUtils code in other projects is permitted provided this // notice is retained and the HEPUtils namespace and include path are changed. diff --git a/contrib/heputils/include/HEPUtils/Vectors.h b/contrib/heputils/include/HEPUtils/Vectors.h index e4264f32a2..58db41bc02 100644 --- a/contrib/heputils/include/HEPUtils/Vectors.h +++ b/contrib/heputils/include/HEPUtils/Vectors.h @@ -1,7 +1,7 @@ // -*- C++ -*- // -// This file is part of HEPUtils -- https://bitbucket.org/andybuckley/heputils -// Copyright (C) 2013-2021 Andy Buckley +// This file is part of HEPUtils -- https://gitlab.com/hepcedar/heputils/ +// Copyright (C) 2013-2022 Andy Buckley // // Embedding of HEPUtils code in other projects is permitted provided this // notice is retained and the HEPUtils namespace and include path are changed. @@ -370,13 +370,13 @@ namespace HEPUtils { /// Get the spatial phi double phi() const { if (rho2() == 0) return 0; else return atan2(py(),px()); } /// Get the spatial theta - double theta() const { if (p2() == 0) return 0; else if (pz() == 0) return M_PI; else return atan2(rho(),pz()); } + double theta() const { if (p2() == 0) return 0; else if (pz() == 0) return 0.5 * M_PI; else return atan2(rho(),pz()); } /// Get the spatial-vector pseudorapidity double eta() const { return std::copysign(log((p() + fabs(pz())) / pT()), pz()); } /// Get the spatial-vector absolute pseudorapidity double abseta() const { return fabs(eta()); } /// Get the 4-momentum rapidity - double rap() const { return 0.5 * (E() + pz()) / (E() - pz()); } + double rap() const { return 0.5 * log((E() + pz()) / (E() - pz())); } /// Get the 4-momentum absolute rapidity double absrap() const { return fabs(rap()); } From bf4c0a2e9a5ce18a09d7cee89fbb9294cd9b8a47 Mon Sep 17 00:00:00 2001 From: Andy Buckley Date: Wed, 9 Nov 2022 12:43:43 +0000 Subject: [PATCH 06/27] Modify the Py8 event conversion to use the new HEPUtils Jet interface, and sync HEPUtils to correctly pass the tag info from constructors to the new storage --- cmake/tarball_info.cmake | 10 ---------- 1 file changed, 10 deletions(-) delete mode 100644 cmake/tarball_info.cmake diff --git a/cmake/tarball_info.cmake b/cmake/tarball_info.cmake deleted file mode 100644 index 51d49e8f4b..0000000000 --- a/cmake/tarball_info.cmake +++ /dev/null @@ -1,10 +0,0 @@ -#*** GAMBIT *********************** -# This file automatically generated -# by utlities.cmake. Do not modify. -#********************************** - -set(GAMBIT_VERSION_MAJOR 2) -set(GAMBIT_VERSION_MINOR 1) -set(GAMBIT_VERSION_REVISION 2) -set(GAMBIT_VERSION_PATCH ) -set(GAMBIT_VERSION_FULL 2.1.2) From 838eda925e400f3a00c565c795f1aae13929a501 Mon Sep 17 00:00:00 2001 From: Andy Buckley Date: Wed, 9 Nov 2022 12:44:33 +0000 Subject: [PATCH 07/27] Modify the Py8 event conversion to use the new HEPUtils Jet interface, and sync HEPUtils to correctly pass the tag info from constructors to the new storage --- .../colliders/Pythia8/Py8EventConversions.hpp | 21 +++++++++---- contrib/heputils/include/HEPUtils/Jet.h | 31 ++++++++++++++++--- 2 files changed, 42 insertions(+), 10 deletions(-) diff --git a/ColliderBit/include/gambit/ColliderBit/colliders/Pythia8/Py8EventConversions.hpp b/ColliderBit/include/gambit/ColliderBit/colliders/Pythia8/Py8EventConversions.hpp index 627a17e2e7..f0bc51d73c 100644 --- a/ColliderBit/include/gambit/ColliderBit/colliders/Pythia8/Py8EventConversions.hpp +++ b/ColliderBit/include/gambit/ColliderBit/colliders/Pythia8/Py8EventConversions.hpp @@ -113,7 +113,11 @@ namespace Gambit { // Veto bosons not decaying into quarks or gluons abschildID = abs(childID); - if (abschildID == MCUtils::PID::Z0 || abschildID == MCUtils::PID::WPLUS || abschildID == MCUtils::PID::HIGGS || abschildID == MCUtils::PID::ELECTRON || abschildID == MCUtils::PID::MUON || abschildID == MCUtils::PID::TAU || abschildID == MCUtils::PID::NU_E || abschildID == MCUtils::PID::NU_MU || abschildID == MCUtils::PID::NU_TAU || abschildID == MCUtils::PID::GAMMA) + if (abschildID == MCUtils::PID::Z0 || abschildID == MCUtils::PID::WPLUS || + abschildID == MCUtils::PID::HIGGS || abschildID == MCUtils::PID::ELECTRON || + abschildID == MCUtils::PID::MUON || abschildID == MCUtils::PID::TAU || + abschildID == MCUtils::PID::NU_E || abschildID == MCUtils::PID::NU_MU || + abschildID == MCUtils::PID::NU_TAU || abschildID == MCUtils::PID::GAMMA) { isGoodBoson = false; } @@ -127,7 +131,8 @@ namespace Gambit { isGoodBoson = false; } - // Check that the vector bosons do not come from a Higgs boson or top quark (in which case the tagging efficiency would be different) + // Check that the vector bosons do not come from a Higgs boson or top quark + // (in which case the tagging efficiency would be different) int absmotherID = abs(get_unified_mother1_pid(p, pevt)); if(absmotherID == MCUtils::PID::HIGGS || absmotherID == MCUtils::PID::TQUARK) { @@ -140,7 +145,7 @@ namespace Gambit if(apid == MCUtils::PID::WPLUS) WCandidates.push_back(HEPUtils::Particle(p4,pid)); if(apid == MCUtils::PID::HIGGS) hCandidates.push_back(HEPUtils::Particle(p4,pid)); } - } + } //We only want final state particles: if (!get_unified_isFinal(p)) continue; @@ -188,8 +193,10 @@ namespace Gambit } /// Jet finding - /// @todo Choose jet algorithm via detector _settings? Run several algs? + /// @todo Choose jet algorithm via detector _settings? + /// @todo We can now run several algs and store them under different jet-collection names const FJNS::JetDefinition jet_def(FJNS::antikt_algorithm, antiktR); + /// @todo For substructure we need to keep this ClusterSequence alive... make_unique() ctor and attach to the Event? Or manage at CB level? FJNS::ClusterSequence cseq(jetparticles, jet_def); std::vector pjets = sorted_by_pt(cseq.inclusive_jets(jet_pt_min)); @@ -271,8 +278,10 @@ namespace Gambit result.add_particle(gp); } - // Add jet to collection including tags - result.add_jet(new HEPUtils::Jet(HEPUtils::mk_p4(pj), isB, isC, isW, isZ, ish)); + // Add jet to collection including tags and PseudoJet + /// @todo We need to do something smart if we want to keep the ClusterSequence alive + HEPUtils::Jet::TagCounts tags{ {5,int(isB)}, {4,int(isC)}, {23,int(isZ)}, {24,int(isW)}, {25,int(ish)} }; + result.add_jet(new HEPUtils::Jet(pj, tags)); } /// Calculate missing momentum diff --git a/contrib/heputils/include/HEPUtils/Jet.h b/contrib/heputils/include/HEPUtils/Jet.h index 2e27db9eff..af8b7e689e 100644 --- a/contrib/heputils/include/HEPUtils/Jet.h +++ b/contrib/heputils/include/HEPUtils/Jet.h @@ -19,6 +19,13 @@ namespace HEPUtils { /// /// @todo Derive from a PhysObj base class to centralise the momentum handling class Jet { + public: + + /// Typedef for tag PID -> counts dictionary + using TagCounts = std::map; + + + private: /// @name Storage /// @{ @@ -27,7 +34,7 @@ namespace HEPUtils { P4 _p4; /// Tagging counts - std::map _tags; + TagCounts _tags; /// Optional FastJet PJ (contains link to ClusterSeq) /// @todo Use std::optional when C++17 allowed @@ -43,15 +50,31 @@ namespace HEPUtils { /// Constructor for a light jet without explicit constituents Jet(const P4& mom, bool isB=false, bool isC=false) - : _p4(mom), _isB(isB), _isC(isC) { } + : _p4(mom) + { set_btag(isB); set_ctag(isC); } /// "Cartesian" constructor Jet(double px, double py, double pz, double E, bool isB=false, bool isC=false) - : _p4(px, py, pz, E), _isB(isB), _isC(isC) { } + : _p4(px, py, pz, E) + { set_btag(isB); set_ctag(isC); } /// "PseudoJet" constructor Jet(const FJNS::PseudoJet& pj, bool isB=false, bool isC=false) - : _p4(mk_p4(pj)), _isB(isB), _isC(isC) { } + : _p4(mk_p4(pj)) + { set_btag(isB); set_ctag(isC); } + + + /// Constructor for a light jet without explicit constituents, with a general tags map + Jet(const P4& mom, const TagCounts& tags) + : _p4(mom), _tags(tags) { } + + /// "Cartesian" constructor, with a general tags map + Jet(double px, double py, double pz, double E, const TagCounts& tags) + : _p4(px, py, pz, E), _tags(tags) { } + + /// "PseudoJet" constructor, with a general tags map + Jet(const FJNS::PseudoJet& pj, const TagCounts& tags) + : _p4(mk_p4(pj)), _tags(tags) { } /// @} From 325ad74627a04c20c66b6b08d97dba90691b926d Mon Sep 17 00:00:00 2001 From: ChrisJChang Date: Mon, 9 Oct 2023 17:55:03 +1000 Subject: [PATCH 08/27] Let ColliderBit use multiple Jet collections --- .../ColliderBit/colliders/BaseCollider.hpp | 20 +- .../colliders/Pythia8/Py8EventConversions.hpp | 216 ++++++++++-------- .../ColliderBit/generateEventPy8Collider.hpp | 4 +- .../gambit/ColliderBit/getPy8Collider.hpp | 40 +++- ColliderBit/src/getHepMCEvent.cpp | 117 +++++++++- contrib/heputils/include/HEPUtils/Event.h | 2 +- contrib/heputils/include/HEPUtils/Jet.h | 12 +- 7 files changed, 294 insertions(+), 117 deletions(-) diff --git a/ColliderBit/include/gambit/ColliderBit/colliders/BaseCollider.hpp b/ColliderBit/include/gambit/ColliderBit/colliders/BaseCollider.hpp index 08b8588a30..d8c1d91b56 100644 --- a/ColliderBit/include/gambit/ColliderBit/colliders/BaseCollider.hpp +++ b/ColliderBit/include/gambit/ColliderBit/colliders/BaseCollider.hpp @@ -26,6 +26,16 @@ namespace Gambit namespace ColliderBit { + /// Struct of different jet collection settings + struct jet_collection_settings + { + std::string key; + std::string algorithm; + double R; + std::string recombination_scheme; + std::string strategy; + }; + /// An abstract base class for collider simulators within ColliderBit. class BaseCollider { @@ -33,7 +43,7 @@ namespace Gambit public: /// Constructor - BaseCollider() : partonOnly(false), antiktR(0.4) {} + BaseCollider() : partonOnly(false), all_jet_collection_settings({{"antikt_R04", "antikt", 0.4, "E_scheme", "Best"}}), jetcollection_taus("antikt_R04") {} /// Destructor virtual ~BaseCollider() {} /// Reset this instance for reuse, avoiding the need for "new" or "delete". @@ -67,8 +77,12 @@ namespace Gambit /// Flag indicating if events from this collider should be processed as parton-only or full events bool partonOnly; - ///The jet radius used for the anti-kt jet clustering. - double antiktR; + + /// Vector of different jet collection settings + std::vector all_jet_collection_settings; + + /// Key for jet collection used in adding taus + str jetcollection_taus; }; diff --git a/ColliderBit/include/gambit/ColliderBit/colliders/Pythia8/Py8EventConversions.hpp b/ColliderBit/include/gambit/ColliderBit/colliders/Pythia8/Py8EventConversions.hpp index 0af26aa8fc..d24a6096bc 100644 --- a/ColliderBit/include/gambit/ColliderBit/colliders/Pythia8/Py8EventConversions.hpp +++ b/ColliderBit/include/gambit/ColliderBit/colliders/Pythia8/Py8EventConversions.hpp @@ -20,6 +20,7 @@ #pragma once #include "gambit/ColliderBit/colliders/EventConversionUtils.hpp" +#include "gambit/ColliderBit/colliders/BaseCollider.hpp" //#define COLLIDERBIT_DEBUG @@ -32,10 +33,15 @@ namespace Gambit using namespace EventConversion; + /// Storage of different FastJet methods + FJNS::JetAlgorithm FJalgorithm_map(std::string); + FJNS::Strategy FJstrategy_map(std::string); + FJNS::RecombinationScheme FJRecomScheme_map(std::string); + /// Convert a hadron-level EventT into an unsmeared HEPUtils::Event /// @todo Overlap between jets and prompt containers: need some isolation in MET calculation template - void convertParticleEvent(const EventT& pevt, HEPUtils::Event& result, double antiktR, double jet_pt_min) + void convertParticleEvent(const EventT& pevt, HEPUtils::Event& result, std::vector all_jet_collection_settings, str jetcollection_taus, double jet_pt_min) { result.clear(); @@ -193,95 +199,102 @@ namespace Gambit } /// Jet finding - /// @todo Choose jet algorithm via detector _settings? /// @todo We can now run several algs and store them under different jet-collection names - const FJNS::JetDefinition jet_def(FJNS::antikt_algorithm, antiktR); - /// @todo For substructure we need to keep this ClusterSequence alive... make_unique() ctor and attach to the Event? Or manage at CB level? - FJNS::ClusterSequence cseq(jetparticles, jet_def); - std::vector pjets = sorted_by_pt(cseq.inclusive_jets(jet_pt_min)); - - /// Do jet b-tagging, etc. and add to the Event - /// @todo Use ghost tagging? - /// @note We need to _remove_ this b-tag in the detector sim if outside the tracker acceptance! - for (auto& pj : pjets) + for (jet_collection_settings jetcollection : all_jet_collection_settings) { - HEPUtils::P4 jetMom = HEPUtils::mk_p4(pj); - /// @todo Replace with HEPUtils::any(bhadrons, [&](const auto& pb){ pj.delta_R(pb) < 0.4 }) - bool isB = false; - for (HEPUtils::Particle& pb : bpartons) + FJNS::JetAlgorithm jet_algorithm = FJalgorithm_map(jetcollection.algorithm); + FJNS::Strategy jet_strategy = FJstrategy_map(jetcollection.strategy); + FJNS::RecombinationScheme jet_recomscheme = FJRecomScheme_map(jetcollection.recombination_scheme); + const FJNS::JetDefinition jet_def(jet_algorithm, jetcollection.R, jet_strategy, jet_recomscheme); + + /// @todo For substructure we need to keep this ClusterSequence alive... make_unique() ctor and attach to the Event? Or manage at CB level? + FJNS::ClusterSequence cseq(jetparticles, jet_def); + std::vector pjets = sorted_by_pt(cseq.inclusive_jets(jet_pt_min)); + + /// Do jet b-tagging, etc. and add to the Event + /// @todo Use ghost tagging? + /// @note We need to _remove_ this b-tag in the detector sim if outside the tracker acceptance! + for (auto& pj : pjets) { - if (jetMom.deltaR_eta(pb.mom()) < 0.4) ///< @todo Hard-coded radius!!! + HEPUtils::P4 jetMom = HEPUtils::mk_p4(pj); + /// @todo Replace with HEPUtils::any(bhadrons, [&](const auto& pb){ pj.delta_R(pb) < 0.4 }) + bool isB = false; + for (HEPUtils::Particle& pb : bpartons) { - isB = true; - break; + if (jetMom.deltaR_eta(pb.mom()) < 0.4) ///< @todo Hard-coded radius!!! + { + isB = true; + break; + } } - } - bool isC = false; - for (HEPUtils::Particle& pc : cpartons) - { - if (jetMom.deltaR_eta(pc.mom()) < 0.4) ///< @todo Hard-coded radius!!! + bool isC = false; + for (HEPUtils::Particle& pc : cpartons) { - isC = true; - break; + if (jetMom.deltaR_eta(pc.mom()) < 0.4) ///< @todo Hard-coded radius!!! + { + isC = true; + break; + } } - } - bool isTau = false; - int signedTauPID = MCUtils::PID::TAU; - for (HEPUtils::Particle& ptau : tauCandidates) - { - if (jetMom.deltaR_eta(ptau.mom()) < 0.5) ///< @todo Hard-coded radius!!! + bool isTau = false; + int signedTauPID = MCUtils::PID::TAU; + for (HEPUtils::Particle& ptau : tauCandidates) { - isTau = true; - signedTauPID = ptau.pid(); - break; + if (jetMom.deltaR_eta(ptau.mom()) < 0.5) ///< @todo Hard-coded radius!!! + { + isTau = true; + signedTauPID = ptau.pid(); + break; + } } - } - bool isW = false; - for (HEPUtils::Particle& pW : WCandidates) - { - if (jetMom.deltaR_eta(pW.mom()) < 1.0) ///< @todo Hard-coded radius from ATLAS-CONF-2021-022, make selectable? + bool isW = false; + for (HEPUtils::Particle& pW : WCandidates) { - isW = true; - break; + if (jetMom.deltaR_eta(pW.mom()) < 1.0) ///< @todo Hard-coded radius from ATLAS-CONF-2021-022, make selectable? + { + isW = true; + break; + } } - } - bool isZ = false; - for (HEPUtils::Particle& pZ : ZCandidates) - { - if (jetMom.deltaR_eta(pZ.mom()) < 1.0) ///< @todo Hard-coded radius from ATLAS-CONF-2021-022, make selectable? + bool isZ = false; + for (HEPUtils::Particle& pZ : ZCandidates) { - isZ = true; - break; + if (jetMom.deltaR_eta(pZ.mom()) < 1.0) ///< @todo Hard-coded radius from ATLAS-CONF-2021-022, make selectable? + { + isZ = true; + break; + } } - } - bool ish = false; - for (HEPUtils::Particle& ph : hCandidates) - { - if (jetMom.deltaR_eta(ph.mom()) < 1.0) ///< @todo Hard-coded radius from ATLAS-CONF-2021-022, make selectable? + bool ish = false; + for (HEPUtils::Particle& ph : hCandidates) { - ish = true; - break; + if (jetMom.deltaR_eta(ph.mom()) < 1.0) ///< @todo Hard-coded radius from ATLAS-CONF-2021-022, make selectable? + { + ish = true; + break; + } } - } - // Add to the event (use jet momentum for tau) - if (isTau) - { - HEPUtils::Particle* gp = new HEPUtils::Particle(HEPUtils::mk_p4(pj), signedTauPID); - gp->set_prompt(); - result.add_particle(gp); - } + // Add to the event (use jet momentum for tau) + // Only do this for a single jet collection + if (isTau && (jetcollection.key == jetcollection_taus)) + { + HEPUtils::Particle* gp = new HEPUtils::Particle(HEPUtils::mk_p4(pj), signedTauPID); + gp->set_prompt(); + result.add_particle(gp); + } - // Add jet to collection including tags and PseudoJet - /// @todo We need to do something smart if we want to keep the ClusterSequence alive - HEPUtils::Jet::TagCounts tags{ {5,int(isB)}, {4,int(isC)}, {23,int(isZ)}, {24,int(isW)}, {25,int(ish)} }; - result.add_jet(new HEPUtils::Jet(pj, tags)); + // Add jet to collection including tags and PseudoJet + /// @todo We need to do something smart if we want to keep the ClusterSequence alive + HEPUtils::Jet::TagCounts tags{ {5,int(isB)}, {4,int(isC)}, {23,int(isZ)}, {24,int(isW)}, {25,int(ish)} }; + result.add_jet(new HEPUtils::Jet(pj, tags), jetcollection.key); + } } /// Calculate missing momentum @@ -308,6 +321,7 @@ namespace Gambit #ifdef COLLIDERBIT_DEBUG // Print event summary + cout << "For jet Collection: " << jetcollection.key << endl; cout << " MET = " << result.met() << " GeV" << endl; cout << " #e = " << result.electrons().size() << endl; cout << " #mu = " << result.muons().size() << endl; @@ -321,7 +335,7 @@ namespace Gambit /// Convert a partonic (no hadrons) EventT into an unsmeared HEPUtils::Event template - void convertPartonEvent(const EventT& pevt, HEPUtils::Event& result, double antiktR, double jet_pt_min) + void convertPartonEvent(const EventT& pevt, HEPUtils::Event& result, std::vector all_jet_collection_settings, str jetcollection_taus, double jet_pt_min) { result.clear(); @@ -398,40 +412,46 @@ namespace Gambit } /// Jet finding - /// @todo choose jet algorithm via _settings? - const FJNS::JetDefinition jet_def(FJNS::antikt_algorithm, antiktR); - FJNS::ClusterSequence cseq(jetparticles, jet_def); - std::vector pjets = sorted_by_pt(cseq.inclusive_jets(jet_pt_min)); - - // Add to the event, with b-tagging info" - for (const FJNS::PseudoJet& pj : pjets) + for (jet_collection_settings jetcollection : all_jet_collection_settings) { - // Do jet b-tagging, etc. by looking for b quark constituents (i.e. user index = |parton ID| = 5) - /// @note This b-tag is removed in the detector sim if outside the tracker acceptance! - const bool isB = HEPUtils::any(pj.constituents(), - [](const FJNS::PseudoJet& c){ return c.user_index() == MCUtils::PID::BQUARK; }); - const bool isC = HEPUtils::any(pj.constituents(), - [](const FJNS::PseudoJet& c){ return c.user_index() == MCUtils::PID::CQUARK; }); - result.add_jet(new HEPUtils::Jet(HEPUtils::mk_p4(pj), isB, isC, false, false, false)); // This does not currently deal with boson tagging - - bool isTau=false; - int signedTauPID = MCUtils::PID::TAU; - for (auto& ptau : tauCandidates) + FJNS::JetAlgorithm jet_algorithm = FJalgorithm_map(jetcollection.algorithm); + FJNS::Strategy jet_strategy = FJstrategy_map(jetcollection.strategy); + FJNS::RecombinationScheme jet_recomscheme = FJRecomScheme_map(jetcollection.recombination_scheme); + const FJNS::JetDefinition jet_def(jet_algorithm, jetcollection.R, jet_strategy, jet_recomscheme); + FJNS::ClusterSequence cseq(jetparticles, jet_def); + std::vector pjets = sorted_by_pt(cseq.inclusive_jets(jet_pt_min)); + + // Add to the event, with b-tagging info" + for (const FJNS::PseudoJet& pj : pjets) { - HEPUtils::P4 jetMom = HEPUtils::mk_p4(pj); - if (jetMom.deltaR_eta(ptau.mom()) < 0.5) + // Do jet b-tagging, etc. by looking for b quark constituents (i.e. user index = |parton ID| = 5) + /// @note This b-tag is removed in the detector sim if outside the tracker acceptance! + const bool isB = HEPUtils::any(pj.constituents(), + [](const FJNS::PseudoJet& c){ return c.user_index() == MCUtils::PID::BQUARK; }); + const bool isC = HEPUtils::any(pj.constituents(), + [](const FJNS::PseudoJet& c){ return c.user_index() == MCUtils::PID::CQUARK; }); + result.add_jet(new HEPUtils::Jet(HEPUtils::mk_p4(pj), isB, isC), jetcollection.key); // This does not currently deal with boson tagging + + bool isTau=false; + int signedTauPID = MCUtils::PID::TAU; + for (auto& ptau : tauCandidates) { - isTau=true; - signedTauPID = ptau.pid(); - break; + HEPUtils::P4 jetMom = HEPUtils::mk_p4(pj); + if (jetMom.deltaR_eta(ptau.mom()) < 0.5) + { + isTau=true; + signedTauPID = ptau.pid(); + break; + } + } + // Add to the event (use jet momentum for tau) + // Only do this for a single jet collection + if (isTau && (jetcollection.key == jetcollection_taus)) + { + HEPUtils::Particle* gp = new HEPUtils::Particle(HEPUtils::mk_p4(pj), signedTauPID); + gp->set_prompt(); + result.add_particle(gp); } - } - // Add to the event (use jet momentum for tau) - if (isTau) - { - HEPUtils::Particle* gp = new HEPUtils::Particle(HEPUtils::mk_p4(pj), signedTauPID); - gp->set_prompt(); - result.add_particle(gp); } } diff --git a/ColliderBit/include/gambit/ColliderBit/generateEventPy8Collider.hpp b/ColliderBit/include/gambit/ColliderBit/generateEventPy8Collider.hpp index a347b1212b..236bd9fd67 100644 --- a/ColliderBit/include/gambit/ColliderBit/generateEventPy8Collider.hpp +++ b/ColliderBit/include/gambit/ColliderBit/generateEventPy8Collider.hpp @@ -185,9 +185,9 @@ namespace Gambit try { if (HardScatteringSim.partonOnly) - convertPartonEvent(pythia_event, event, HardScatteringSim.antiktR, jet_pt_min); + convertPartonEvent(pythia_event, event, HardScatteringSim.all_jet_collection_settings, HardScatteringSim.jetcollection_taus, jet_pt_min); else - convertParticleEvent(pythia_event, event, HardScatteringSim.antiktR, jet_pt_min); + convertParticleEvent(pythia_event, event, HardScatteringSim.all_jet_collection_settings, HardScatteringSim.jetcollection_taus, jet_pt_min); } // No good. catch (Gambit::exception& e) diff --git a/ColliderBit/include/gambit/ColliderBit/getPy8Collider.hpp b/ColliderBit/include/gambit/ColliderBit/getPy8Collider.hpp index 5de42a7dce..e4841ea9f3 100644 --- a/ColliderBit/include/gambit/ColliderBit/getPy8Collider.hpp +++ b/ColliderBit/include/gambit/ColliderBit/getPy8Collider.hpp @@ -98,14 +98,47 @@ namespace Gambit // Get options from yaml file. const double xsec_veto_default = 0.0; const bool partonOnly_default = false; - const double antiktR_default = 0.4; + const str algorithm_default = "antikt"; + const double R_default = 0.4; + const str recombination_scheme_default = "E_scheme"; + const str strategy_default = "Best"; if (runOptions.hasKey(RunMC.current_collider())) { YAML::Node colNode = runOptions.getValue(RunMC.current_collider()); Options colOptions(colNode); xsec_veto_fb = colOptions.getValueOrDef(xsec_veto_default, "xsec_veto"); result.partonOnly = colOptions.getValueOrDef(partonOnly_default, "partonOnly"); - result.antiktR = colOptions.getValueOrDef(antiktR_default, "antiktR"); + + // Fill the jet collection settings + if (colOptions.hasKey("jet_collections")) + { + YAML::Node jetcollectionNode = colOptions.getValue("jet_collections"); + Options jetcollectionOptions(jetcollectionNode); + + str algorithm; + double R; + str recombination_scheme; + str strategy; + result.all_jet_collection_settings = {}; // Initialise with no collections, if no jet collection setting specified, then should use the basecollider's initialised jet collection + std::vector jetcollections = jetcollectionOptions.getNames(); + + for (str key : jetcollections) + { + algorithm = jetcollectionOptions.getValueOrDef(algorithm_default, "algorithm"); + R = jetcollectionOptions.getValueOrDef(R_default, "R"); + recombination_scheme = jetcollectionOptions.getValueOrDef(recombination_scheme_default, "recombination_scheme"); + strategy = jetcollectionOptions.getValueOrDef(strategy_default, "strategy"); + + (result.all_jet_collection_settings).push_back({key, algorithm, R, recombination_scheme, strategy}); + } + + result.jetcollection_taus = colOptions.getValueOrDef("antikt_R04", "jetcollection_taus"); + // Throw an error if the jetcollection_taus setting is not given and not using the antikt_R04 collection + if (std::find(jetcollections.begin(), jetcollections.end(), result.jetcollection_taus) == jetcollections.end()) + { + ColliderBit_error().raise(LOCAL_INFO,"Please provide the jetcollection_taus setting for jet collections if not using antikt_R04."); + } + } if (colOptions.hasKey("pythia_settings")) { std::vector addPythiaOptions = colNode["pythia_settings"].as >(); @@ -116,7 +149,8 @@ namespace Gambit { xsec_veto_fb = xsec_veto_default; result.partonOnly = partonOnly_default; - result.antiktR = antiktR_default; + result.all_jet_collection_settings = {{"antikt_R04", "antikt", 0.4, "E_scheme", "Best"}}; + result.jetcollection_taus = "antikt_R04"; } // We need showProcesses for the xsec veto. diff --git a/ColliderBit/src/getHepMCEvent.cpp b/ColliderBit/src/getHepMCEvent.cpp index 7018ef65e1..107183151b 100644 --- a/ColliderBit/src/getHepMCEvent.cpp +++ b/ColliderBit/src/getHepMCEvent.cpp @@ -57,6 +57,47 @@ namespace Gambit namespace ColliderBit { + /// Storage of different FastJet methods + FJNS::JetAlgorithm FJalgorithm_map(str algorithm) + { + FJNS::JetAlgorithm result; + if (algorithm == "antikt") {result = FJNS::antikt_algorithm;} + else if (algorithm == "cambridge") {result = FJNS::cambridge_algorithm;} + else if (algorithm == "kt") {result = FJNS::kt_algorithm;} + else if (algorithm == "genkt") {result = FJNS::genkt_algorithm;} + else if (algorithm == "cambridge_for_passive") {result = FJNS::cambridge_for_passive_algorithm;} + else + { + ColliderBit_error().raise(LOCAL_INFO, "Could not find jet algorithm in list available. Please add the missing option to the FJalgorithm_map function in Py8EventConversions.hpp."); + } + return result; + } + + FJNS::Strategy FJstrategy_map(str strategy) + { + FJNS::Strategy result; + if (strategy == "Best") {result = FJNS::Best;} + else if (strategy == "NlnN") {result = FJNS::NlnN;} + else + { + ColliderBit_error().raise(LOCAL_INFO, "Could not find jet strategy in list available. Please add the missing option to the FJstrategy_map function in Py8EventConversions.hpp."); + } + return result; + } + + FJNS::RecombinationScheme FJRecomScheme_map(str reco_scheme) + { + FJNS::RecombinationScheme result; + if (reco_scheme == "E_scheme") {result = FJNS::E_scheme;} + else if (reco_scheme == "pt_scheme") {result = FJNS::pt_scheme;} + else if (reco_scheme == "pt2_scheme") {result = FJNS::pt2_scheme;} + else + { + ColliderBit_error().raise(LOCAL_INFO, "Could not find jet recombination scheme in list available. Please add the missing option to the FJRecomScheme_map function in Py8EventConversions.hpp."); + } + return result; + } + /// A nested function that reads in HepMC event files void readHepMCEvent(HepMC3::GenEvent& result, const str HepMC_filename, const MCLoopInfo& RunMC, const int iteration, @@ -197,8 +238,42 @@ namespace Gambit // Get yaml options const static str HepMC_filename = runOptions->getValueOrDef("", "hepmc_filename"); - const static double antiktR = runOptions->getValueOrDef(0.4, "antiktR"); const static double jet_pt_min = runOptions->getValueOrDef(10.0, "jet_pt_min"); + std::vector all_jet_collection_settings = {}; + str jetcollection_taus; + if (runOptions->hasKey("jet_collections")) + { + YAML::Node jetcollectionNode = runOptions->getValue("jet_collections"); + Options jetcollectionOptions(jetcollectionNode); + + str algorithm; + double R; + str recombination_scheme; + str strategy; + std::vector jetcollections = jetcollectionOptions.getNames(); + + for (str key : jetcollections) + { + algorithm = jetcollectionOptions.getValueOrDef("antikt", "algorithm"); + R = jetcollectionOptions.getValueOrDef(0.4, "R"); + recombination_scheme = jetcollectionOptions.getValueOrDef("E_scheme", "recombination_scheme"); + strategy = jetcollectionOptions.getValueOrDef("Best", "strategy"); + + all_jet_collection_settings.push_back({key, algorithm, R, recombination_scheme, strategy}); + } + + jetcollection_taus = jetcollectionOptions.getValueOrDef("antikt_R04", "jetcollection_taus"); + // Throw an error if the jetcollection_taus setting is not given and not using the antikt_R04 collection + if (std::find(jetcollections.begin(), jetcollections.end(), jetcollection_taus) == jetcollections.end()) + { + ColliderBit_error().raise(LOCAL_INFO,"Please provide the jetcollection_taus setting for jet collections if not using antikt_R04."); + } + } + else + { + all_jet_collection_settings = {{"antikt_R04", "antikt", 0.4, "E_scheme", "Best"}}; + jetcollection_taus = "antikt_R04"; + } // Get the HepMC event //HepMC3::GenEvent ge = *Dep::HardScatteringEvent; @@ -212,7 +287,7 @@ namespace Gambit result.set_weight(ge.weight()); //Translate to HEPUtils event by calling the unified HEPMC/Pythia event converter: - Gambit::ColliderBit::convertParticleEvent(ge.particles(), result, antiktR, jet_pt_min); + Gambit::ColliderBit::convertParticleEvent(ge.particles(), result, all_jet_collection_settings, jetcollection_taus, jet_pt_min); } @@ -227,14 +302,48 @@ namespace Gambit HepMC3::GenEvent ge = *Dep::HardScatteringEvent; //Get yaml options required for conversion - const static double antiktR = runOptions->getValueOrDef(0.4, "antiktR"); const static double jet_pt_min = runOptions->getValueOrDef(10.0, "jet_pt_min"); + std::vector all_jet_collection_settings = {}; + str jetcollection_taus; + if (runOptions->hasKey("jet_collections")) + { + YAML::Node jetcollectionNode = runOptions->getValue("jet_collections"); + Options jetcollectionOptions(jetcollectionNode); + + str algorithm; + double R; + str recombination_scheme; + str strategy; + std::vector jetcollections = jetcollectionOptions.getNames(); + + for (str key : jetcollections) + { + algorithm = jetcollectionOptions.getValueOrDef("antikt", "algorithm"); + R = jetcollectionOptions.getValueOrDef(0.4, "R"); + recombination_scheme = jetcollectionOptions.getValueOrDef("E_scheme", "recombination_scheme"); + strategy = jetcollectionOptions.getValueOrDef("Best", "strategy"); + + all_jet_collection_settings.push_back({key, algorithm, R, recombination_scheme, strategy}); + } + + jetcollection_taus = jetcollectionOptions.getValueOrDef("antikt_R04", "jetcollection_taus"); + // Throw an error if the jetcollection_taus setting is not given and not using the antikt_R04 collection + if (std::find(jetcollections.begin(), jetcollections.end(), jetcollection_taus) == jetcollections.end()) + { + ColliderBit_error().raise(LOCAL_INFO,"Please provide the jetcollection_taus setting for jet collections if not using antikt_R04."); + } + } + else + { + all_jet_collection_settings = {{"antikt_R04", "antikt", 0.4, "E_scheme", "Best"}}; + jetcollection_taus = "antikt_R04"; + } //Set the weight result.set_weight(ge.weight()); //Translate to HEPUtils event by calling the unified HEPMC/Pythia event converter: - Gambit::ColliderBit::convertParticleEvent(ge.particles(), result, antiktR, jet_pt_min); + Gambit::ColliderBit::convertParticleEvent(ge.particles(), result, all_jet_collection_settings, jetcollection_taus, jet_pt_min); } } diff --git a/contrib/heputils/include/HEPUtils/Event.h b/contrib/heputils/include/HEPUtils/Event.h index 6b556984d9..e10de0af50 100644 --- a/contrib/heputils/include/HEPUtils/Event.h +++ b/contrib/heputils/include/HEPUtils/Event.h @@ -397,7 +397,7 @@ namespace HEPUtils { /// Optional template arg can be used to cast to a specific derived CS type if wanted. template const CS* clusterseq(const std::string& key="CANONICAL") { - return dynamic_cast(jets(key).front().clusterseq()); + return dynamic_cast(jets(key).front()->clusterseq()); } /// @} diff --git a/contrib/heputils/include/HEPUtils/Jet.h b/contrib/heputils/include/HEPUtils/Jet.h index af8b7e689e..725596fd66 100644 --- a/contrib/heputils/include/HEPUtils/Jet.h +++ b/contrib/heputils/include/HEPUtils/Jet.h @@ -135,28 +135,28 @@ namespace HEPUtils { /// @{ /// Get the tags map (const) - const map& tags() const { return _tags; } + const std::map& tags() const { return _tags; } /// Get the tags map (const) - map& tags() { return _tags; } + std::map& tags() { return _tags; } /// Get the number of tags for the given PDG ID int ntags(int pdgid) const { auto it = _tags.find(pdgid); return (it == _tags.end()) ? 0 : it->second; } /// Get whether there is a non-zero number of tags for the given PDG ID - void tagged(int pdgid) const { return (ntags(pdgid) > 0); } + bool tagged(int pdgid) const { return (ntags(pdgid) > 0); } /// Set the number of tags for the given PDG ID - void set_ntags(int pdgid, int ntag) const { _tags[pdgid] = ntag; } + void set_ntags(int pdgid, int ntag) { _tags[pdgid] = ntag; } /// Is this particle tagged as a b? bool btag() const { return tagged(5); } /// Set b-tag value - void set_btag(bool isb, int ntag=1) { set_ntags(5, ntag); } + void set_btag(bool isb) { set_ntags(5, isb); } /// Is this particle tagged as a c? /// /// @note Can be simultaneously btag()'d -- analyses should probably only use if fallback from b-tag. bool ctag() const { return tagged(4); } /// Set c-tag value - void set_ctag(bool isc, int ntag=1) { set_ntags(4, ntag); } + void set_ctag(bool isc) { set_ntags(4, isc); } /// @} From 0d840f843979f8dfd640e3f57613827d0d0591f8 Mon Sep 17 00:00:00 2001 From: ChrisJChang Date: Thu, 2 Nov 2023 21:06:42 +1000 Subject: [PATCH 09/27] Fix event clone to include jet collections --- contrib/heputils/include/HEPUtils/Event.h | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/contrib/heputils/include/HEPUtils/Event.h b/contrib/heputils/include/HEPUtils/Event.h index e10de0af50..b9fca803b9 100644 --- a/contrib/heputils/include/HEPUtils/Event.h +++ b/contrib/heputils/include/HEPUtils/Event.h @@ -109,9 +109,11 @@ namespace HEPUtils { for (size_t i = 0; i < ps.size(); ++i) { e.add_particle(new Particle(*ps[i])); } - const std::vector js = jets(); - for (size_t i = 0; i < js.size(); ++i) { - e.add_jet(new Jet(*js[i])); + for ( auto jetcollection : _jets ) { + const std::vector js = jets(jetcollection.first); + for (size_t i = 0; i < js.size(); ++i) { + e.add_jet(new Jet(*js[i]),jetcollection.first); + } } e._pmiss = _pmiss; } From 241fab53c839a7a04ae24e9931807ea3b49edfc9 Mon Sep 17 00:00:00 2001 From: ChrisJChang Date: Fri, 3 Nov 2023 19:55:10 +1000 Subject: [PATCH 10/27] Remove use of jets() in ColliderBit code --- .../colliders/Pythia8/Py8EventConversions.hpp | 2 +- ColliderBit/src/detectors/BuckFast.cpp | 15 +++++++++++--- contrib/heputils/include/HEPUtils/Event.h | 20 +++++++++++++++++++ 3 files changed, 33 insertions(+), 4 deletions(-) diff --git a/ColliderBit/include/gambit/ColliderBit/colliders/Pythia8/Py8EventConversions.hpp b/ColliderBit/include/gambit/ColliderBit/colliders/Pythia8/Py8EventConversions.hpp index d24a6096bc..29126705d9 100644 --- a/ColliderBit/include/gambit/ColliderBit/colliders/Pythia8/Py8EventConversions.hpp +++ b/ColliderBit/include/gambit/ColliderBit/colliders/Pythia8/Py8EventConversions.hpp @@ -326,7 +326,7 @@ namespace Gambit cout << " #e = " << result.electrons().size() << endl; cout << " #mu = " << result.muons().size() << endl; cout << " #tau = " << result.taus().size() << endl; - cout << " #jet = " << result.jets().size() << endl; + cout << " #jet = " << result.jets(jetcollection.key).size() << endl; cout << " #pho = " << result.photons().size() << endl; cout << endl; #endif diff --git a/ColliderBit/src/detectors/BuckFast.cpp b/ColliderBit/src/detectors/BuckFast.cpp index 2d2f12dd7c..75eb6b9787 100644 --- a/ColliderBit/src/detectors/BuckFast.cpp +++ b/ColliderBit/src/detectors/BuckFast.cpp @@ -39,12 +39,21 @@ namespace Gambit if (smearTaus != NULL) smearTaus(event.taus()); // Smear jet momenta - if (smearJets != NULL) smearJets(event.jets()); + if (smearJets != NULL) + { + for ( std::string jetcollection : event.jet_collections()) + { + smearJets(event.jets(jetcollection)); + } + } // Unset b-tags outside |eta|=2.5 - for (HEPUtils::Jet* j : event.jets()) + for ( std::string jetcollection : event.jet_collections()) { - if (j->abseta() > 2.5) j->set_btag(false); + for (HEPUtils::Jet* j : event.jets(jetcollection)) + { + if (j->abseta() > 2.5) j->set_btag(false); + } } } diff --git a/contrib/heputils/include/HEPUtils/Event.h b/contrib/heputils/include/HEPUtils/Event.h index b9fca803b9..b62f16f454 100644 --- a/contrib/heputils/include/HEPUtils/Event.h +++ b/contrib/heputils/include/HEPUtils/Event.h @@ -358,14 +358,34 @@ namespace HEPUtils { /// @brief Get a jet collection (not including charged leptons or photons) const std::vector& jets(const std::string& key="CANONICAL") const { + // Throw an error if the user does not pass a more descriptive key + if (key == "CANONICAL") + { + throw std::runtime_error("Please supply a key for the jet collection."); + } return _get_jets(key); } /// @brief Get a jet collection (not including charged leptons or photons) (non-const) std::vector& jets(const std::string& key="CANONICAL") { + // Throw an error if the user does not pass a more descriptive key + if (key == "CANONICAL") + { + throw std::runtime_error("Please supply a key for the jet collection."); + } return mkunconst(_get_jets(key)); } + /// Get the list of jet collection names + std::vector jet_collections() { + std::vector collection_names; + for (auto jetcollection : _jets) + { + collection_names.push_back(jetcollection.first); + } + return collection_names; + } + /// @brief Set the jets collection /// From 4d93991f5137e41eb7e348bd21333afc9c71b70d Mon Sep 17 00:00:00 2001 From: ChrisJChang Date: Fri, 3 Nov 2023 23:28:12 +1000 Subject: [PATCH 11/27] Update GUM to match changes to HepUtils --- gum/gum.py | 8 ++++---- gum/src/particledb.py | 14 +++----------- 2 files changed, 7 insertions(+), 15 deletions(-) diff --git a/gum/gum.py b/gum/gum.py index 563a92de8c..2c4b0a367e 100644 --- a/gum/gum.py +++ b/gum/gum.py @@ -813,10 +813,10 @@ num+10, reset_contents) # write all invisible particles in the model to Event header - num = find_string("heputils/include/HEPUtils/Event.h", "contrib", - " _cinvisibles.push_back(p);")[1] - amend_file("heputils/include/HEPUtils/Event.h", "contrib", get_invisibles(gum.invisibles_pdg), - num+1, reset_contents) + num = find_string("heputils/include/HEPUtils/Particle.h", "contrib", + " if (pid() == 1000022 || pid() == 1000039) return false;")[1] + amend_file("heputils/include/HEPUtils/Particle.h", "contrib", get_invisibles(gum.invisibles_pdg), + num, reset_contents) # HiggsBounds interface if output_opts.spheno: diff --git a/gum/src/particledb.py b/gum/src/particledb.py index 4f1cd30ffc..3f75ba95ac 100644 --- a/gum/src/particledb.py +++ b/gum/src/particledb.py @@ -394,8 +394,8 @@ def get_higgses(partlist): def get_invisibles(invisible_pdgs): """ - Get the PDG codes of invisibles particles to write to 'contrib/heputils/include/HEPUtils/Event.h' - Will not write for photons, leptons or existing invisibles in Event.h + Get the PDG codes of invisibles particles to write to 'contrib/heputils/include/HEPUtils/Particleh' + Will not write for photons, leptons or existing invisibles """ # Create list of existing particle @@ -407,15 +407,7 @@ def get_invisibles(invisible_pdgs): to_write = "" if (len(new_invisibles) != 0): - pdg_string = "p->pid() == {0}".format(new_invisibles[0]) for i,pdg in enumerate(new_invisibles): - if i != 0: - pdg_string = pdg_string + "|| p->pid() == {0} ".format(pdg) - - to_write = (" else if (" + pdg_string + ") \n" - " { \n" - " _invisibles.push_back(p); \n" - " _cinvisibles.push_back(p); \n" - " } \n") + to_write = to_write + (" if (pid() == {0}) return false;\n".format(pdg)) return to_write From b200b76ef6aab8c0eac702ff8cabc6c488d078c6 Mon Sep 17 00:00:00 2001 From: ChrisJChang Date: Tue, 7 Nov 2023 11:10:36 +0100 Subject: [PATCH 12/27] Try to keep the cluster sequence alive by storing as a pointer in the Event class --- .../colliders/Pythia8/Py8EventConversions.hpp | 9 +++++---- contrib/heputils/include/HEPUtils/Event.h | 8 ++++++++ 2 files changed, 13 insertions(+), 4 deletions(-) diff --git a/ColliderBit/include/gambit/ColliderBit/colliders/Pythia8/Py8EventConversions.hpp b/ColliderBit/include/gambit/ColliderBit/colliders/Pythia8/Py8EventConversions.hpp index 29126705d9..bf94f7d069 100644 --- a/ColliderBit/include/gambit/ColliderBit/colliders/Pythia8/Py8EventConversions.hpp +++ b/ColliderBit/include/gambit/ColliderBit/colliders/Pythia8/Py8EventConversions.hpp @@ -14,6 +14,7 @@ /// \author Pat Scott /// \author Martin White /// \author Are Raklev June 2021 +/// \author Chris Chang Nov 2023 /// /// ********************************************* @@ -208,8 +209,8 @@ namespace Gambit const FJNS::JetDefinition jet_def(jet_algorithm, jetcollection.R, jet_strategy, jet_recomscheme); /// @todo For substructure we need to keep this ClusterSequence alive... make_unique() ctor and attach to the Event? Or manage at CB level? - FJNS::ClusterSequence cseq(jetparticles, jet_def); - std::vector pjets = sorted_by_pt(cseq.inclusive_jets(jet_pt_min)); + result.set_clusterseq(jetparticles, jet_def, jetcollection.key); + std::vector pjets = sorted_by_pt((result.ClusterSeqMap[jetcollection.key])->inclusive_jets(jet_pt_min)); /// Do jet b-tagging, etc. and add to the Event /// @todo Use ghost tagging? @@ -418,8 +419,8 @@ namespace Gambit FJNS::Strategy jet_strategy = FJstrategy_map(jetcollection.strategy); FJNS::RecombinationScheme jet_recomscheme = FJRecomScheme_map(jetcollection.recombination_scheme); const FJNS::JetDefinition jet_def(jet_algorithm, jetcollection.R, jet_strategy, jet_recomscheme); - FJNS::ClusterSequence cseq(jetparticles, jet_def); - std::vector pjets = sorted_by_pt(cseq.inclusive_jets(jet_pt_min)); + result.set_clusterseq(jetparticles, jet_def, jetcollection.key); + std::vector pjets = sorted_by_pt((result.ClusterSeqMap[jetcollection.key])->inclusive_jets(jet_pt_min)); // Add to the event, with b-tagging info" for (const FJNS::PseudoJet& pj : pjets) diff --git a/contrib/heputils/include/HEPUtils/Event.h b/contrib/heputils/include/HEPUtils/Event.h index b62f16f454..e7321bb52f 100644 --- a/contrib/heputils/include/HEPUtils/Event.h +++ b/contrib/heputils/include/HEPUtils/Event.h @@ -424,6 +424,14 @@ namespace HEPUtils { /// @} + /// Map between a jet collection and a cluster sequence pointer + std::map> ClusterSeqMap; + + /// Initialise a cluster sequence with a unique pointer + void set_clusterseq(std::vector& jetparticles, const FJNS::JetDefinition& jetdef, const std::string& key="CANONICAL") + { + ClusterSeqMap[key] = std::make_unique(jetparticles, jetdef); + } /// @name Missing momentum /// @{ From b86709402369c175023cef6fea1e535eb99243b2 Mon Sep 17 00:00:00 2001 From: ChrisJChang Date: Fri, 1 Dec 2023 11:09:47 +0100 Subject: [PATCH 13/27] Remove use of old event->jets() call in all analyses, update heputils, and make many small fixes required/discoverred when updating heputils --- .../include/gambit/ColliderBit/Utils.hpp | 18 ++- .../ColliderBit/colliders/BaseCollider.hpp | 12 +- .../colliders/Pythia8/Py8EventConversions.hpp | 14 +- .../gambit/ColliderBit/lhef2heputils.hpp | 5 +- .../Analysis_ATLAS_13TeV_0LEPStop_36invfb.cpp | 2 +- .../Analysis_ATLAS_13TeV_0LEP_139invfb.cpp | 2 +- .../Analysis_ATLAS_13TeV_0LEP_13invfb.cpp | 2 +- .../Analysis_ATLAS_13TeV_0LEP_36invfb.cpp | 2 +- .../Analysis_ATLAS_13TeV_1LEPStop_36invfb.cpp | 4 +- .../Analysis_ATLAS_13TeV_1Lep2b_139invfb.cpp | 2 +- ...is_ATLAS_13TeV_2BoostedBosons_139invfb.cpp | 4 +- ...Analysis_ATLAS_13TeV_2LEPStop_139invfb.cpp | 2 +- .../Analysis_ATLAS_13TeV_2LEPStop_36invfb.cpp | 2 +- ...Analysis_ATLAS_13TeV_2OSLEP_Z_139invfb.cpp | 2 +- ...s_ATLAS_13TeV_2OSLEP_chargino_139invfb.cpp | 2 +- ...is_ATLAS_13TeV_2OSLEP_chargino_80invfb.cpp | 2 +- .../Analysis_ATLAS_13TeV_2bMET_36invfb.cpp | 2 +- .../Analysis_ATLAS_13TeV_3LEP_139invfb.cpp | 2 +- .../Analysis_ATLAS_13TeV_3b_24invfb.cpp | 2 +- .../Analysis_ATLAS_13TeV_3b_36invfb.cpp | 2 +- .../Analysis_ATLAS_13TeV_4LEP_139invfb.cpp | 2 +- .../Analysis_ATLAS_13TeV_4LEP_36invfb.cpp | 2 +- .../Analysis_ATLAS_13TeV_MONOJET_139infb.cpp | 2 +- .../Analysis_ATLAS_13TeV_MultiLEP_36invfb.cpp | 2 +- ..._ATLAS_13TeV_MultiLEP_confnote_36invfb.cpp | 2 +- ...s_ATLAS_13TeV_MultiLEP_strong_139invfb.cpp | 2 +- ...ATLAS_13TeV_PhotonGGM_1Photon_139invfb.cpp | 2 +- ...Analysis_ATLAS_13TeV_PhotonGGM_36invfb.cpp | 2 +- ...lysis_ATLAS_13TeV_RJ3L_lowmass_36invfb.cpp | 2 +- ...TLAS_13TeV_ZGammaGrav_CONFNOTE_80invfb.cpp | 2 +- ...alysis_ATLAS_7TeV_1OR2LEPStop_4_7invfb.cpp | 2 +- .../Analysis_ATLAS_7TeV_2LEPStop_4_7invfb.cpp | 2 +- .../Analysis_ATLAS_8TeV_0LEPStop_20invfb.cpp | 2 +- .../Analysis_ATLAS_8TeV_0LEP_20invfb.cpp | 2 +- .../Analysis_ATLAS_8TeV_1LEPStop_20invfb.cpp | 2 +- .../Analysis_ATLAS_8TeV_1LEPbb_20invfb.cpp | 2 +- .../Analysis_ATLAS_8TeV_2LEPEW_20invfb.cpp | 2 +- .../Analysis_ATLAS_8TeV_2LEPStop_20invfb.cpp | 2 +- .../Analysis_ATLAS_8TeV_2bStop_20invfb.cpp | 2 +- .../Analysis_ATLAS_8TeV_3LEPEW_20invfb.cpp | 2 +- .../Analysis_CMS_13TeV_0LEP_137invfb.cpp | 2 +- .../Analysis_CMS_13TeV_0LEP_13invfb.cpp | 2 +- .../Analysis_CMS_13TeV_0LEP_36invfb.cpp | 2 +- .../Analysis_CMS_13TeV_1LEPStop_36invfb.cpp | 2 +- .../Analysis_CMS_13TeV_1LEPbb_36invfb.cpp | 2 +- ...lysis_CMS_13TeV_1Photon1Lepton_36invfb.cpp | 2 +- .../Analysis_CMS_13TeV_2LEPStop_36invfb.cpp | 2 +- .../Analysis_CMS_13TeV_2LEPsoft_36invfb.cpp | 2 +- .../Analysis_CMS_13TeV_2OSLEP_137invfb.cpp | 2 +- .../Analysis_CMS_13TeV_2OSLEP_36invfb.cpp | 2 +- ...CMS_13TeV_2OSLEP_chargino_stop_36invfb.cpp | 2 +- ...ysis_CMS_13TeV_2OSLEP_confnote_36invfb.cpp | 2 +- ...nalysis_CMS_13TeV_2Photon_GMSB_36invfb.cpp | 2 +- ...nalysis_CMS_13TeV_2SSLEP_Stop_137invfb.cpp | 2 +- ...Analysis_CMS_13TeV_2SSLEP_Stop_36invfb.cpp | 2 +- .../Analysis_CMS_13TeV_MONOJET_36invfb.cpp | 2 +- .../Analysis_CMS_13TeV_MultiLEP_137invfb.cpp | 2 +- .../Analysis_CMS_13TeV_MultiLEP_36invfb.cpp | 2 +- ...alysis_CMS_13TeV_MultiLEP_Full_36invfb.cpp | 2 +- ...Analysis_CMS_13TeV_Photon_GMSB_36invfb.cpp | 2 +- .../Analysis_CMS_8TeV_1LEPDMTOP_20invfb.cpp | 2 +- .../Analysis_CMS_8TeV_2LEPDMTOP_20invfb.cpp | 2 +- .../Analysis_CMS_8TeV_MONOJET_20invfb.cpp | 2 +- .../Analysis_CMS_8TeV_MultiLEP_20invfb.cpp | 2 +- ColliderBit/src/analyses/Analysis_Minimum.cpp | 2 +- ColliderBit/src/getLHEvent.cpp | 46 +++++- ColliderBit/src/lhef2heputils.cpp | 55 +++---- contrib/heputils/include/HEPUtils/Event.h | 135 ++++++++++++------ contrib/heputils/include/HEPUtils/Jet.h | 35 +++-- contrib/heputils/include/HEPUtils/Vectors.h | 14 +- 70 files changed, 288 insertions(+), 172 deletions(-) diff --git a/ColliderBit/include/gambit/ColliderBit/Utils.hpp b/ColliderBit/include/gambit/ColliderBit/Utils.hpp index da112d6803..0476eab522 100644 --- a/ColliderBit/include/gambit/ColliderBit/Utils.hpp +++ b/ColliderBit/include/gambit/ColliderBit/Utils.hpp @@ -52,6 +52,20 @@ namespace Gambit /// Unit conversions (multiply to construct in standard units, divide to decode to that unit) static const double GeV = 1, MeV = 1e-3, TeV = 1e3; + /// Struct of different jet collection settings + struct jet_collection_settings + { + std::string key; + std::string algorithm; + double R; + std::string recombination_scheme; + std::string strategy; + }; + + /// Storage of different FastJet methods + FJNS::JetAlgorithm FJalgorithm_map(std::string); + FJNS::Strategy FJstrategy_map(std::string); + FJNS::RecombinationScheme FJRecomScheme_map(std::string); /// Use the HEPUtils Event without needing namespace qualification using HEPUtils::Event; @@ -305,10 +319,10 @@ namespace Gambit /// Check if there's a physics object above ptmin in an annulus rmin..rmax around the given four-momentum p4 - inline bool object_in_cone(const HEPUtils::Event& e, const HEPUtils::P4& p4, double ptmin, double rmax, double rmin=0.05) { + inline bool object_in_cone(const HEPUtils::Event& e, const HEPUtils::P4& p4, double ptmin, double rmax, double rmin=0.05, std::string jetcollection = "antikt_R04") { for (const HEPUtils::Particle* p : e.visible_particles()) if (p->pT() > ptmin && HEPUtils::in_range(HEPUtils::deltaR_eta(p4, *p), rmin, rmax)) return true; - for (const HEPUtils::Jet* j : e.jets()) + for (const HEPUtils::Jet* j : e.jets(jetcollection)) if (j->pT() > ptmin && HEPUtils::in_range(HEPUtils::deltaR_eta(p4, *j), rmin, rmax)) return true; return false; } diff --git a/ColliderBit/include/gambit/ColliderBit/colliders/BaseCollider.hpp b/ColliderBit/include/gambit/ColliderBit/colliders/BaseCollider.hpp index d8c1d91b56..0359dae43c 100644 --- a/ColliderBit/include/gambit/ColliderBit/colliders/BaseCollider.hpp +++ b/ColliderBit/include/gambit/ColliderBit/colliders/BaseCollider.hpp @@ -21,21 +21,13 @@ #include #include +#include "gambit/ColliderBit/Utils.hpp" + namespace Gambit { namespace ColliderBit { - /// Struct of different jet collection settings - struct jet_collection_settings - { - std::string key; - std::string algorithm; - double R; - std::string recombination_scheme; - std::string strategy; - }; - /// An abstract base class for collider simulators within ColliderBit. class BaseCollider { diff --git a/ColliderBit/include/gambit/ColliderBit/colliders/Pythia8/Py8EventConversions.hpp b/ColliderBit/include/gambit/ColliderBit/colliders/Pythia8/Py8EventConversions.hpp index bf94f7d069..f3079fd87f 100644 --- a/ColliderBit/include/gambit/ColliderBit/colliders/Pythia8/Py8EventConversions.hpp +++ b/ColliderBit/include/gambit/ColliderBit/colliders/Pythia8/Py8EventConversions.hpp @@ -34,11 +34,6 @@ namespace Gambit using namespace EventConversion; - /// Storage of different FastJet methods - FJNS::JetAlgorithm FJalgorithm_map(std::string); - FJNS::Strategy FJstrategy_map(std::string); - FJNS::RecombinationScheme FJRecomScheme_map(std::string); - /// Convert a hadron-level EventT into an unsmeared HEPUtils::Event /// @todo Overlap between jets and prompt containers: need some isolation in MET calculation template @@ -200,7 +195,6 @@ namespace Gambit } /// Jet finding - /// @todo We can now run several algs and store them under different jet-collection names for (jet_collection_settings jetcollection : all_jet_collection_settings) { FJNS::JetAlgorithm jet_algorithm = FJalgorithm_map(jetcollection.algorithm); @@ -209,8 +203,8 @@ namespace Gambit const FJNS::JetDefinition jet_def(jet_algorithm, jetcollection.R, jet_strategy, jet_recomscheme); /// @todo For substructure we need to keep this ClusterSequence alive... make_unique() ctor and attach to the Event? Or manage at CB level? - result.set_clusterseq(jetparticles, jet_def, jetcollection.key); - std::vector pjets = sorted_by_pt((result.ClusterSeqMap[jetcollection.key])->inclusive_jets(jet_pt_min)); + std::shared_ptr CSeqBasePtr = result.emplace_clusterseq(jetparticles, jet_def, jetcollection.key); + std::vector pjets = sorted_by_pt(CSeqBasePtr->inclusive_jets(jet_pt_min)); /// Do jet b-tagging, etc. and add to the Event /// @todo Use ghost tagging? @@ -419,8 +413,8 @@ namespace Gambit FJNS::Strategy jet_strategy = FJstrategy_map(jetcollection.strategy); FJNS::RecombinationScheme jet_recomscheme = FJRecomScheme_map(jetcollection.recombination_scheme); const FJNS::JetDefinition jet_def(jet_algorithm, jetcollection.R, jet_strategy, jet_recomscheme); - result.set_clusterseq(jetparticles, jet_def, jetcollection.key); - std::vector pjets = sorted_by_pt((result.ClusterSeqMap[jetcollection.key])->inclusive_jets(jet_pt_min)); + std::shared_ptr CSeqBasePtr = result.emplace_clusterseq(jetparticles, jet_def, jetcollection.key); + std::vector pjets = sorted_by_pt(CSeqBasePtr->inclusive_jets(jet_pt_min)); // Add to the event, with b-tagging info" for (const FJNS::PseudoJet& pj : pjets) diff --git a/ColliderBit/include/gambit/ColliderBit/lhef2heputils.hpp b/ColliderBit/include/gambit/ColliderBit/lhef2heputils.hpp index 82d2dd98b0..95a28a1415 100644 --- a/ColliderBit/include/gambit/ColliderBit/lhef2heputils.hpp +++ b/ColliderBit/include/gambit/ColliderBit/lhef2heputils.hpp @@ -24,6 +24,7 @@ /// ********************************************* #include "gambit/cmake/cmake_variables.hpp" +#include "gambit/ColliderBit/Utils.hpp" #ifndef EXCLUDE_HEPMC @@ -33,6 +34,6 @@ namespace LHEF { class Reader; } /// Extract an LHE event as a HEPUtils::Event -void get_HEPUtils_event(const LHEF::Reader&, HEPUtils::Event&, double); +void get_HEPUtils_event(const LHEF::Reader&, HEPUtils::Event&, double, std::vector); -#endif \ No newline at end of file +#endif diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_0LEPStop_36invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_0LEPStop_36invfb.cpp index d0b6028697..a870ad6e9c 100644 --- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_0LEPStop_36invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_0LEPStop_36invfb.cpp @@ -279,7 +279,7 @@ namespace Gambit { const std::vector b = {0,10000.}; const std::vector c = {0.77}; // set b-tag efficiency to 77% HEPUtils::BinnedFn2D _eff2d(a,b,c); - for (const HEPUtils::Jet* jet : event->jets()) + for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) { bool hasTag=has_tag(_eff2d, fabs(jet->eta()), jet->pT()); if (jet->pT() > 20. && fabs(jet->eta()) < 2.8) diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_0LEP_139invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_0LEP_139invfb.cpp index 2bb39009e9..eb91ebffd0 100644 --- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_0LEP_139invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_0LEP_139invfb.cpp @@ -101,7 +101,7 @@ namespace Gambit // Get baseline jets /// @todo Drop b-tag if pT < 50 GeV or |eta| > 2.5? vector baselineJets; - for (const Jet* jet : event->jets()) + for (const Jet* jet : event->jets("antikt_R04")) { if (jet->pT() > 20. && jet->abseta() < 2.8) { diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_0LEP_13invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_0LEP_13invfb.cpp index 0b4206c9c0..61896dd221 100644 --- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_0LEP_13invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_0LEP_13invfb.cpp @@ -88,7 +88,7 @@ namespace Gambit // Get baseline jets /// @todo Drop b-tag if pT < 50 GeV or |eta| > 2.5? vector baselineJets; - for (const Jet* jet : event->jets()) + for (const Jet* jet : event->jets("antikt_R04")) if (jet->pT() > 20. && jet->abseta() < 2.8) { baselineJets.push_back(jet); diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_0LEP_36invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_0LEP_36invfb.cpp index 406a083714..dfee9f7788 100644 --- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_0LEP_36invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_0LEP_36invfb.cpp @@ -115,7 +115,7 @@ namespace Gambit { // Get baseline jets /// @todo Drop b-tag if pT < 50 GeV or |eta| > 2.5? vector baselineJets; - for (const Jet* jet : event->jets()) + for (const Jet* jet : event->jets("antikt_R04")) if (jet->pT() > 20. && jet->abseta() < 2.8) { baselineJets.push_back(jet); } diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_1LEPStop_36invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_1LEPStop_36invfb.cpp index 20b46b48e1..d40bda7b5c 100644 --- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_1LEPStop_36invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_1LEPStop_36invfb.cpp @@ -437,7 +437,7 @@ namespace Gambit { const std::vector b = {0,10000.}; const std::vector c = {0.77}; // set b-tag efficiency to 77% HEPUtils::BinnedFn2D _eff2d(a,b,c); - for (const HEPUtils::Jet* jet : event->jets()) + for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) { bool hasTag=has_tag(_eff2d, fabs(jet->eta()), jet->pT()); if (jet->pT() > 20. && fabs(jet->eta()) < 4.9) @@ -1423,4 +1423,4 @@ namespace Gambit { } } -#endif \ No newline at end of file +#endif diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_1Lep2b_139invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_1Lep2b_139invfb.cpp index 468b221f80..70227f866b 100644 --- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_1Lep2b_139invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_1Lep2b_139invfb.cpp @@ -180,7 +180,7 @@ namespace Gambit vector baselineJets; - for (const HEPUtils::Jet* jet : event->jets()) + for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) { if (jet->pT() > 20.0 && jet->abseta() < 4.5) { diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2BoostedBosons_139invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2BoostedBosons_139invfb.cpp index 2180e28395..5e07a98afe 100644 --- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2BoostedBosons_139invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2BoostedBosons_139invfb.cpp @@ -167,7 +167,7 @@ namespace Gambit // Look at jets to see if they fulfil criteria for fat jets vector fatJets; - for (const HEPUtils::Jet* jet : event->jets()) + for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) { // cout << jet->pT() << " " << jet->mass() << " Z-tag " << jet-tagged(23) << " W-tag " << jet->tagged(24) << " " << endl; if (jet->pT() > 200. && fabs(jet->eta()) < 2.0 && jet->mass() > 40.) @@ -223,7 +223,7 @@ namespace Gambit // double btag = 0.83; double cmisstag = 1/3.; double misstag = 1./33.; int nb = 0; - for ( const HEPUtils::Jet* jet : event->jets() ) + for ( const HEPUtils::Jet* jet : event->jets("antikt_R04") ) { // Tag b-jet if( jet->btag() ) nb++; diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2LEPStop_139invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2LEPStop_139invfb.cpp index c2b0f1f15d..201f9b0f24 100644 --- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2LEPStop_139invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2LEPStop_139invfb.cpp @@ -163,7 +163,7 @@ namespace Gambit // primary vertex (see paper) vector baselineJets; - for (const HEPUtils::Jet* jet : event->jets()) + for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) { if (jet->pT() > 20. && jet->abseta() < 2.8) { diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2LEPStop_36invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2LEPStop_36invfb.cpp index f309b7dacf..974df46699 100644 --- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2LEPStop_36invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2LEPStop_36invfb.cpp @@ -175,7 +175,7 @@ namespace Gambit { // Jets vector blJets; // Used for SR-2body and SR-3body vector baselineJets; // Used for SR-4body - for (const HEPUtils::Jet* jet : event->jets()) { + for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) { if (jet->pT() > 20. && fabs(jet->eta()) < 2.8) blJets.push_back(jet); if (jet->pT() > 20. && fabs(jet->eta()) < 2.8) baselineJets.push_back(jet); } diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2OSLEP_Z_139invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2OSLEP_Z_139invfb.cpp index a11ff77657..f275327166 100644 --- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2OSLEP_Z_139invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2OSLEP_Z_139invfb.cpp @@ -148,7 +148,7 @@ namespace Gambit // Jets with pT < 120 GeV and |η| < 2.8 have an efficiency of 90% // Mising: cut based on detector noise and non-collision backgrounds double jet_eff = 0.9; - for (const HEPUtils::Jet* jet : event->jets()) + for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) { if (jet->pT()>20. && jet->abseta()<2.8) if( (jet->pT() >= 120. || jet->abseta() >= 2.5) || random_bool(jet_eff) ) baselineJets.push_back(jet); diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2OSLEP_chargino_139invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2OSLEP_chargino_139invfb.cpp index bee304bfe5..af0c953e8b 100644 --- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2OSLEP_chargino_139invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2OSLEP_chargino_139invfb.cpp @@ -213,7 +213,7 @@ namespace Gambit // Jets vector candJets; - for (const HEPUtils::Jet* jet : event->jets()) { + for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) { if (jet->pT() > 20. && fabs(jet->eta()) < 2.4) candJets.push_back(jet); } diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2OSLEP_chargino_80invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2OSLEP_chargino_80invfb.cpp index 6426698cf9..f7cb98872d 100644 --- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2OSLEP_chargino_80invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2OSLEP_chargino_80invfb.cpp @@ -199,7 +199,7 @@ namespace Gambit // Jets vector candJets; - for (const HEPUtils::Jet* jet : event->jets()) { + for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) { if (jet->pT() > 20. && fabs(jet->eta()) < 2.5) candJets.push_back(jet); } diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2bMET_36invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2bMET_36invfb.cpp index 93d5939deb..63e54cfcb7 100644 --- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2bMET_36invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2bMET_36invfb.cpp @@ -357,7 +357,7 @@ namespace Gambit const std::vector b = {0,10000.}; const std::vector c = {0.77}; // set b-tag efficiency to 77% HEPUtils::BinnedFn2D _eff2d(a,b,c); - for (const HEPUtils::Jet* jet : event->jets()) + for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) { bool hasTag=has_tag(_eff2d, fabs(jet->eta()), jet->pT()); if (jet->pT() > 20. && fabs(jet->eta()) < 4.8) diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_3LEP_139invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_3LEP_139invfb.cpp index e68af59ba9..d5368bc4d4 100644 --- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_3LEP_139invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_3LEP_139invfb.cpp @@ -219,7 +219,7 @@ namespace Gambit // Baseline jets vector baselineJets; - for (const HEPUtils::Jet* jet : event->jets()) + for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) { if (jet->pT() > 20. && fabs(jet->eta()) < 4.5 ) baselineJets.push_back(jet); diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_3b_24invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_3b_24invfb.cpp index 957738b8d2..4fc7e6def0 100644 --- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_3b_24invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_3b_24invfb.cpp @@ -197,7 +197,7 @@ namespace Gambit { ATLAS::applyMuonEff(muons); vector candJets; - for (const HEPUtils::Jet* jet : event->jets()) { + for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) { if (jet->pT() > 25. && fabs(jet->eta()) < 2.5) candJets.push_back(jet); } diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_3b_36invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_3b_36invfb.cpp index d82199d666..43b157e0a0 100644 --- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_3b_36invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_3b_36invfb.cpp @@ -162,7 +162,7 @@ namespace Gambit { ATLAS::applyMuonEff(muons); vector candJets; - for (const HEPUtils::Jet* jet : event->jets()) { + for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) { if (jet->pT() > 20. && fabs(jet->eta()) < 2.8) candJets.push_back(jet); } diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_4LEP_139invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_4LEP_139invfb.cpp index 54e610942a..9950a14d55 100644 --- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_4LEP_139invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_4LEP_139invfb.cpp @@ -326,7 +326,7 @@ namespace Gambit // Since tau efficiencies are not applied as part of the BuckFast ATLAS sim we apply it here ATLAS::applyTauEfficiencyR2(baselineTaus); - for (const HEPUtils::Jet* jet : event->jets()) + for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) { if (jet->pT()>20. && jet->abseta()<2.8) baselineJets.push_back(jet); } diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_4LEP_36invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_4LEP_36invfb.cpp index a40e20e089..4150d66976 100644 --- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_4LEP_36invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_4LEP_36invfb.cpp @@ -224,7 +224,7 @@ namespace Gambit // Since tau efficiencies are not applied as part of the BuckFast ATLAS sim we apply it here ATLAS::applyTauEfficiencyR2(baselineTaus); - for (const HEPUtils::Jet* jet : event->jets()) + for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) { if (jet->pT()>20. && jet->abseta()<2.8) baselineJets.push_back(jet); } diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_MONOJET_139infb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_MONOJET_139infb.cpp index 0a1750e5c9..46e599ea3c 100644 --- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_MONOJET_139infb.cpp +++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_MONOJET_139infb.cpp @@ -141,7 +141,7 @@ namespace Gambit // Get jets (0.9 is to emulate the requirement of coming from a primary vertex) vector baselineJets; - for (const Jet* jet : event->jets()) + for (const Jet* jet : event->jets("antikt_R04")) { if ((jet->pT() > 30) && (jet->abseta() < 2.8)) { diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_MultiLEP_36invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_MultiLEP_36invfb.cpp index 1e6d422dd1..b4acd611f9 100644 --- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_MultiLEP_36invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_MultiLEP_36invfb.cpp @@ -187,7 +187,7 @@ namespace Gambit { ATLAS::applyMuonEff(baselineMuons); vector baselineJets; - for (const HEPUtils::Jet* jet : event->jets()) { + for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) { if (jet->pT()>20. && jet->abseta()<4.5)baselineJets.push_back(jet); } diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_MultiLEP_confnote_36invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_MultiLEP_confnote_36invfb.cpp index 261b9b45fa..5326185a2e 100644 --- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_MultiLEP_confnote_36invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_MultiLEP_confnote_36invfb.cpp @@ -176,7 +176,7 @@ namespace Gambit { ATLAS::applyMuonEff(baselineMuons); vector baselineJets; - for (const HEPUtils::Jet* jet : event->jets()) { + for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) { if (jet->pT()>20. && jet->abseta()<4.5)baselineJets.push_back(jet); } diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_MultiLEP_strong_139invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_MultiLEP_strong_139invfb.cpp index 7644e2fdbd..84c9489362 100644 --- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_MultiLEP_strong_139invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_MultiLEP_strong_139invfb.cpp @@ -179,7 +179,7 @@ namespace Gambit { // Get baseline jets /// @todo Drop b-tag if |eta| > 2.5? - for (const Jet* jet : event->jets()) + for (const Jet* jet : event->jets("antikt_R04")) { if (jet->pT() > 20. && jet->abseta() < 2.8) { diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_PhotonGGM_1Photon_139invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_PhotonGGM_1Photon_139invfb.cpp index e75cb51172..648f0d7691 100644 --- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_PhotonGGM_1Photon_139invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_PhotonGGM_1Photon_139invfb.cpp @@ -120,7 +120,7 @@ namespace Gambit // - pT > 30 // - |eta| < 2.5 vector signalJets; - for (const HEPUtils::Jet* jet : event->jets()) + for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) { if (jet->pT() > 30. && fabs(jet->eta()) < 2.5) { diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_PhotonGGM_36invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_PhotonGGM_36invfb.cpp index 5f54059dc5..47de7698a0 100644 --- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_PhotonGGM_36invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_PhotonGGM_36invfb.cpp @@ -183,7 +183,7 @@ namespace Gambit // Jets vector jets28; vector jets28_nophooverlap; - for (const HEPUtils::Jet* jet : event->jets()) + for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) { if (jet->pT() > 30. && fabs(jet->eta()) < 2.8) { diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_RJ3L_lowmass_36invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_RJ3L_lowmass_36invfb.cpp index 14b3c7f742..5dd2c42d55 100644 --- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_RJ3L_lowmass_36invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_RJ3L_lowmass_36invfb.cpp @@ -645,7 +645,7 @@ namespace Gambit { const std::vector b = {0,10000.}; const std::vector c = {0.77}; // set b-tag efficiency to 77% HEPUtils::BinnedFn2D _eff2d(a,b,c); - for (const HEPUtils::Jet* jet : event->jets()) { + for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) { bool hasTag=has_tag(_eff2d, jet->abseta(), jet->pT()); if (jet->pT() > 20. && fabs(jet->eta()) < 2.4) { if(jet->btag() && hasTag){ diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_ZGammaGrav_CONFNOTE_80invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_ZGammaGrav_CONFNOTE_80invfb.cpp index d5a0dce055..f521963413 100644 --- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_ZGammaGrav_CONFNOTE_80invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_ZGammaGrav_CONFNOTE_80invfb.cpp @@ -67,7 +67,7 @@ namespace Gambit { // Jets JetPtrs jets; - for (const Jet* j : event->jets()) + for (const Jet* j : event->jets("antikt_R04")) if (j->pT() > 20. && j->absrap() < 4.4) jets.push_back(j); diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_7TeV_1OR2LEPStop_4_7invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_7TeV_1OR2LEPStop_4_7invfb.cpp index 0d8dbe10f6..ec043c17ee 100644 --- a/ColliderBit/src/analyses/Analysis_ATLAS_7TeV_1OR2LEPStop_4_7invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_ATLAS_7TeV_1OR2LEPStop_4_7invfb.cpp @@ -211,7 +211,7 @@ Analysis_ATLAS_7TeV_1OR2LEPStop_4_7invfb() incrementCut(Total_events); std::vector electrons = event->electrons(); std::vector muons = event->muons(); - std::vector jets = event->jets(); + std::vector jets = event->jets("antikt_R04"); electrons = AnalysisUtil::filterPtEta(electrons, 20, 2.47); muons = AnalysisUtil::filterPtEta(muons, 10, 2.4); diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_7TeV_2LEPStop_4_7invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_7TeV_2LEPStop_4_7invfb.cpp index f48ed86d47..fc63bbbeb9 100644 --- a/ColliderBit/src/analyses/Analysis_ATLAS_7TeV_2LEPStop_4_7invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_ATLAS_7TeV_2LEPStop_4_7invfb.cpp @@ -49,7 +49,7 @@ namespace Gambit // get the jets and leptons filtered by their pt and eta requirements electrons = AnalysisUtil::filterPtEta(electrons, 10, 2.47); muons = AnalysisUtil::filterPtEta(muons, 10, 2.4); - std::vector jets = AnalysisUtil::filterPtEta(event->jets(), 20, 4.5); + std::vector jets = AnalysisUtil::filterPtEta(event->jets("antikt_R04"), 20, 4.5); // check if any of the triggers were triggered bool eeTrigger = AnalysisUtil::isMultipleParticleTriggered(electrons, {17, 17}); diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_0LEPStop_20invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_0LEPStop_20invfb.cpp index d8de98095c..83bb54ca65 100644 --- a/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_0LEPStop_20invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_0LEPStop_20invfb.cpp @@ -157,7 +157,7 @@ namespace Gambit { const std::vector b = {0,10000.}; const std::vector c = {0.7}; HEPUtils::BinnedFn2D _eff2d(a,b,c); - for (const HEPUtils::Jet* jet : event->jets()) { + for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) { if (jet->pT() > 20. && fabs(jet->eta()) < 4.5) baselineJets.push_back(jet); bool hasTag=has_tag(_eff2d, fabs(jet->eta()), jet->pT()); if (jet->pT() > 20. && fabs(jet->eta()) < 4.5) { diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_0LEP_20invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_0LEP_20invfb.cpp index 0192565cd6..2109e92358 100644 --- a/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_0LEP_20invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_0LEP_20invfb.cpp @@ -77,7 +77,7 @@ namespace Gambit { ATLAS::applyMuonEff(baselineMuons); vector baselineJets; - for (const HEPUtils::Jet* jet : event->jets()) { + for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) { if (jet->pT() > 20. && fabs(jet->eta()) < 4.5) baselineJets.push_back(jet); } diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_1LEPStop_20invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_1LEPStop_20invfb.cpp index 7ceba10beb..45c960527f 100644 --- a/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_1LEPStop_20invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_1LEPStop_20invfb.cpp @@ -232,7 +232,7 @@ namespace Gambit { // Get b jets with efficiency and mistag (fake) rates vector baselineJets, bJets; // trueBJets; //for debugging - for (const HEPUtils::Jet* jet : event->jets()) { + for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) { if (jet->pT() > 20. && jet->abseta() < 10.0) baselineJets.push_back(jet); if (jet->abseta() < 2.5 && jet->pT() > 25.) { if ((jet->btag() && Random::draw() < 0.75) || (!jet->btag() && Random::draw() < 0.02)) bJets.push_back(jet); diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_1LEPbb_20invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_1LEPbb_20invfb.cpp index 0c593b452a..dcb7d864b7 100644 --- a/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_1LEPbb_20invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_1LEPbb_20invfb.cpp @@ -109,7 +109,7 @@ namespace Gambit { ATLAS::applyMuonEff(baselineMuons); vector baselineJets; - for (const HEPUtils::Jet* jet : event->jets()) { + for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) { if (jet->pT()>20. && fabs(jet->eta())<4.5) baselineJets.push_back(jet); } diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_2LEPEW_20invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_2LEPEW_20invfb.cpp index dbe5c070b4..cb9c2e0ad2 100644 --- a/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_2LEPEW_20invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_2LEPEW_20invfb.cpp @@ -287,7 +287,7 @@ namespace Gambit ATLAS::applyMuonEff(signalMuons); vector signalJets; - for (const HEPUtils::Jet* jet : event->jets()) + for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) { if (jet->pT() > 20. && fabs(jet->eta()) < 4.5) signalJets.push_back(jet); //if(jet->btag() && fabs(jet->eta()) < 2.5 && jet->pT() > 20.) bJets.push_back(jet); diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_2LEPStop_20invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_2LEPStop_20invfb.cpp index 2c681d4b2d..886bde501b 100644 --- a/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_2LEPStop_20invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_2LEPStop_20invfb.cpp @@ -100,7 +100,7 @@ namespace Gambit { vector baselineJets; vector bJets; vector trueBJets; //for debugging - for (const HEPUtils::Jet* jet : event->jets()) { + for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) { if (jet->pT() > 20. && fabs(jet->eta()) < 2.5) baselineJets.push_back(jet); if(jet->btag() && fabs(jet->eta()) < 2.5 && jet->pT() > 20.) bJets.push_back(jet); } diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_2bStop_20invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_2bStop_20invfb.cpp index 18b5eaba81..367131cf37 100644 --- a/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_2bStop_20invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_2bStop_20invfb.cpp @@ -123,7 +123,7 @@ namespace Gambit vector bJets; vector trueBJets; //for debugging - for (const HEPUtils::Jet* jet : event->jets()) + for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) { if (jet->pT() > 20. && fabs(jet->eta()) < 4.9) baselineJets.push_back(jet); } diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_3LEPEW_20invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_3LEPEW_20invfb.cpp index 12f006036e..58192f8470 100644 --- a/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_3LEPEW_20invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_3LEPEW_20invfb.cpp @@ -215,7 +215,7 @@ namespace Gambit { vector signalJets; vector bJets; - for (const HEPUtils::Jet* jet : event->jets()) { + for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) { if (jet->pT() > 20. && fabs(jet->eta()) < 2.5) signalJets.push_back(jet); //if(jet->btag() && fabs(jet->eta()) < 2.5 && jet->pT() > 20.) bJets.push_back(jet); } diff --git a/ColliderBit/src/analyses/Analysis_CMS_13TeV_0LEP_137invfb.cpp b/ColliderBit/src/analyses/Analysis_CMS_13TeV_0LEP_137invfb.cpp index 265258625e..9b150cd6d9 100644 --- a/ColliderBit/src/analyses/Analysis_CMS_13TeV_0LEP_137invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_CMS_13TeV_0LEP_137invfb.cpp @@ -65,7 +65,7 @@ namespace Gambit // Get jets vector jets24, jets50; - for (const Jet* jet : event->jets()) + for (const Jet* jet : event->jets("antikt_R04")) { if (jet->pT() < 30) continue; if (jet->abseta() < 2.4) jets24.push_back(jet); diff --git a/ColliderBit/src/analyses/Analysis_CMS_13TeV_0LEP_13invfb.cpp b/ColliderBit/src/analyses/Analysis_CMS_13TeV_0LEP_13invfb.cpp index 985c2a5bac..6f0d7ca8ef 100644 --- a/ColliderBit/src/analyses/Analysis_CMS_13TeV_0LEP_13invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_CMS_13TeV_0LEP_13invfb.cpp @@ -57,7 +57,7 @@ namespace Gambit { // Get baseline jets vector jets24, jets50; - for (const Jet* jet : event->jets()) { + for (const Jet* jet : event->jets("antikt_R04")) { if (jet->pT() < 30) continue; if (jet->abseta() < 2.4) jets24.push_back(jet); if (jet->abseta() < 5.0) jets50.push_back(jet); diff --git a/ColliderBit/src/analyses/Analysis_CMS_13TeV_0LEP_36invfb.cpp b/ColliderBit/src/analyses/Analysis_CMS_13TeV_0LEP_36invfb.cpp index 2f79a8fd7a..7d0b853edf 100644 --- a/ColliderBit/src/analyses/Analysis_CMS_13TeV_0LEP_36invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_CMS_13TeV_0LEP_36invfb.cpp @@ -57,7 +57,7 @@ namespace Gambit { // Get baseline jets vector jets24, jets50; - for (const Jet* jet : event->jets()) { + for (const Jet* jet : event->jets("antikt_R04")) { if (jet->pT() < 30) continue; if (jet->abseta() < 2.4) jets24.push_back(jet); if (jet->abseta() < 5.0) jets50.push_back(jet); diff --git a/ColliderBit/src/analyses/Analysis_CMS_13TeV_1LEPStop_36invfb.cpp b/ColliderBit/src/analyses/Analysis_CMS_13TeV_1LEPStop_36invfb.cpp index eea60e0dac..b097fe873b 100644 --- a/ColliderBit/src/analyses/Analysis_CMS_13TeV_1LEPStop_36invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_CMS_13TeV_1LEPStop_36invfb.cpp @@ -154,7 +154,7 @@ namespace Gambit { // Jets vector baselineJets; vector fullJets; - for (const HEPUtils::Jet* jet : event->jets()) { + for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) { if (jet->pT() > 30. && jet->abseta() < 2.4) baselineJets.push_back(jet); if (jet->abseta() < 5.0) fullJets.push_back(jet); } diff --git a/ColliderBit/src/analyses/Analysis_CMS_13TeV_1LEPbb_36invfb.cpp b/ColliderBit/src/analyses/Analysis_CMS_13TeV_1LEPbb_36invfb.cpp index 462e5af12a..4f1e309e23 100644 --- a/ColliderBit/src/analyses/Analysis_CMS_13TeV_1LEPbb_36invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_CMS_13TeV_1LEPbb_36invfb.cpp @@ -119,7 +119,7 @@ namespace Gambit } vector baselineJets; - for (const HEPUtils::Jet* jet : event->jets()) + for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) { if (jet->pT()>25. &&fabs(jet->eta())<2.4)baselineJets.push_back(jet); } diff --git a/ColliderBit/src/analyses/Analysis_CMS_13TeV_1Photon1Lepton_36invfb.cpp b/ColliderBit/src/analyses/Analysis_CMS_13TeV_1Photon1Lepton_36invfb.cpp index 614fd45578..922e185e62 100644 --- a/ColliderBit/src/analyses/Analysis_CMS_13TeV_1Photon1Lepton_36invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_CMS_13TeV_1Photon1Lepton_36invfb.cpp @@ -270,7 +270,7 @@ namespace Gambit { // Jets vector jets; - for (const HEPUtils::Jet* jet : event->jets()) + for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) { if (jet->pT()>30. && jet->abseta()<2.5) jets.push_back(jet); } diff --git a/ColliderBit/src/analyses/Analysis_CMS_13TeV_2LEPStop_36invfb.cpp b/ColliderBit/src/analyses/Analysis_CMS_13TeV_2LEPStop_36invfb.cpp index a5062c83fa..4f295d8588 100644 --- a/ColliderBit/src/analyses/Analysis_CMS_13TeV_2LEPStop_36invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_CMS_13TeV_2LEPStop_36invfb.cpp @@ -178,7 +178,7 @@ namespace Gambit { ATLAS::applyLooseIDElectronSelectionR2(baselineElectrons); // Jets vector baselineJets; - for (const HEPUtils::Jet* jet : event->jets()) { + for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) { if (jet->pT() > 30. && fabs(jet->eta()) < 2.4) baselineJets.push_back(jet); } diff --git a/ColliderBit/src/analyses/Analysis_CMS_13TeV_2LEPsoft_36invfb.cpp b/ColliderBit/src/analyses/Analysis_CMS_13TeV_2LEPsoft_36invfb.cpp index 011e00d72d..57705b9ba8 100644 --- a/ColliderBit/src/analyses/Analysis_CMS_13TeV_2LEPsoft_36invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_CMS_13TeV_2LEPsoft_36invfb.cpp @@ -168,7 +168,7 @@ namespace Gambit { else if (met > 300. && muon->pT()>3.5 && muon->pT()<30. && fabs(muon->eta())<2.4 && isMu) signalMuons.push_back(muon); } - for (const HEPUtils::Jet* jet : event->jets()) { + for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) { if (jet->pT()>25. && fabs(jet->eta())<2.4) { signalJets.push_back(jet); if (jet->btag())signalBJets.push_back(jet); diff --git a/ColliderBit/src/analyses/Analysis_CMS_13TeV_2OSLEP_137invfb.cpp b/ColliderBit/src/analyses/Analysis_CMS_13TeV_2OSLEP_137invfb.cpp index 826caeaed0..8875551ee2 100644 --- a/ColliderBit/src/analyses/Analysis_CMS_13TeV_2OSLEP_137invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_CMS_13TeV_2OSLEP_137invfb.cpp @@ -200,7 +200,7 @@ namespace Gambit // Baseline jets vector baselineJets; - for (const HEPUtils::Jet* jet : event->jets()) + for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) { if (jet->pT()>25. && fabs(jet->eta())<2.4) baselineJets.push_back(jet); } diff --git a/ColliderBit/src/analyses/Analysis_CMS_13TeV_2OSLEP_36invfb.cpp b/ColliderBit/src/analyses/Analysis_CMS_13TeV_2OSLEP_36invfb.cpp index 87181c3812..753a2541b7 100644 --- a/ColliderBit/src/analyses/Analysis_CMS_13TeV_2OSLEP_36invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_CMS_13TeV_2OSLEP_36invfb.cpp @@ -139,7 +139,7 @@ namespace Gambit { // Baseline jets vector baselineJets; - for (const HEPUtils::Jet* jet : event->jets()) + for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) { // We use 25 GeV rather than 35 GeV // if (jet->pT()>35. &&fabs(jet->eta())<2.4) baselineJets.push_back(jet); diff --git a/ColliderBit/src/analyses/Analysis_CMS_13TeV_2OSLEP_chargino_stop_36invfb.cpp b/ColliderBit/src/analyses/Analysis_CMS_13TeV_2OSLEP_chargino_stop_36invfb.cpp index 66aa38e9e5..604b4d6274 100644 --- a/ColliderBit/src/analyses/Analysis_CMS_13TeV_2OSLEP_chargino_stop_36invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_CMS_13TeV_2OSLEP_chargino_stop_36invfb.cpp @@ -265,7 +265,7 @@ namespace Gambit { // Baseline jets vector baselineJets; - for (const HEPUtils::Jet* jet : event->jets()) + for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) { if (jet->pT()>20. &&fabs(jet->eta())<2.4) baselineJets.push_back(jet); } diff --git a/ColliderBit/src/analyses/Analysis_CMS_13TeV_2OSLEP_confnote_36invfb.cpp b/ColliderBit/src/analyses/Analysis_CMS_13TeV_2OSLEP_confnote_36invfb.cpp index fb3e049aaf..62a9a7a4b7 100644 --- a/ColliderBit/src/analyses/Analysis_CMS_13TeV_2OSLEP_confnote_36invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_CMS_13TeV_2OSLEP_confnote_36invfb.cpp @@ -97,7 +97,7 @@ namespace Gambit { } vector baselineJets; - for (const HEPUtils::Jet* jet : event->jets()) { + for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) { if (jet->pT()>35. &&fabs(jet->eta())<2.4)baselineJets.push_back(jet); } diff --git a/ColliderBit/src/analyses/Analysis_CMS_13TeV_2Photon_GMSB_36invfb.cpp b/ColliderBit/src/analyses/Analysis_CMS_13TeV_2Photon_GMSB_36invfb.cpp index 1f4a282f56..4549ce53c5 100644 --- a/ColliderBit/src/analyses/Analysis_CMS_13TeV_2Photon_GMSB_36invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_CMS_13TeV_2Photon_GMSB_36invfb.cpp @@ -207,7 +207,7 @@ namespace Gambit { // Jets vector jets; - for (const HEPUtils::Jet* jet : event->jets()) { + for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) { // No info on baseline jet cuts in the paper, so for now we'll // apply an|eta| cut for HCAL coverage and a loose jet pT cut if (jet->pT()>10. && jet->abseta()<3.0) jets.push_back(jet); diff --git a/ColliderBit/src/analyses/Analysis_CMS_13TeV_2SSLEP_Stop_137invfb.cpp b/ColliderBit/src/analyses/Analysis_CMS_13TeV_2SSLEP_Stop_137invfb.cpp index 97f875a317..e74cab7cbb 100644 --- a/ColliderBit/src/analyses/Analysis_CMS_13TeV_2SSLEP_Stop_137invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_CMS_13TeV_2SSLEP_Stop_137invfb.cpp @@ -296,7 +296,7 @@ namespace Gambit { // Jets double HT = 0.; vector candJets; - for (const HEPUtils::Jet* jet : event->jets()) { + for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) { if (jet->pT() > 25. && fabs(jet->eta()) < 2.4){ HT += jet->pT(); candJets.push_back(jet); diff --git a/ColliderBit/src/analyses/Analysis_CMS_13TeV_2SSLEP_Stop_36invfb.cpp b/ColliderBit/src/analyses/Analysis_CMS_13TeV_2SSLEP_Stop_36invfb.cpp index 5418297800..a7e4f8ce1f 100644 --- a/ColliderBit/src/analyses/Analysis_CMS_13TeV_2SSLEP_Stop_36invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_CMS_13TeV_2SSLEP_Stop_36invfb.cpp @@ -293,7 +293,7 @@ namespace Gambit { // Jets double HT = 0.; vector candJets; - for (const HEPUtils::Jet* jet : event->jets()) { + for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) { if (jet->pT() > 25. && fabs(jet->eta()) < 2.4){ HT += jet->pT(); candJets.push_back(jet); diff --git a/ColliderBit/src/analyses/Analysis_CMS_13TeV_MONOJET_36invfb.cpp b/ColliderBit/src/analyses/Analysis_CMS_13TeV_MONOJET_36invfb.cpp index 1a9388fde1..650820373a 100644 --- a/ColliderBit/src/analyses/Analysis_CMS_13TeV_MONOJET_36invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_CMS_13TeV_MONOJET_36invfb.cpp @@ -91,7 +91,7 @@ namespace Gambit { // Get jets vector jets4; - for (const Jet* jet : event->jets()) + for (const Jet* jet : event->jets("antikt_R04")) if (jet->pT() > 20) jets4.push_back(jet); // Veto if there are any b-tagged jets (reduce top background) diff --git a/ColliderBit/src/analyses/Analysis_CMS_13TeV_MultiLEP_137invfb.cpp b/ColliderBit/src/analyses/Analysis_CMS_13TeV_MultiLEP_137invfb.cpp index fa8d6a9b73..ab3d3f7bc7 100644 --- a/ColliderBit/src/analyses/Analysis_CMS_13TeV_MultiLEP_137invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_CMS_13TeV_MultiLEP_137invfb.cpp @@ -184,7 +184,7 @@ namespace Gambit // Baseline jets std::vector baselineJets; - for (const HEPUtils::Jet* jet : event->jets()) + for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) { if (jet->pT()>25. &&fabs(jet->eta())<2.4) baselineJets.push_back(jet); diff --git a/ColliderBit/src/analyses/Analysis_CMS_13TeV_MultiLEP_36invfb.cpp b/ColliderBit/src/analyses/Analysis_CMS_13TeV_MultiLEP_36invfb.cpp index 1c9351564a..50eb14f3a0 100644 --- a/ColliderBit/src/analyses/Analysis_CMS_13TeV_MultiLEP_36invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_CMS_13TeV_MultiLEP_36invfb.cpp @@ -193,7 +193,7 @@ namespace Gambit } vector baselineJets; - for (const HEPUtils::Jet* jet : event->jets()) + for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) { if (jet->pT()>25. &&fabs(jet->eta())<2.4)baselineJets.push_back(jet); } diff --git a/ColliderBit/src/analyses/Analysis_CMS_13TeV_MultiLEP_Full_36invfb.cpp b/ColliderBit/src/analyses/Analysis_CMS_13TeV_MultiLEP_Full_36invfb.cpp index 24b5f195ef..1885609fc5 100644 --- a/ColliderBit/src/analyses/Analysis_CMS_13TeV_MultiLEP_Full_36invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_CMS_13TeV_MultiLEP_Full_36invfb.cpp @@ -312,7 +312,7 @@ namespace Gambit { } vector baselineJets; - for (const HEPUtils::Jet* jet : event->jets()) { + for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) { if (jet->pT()>25. &&fabs(jet->eta())<2.4)baselineJets.push_back(jet); } diff --git a/ColliderBit/src/analyses/Analysis_CMS_13TeV_Photon_GMSB_36invfb.cpp b/ColliderBit/src/analyses/Analysis_CMS_13TeV_Photon_GMSB_36invfb.cpp index cd1cfcf68a..3b5e209a19 100644 --- a/ColliderBit/src/analyses/Analysis_CMS_13TeV_Photon_GMSB_36invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_CMS_13TeV_Photon_GMSB_36invfb.cpp @@ -87,7 +87,7 @@ namespace Gambit { // jets vector Jets; - for (const HEPUtils::Jet* jet : event->jets()) + for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) { if (jet->pT()>30. &&fabs(jet->eta())<3.0) Jets.push_back(jet); } diff --git a/ColliderBit/src/analyses/Analysis_CMS_8TeV_1LEPDMTOP_20invfb.cpp b/ColliderBit/src/analyses/Analysis_CMS_8TeV_1LEPDMTOP_20invfb.cpp index b10a199891..92fbd513f9 100644 --- a/ColliderBit/src/analyses/Analysis_CMS_8TeV_1LEPDMTOP_20invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_CMS_8TeV_1LEPDMTOP_20invfb.cpp @@ -114,7 +114,7 @@ namespace Gambit { const std::vector c = {0.60}; HEPUtils::BinnedFn2D _eff2d(a,b,c); - for (const HEPUtils::Jet* jet : event->jets()) { + for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) { if (jet->pT() > 30. && fabs(jet->eta()) < 4.0) { baselineJets.push_back(jet); //LorentzVector j1 (jet->mom().px(),jet->mom().py(),jet->mom().pz(),jet->mom().E()) ; diff --git a/ColliderBit/src/analyses/Analysis_CMS_8TeV_2LEPDMTOP_20invfb.cpp b/ColliderBit/src/analyses/Analysis_CMS_8TeV_2LEPDMTOP_20invfb.cpp index 56871097d5..77619a6910 100644 --- a/ColliderBit/src/analyses/Analysis_CMS_8TeV_2LEPDMTOP_20invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_CMS_8TeV_2LEPDMTOP_20invfb.cpp @@ -111,7 +111,7 @@ namespace Gambit { const std::vector c = {0.60}; HEPUtils::BinnedFn2D _eff2d(a,b,c); - for (const HEPUtils::Jet* jet : event->jets()) { + for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) { if (jet->pT() > 30. && fabs(jet->eta()) < 5.0) { baselineJets.push_back(jet); //LorentzVector j1 (jet->mom().px(),jet->mom().py(),jet->mom().pz(),jet->mom().E()) ; diff --git a/ColliderBit/src/analyses/Analysis_CMS_8TeV_MONOJET_20invfb.cpp b/ColliderBit/src/analyses/Analysis_CMS_8TeV_MONOJET_20invfb.cpp index 4b3fe973ca..b577d8d2b7 100644 --- a/ColliderBit/src/analyses/Analysis_CMS_8TeV_MONOJET_20invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_CMS_8TeV_MONOJET_20invfb.cpp @@ -110,7 +110,7 @@ namespace Gambit { vector baselineJets; vector jets; - for (const HEPUtils::Jet* jet : event->jets()) { + for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) { if (jet->pT() > 30. && fabs(jet->eta()) < 4.5) { baselineJets.push_back(jet); } diff --git a/ColliderBit/src/analyses/Analysis_CMS_8TeV_MultiLEP_20invfb.cpp b/ColliderBit/src/analyses/Analysis_CMS_8TeV_MultiLEP_20invfb.cpp index 3c1e31b638..0780ee2dd2 100644 --- a/ColliderBit/src/analyses/Analysis_CMS_8TeV_MultiLEP_20invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_CMS_8TeV_MultiLEP_20invfb.cpp @@ -477,7 +477,7 @@ namespace Gambit { // - jets vector signalJets; vector signalBjets; - for (const HEPUtils::Jet* jet : event->jets()) { + for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) { if (jet->pT() > 30. && fabs(jet->eta()) < 2.5) signalJets.push_back(jet); if(jet->btag() && fabs(jet->eta()) < 2.5 && jet->pT() > 30.) signalBjets.push_back(jet); } diff --git a/ColliderBit/src/analyses/Analysis_Minimum.cpp b/ColliderBit/src/analyses/Analysis_Minimum.cpp index 4044c022a3..683332bfa4 100644 --- a/ColliderBit/src/analyses/Analysis_Minimum.cpp +++ b/ColliderBit/src/analyses/Analysis_Minimum.cpp @@ -56,7 +56,7 @@ namespace Gambit { // Baseline jets vector baselineJets; - for (const HEPUtils::Jet* jet : event->jets()) { + for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) { if (jet->pT() > 20. && fabs(jet->eta()) < 4.5) baselineJets.push_back(jet); } diff --git a/ColliderBit/src/getLHEvent.cpp b/ColliderBit/src/getLHEvent.cpp index c21dc95d30..535649e97d 100644 --- a/ColliderBit/src/getLHEvent.cpp +++ b/ColliderBit/src/getLHEvent.cpp @@ -55,6 +55,50 @@ namespace Gambit } static LHEF::Reader lhe(lhef_filename); + // Get all jet collection settings + str jetcollection_taus; + std::vector all_jet_collection_settings = {}; // Initialise with no collections, if no jet collection setting specified, then should use the basecollider's initialised jet collection + if (runOptions->hasKey((*Dep::RunMC).current_collider())) + { + YAML::Node colNode = runOptions->getValue((*Dep::RunMC).current_collider()); + Options colOptions(colNode); + + // Fill the jet collection settings + if (colOptions.hasKey("jet_collections")) + { + YAML::Node jetcollectionNode = colOptions.getValue("jet_collections"); + Options jetcollectionOptions(jetcollectionNode); + + str algorithm; + double R; + str recombination_scheme; + str strategy; + std::vector jetcollections = jetcollectionOptions.getNames(); + + for (str key : jetcollections) + { + algorithm = jetcollectionOptions.getValueOrDef("antikt", "algorithm"); + R = jetcollectionOptions.getValueOrDef(0.4, "R"); + recombination_scheme = jetcollectionOptions.getValueOrDef("E_scheme", "recombination_scheme"); + strategy = jetcollectionOptions.getValueOrDef("Best", "strategy"); + + all_jet_collection_settings.push_back({key, algorithm, R, recombination_scheme, strategy}); + } + + jetcollection_taus = colOptions.getValueOrDef("antikt_R04", "jetcollection_taus"); + // Throw an error if the jetcollection_taus setting is not given and not using the antikt_R04 collection + if (std::find(jetcollections.begin(), jetcollections.end(), jetcollection_taus) == jetcollections.end()) + { + ColliderBit_error().raise(LOCAL_INFO,"Please provide the jetcollection_taus setting for jet collections if not using antikt_R04."); + } + } + } + else + { + all_jet_collection_settings = {{"antikt_R04", "antikt", 0.4, "E_scheme", "Best"}}; + jetcollection_taus = "antikt_R04"; + } + // Don't do anything during special iterations if (*Loop::iteration < 0) return; @@ -62,7 +106,7 @@ namespace Gambit bool event_retrieved = true; #pragma omp critical (reading_LHEvent) { - if (lhe.readEvent()) get_HEPUtils_event(lhe, result, jet_pt_min); + if (lhe.readEvent()) get_HEPUtils_event(lhe, result, jet_pt_min, all_jet_collection_settings); else event_retrieved = false; } if (not event_retrieved) diff --git a/ColliderBit/src/lhef2heputils.cpp b/ColliderBit/src/lhef2heputils.cpp index b267d279db..ffa8347fe0 100644 --- a/ColliderBit/src/lhef2heputils.cpp +++ b/ColliderBit/src/lhef2heputils.cpp @@ -44,8 +44,10 @@ using namespace std; using namespace HEPUtils; using namespace FJNS; +using namespace Gambit::ColliderBit; + /// Extract an LHE event as a HEPUtils::Event -void get_HEPUtils_event(const LHEF::Reader& lhe, Event& evt, double jet_pt_min) +void get_HEPUtils_event(const LHEF::Reader& lhe, Event& evt, double jet_pt_min, std::vector all_jet_collection_settings) { P4 vmet; @@ -99,31 +101,36 @@ void get_HEPUtils_event(const LHEF::Reader& lhe, Event& evt, double jet_pt_min) vmet.setM(0); evt.set_missingmom(vmet); - // Jets - vector jets = get_jets(jetparticles, 0.4, jet_pt_min); - - for (const PseudoJet& pj : jets) + // Jet Finding + for (jet_collection_settings jetcollection : all_jet_collection_settings) { - bool hasC = false, hasB = false; - /// @todo Bug in HEPUtils::get_jets means that constituent info is lost for now... - // for (const PseudoJet& c : pj.constituents()) { - // if (c.user_index() == 4) hasC = true; - // if (c.user_index() == 5) hasB = true; - // } - evt.add_jet(new Jet(mk_p4(pj), hasB, hasC)); - } - #ifdef COLLIDERBIT_DEBUG - // Print event summary - cout << " MET = " << evt.met() << " GeV" << endl; - cout << " #e = " << evt.electrons().size() << endl; - cout << " #mu = " << evt.muons().size() << endl; - cout << " #tau = " << evt.taus().size() << endl; - cout << " #jet = " << evt.jets().size() << endl; - cout << " #pho = " << evt.photons().size() << endl; - cout << endl; - #endif + // @todo get_jets function could accept a more general jet definition + vector jets = get_jets(jetparticles, 0.4, jet_pt_min, FJalgorithm_map(jetcollection.algorithm)); + + for (const PseudoJet& pj : jets) + { + bool hasC = false, hasB = false; + /// @todo Bug in HEPUtils::get_jets means that constituent info is lost for now... + // for (const PseudoJet& c : pj.constituents()) { + // if (c.user_index() == 4) hasC = true; + // if (c.user_index() == 5) hasB = true; + // } + evt.add_jet(new Jet(mk_p4(pj), hasB, hasC), jetcollection.key); + } + + #ifdef COLLIDERBIT_DEBUG + // Print event summary + cout << " MET = " << evt.met() << " GeV" << endl; + cout << " #e = " << evt.electrons().size() << endl; + cout << " #mu = " << evt.muons().size() << endl; + cout << " #tau = " << evt.taus().size() << endl; + cout << " #jet = " << evt.jets().size() << endl; + cout << " #pho = " << evt.photons().size() << endl; + cout << endl; + #endif + } } -#endif \ No newline at end of file +#endif diff --git a/contrib/heputils/include/HEPUtils/Event.h b/contrib/heputils/include/HEPUtils/Event.h index e7321bb52f..9b33158f99 100644 --- a/contrib/heputils/include/HEPUtils/Event.h +++ b/contrib/heputils/include/HEPUtils/Event.h @@ -1,7 +1,7 @@ // -*- C++ -*- // // This file is part of HEPUtils -- https://gitlab.com/hepcedar/heputils/ -// Copyright (C) 2013-2022 Andy Buckley +// Copyright (C) 2013-2023 Andy Buckley // // Embedding of HEPUtils code in other projects is permitted provided this // notice is retained and the HEPUtils namespace and include path are changed. @@ -11,6 +11,7 @@ #include "HEPUtils/Particle.h" #include "HEPUtils/Jet.h" #include +#include #include namespace HEPUtils { @@ -36,6 +37,14 @@ namespace HEPUtils { /// Jets collection(s) (mutable to allow sorting) mutable std::map> _jets; + /// Typedef for the generic cluster-sequence type + using CSeqBase = FJNS::ClusterSequence; + /// Typedef for a smart ptr to the generic cluster-sequence type + using CSeqBasePtr = std::shared_ptr; + + /// Hold the cluster sequences corresponding to jets, to keep them alive + std::map _cseqs; + /// Missing momentum vector P4 _pmiss; @@ -60,6 +69,7 @@ namespace HEPUtils { _taus = e._taus; // _jets = e._jets; + _cseqs = e._cseqs; _pmiss = e._pmiss; } @@ -109,10 +119,10 @@ namespace HEPUtils { for (size_t i = 0; i < ps.size(); ++i) { e.add_particle(new Particle(*ps[i])); } - for ( auto jetcollection : _jets ) { - const std::vector js = jets(jetcollection.first); + for (const auto& kv : _jets ) { + const std::vector js = jets(kv.first); for (size_t i = 0; i < js.size(); ++i) { - e.add_jet(new Jet(*js[i]),jetcollection.first); + e.add_jet(new Jet(*js[i]), kv.first); } } e._pmiss = _pmiss; @@ -139,17 +149,15 @@ namespace HEPUtils { _taus.clear(); // Jets - for (auto& js : _jets) { - for (const Jet* j : js.second) delete j; - } + for (const std::string& jc : jet_collections()) clear_jets(jc); _jets.clear(); + _cseqs.clear(); // MET _pmiss.clear(); } - /// @name Weights /// @{ @@ -357,81 +365,122 @@ namespace HEPUtils { } /// @brief Get a jet collection (not including charged leptons or photons) - const std::vector& jets(const std::string& key="CANONICAL") const { - // Throw an error if the user does not pass a more descriptive key - if (key == "CANONICAL") - { - throw std::runtime_error("Please supply a key for the jet collection."); - } + const std::vector& jets(const std::string& key) const { return _get_jets(key); } /// @brief Get a jet collection (not including charged leptons or photons) (non-const) - std::vector& jets(const std::string& key="CANONICAL") { - // Throw an error if the user does not pass a more descriptive key - if (key == "CANONICAL") - { - throw std::runtime_error("Please supply a key for the jet collection."); - } + std::vector& jets(const std::string& key) { return mkunconst(_get_jets(key)); } - /// Get the list of jet collection names + + /// Get the list of jet-collection names std::vector jet_collections() { - std::vector collection_names; - for (auto jetcollection : _jets) - { - collection_names.push_back(jetcollection.first); - } - return collection_names; + std::vector rtn; + for (const auto& kv : _jets) rtn.push_back(kv.first); + return rtn; } - - /// @brief Set the jets collection + + /// @brief Set a jet collection /// /// @warning The Jets should be new'd; Event will take ownership. + /// /// @todo "Lock" at some point so that jet finding etc. only get done once - void set_jets(const std::vector& jets, const std::string& key="CANONICAL") { + /// + /// @note This resets the cluster sequence, but as a shared_ptr is used for storage, + /// any existing shared_ptr links to the previous one will keep it alive. + void set_jets(const std::vector& jets, const std::string& key) { _jets[key] = jets; std::sort(_jets[key].begin(), _jets[key].end(), _cmpPtDescPtr); + set_clusterseq(jets.front()->clusterseq(), key); } // /// @brief Set the jets collection (non-const input) - // void set_jets(const std::vector& jets, const std::string& key="CANONICAL") { + // void set_jets(const std::vector& jets, const std::string& key) { // set_jets(mkconst(jets), key); // } - /// @brief Add a jet to the jets collection + /// @brief Clear a jet collection + /// + /// @note This resets the cluster sequence, but as a shared_ptr is used for storage, + /// any existing shared_ptr links to the previous one will keep it alive. + void clear_jets(const std::string& key) { + for (const Jet* j : jets(key)) delete j; + _jets.erase(key); + _cseqs.erase(key); + } + + + /// @brief Add a jet to a jet collection /// /// @warning The Jet should be new'd; Event will take ownership. + /// /// @todo "Lock" at some point so that jet finding etc. only get done once - void add_jet(const Jet* j, const std::string& key="CANONICAL") { + void add_jet(const Jet* j, const std::string& key) { _jets[key].push_back(j); std::sort(_jets[key].begin(), _jets[key].end(), _cmpPtDescPtr); + // Check that the CSeq on the added jet is consistent with this collection + /// @todo Needs more care that the cseq pointer is live + // if (_cseqs.find(key) != _cseqs.end() && _cseqs.at(key) && + // !_jets.at(key).empty() && !_jets.at(key).front()->clusterseq() && + // _cseqs.at(key) != _jets.at(key).front()->clusterseq()) { + // throw std::runtime_error("Event::add_jet() received a jet whose cluster sequence mismatched the active one for that collection"); + // } } // /// @brief Add a jet to the jets collection (non-const input) - // void add_jet(Jet* j, const std::string& key="CANONICAL") { + // void add_jet(Jet* j, const std::string& key) { // add_jet(mkconst(j), key); // } + /// @brief Access the jets' ClusterSequence object if possible (can be null) /// /// Optional template arg can be used to cast to a specific derived CS type if wanted. template - const CS* clusterseq(const std::string& key="CANONICAL") { - return dynamic_cast(jets(key).front()->clusterseq()); + typename std::shared_ptr clusterseq(const std::string& key) const { + return std::dynamic_pointer_cast(_cseqs.find(key)->second); } - /// @} + // /// @brief Non-const access to the jets' ClusterSequence object if possible (can be null) + // /// + // /// Optional template arg can be used to cast to a specific derived CS type if wanted. + // template + // std::shared_ptr clusterseq(const std::string& key) { + // return std::dynamic_pointer_cast(_cseq->find(key).second); + // } - /// Map between a jet collection and a cluster sequence pointer - std::map> ClusterSeqMap; + /// Set a specific cluster sequence for the named jet collection (which must currently be empty) + /// + /// @warning The CS should be new'd; Event will take ownership via a shared_ptr + template + void set_clusterseq(std::shared_ptr cseq, const std::string& key) { + if (_cseqs.find(key) != _cseqs.end() && !_cseqs.empty()) { + throw std::runtime_error("Event::set_clusterseq() called for a non-empty jet collection"); + } + _cseqs[key] = cseq; + } - /// Initialise a cluster sequence with a unique pointer - void set_clusterseq(std::vector& jetparticles, const FJNS::JetDefinition& jetdef, const std::string& key="CANONICAL") - { - ClusterSeqMap[key] = std::make_unique(jetparticles, jetdef); + + /// Create and run a new cluster sequence for the named jet collection + /// + /// @warning The CS will be new'd, and Event will take ownership via a shared_ptr + /// + /// @note The resulting pseudojets from this still need to be manually set as Jets. + /// + /// @todo How to run a more advanced CS like the active- or Voronoi-area ones? + template + CSeqBasePtr emplace_clusterseq(std::vector& jetparticles, const FJNS::JetDefinition& jetdef, const std::string& key) { + if (_cseqs.find(key) != _cseqs.end() && !_cseqs.empty()) { + throw std::runtime_error("Event::emplace_clusterseq() called for a non-empty jet collection"); + } + _cseqs[key] = std::make_shared(jetparticles, jetdef); + return _cseqs[key]; } + + /// @} + /// @name Missing momentum /// @{ diff --git a/contrib/heputils/include/HEPUtils/Jet.h b/contrib/heputils/include/HEPUtils/Jet.h index 725596fd66..9255f848e8 100644 --- a/contrib/heputils/include/HEPUtils/Jet.h +++ b/contrib/heputils/include/HEPUtils/Jet.h @@ -1,7 +1,7 @@ // -*- C++ -*- // // This file is part of HEPUtils -- https://gitlab.com/hepcedar/heputils/ -// Copyright (C) 2013-2022 Andy Buckley +// Copyright (C) 2013-2023 Andy Buckley // // Embedding of HEPUtils code in other projects is permitted provided this // notice is retained and the HEPUtils namespace and include path are changed. @@ -51,17 +51,17 @@ namespace HEPUtils { /// Constructor for a light jet without explicit constituents Jet(const P4& mom, bool isB=false, bool isC=false) : _p4(mom) - { set_btag(isB); set_ctag(isC); } + { set_btag(int(isB)); set_ctag(int(isC)); } /// "Cartesian" constructor Jet(double px, double py, double pz, double E, bool isB=false, bool isC=false) : _p4(px, py, pz, E) - { set_btag(isB); set_ctag(isC); } + { set_btag(int(isB)); set_ctag(int(isC)); } /// "PseudoJet" constructor Jet(const FJNS::PseudoJet& pj, bool isB=false, bool isC=false) - : _p4(mk_p4(pj)) - { set_btag(isB); set_ctag(isC); } + : _p4(mk_p4(pj)), _pj(pj) + { set_btag(int(isB)); set_ctag(int(isC)); } /// Constructor for a light jet without explicit constituents, with a general tags map @@ -74,7 +74,7 @@ namespace HEPUtils { /// "PseudoJet" constructor, with a general tags map Jet(const FJNS::PseudoJet& pj, const TagCounts& tags) - : _p4(mk_p4(pj)), _tags(tags) { } + : _p4(mk_p4(pj)), _tags(tags), _pj(pj) { } /// @} @@ -149,14 +149,17 @@ namespace HEPUtils { /// Is this particle tagged as a b? bool btag() const { return tagged(5); } /// Set b-tag value - void set_btag(bool isb) { set_ntags(5, isb); } + void set_btag(int ntag=1) { set_ntags(5, ntag); } /// Is this particle tagged as a c? /// - /// @note Can be simultaneously btag()'d -- analyses should probably only use if fallback from b-tag. - bool ctag() const { return tagged(4); } + /// @note Can be simultaneously b-tagged unless the optional bool is set true + bool ctag(bool not_if_b=false) const { + if (not_if_b) return tagged(5) ? false : tagged(4); + return tagged(4); + } /// Set c-tag value - void set_ctag(bool isc) { set_ntags(4, isc); } + void set_ctag(int ntag=1) { set_ntags(4, ntag); } /// @} @@ -177,8 +180,16 @@ namespace HEPUtils { /// /// Optional template arg can be used to cast to a specific derived CS type if wanted. template - const CS* clusterseq() { - return dynamic_cast(_pj.associated_cs()); + typename std::shared_ptr clusterseq() const { + return std::shared_ptr(clusterseq_raw()); + } + + /// @brief Access the ClusterSequence object if possible (can be null) + /// + /// Optional template arg can be used to cast to a specific derived CS type if wanted. + template + const CS* clusterseq_raw() const { + return dynamic_cast(_pj.associated_cs()); } /// @} diff --git a/contrib/heputils/include/HEPUtils/Vectors.h b/contrib/heputils/include/HEPUtils/Vectors.h index 58db41bc02..f632aea81f 100644 --- a/contrib/heputils/include/HEPUtils/Vectors.h +++ b/contrib/heputils/include/HEPUtils/Vectors.h @@ -1,7 +1,7 @@ // -*- C++ -*- // // This file is part of HEPUtils -- https://gitlab.com/hepcedar/heputils/ -// Copyright (C) 2013-2022 Andy Buckley +// Copyright (C) 2013-2023 Andy Buckley // // Embedding of HEPUtils code in other projects is permitted provided this // notice is retained and the HEPUtils namespace and include path are changed. @@ -367,10 +367,14 @@ namespace HEPUtils { /// Get the transverse momentum (same as rho) double pT() const { return rho(); } - /// Get the spatial phi - double phi() const { if (rho2() == 0) return 0; else return atan2(py(),px()); } - /// Get the spatial theta - double theta() const { if (p2() == 0) return 0; else if (pz() == 0) return 0.5 * M_PI; else return atan2(rho(),pz()); } + /// Get the spatial phi (in the range -pi .. pi) + double phi() const { if (rho2() == 0) return 0; else return atan2(py(), px()); } + /// Get the spatial phi (in the range 0 .. 2pi) + double phi_02pi() const { if (rho2() == 0) return 0; else return phi() + M_PI; } + + /// Get the spatial theta (in the range 0 .. pi) + double theta() const { if (p2() == 0) return 0; else + if (pz() == 0) return M_PI/2; else return atan2(rho(), pz()); } //< atan2(+ve, z) is +ve /// Get the spatial-vector pseudorapidity double eta() const { return std::copysign(log((p() + fabs(pz())) / pT()), pz()); } /// Get the spatial-vector absolute pseudorapidity From 2bd90cd65bd386e7855e3774a66ff5e3fb2c5d88 Mon Sep 17 00:00:00 2001 From: ChrisJChang Date: Fri, 1 Dec 2023 14:59:28 +0100 Subject: [PATCH 14/27] One last update to heputils to make sure the cluster sequence map is cloned --- contrib/heputils/include/HEPUtils/Event.h | 1 + 1 file changed, 1 insertion(+) diff --git a/contrib/heputils/include/HEPUtils/Event.h b/contrib/heputils/include/HEPUtils/Event.h index 9b33158f99..acc29dfa23 100644 --- a/contrib/heputils/include/HEPUtils/Event.h +++ b/contrib/heputils/include/HEPUtils/Event.h @@ -126,6 +126,7 @@ namespace HEPUtils { } } e._pmiss = _pmiss; + e._cseqs = _cseqs; } /// @} From 193c4b015f38103f4075918fdf3e2b63d829a831 Mon Sep 17 00:00:00 2001 From: Anders Kvellestad Date: Sun, 3 Dec 2023 23:01:13 +0100 Subject: [PATCH 15/27] Changed the object_in_cone function to not have a default value for the jet collection name. --- ColliderBit/include/gambit/ColliderBit/Utils.hpp | 2 +- .../src/analyses/Analysis_ATLAS_8TeV_1LEPStop_20invfb.cpp | 8 +++++--- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/ColliderBit/include/gambit/ColliderBit/Utils.hpp b/ColliderBit/include/gambit/ColliderBit/Utils.hpp index 0476eab522..7b0341ca7c 100644 --- a/ColliderBit/include/gambit/ColliderBit/Utils.hpp +++ b/ColliderBit/include/gambit/ColliderBit/Utils.hpp @@ -319,7 +319,7 @@ namespace Gambit /// Check if there's a physics object above ptmin in an annulus rmin..rmax around the given four-momentum p4 - inline bool object_in_cone(const HEPUtils::Event& e, const HEPUtils::P4& p4, double ptmin, double rmax, double rmin=0.05, std::string jetcollection = "antikt_R04") { + inline bool object_in_cone(const HEPUtils::Event& e, std::string jetcollection, const HEPUtils::P4& p4, double ptmin, double rmax, double rmin=0.05) { for (const HEPUtils::Particle* p : e.visible_particles()) if (p->pT() > ptmin && HEPUtils::in_range(HEPUtils::deltaR_eta(p4, *p), rmin, rmax)) return true; for (const HEPUtils::Jet* j : e.jets(jetcollection)) diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_1LEPStop_20invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_1LEPStop_20invfb.cpp index 45c960527f..b74a546d35 100644 --- a/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_1LEPStop_20invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_1LEPStop_20invfb.cpp @@ -206,6 +206,8 @@ namespace Gambit { void run(const HEPUtils::Event* event) { + static const std::string jet_collection_name = "antikt_R04"; + // Missing energy HEPUtils::P4 ptot = event->missingmom(); double met = event->met(); @@ -214,7 +216,7 @@ namespace Gambit { vector baselineElectrons; for (const HEPUtils::Particle* electron : event->electrons()) { if (electron->pT() > 10. && electron->abseta() < 2.47 && - !object_in_cone(*event, *electron, 0.1*electron->pT(), 0.2)) baselineElectrons.push_back(electron); + !object_in_cone(*event, jet_collection_name, *electron, 0.1*electron->pT(), 0.2)) baselineElectrons.push_back(electron); } // Apply electron efficiency @@ -224,7 +226,7 @@ namespace Gambit { vector baselineMuons; for (const HEPUtils::Particle* muon : event->muons()) { if (muon->pT() > 10. && muon->abseta() < 2.4 && - !object_in_cone(*event, *muon, 1.8, 0.2)) baselineMuons.push_back(muon); + !object_in_cone(*event, jet_collection_name, *muon, 1.8, 0.2)) baselineMuons.push_back(muon); } // Apply muon efficiency @@ -232,7 +234,7 @@ namespace Gambit { // Get b jets with efficiency and mistag (fake) rates vector baselineJets, bJets; // trueBJets; //for debugging - for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) { + for (const HEPUtils::Jet* jet : event->jets(jet_collection_name)) { if (jet->pT() > 20. && jet->abseta() < 10.0) baselineJets.push_back(jet); if (jet->abseta() < 2.5 && jet->pT() > 25.) { if ((jet->btag() && Random::draw() < 0.75) || (!jet->btag() && Random::draw() < 0.02)) bJets.push_back(jet); From f51aba1f02cdba6aea5858c16946879bdb0d179e Mon Sep 17 00:00:00 2001 From: Anders Kvellestad Date: Sun, 3 Dec 2023 23:29:55 +0100 Subject: [PATCH 16/27] Fixed typo. Changed a str --> std::string from in-file consistency. --- .../include/gambit/ColliderBit/colliders/BaseCollider.hpp | 2 +- .../analyses/Analysis_ATLAS_13TeV_2BoostedBosons_139invfb.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/ColliderBit/include/gambit/ColliderBit/colliders/BaseCollider.hpp b/ColliderBit/include/gambit/ColliderBit/colliders/BaseCollider.hpp index 0359dae43c..ca37e14bea 100644 --- a/ColliderBit/include/gambit/ColliderBit/colliders/BaseCollider.hpp +++ b/ColliderBit/include/gambit/ColliderBit/colliders/BaseCollider.hpp @@ -74,7 +74,7 @@ namespace Gambit std::vector all_jet_collection_settings; /// Key for jet collection used in adding taus - str jetcollection_taus; + std::string jetcollection_taus; }; diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2BoostedBosons_139invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2BoostedBosons_139invfb.cpp index 5e07a98afe..8cea825e2b 100644 --- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2BoostedBosons_139invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2BoostedBosons_139invfb.cpp @@ -169,7 +169,7 @@ namespace Gambit vector fatJets; for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) { - // cout << jet->pT() << " " << jet->mass() << " Z-tag " << jet-tagged(23) << " W-tag " << jet->tagged(24) << " " << endl; + // cout << jet->pT() << " " << jet->mass() << " Z-tag " << jet->tagged(23) << " W-tag " << jet->tagged(24) << " " << endl; if (jet->pT() > 200. && fabs(jet->eta()) < 2.0 && jet->mass() > 40.) { fatJets.push_back(jet); From fe3b847cca2fda9af73ac24999df20b673ac6e65 Mon Sep 17 00:00:00 2001 From: Anders Kvellestad Date: Sun, 3 Dec 2023 23:57:49 +0100 Subject: [PATCH 17/27] Some code formatting corrections --- ColliderBit/src/ColliderBit_measurements.cpp | 8 ++++---- .../Analysis_ATLAS_13TeV_2BoostedBosons_139invfb.cpp | 2 +- ColliderBit/src/detectors/BuckFast.cpp | 4 ++-- 3 files changed, 7 insertions(+), 7 deletions(-) diff --git a/ColliderBit/src/ColliderBit_measurements.cpp b/ColliderBit/src/ColliderBit_measurements.cpp index 9e632c5101..1964fc12f2 100644 --- a/ColliderBit/src/ColliderBit_measurements.cpp +++ b/ColliderBit/src/ColliderBit_measurements.cpp @@ -484,7 +484,7 @@ namespace Gambit summary_line << "LHC Contur LogLikes per pool: "; result = (*Dep::LHC_measurements).pool_LLR; - for( auto const& entry : result) + for (auto const& entry : result) { summary_line << entry.first << ":" << entry.second << ", "; } @@ -506,7 +506,7 @@ namespace Gambit result[pool_LLR_entry.first + "_" + contur_output_instance.first] = pool_LLR_entry.second; } } - for( auto const& entry : result) + for (auto const& entry : result) { summary_line << entry.first << ":" << entry.second << ", "; } @@ -522,7 +522,7 @@ namespace Gambit summary_line << "LHC Contur LogLikes per pool: "; result = (*Dep::LHC_measurements).pool_tags; - for( auto const& entry : result) + for (auto const& entry : result) { summary_line << entry.first << ":" << entry.second << ", "; } @@ -547,7 +547,7 @@ namespace Gambit result[pool_LLR_entry.first + "_" + contur_output_instance.first] = pool_LLR_entry.second; } } - for( auto const& entry : result) + for (auto const& entry : result) { summary_line << entry.first << ":" << entry.second << ", "; } diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2BoostedBosons_139invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2BoostedBosons_139invfb.cpp index 8cea825e2b..9e490760cc 100644 --- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2BoostedBosons_139invfb.cpp +++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2BoostedBosons_139invfb.cpp @@ -371,7 +371,7 @@ namespace Gambit // Compare meff spectrum cout << "Meff SR-4Q-VV\t" << "GAMBIT\t" << "ATLAS " << endl; - for( size_t j = 0; j < _meff_4QVV.size(); j++){ + for (size_t j = 0; j < _meff_4QVV.size(); j++){ cout << "[" << _meff_bins[j] << ", " << _meff_bins[j+1] << "]\t" << _meff_4QVV[j]*_scale << "\t" << _meff_4QVV_model[j] << endl; } diff --git a/ColliderBit/src/detectors/BuckFast.cpp b/ColliderBit/src/detectors/BuckFast.cpp index 75eb6b9787..70a945b2c7 100644 --- a/ColliderBit/src/detectors/BuckFast.cpp +++ b/ColliderBit/src/detectors/BuckFast.cpp @@ -41,14 +41,14 @@ namespace Gambit // Smear jet momenta if (smearJets != NULL) { - for ( std::string jetcollection : event.jet_collections()) + for (std::string jetcollection : event.jet_collections()) { smearJets(event.jets(jetcollection)); } } // Unset b-tags outside |eta|=2.5 - for ( std::string jetcollection : event.jet_collections()) + for (std::string jetcollection : event.jet_collections()) { for (HEPUtils::Jet* j : event.jets(jetcollection)) { From 7864f90e4d6a9dd08e0f09812648fabc40ad209d Mon Sep 17 00:00:00 2001 From: Anders Kvellestad Date: Mon, 4 Dec 2023 01:03:50 +0100 Subject: [PATCH 18/27] Added a helper function in getHepMCEvent.cpp to reduce code duplication. --- ColliderBit/src/getHepMCEvent.cpp | 65 ++++++++++--------------------- 1 file changed, 21 insertions(+), 44 deletions(-) diff --git a/ColliderBit/src/getHepMCEvent.cpp b/ColliderBit/src/getHepMCEvent.cpp index 107183151b..8ff3f808ac 100644 --- a/ColliderBit/src/getHepMCEvent.cpp +++ b/ColliderBit/src/getHepMCEvent.cpp @@ -231,19 +231,14 @@ namespace Gambit } - /// A nested function that reads in HepMC event files and converts them to HEPUtils::Event format - void getHepMCEvent_HEPUtils(HEPUtils::Event &result) - { - using namespace Pipes::getHepMCEvent_HEPUtils; - // Get yaml options - const static str HepMC_filename = runOptions->getValueOrDef("", "hepmc_filename"); - const static double jet_pt_min = runOptions->getValueOrDef(10.0, "jet_pt_min"); - std::vector all_jet_collection_settings = {}; - str jetcollection_taus; - if (runOptions->hasKey("jet_collections")) + /// A helper function for collecting the jet_collections yaml settings. + /// Used by functions getHepMCEvent_HEPUtils and convertHepMCEvent_HEPUtils. + void read_jet_collections_settings(const Options& runOptions, std::vector& all_jet_collection_settings, str& jetcollection_taus) + { + if (runOptions.hasKey("jet_collections")) { - YAML::Node jetcollectionNode = runOptions->getValue("jet_collections"); + YAML::Node jetcollectionNode = runOptions.getValue("jet_collections"); Options jetcollectionOptions(jetcollectionNode); str algorithm; @@ -274,6 +269,20 @@ namespace Gambit all_jet_collection_settings = {{"antikt_R04", "antikt", 0.4, "E_scheme", "Best"}}; jetcollection_taus = "antikt_R04"; } + } + + + /// A nested function that reads in HepMC event files and converts them to HEPUtils::Event format + void getHepMCEvent_HEPUtils(HEPUtils::Event &result) + { + using namespace Pipes::getHepMCEvent_HEPUtils; + + // Get yaml options + const static str HepMC_filename = runOptions->getValueOrDef("", "hepmc_filename"); + const static double jet_pt_min = runOptions->getValueOrDef(10.0, "jet_pt_min"); + std::vector all_jet_collection_settings = {}; + str jetcollection_taus; + read_jet_collections_settings(*runOptions, all_jet_collection_settings, jetcollection_taus); // Get the HepMC event //HepMC3::GenEvent ge = *Dep::HardScatteringEvent; @@ -305,39 +314,7 @@ namespace Gambit const static double jet_pt_min = runOptions->getValueOrDef(10.0, "jet_pt_min"); std::vector all_jet_collection_settings = {}; str jetcollection_taus; - if (runOptions->hasKey("jet_collections")) - { - YAML::Node jetcollectionNode = runOptions->getValue("jet_collections"); - Options jetcollectionOptions(jetcollectionNode); - - str algorithm; - double R; - str recombination_scheme; - str strategy; - std::vector jetcollections = jetcollectionOptions.getNames(); - - for (str key : jetcollections) - { - algorithm = jetcollectionOptions.getValueOrDef("antikt", "algorithm"); - R = jetcollectionOptions.getValueOrDef(0.4, "R"); - recombination_scheme = jetcollectionOptions.getValueOrDef("E_scheme", "recombination_scheme"); - strategy = jetcollectionOptions.getValueOrDef("Best", "strategy"); - - all_jet_collection_settings.push_back({key, algorithm, R, recombination_scheme, strategy}); - } - - jetcollection_taus = jetcollectionOptions.getValueOrDef("antikt_R04", "jetcollection_taus"); - // Throw an error if the jetcollection_taus setting is not given and not using the antikt_R04 collection - if (std::find(jetcollections.begin(), jetcollections.end(), jetcollection_taus) == jetcollections.end()) - { - ColliderBit_error().raise(LOCAL_INFO,"Please provide the jetcollection_taus setting for jet collections if not using antikt_R04."); - } - } - else - { - all_jet_collection_settings = {{"antikt_R04", "antikt", 0.4, "E_scheme", "Best"}}; - jetcollection_taus = "antikt_R04"; - } + read_jet_collections_settings(*runOptions, all_jet_collection_settings, jetcollection_taus); //Set the weight result.set_weight(ge.weight()); From dc31720510805629c732bdb8620581c456bab0b2 Mon Sep 17 00:00:00 2001 From: Anders Kvellestad Date: Mon, 4 Dec 2023 01:22:23 +0100 Subject: [PATCH 19/27] Moved the function get_HEPUtils_event (for LHE event --> HEPUtils event conversion) into the Gambit::ColliderBit namespace. --- .../gambit/ColliderBit/lhef2heputils.hpp | 14 +- ColliderBit/src/lhef2heputils.cpp | 160 +++++++++--------- 2 files changed, 96 insertions(+), 78 deletions(-) diff --git a/ColliderBit/include/gambit/ColliderBit/lhef2heputils.hpp b/ColliderBit/include/gambit/ColliderBit/lhef2heputils.hpp index 95a28a1415..934241ec5b 100644 --- a/ColliderBit/include/gambit/ColliderBit/lhef2heputils.hpp +++ b/ColliderBit/include/gambit/ColliderBit/lhef2heputils.hpp @@ -33,7 +33,17 @@ /// Forward declaration to cut down on includes namespace LHEF { class Reader; } -/// Extract an LHE event as a HEPUtils::Event -void get_HEPUtils_event(const LHEF::Reader&, HEPUtils::Event&, double, std::vector); +namespace Gambit +{ + + namespace ColliderBit + { + + /// Extract an LHE event as a HEPUtils::Event + void get_HEPUtils_event(const LHEF::Reader&, HEPUtils::Event&, double, std::vector); + + } + +} #endif diff --git a/ColliderBit/src/lhef2heputils.cpp b/ColliderBit/src/lhef2heputils.cpp index ffa8347fe0..cbbd44d6cb 100644 --- a/ColliderBit/src/lhef2heputils.cpp +++ b/ColliderBit/src/lhef2heputils.cpp @@ -44,91 +44,99 @@ using namespace std; using namespace HEPUtils; using namespace FJNS; -using namespace Gambit::ColliderBit; - -/// Extract an LHE event as a HEPUtils::Event -void get_HEPUtils_event(const LHEF::Reader& lhe, Event& evt, double jet_pt_min, std::vector all_jet_collection_settings) +namespace Gambit { - P4 vmet; - vector jetparticles; - - evt.set_weight(lhe.hepeup.weight()); - - // Loop over all particles in the event - for (int i = 0; i < lhe.hepeup.NUP; ++i) + namespace ColliderBit { - // Get status and PID code - const int st = lhe.hepeup.ISTUP[i]; - const int pid = lhe.hepeup.IDUP[i]; - const int apid = fabs(pid); - - // Use LHE-stable particles only - if (st != 1) continue; - - // Get 4-momentum - const P4 p4 = P4::mkXYZM(lhe.hepeup.PUP[i][0], lhe.hepeup.PUP[i][1], lhe.hepeup.PUP[i][2], lhe.hepeup.PUP[i][4]); - // Store interacting prompt particles - /// @todo Dress leptons? - if (apid == 22 || apid == 11 || apid == 13 || apid == 15) + /// Extract an LHE event as a HEPUtils::Event + void get_HEPUtils_event(const LHEF::Reader& lhe, Event& evt, double jet_pt_min, std::vector all_jet_collection_settings) { - Particle* p = new Particle(p4, pid); // the event will take ownership of this pointer - p->set_prompt(true); - evt.add_particle(p); - } - // Aggregate missing ET - else if (apid == 12 || apid == 14 || apid == 16 || apid == 1000022 || apid == 1000039) - { - Particle* p = new Particle(p4, pid); // the event will take ownership of this pointer - p->set_prompt(true); - evt.add_particle(p); - vmet += p4; - } + P4 vmet; + vector jetparticles; + + evt.set_weight(lhe.hepeup.weight()); + + // Loop over all particles in the event + for (int i = 0; i < lhe.hepeup.NUP; ++i) + { + // Get status and PID code + const int st = lhe.hepeup.ISTUP[i]; + const int pid = lhe.hepeup.IDUP[i]; + const int apid = fabs(pid); + + // Use LHE-stable particles only + if (st != 1) continue; + + // Get 4-momentum + const P4 p4 = P4::mkXYZM(lhe.hepeup.PUP[i][0], lhe.hepeup.PUP[i][1], lhe.hepeup.PUP[i][2], lhe.hepeup.PUP[i][4]); + + // Store interacting prompt particles + /// @todo Dress leptons? + if (apid == 22 || apid == 11 || apid == 13 || apid == 15) + { + Particle* p = new Particle(p4, pid); // the event will take ownership of this pointer + p->set_prompt(true); + evt.add_particle(p); + } + + // Aggregate missing ET + else if (apid == 12 || apid == 14 || apid == 16 || apid == 1000022 || apid == 1000039) + { + Particle* p = new Particle(p4, pid); // the event will take ownership of this pointer + p->set_prompt(true); + evt.add_particle(p); + vmet += p4; + } + + // Store non-prompt momenta for Jet building + else + { + PseudoJet pj = mk_pj(p4); + pj.set_user_index(apid); + jetparticles.push_back(pj); + } + } + + // MET + vmet.setPz(0); + vmet.setM(0); + evt.set_missingmom(vmet); + + // Jet Finding + for (jet_collection_settings jetcollection : all_jet_collection_settings) + { + + // @todo get_jets function could accept a more general jet definition + vector jets = HEPUtils::get_jets(jetparticles, 0.4, jet_pt_min, FJalgorithm_map(jetcollection.algorithm)); + + for (const PseudoJet& pj : jets) + { + bool hasC = false, hasB = false; + /// @todo Bug in HEPUtils::get_jets means that constituent info is lost for now... + // for (const PseudoJet& c : pj.constituents()) { + // if (c.user_index() == 4) hasC = true; + // if (c.user_index() == 5) hasB = true; + // } + evt.add_jet(new Jet(mk_p4(pj), hasB, hasC), jetcollection.key); + } + + #ifdef COLLIDERBIT_DEBUG + // Print event summary + cout << " MET = " << evt.met() << " GeV" << endl; + cout << " #e = " << evt.electrons().size() << endl; + cout << " #mu = " << evt.muons().size() << endl; + cout << " #tau = " << evt.taus().size() << endl; + cout << " #jet = " << evt.jets().size() << endl; + cout << " #pho = " << evt.photons().size() << endl; + cout << endl; + #endif + } - // Store non-prompt momenta for Jet building - else - { - PseudoJet pj = mk_pj(p4); - pj.set_user_index(apid); - jetparticles.push_back(pj); - } - } - - // MET - vmet.setPz(0); - vmet.setM(0); - evt.set_missingmom(vmet); - - // Jet Finding - for (jet_collection_settings jetcollection : all_jet_collection_settings) - { - - // @todo get_jets function could accept a more general jet definition - vector jets = get_jets(jetparticles, 0.4, jet_pt_min, FJalgorithm_map(jetcollection.algorithm)); - - for (const PseudoJet& pj : jets) - { - bool hasC = false, hasB = false; - /// @todo Bug in HEPUtils::get_jets means that constituent info is lost for now... - // for (const PseudoJet& c : pj.constituents()) { - // if (c.user_index() == 4) hasC = true; - // if (c.user_index() == 5) hasB = true; - // } - evt.add_jet(new Jet(mk_p4(pj), hasB, hasC), jetcollection.key); } - #ifdef COLLIDERBIT_DEBUG - // Print event summary - cout << " MET = " << evt.met() << " GeV" << endl; - cout << " #e = " << evt.electrons().size() << endl; - cout << " #mu = " << evt.muons().size() << endl; - cout << " #tau = " << evt.taus().size() << endl; - cout << " #jet = " << evt.jets().size() << endl; - cout << " #pho = " << evt.photons().size() << endl; - cout << endl; - #endif } } From f64bc2a0b6cdf0d0c9583fb50d081e41b095ce1a Mon Sep 17 00:00:00 2001 From: Anders Kvellestad Date: Mon, 4 Dec 2023 01:30:50 +0100 Subject: [PATCH 20/27] Fixed typo --- gum/src/particledb.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gum/src/particledb.py b/gum/src/particledb.py index 3f75ba95ac..ade2ef16ae 100644 --- a/gum/src/particledb.py +++ b/gum/src/particledb.py @@ -394,7 +394,7 @@ def get_higgses(partlist): def get_invisibles(invisible_pdgs): """ - Get the PDG codes of invisibles particles to write to 'contrib/heputils/include/HEPUtils/Particleh' + Get the PDG codes of invisibles particles to write to 'contrib/heputils/include/HEPUtils/Particle.h' Will not write for photons, leptons or existing invisibles """ From 5043efc30c0c02883c69cda46c7968b164811dcf Mon Sep 17 00:00:00 2001 From: Anders Kvellestad Date: Mon, 4 Dec 2023 02:06:20 +0100 Subject: [PATCH 21/27] Fixed some old todo comments. --- .../ColliderBit/colliders/Pythia8/Py8EventConversions.hpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/ColliderBit/include/gambit/ColliderBit/colliders/Pythia8/Py8EventConversions.hpp b/ColliderBit/include/gambit/ColliderBit/colliders/Pythia8/Py8EventConversions.hpp index f3079fd87f..3a07adfd74 100644 --- a/ColliderBit/include/gambit/ColliderBit/colliders/Pythia8/Py8EventConversions.hpp +++ b/ColliderBit/include/gambit/ColliderBit/colliders/Pythia8/Py8EventConversions.hpp @@ -202,8 +202,11 @@ namespace Gambit FJNS::RecombinationScheme jet_recomscheme = FJRecomScheme_map(jetcollection.recombination_scheme); const FJNS::JetDefinition jet_def(jet_algorithm, jetcollection.R, jet_strategy, jet_recomscheme); - /// @todo For substructure we need to keep this ClusterSequence alive... make_unique() ctor and attach to the Event? Or manage at CB level? + /// Create and run a new cluster sequence for the given jet collection. + /// The HEPUtils::Event instance ('result') takes ownership of the + /// cluster sequence and a shared_ptr is returned here. std::shared_ptr CSeqBasePtr = result.emplace_clusterseq(jetparticles, jet_def, jetcollection.key); + /// Get the resulting pseudojets std::vector pjets = sorted_by_pt(CSeqBasePtr->inclusive_jets(jet_pt_min)); /// Do jet b-tagging, etc. and add to the Event @@ -286,7 +289,6 @@ namespace Gambit } // Add jet to collection including tags and PseudoJet - /// @todo We need to do something smart if we want to keep the ClusterSequence alive HEPUtils::Jet::TagCounts tags{ {5,int(isB)}, {4,int(isC)}, {23,int(isZ)}, {24,int(isW)}, {25,int(ish)} }; result.add_jet(new HEPUtils::Jet(pj, tags), jetcollection.key); } From 2f96ed2584f2f5b90e6519be5f4daf1573c566b1 Mon Sep 17 00:00:00 2001 From: Anders Kvellestad Date: Mon, 4 Dec 2023 03:20:32 +0100 Subject: [PATCH 22/27] Fixed a small bug in the parsing of jet_collections yaml options. Now the yaml options are actually used. :) --- .../gambit/ColliderBit/getPy8Collider.hpp | 25 +++++++++-------- ColliderBit/src/getHepMCEvent.cpp | 27 ++++++++++--------- ColliderBit/src/getLHEvent.cpp | 25 +++++++++-------- yaml_files/ColliderBit_CMSSM.yaml | 21 +++++++++++++-- 4 files changed, 62 insertions(+), 36 deletions(-) diff --git a/ColliderBit/include/gambit/ColliderBit/getPy8Collider.hpp b/ColliderBit/include/gambit/ColliderBit/getPy8Collider.hpp index e4841ea9f3..b59a8c35a3 100644 --- a/ColliderBit/include/gambit/ColliderBit/getPy8Collider.hpp +++ b/ColliderBit/include/gambit/ColliderBit/getPy8Collider.hpp @@ -112,31 +112,34 @@ namespace Gambit // Fill the jet collection settings if (colOptions.hasKey("jet_collections")) { - YAML::Node jetcollectionNode = colOptions.getValue("jet_collections"); - Options jetcollectionOptions(jetcollectionNode); + YAML::Node all_jetcollections_node = colOptions.getValue("jet_collections"); + Options all_jetcollection_options(all_jetcollections_node); str algorithm; double R; str recombination_scheme; str strategy; result.all_jet_collection_settings = {}; // Initialise with no collections, if no jet collection setting specified, then should use the basecollider's initialised jet collection - std::vector jetcollections = jetcollectionOptions.getNames(); + std::vector jetcollection_names = all_jetcollection_options.getNames(); - for (str key : jetcollections) + for (str key : jetcollection_names) { - algorithm = jetcollectionOptions.getValueOrDef(algorithm_default, "algorithm"); - R = jetcollectionOptions.getValueOrDef(R_default, "R"); - recombination_scheme = jetcollectionOptions.getValueOrDef(recombination_scheme_default, "recombination_scheme"); - strategy = jetcollectionOptions.getValueOrDef(strategy_default, "strategy"); + YAML::Node current_jc_node = all_jetcollection_options.getValue(key); + Options current_jc_options(current_jc_node); + + algorithm = current_jc_options.getValueOrDef(algorithm_default, "algorithm"); + R = current_jc_options.getValueOrDef(R_default, "R"); + recombination_scheme = current_jc_options.getValueOrDef(recombination_scheme_default, "recombination_scheme"); + strategy = current_jc_options.getValueOrDef(strategy_default, "strategy"); (result.all_jet_collection_settings).push_back({key, algorithm, R, recombination_scheme, strategy}); } - result.jetcollection_taus = colOptions.getValueOrDef("antikt_R04", "jetcollection_taus"); + result.jetcollection_taus = colOptions.getValueOrDef("antikt_R04", "jet_collection_taus"); // Throw an error if the jetcollection_taus setting is not given and not using the antikt_R04 collection - if (std::find(jetcollections.begin(), jetcollections.end(), result.jetcollection_taus) == jetcollections.end()) + if (std::find(jetcollection_names.begin(), jetcollection_names.end(), result.jetcollection_taus) == jetcollection_names.end()) { - ColliderBit_error().raise(LOCAL_INFO,"Please provide the jetcollection_taus setting for jet collections if not using antikt_R04."); + ColliderBit_error().raise(LOCAL_INFO,"Please provide the jet_collection_taus setting for jet collections if not using antikt_R04."); } } if (colOptions.hasKey("pythia_settings")) diff --git a/ColliderBit/src/getHepMCEvent.cpp b/ColliderBit/src/getHepMCEvent.cpp index 8ff3f808ac..e3f0580244 100644 --- a/ColliderBit/src/getHepMCEvent.cpp +++ b/ColliderBit/src/getHepMCEvent.cpp @@ -238,30 +238,33 @@ namespace Gambit { if (runOptions.hasKey("jet_collections")) { - YAML::Node jetcollectionNode = runOptions.getValue("jet_collections"); - Options jetcollectionOptions(jetcollectionNode); - + YAML::Node all_jetcollections_node = runOptions.getValue("jet_collections"); + Options all_jetcollection_options(all_jetcollections_node); + str algorithm; double R; str recombination_scheme; str strategy; - std::vector jetcollections = jetcollectionOptions.getNames(); + std::vector jetcollection_names = all_jetcollection_options.getNames(); - for (str key : jetcollections) + for (str key : jetcollection_names) { - algorithm = jetcollectionOptions.getValueOrDef("antikt", "algorithm"); - R = jetcollectionOptions.getValueOrDef(0.4, "R"); - recombination_scheme = jetcollectionOptions.getValueOrDef("E_scheme", "recombination_scheme"); - strategy = jetcollectionOptions.getValueOrDef("Best", "strategy"); + YAML::Node current_jc_node = all_jetcollection_options.getValue(key); + Options current_jc_options(current_jc_node); + + algorithm = current_jc_options.getValueOrDef("antikt", "algorithm"); + R = current_jc_options.getValueOrDef(0.4, "R"); + recombination_scheme = current_jc_options.getValueOrDef("E_scheme", "recombination_scheme"); + strategy = current_jc_options.getValueOrDef("Best", "strategy"); all_jet_collection_settings.push_back({key, algorithm, R, recombination_scheme, strategy}); } - jetcollection_taus = jetcollectionOptions.getValueOrDef("antikt_R04", "jetcollection_taus"); + jetcollection_taus = runOptions.getValueOrDef("antikt_R04", "jet_collection_taus"); // Throw an error if the jetcollection_taus setting is not given and not using the antikt_R04 collection - if (std::find(jetcollections.begin(), jetcollections.end(), jetcollection_taus) == jetcollections.end()) + if (std::find(jetcollection_names.begin(), jetcollection_names.end(), jetcollection_taus) == jetcollection_names.end()) { - ColliderBit_error().raise(LOCAL_INFO,"Please provide the jetcollection_taus setting for jet collections if not using antikt_R04."); + ColliderBit_error().raise(LOCAL_INFO,"Please provide the jet_collection_taus setting for jet collections if not using antikt_R04."); } } else diff --git a/ColliderBit/src/getLHEvent.cpp b/ColliderBit/src/getLHEvent.cpp index 535649e97d..0b84d6c93d 100644 --- a/ColliderBit/src/getLHEvent.cpp +++ b/ColliderBit/src/getLHEvent.cpp @@ -66,30 +66,33 @@ namespace Gambit // Fill the jet collection settings if (colOptions.hasKey("jet_collections")) { - YAML::Node jetcollectionNode = colOptions.getValue("jet_collections"); - Options jetcollectionOptions(jetcollectionNode); + YAML::Node all_jetcollections_node = colOptions.getValue("jet_collections"); + Options all_jetcollection_options(all_jetcollections_node); str algorithm; double R; str recombination_scheme; str strategy; - std::vector jetcollections = jetcollectionOptions.getNames(); + std::vector jetcollection_names = all_jetcollection_options.getNames(); - for (str key : jetcollections) + for (str key : jetcollection_names) { - algorithm = jetcollectionOptions.getValueOrDef("antikt", "algorithm"); - R = jetcollectionOptions.getValueOrDef(0.4, "R"); - recombination_scheme = jetcollectionOptions.getValueOrDef("E_scheme", "recombination_scheme"); - strategy = jetcollectionOptions.getValueOrDef("Best", "strategy"); + YAML::Node current_jc_node = all_jetcollection_options.getValue(key); + Options current_jc_options(current_jc_node); + + algorithm = current_jc_options.getValueOrDef("antikt", "algorithm"); + R = current_jc_options.getValueOrDef(0.4, "R"); + recombination_scheme = current_jc_options.getValueOrDef("E_scheme", "recombination_scheme"); + strategy = current_jc_options.getValueOrDef("Best", "strategy"); all_jet_collection_settings.push_back({key, algorithm, R, recombination_scheme, strategy}); } - jetcollection_taus = colOptions.getValueOrDef("antikt_R04", "jetcollection_taus"); + jetcollection_taus = colOptions.getValueOrDef("antikt_R04", "jet_collection_taus"); // Throw an error if the jetcollection_taus setting is not given and not using the antikt_R04 collection - if (std::find(jetcollections.begin(), jetcollections.end(), jetcollection_taus) == jetcollections.end()) + if (std::find(jetcollection_names.begin(), jetcollection_names.end(), jetcollection_taus) == jetcollection_names.end()) { - ColliderBit_error().raise(LOCAL_INFO,"Please provide the jetcollection_taus setting for jet collections if not using antikt_R04."); + ColliderBit_error().raise(LOCAL_INFO,"Please provide the jet_collection_taus setting for jet collections if not using antikt_R04."); } } } diff --git a/yaml_files/ColliderBit_CMSSM.yaml b/yaml_files/ColliderBit_CMSSM.yaml index 313f2a1fc1..d13f943650 100644 --- a/yaml_files/ColliderBit_CMSSM.yaml +++ b/yaml_files/ColliderBit_CMSSM.yaml @@ -310,7 +310,13 @@ Rules: LHC_8TeV: xsec_veto: 0.00 partonOnly: false - antiktR: 0.4 + jet_collections: + antikt_R04: + algorithm: antikt + R: 0.4 + recombination_scheme: E_scheme + strategy: Best + jet_collection_taus: antikt_R04 pythia_settings: - Beams:eCM = 8000 - Next:numberShowProcess = 1 @@ -329,7 +335,18 @@ Rules: # 0.028 fb corresponds to ~1 expected event at L = 36 fb^-1. xsec_veto: 0.028 partonOnly: false - antiktR: 0.4 + jet_collections: + antikt_R04: + algorithm: antikt + R: 0.4 + recombination_scheme: E_scheme + strategy: Best + # antikt_R08: + # algorithm: antikt + # R: 0.8 + # recombination_scheme: E_scheme + # strategy: Best + jet_collection_taus: antikt_R04 pythia_settings: - Beams:eCM = 13000 - Next:numberShowProcess = 1 From 0762388f4005e259dedf9da08965897e1c97d6fd Mon Sep 17 00:00:00 2001 From: Anders Kvellestad Date: Mon, 4 Dec 2023 11:34:24 +0100 Subject: [PATCH 23/27] Updated the jet settings in YAML files to match the new jet_collections system. --- yaml_files/CMSSM.yaml | 8 +++++++- yaml_files/ColliderBit_CMSSM.yaml | 6 +++--- yaml_files/MSSM7.yaml | 8 +++++++- yaml_files/MSSM9.yaml | 8 +++++++- yaml_files/NUHM1.yaml | 8 +++++++- yaml_files/NUHM2.yaml | 8 +++++++- 6 files changed, 38 insertions(+), 8 deletions(-) diff --git a/yaml_files/CMSSM.yaml b/yaml_files/CMSSM.yaml index 46f2fc16f6..9f34496e15 100644 --- a/yaml_files/CMSSM.yaml +++ b/yaml_files/CMSSM.yaml @@ -708,7 +708,13 @@ Rules: LHC_13TeV: # 0.028 fb corresponds to ~1 expected event at L = 36 fb^-1. xsec_veto: 0.028 - antiktR: 0.4 + jet_collections: + antikt_R04: + algorithm: antikt + R: 0.4 + recombination_scheme: E_scheme + strategy: Best + jet_collection_taus: antikt_R04 partonOnly: false pythia_settings: - Print:quiet = on diff --git a/yaml_files/ColliderBit_CMSSM.yaml b/yaml_files/ColliderBit_CMSSM.yaml index d13f943650..ef6835e9c0 100644 --- a/yaml_files/ColliderBit_CMSSM.yaml +++ b/yaml_files/ColliderBit_CMSSM.yaml @@ -337,10 +337,10 @@ Rules: partonOnly: false jet_collections: antikt_R04: - algorithm: antikt + algorithm: antikt # antikt, cambridge, kt, genkt, cambridge_for_passive R: 0.4 - recombination_scheme: E_scheme - strategy: Best + recombination_scheme: E_scheme # E_scheme, pt_scheme, pt2_scheme + strategy: Best # Best, NlnN # antikt_R08: # algorithm: antikt # R: 0.8 diff --git a/yaml_files/MSSM7.yaml b/yaml_files/MSSM7.yaml index 2e629db6b1..a4b32aa393 100644 --- a/yaml_files/MSSM7.yaml +++ b/yaml_files/MSSM7.yaml @@ -616,7 +616,13 @@ Rules: # 0.028 fb corresponds to ~1 expected event at L = 36 fb^-1. xsec_veto: 0.028 partonOnly: false - antiktR: 0.4 + jet_collections: + antikt_R04: + algorithm: antikt + R: 0.4 + recombination_scheme: E_scheme + strategy: Best + jet_collection_taus: antikt_R04 pythia_settings: - Print:quiet = on - Next:numberCount = 0 diff --git a/yaml_files/MSSM9.yaml b/yaml_files/MSSM9.yaml index fe6de4ccfd..d18401677d 100644 --- a/yaml_files/MSSM9.yaml +++ b/yaml_files/MSSM9.yaml @@ -624,7 +624,13 @@ Rules: # 0.028 fb corresponds to ~1 expected event at L = 36 fb^-1. xsec_veto: 0.028 partonOnly: false - antiktR: 0.4 + jet_collections: + antikt_R04: + algorithm: antikt + R: 0.4 + recombination_scheme: E_scheme + strategy: Best + jet_collection_taus: antikt_R04 pythia_settings: - Print:quiet = on - Next:numberCount = 0 diff --git a/yaml_files/NUHM1.yaml b/yaml_files/NUHM1.yaml index 8a160a4ddb..1693527f74 100644 --- a/yaml_files/NUHM1.yaml +++ b/yaml_files/NUHM1.yaml @@ -608,7 +608,13 @@ Rules: # 0.028 fb corresponds to ~1 expected event at L = 36 fb^-1. xsec_veto: 0.028 partonOnly: false - antiktR: 0.4 + jet_collections: + antikt_R04: + algorithm: antikt + R: 0.4 + recombination_scheme: E_scheme + strategy: Best + jet_collection_taus: antikt_R04 pythia_settings: - Print:quiet = on - Next:numberCount = 0 diff --git a/yaml_files/NUHM2.yaml b/yaml_files/NUHM2.yaml index f5a53e5440..2d3af567c6 100644 --- a/yaml_files/NUHM2.yaml +++ b/yaml_files/NUHM2.yaml @@ -611,7 +611,13 @@ Rules: # 0.028 fb corresponds to ~1 expected event at L = 36 fb^-1. xsec_veto: 0.028 partonOnly: false - antiktR: 0.4 + jet_collections: + antikt_R04: + algorithm: antikt + R: 0.4 + recombination_scheme: E_scheme + strategy: Best + jet_collection_taus: antikt_R04 pythia_settings: - Print:quiet = on - Next:numberCount = 0 From d593582edb8040b5588b791ce118f0908fe6ac62 Mon Sep 17 00:00:00 2001 From: ChrisJChang Date: Mon, 4 Dec 2023 14:42:42 +0100 Subject: [PATCH 24/27] Remove default jet collections settings. Throw error if none present. Move location of stored fastjet strategies, algorithms, etc to Utils.cpp --- .../ColliderBit/colliders/BaseCollider.hpp | 2 +- .../gambit/ColliderBit/getPy8Collider.hpp | 31 ++++------ ColliderBit/src/Utils.cpp | 42 +++++++++++++ ColliderBit/src/getHepMCEvent.cpp | 59 ++----------------- ColliderBit/src/getLHEvent.cpp | 29 ++++----- 5 files changed, 75 insertions(+), 88 deletions(-) diff --git a/ColliderBit/include/gambit/ColliderBit/colliders/BaseCollider.hpp b/ColliderBit/include/gambit/ColliderBit/colliders/BaseCollider.hpp index ca37e14bea..d29690567a 100644 --- a/ColliderBit/include/gambit/ColliderBit/colliders/BaseCollider.hpp +++ b/ColliderBit/include/gambit/ColliderBit/colliders/BaseCollider.hpp @@ -35,7 +35,7 @@ namespace Gambit public: /// Constructor - BaseCollider() : partonOnly(false), all_jet_collection_settings({{"antikt_R04", "antikt", 0.4, "E_scheme", "Best"}}), jetcollection_taus("antikt_R04") {} + BaseCollider() : partonOnly(false), all_jet_collection_settings({}), jetcollection_taus("") {} /// Destructor virtual ~BaseCollider() {} /// Reset this instance for reuse, avoiding the need for "new" or "delete". diff --git a/ColliderBit/include/gambit/ColliderBit/getPy8Collider.hpp b/ColliderBit/include/gambit/ColliderBit/getPy8Collider.hpp index b59a8c35a3..978c414f59 100644 --- a/ColliderBit/include/gambit/ColliderBit/getPy8Collider.hpp +++ b/ColliderBit/include/gambit/ColliderBit/getPy8Collider.hpp @@ -98,10 +98,6 @@ namespace Gambit // Get options from yaml file. const double xsec_veto_default = 0.0; const bool partonOnly_default = false; - const str algorithm_default = "antikt"; - const double R_default = 0.4; - const str recombination_scheme_default = "E_scheme"; - const str strategy_default = "Best"; if (runOptions.hasKey(RunMC.current_collider())) { YAML::Node colNode = runOptions.getValue(RunMC.current_collider()); @@ -114,12 +110,6 @@ namespace Gambit { YAML::Node all_jetcollections_node = colOptions.getValue("jet_collections"); Options all_jetcollection_options(all_jetcollections_node); - - str algorithm; - double R; - str recombination_scheme; - str strategy; - result.all_jet_collection_settings = {}; // Initialise with no collections, if no jet collection setting specified, then should use the basecollider's initialised jet collection std::vector jetcollection_names = all_jetcollection_options.getNames(); for (str key : jetcollection_names) @@ -127,21 +117,25 @@ namespace Gambit YAML::Node current_jc_node = all_jetcollection_options.getValue(key); Options current_jc_options(current_jc_node); - algorithm = current_jc_options.getValueOrDef(algorithm_default, "algorithm"); - R = current_jc_options.getValueOrDef(R_default, "R"); - recombination_scheme = current_jc_options.getValueOrDef(recombination_scheme_default, "recombination_scheme"); - strategy = current_jc_options.getValueOrDef(strategy_default, "strategy"); + str algorithm = current_jc_options.getValue("algorithm"); + double R = current_jc_options.getValue("R"); + str recombination_scheme = current_jc_options.getValue("recombination_scheme"); + str strategy = current_jc_options.getValue("strategy"); (result.all_jet_collection_settings).push_back({key, algorithm, R, recombination_scheme, strategy}); } - result.jetcollection_taus = colOptions.getValueOrDef("antikt_R04", "jet_collection_taus"); + result.jetcollection_taus = colOptions.getValue("jet_collection_taus"); // Throw an error if the jetcollection_taus setting is not given and not using the antikt_R04 collection if (std::find(jetcollection_names.begin(), jetcollection_names.end(), result.jetcollection_taus) == jetcollection_names.end()) { - ColliderBit_error().raise(LOCAL_INFO,"Please provide the jet_collection_taus setting for jet collections if not using antikt_R04."); + ColliderBit_error().raise(LOCAL_INFO,"Please provide the jet_collection_taus setting for jet collections."); } } + else + { + ColliderBit_error().raise(LOCAL_INFO,"Could not find jet_collections option for collider " + RunMC.current_collider() + ". Please provide this in the YAML file."); + } if (colOptions.hasKey("pythia_settings")) { std::vector addPythiaOptions = colNode["pythia_settings"].as >(); @@ -150,10 +144,7 @@ namespace Gambit } else { - xsec_veto_fb = xsec_veto_default; - result.partonOnly = partonOnly_default; - result.all_jet_collection_settings = {{"antikt_R04", "antikt", 0.4, "E_scheme", "Best"}}; - result.jetcollection_taus = "antikt_R04"; + ColliderBit_error().raise(LOCAL_INFO,"Could not find runOptions for collider " + RunMC.current_collider() + "."); } // We need showProcesses for the xsec veto. diff --git a/ColliderBit/src/Utils.cpp b/ColliderBit/src/Utils.cpp index a4e93a2946..c7fb69d174 100644 --- a/ColliderBit/src/Utils.cpp +++ b/ColliderBit/src/Utils.cpp @@ -31,6 +31,7 @@ /// ********************************************* #include "gambit/ColliderBit/Utils.hpp" +#include "gambit/ColliderBit/ColliderBit_eventloop.hpp" #include "gambit/Utils/threadsafe_rng.hpp" #include using namespace std; @@ -41,6 +42,47 @@ namespace Gambit { + /// Storage of different FastJet methods + FJNS::JetAlgorithm FJalgorithm_map(str algorithm) + { + FJNS::JetAlgorithm result; + if (algorithm == "antikt") {result = FJNS::antikt_algorithm;} + else if (algorithm == "cambridge") {result = FJNS::cambridge_algorithm;} + else if (algorithm == "kt") {result = FJNS::kt_algorithm;} + else if (algorithm == "genkt") {result = FJNS::genkt_algorithm;} + else if (algorithm == "cambridge_for_passive") {result = FJNS::cambridge_for_passive_algorithm;} + else + { + ColliderBit_error().raise(LOCAL_INFO, "Could not find jet algorithm in list available. Please add the missing option to the FJalgorithm_map function in Py8EventConversions.hpp."); + } + return result; + } + + FJNS::Strategy FJstrategy_map(str strategy) + { + FJNS::Strategy result; + if (strategy == "Best") {result = FJNS::Best;} + else if (strategy == "NlnN") {result = FJNS::NlnN;} + else + { + ColliderBit_error().raise(LOCAL_INFO, "Could not find jet strategy in list available. Please add the missing option to the FJstrategy_map function in Py8EventConversions.hpp."); + } + return result; + } + + FJNS::RecombinationScheme FJRecomScheme_map(str reco_scheme) + { + FJNS::RecombinationScheme result; + if (reco_scheme == "E_scheme") {result = FJNS::E_scheme;} + else if (reco_scheme == "pt_scheme") {result = FJNS::pt_scheme;} + else if (reco_scheme == "pt2_scheme") {result = FJNS::pt2_scheme;} + else + { + ColliderBit_error().raise(LOCAL_INFO, "Could not find jet recombination scheme in list available. Please add the missing option to the FJRecomScheme_map function in Py8EventConversions.hpp."); + } + return result; + } + bool random_bool(double eff) { /// @todo Handle out-of-range eff values diff --git a/ColliderBit/src/getHepMCEvent.cpp b/ColliderBit/src/getHepMCEvent.cpp index e3f0580244..21d0af1790 100644 --- a/ColliderBit/src/getHepMCEvent.cpp +++ b/ColliderBit/src/getHepMCEvent.cpp @@ -57,47 +57,6 @@ namespace Gambit namespace ColliderBit { - /// Storage of different FastJet methods - FJNS::JetAlgorithm FJalgorithm_map(str algorithm) - { - FJNS::JetAlgorithm result; - if (algorithm == "antikt") {result = FJNS::antikt_algorithm;} - else if (algorithm == "cambridge") {result = FJNS::cambridge_algorithm;} - else if (algorithm == "kt") {result = FJNS::kt_algorithm;} - else if (algorithm == "genkt") {result = FJNS::genkt_algorithm;} - else if (algorithm == "cambridge_for_passive") {result = FJNS::cambridge_for_passive_algorithm;} - else - { - ColliderBit_error().raise(LOCAL_INFO, "Could not find jet algorithm in list available. Please add the missing option to the FJalgorithm_map function in Py8EventConversions.hpp."); - } - return result; - } - - FJNS::Strategy FJstrategy_map(str strategy) - { - FJNS::Strategy result; - if (strategy == "Best") {result = FJNS::Best;} - else if (strategy == "NlnN") {result = FJNS::NlnN;} - else - { - ColliderBit_error().raise(LOCAL_INFO, "Could not find jet strategy in list available. Please add the missing option to the FJstrategy_map function in Py8EventConversions.hpp."); - } - return result; - } - - FJNS::RecombinationScheme FJRecomScheme_map(str reco_scheme) - { - FJNS::RecombinationScheme result; - if (reco_scheme == "E_scheme") {result = FJNS::E_scheme;} - else if (reco_scheme == "pt_scheme") {result = FJNS::pt_scheme;} - else if (reco_scheme == "pt2_scheme") {result = FJNS::pt2_scheme;} - else - { - ColliderBit_error().raise(LOCAL_INFO, "Could not find jet recombination scheme in list available. Please add the missing option to the FJRecomScheme_map function in Py8EventConversions.hpp."); - } - return result; - } - /// A nested function that reads in HepMC event files void readHepMCEvent(HepMC3::GenEvent& result, const str HepMC_filename, const MCLoopInfo& RunMC, const int iteration, @@ -240,11 +199,6 @@ namespace Gambit { YAML::Node all_jetcollections_node = runOptions.getValue("jet_collections"); Options all_jetcollection_options(all_jetcollections_node); - - str algorithm; - double R; - str recombination_scheme; - str strategy; std::vector jetcollection_names = all_jetcollection_options.getNames(); for (str key : jetcollection_names) @@ -252,10 +206,10 @@ namespace Gambit YAML::Node current_jc_node = all_jetcollection_options.getValue(key); Options current_jc_options(current_jc_node); - algorithm = current_jc_options.getValueOrDef("antikt", "algorithm"); - R = current_jc_options.getValueOrDef(0.4, "R"); - recombination_scheme = current_jc_options.getValueOrDef("E_scheme", "recombination_scheme"); - strategy = current_jc_options.getValueOrDef("Best", "strategy"); + str algorithm = current_jc_options.getValue("algorithm"); + double R = current_jc_options.getValue("R"); + str recombination_scheme = current_jc_options.getValue("recombination_scheme"); + str strategy = current_jc_options.getValue("strategy"); all_jet_collection_settings.push_back({key, algorithm, R, recombination_scheme, strategy}); } @@ -264,13 +218,12 @@ namespace Gambit // Throw an error if the jetcollection_taus setting is not given and not using the antikt_R04 collection if (std::find(jetcollection_names.begin(), jetcollection_names.end(), jetcollection_taus) == jetcollection_names.end()) { - ColliderBit_error().raise(LOCAL_INFO,"Please provide the jet_collection_taus setting for jet collections if not using antikt_R04."); + ColliderBit_error().raise(LOCAL_INFO,"Please provide the jet_collection_taus setting for jet collections."); } } else { - all_jet_collection_settings = {{"antikt_R04", "antikt", 0.4, "E_scheme", "Best"}}; - jetcollection_taus = "antikt_R04"; + ColliderBit_error().raise(LOCAL_INFO,"Could not find jet_collections options. Please provide this in the YAML file."); } } diff --git a/ColliderBit/src/getLHEvent.cpp b/ColliderBit/src/getLHEvent.cpp index 0b84d6c93d..7c8f778129 100644 --- a/ColliderBit/src/getLHEvent.cpp +++ b/ColliderBit/src/getLHEvent.cpp @@ -57,7 +57,7 @@ namespace Gambit // Get all jet collection settings str jetcollection_taus; - std::vector all_jet_collection_settings = {}; // Initialise with no collections, if no jet collection setting specified, then should use the basecollider's initialised jet collection + std::vector all_jet_collection_settings = {}; if (runOptions->hasKey((*Dep::RunMC).current_collider())) { YAML::Node colNode = runOptions->getValue((*Dep::RunMC).current_collider()); @@ -68,11 +68,6 @@ namespace Gambit { YAML::Node all_jetcollections_node = colOptions.getValue("jet_collections"); Options all_jetcollection_options(all_jetcollections_node); - - str algorithm; - double R; - str recombination_scheme; - str strategy; std::vector jetcollection_names = all_jetcollection_options.getNames(); for (str key : jetcollection_names) @@ -80,26 +75,32 @@ namespace Gambit YAML::Node current_jc_node = all_jetcollection_options.getValue(key); Options current_jc_options(current_jc_node); - algorithm = current_jc_options.getValueOrDef("antikt", "algorithm"); - R = current_jc_options.getValueOrDef(0.4, "R"); - recombination_scheme = current_jc_options.getValueOrDef("E_scheme", "recombination_scheme"); - strategy = current_jc_options.getValueOrDef("Best", "strategy"); + str algorithm = current_jc_options.getValue("algorithm"); + double R = current_jc_options.getValue("R"); + str recombination_scheme = current_jc_options.getValue("recombination_scheme"); + str strategy = current_jc_options.getValue("strategy"); all_jet_collection_settings.push_back({key, algorithm, R, recombination_scheme, strategy}); } - jetcollection_taus = colOptions.getValueOrDef("antikt_R04", "jet_collection_taus"); + // Note that jetcollection_taus is not used here as get_HEPUtils_event(...) has a much simpler jet definition than in Py8Conversions.hpp + jetcollection_taus = colOptions.getValue("jet_collection_taus"); // Throw an error if the jetcollection_taus setting is not given and not using the antikt_R04 collection if (std::find(jetcollection_names.begin(), jetcollection_names.end(), jetcollection_taus) == jetcollection_names.end()) { - ColliderBit_error().raise(LOCAL_INFO,"Please provide the jet_collection_taus setting for jet collections if not using antikt_R04."); + ColliderBit_error().raise(LOCAL_INFO,"Please provide the jet_collection_taus setting for jet collections."); } } + else + { + str current_collider = (*Dep::RunMC).current_collider(); + ColliderBit_error().raise(LOCAL_INFO,"Could not find jet_collections option for collider " + current_collider + ". Please provide this in the YAML file."); + } } else { - all_jet_collection_settings = {{"antikt_R04", "antikt", 0.4, "E_scheme", "Best"}}; - jetcollection_taus = "antikt_R04"; + str current_collider = (*Dep::RunMC).current_collider(); + ColliderBit_error().raise(LOCAL_INFO,"Could not find runOptions for collider " + current_collider + "."); } // Don't do anything during special iterations From 2fb9f93b8493df60baa5400c48ccc030005f1733 Mon Sep 17 00:00:00 2001 From: ChrisJChang Date: Mon, 4 Dec 2023 14:45:45 +0100 Subject: [PATCH 25/27] update error message to point to the right location --- ColliderBit/src/Utils.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/ColliderBit/src/Utils.cpp b/ColliderBit/src/Utils.cpp index c7fb69d174..0a51254fa6 100644 --- a/ColliderBit/src/Utils.cpp +++ b/ColliderBit/src/Utils.cpp @@ -53,7 +53,7 @@ namespace Gambit else if (algorithm == "cambridge_for_passive") {result = FJNS::cambridge_for_passive_algorithm;} else { - ColliderBit_error().raise(LOCAL_INFO, "Could not find jet algorithm in list available. Please add the missing option to the FJalgorithm_map function in Py8EventConversions.hpp."); + ColliderBit_error().raise(LOCAL_INFO, "Could not find jet algorithm in list available. Please add the missing option to the FJalgorithm_map function in ColliderBit/src/Utils.hpp."); } return result; } @@ -65,7 +65,7 @@ namespace Gambit else if (strategy == "NlnN") {result = FJNS::NlnN;} else { - ColliderBit_error().raise(LOCAL_INFO, "Could not find jet strategy in list available. Please add the missing option to the FJstrategy_map function in Py8EventConversions.hpp."); + ColliderBit_error().raise(LOCAL_INFO, "Could not find jet strategy in list available. Please add the missing option to the FJstrategy_map function in ColliderBit/src/Utils.hpp."); } return result; } @@ -78,7 +78,7 @@ namespace Gambit else if (reco_scheme == "pt2_scheme") {result = FJNS::pt2_scheme;} else { - ColliderBit_error().raise(LOCAL_INFO, "Could not find jet recombination scheme in list available. Please add the missing option to the FJRecomScheme_map function in Py8EventConversions.hpp."); + ColliderBit_error().raise(LOCAL_INFO, "Could not find jet recombination scheme in list available. Please add the missing option to the FJRecomScheme_map function in ColliderBit/src/Utils.hpp."); } return result; } From 1cc7cf17e283d15fba9dfe4b62a0be57df9e1478 Mon Sep 17 00:00:00 2001 From: ChrisJChang Date: Mon, 4 Dec 2023 14:46:48 +0100 Subject: [PATCH 26/27] fix typo --- ColliderBit/src/Utils.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/ColliderBit/src/Utils.cpp b/ColliderBit/src/Utils.cpp index 0a51254fa6..f8250bc34e 100644 --- a/ColliderBit/src/Utils.cpp +++ b/ColliderBit/src/Utils.cpp @@ -53,7 +53,7 @@ namespace Gambit else if (algorithm == "cambridge_for_passive") {result = FJNS::cambridge_for_passive_algorithm;} else { - ColliderBit_error().raise(LOCAL_INFO, "Could not find jet algorithm in list available. Please add the missing option to the FJalgorithm_map function in ColliderBit/src/Utils.hpp."); + ColliderBit_error().raise(LOCAL_INFO, "Could not find jet algorithm in list available. Please add the missing option to the FJalgorithm_map function in ColliderBit/src/Utils.cpp."); } return result; } @@ -65,7 +65,7 @@ namespace Gambit else if (strategy == "NlnN") {result = FJNS::NlnN;} else { - ColliderBit_error().raise(LOCAL_INFO, "Could not find jet strategy in list available. Please add the missing option to the FJstrategy_map function in ColliderBit/src/Utils.hpp."); + ColliderBit_error().raise(LOCAL_INFO, "Could not find jet strategy in list available. Please add the missing option to the FJstrategy_map function in ColliderBit/src/Utils.cpp."); } return result; } @@ -78,7 +78,7 @@ namespace Gambit else if (reco_scheme == "pt2_scheme") {result = FJNS::pt2_scheme;} else { - ColliderBit_error().raise(LOCAL_INFO, "Could not find jet recombination scheme in list available. Please add the missing option to the FJRecomScheme_map function in ColliderBit/src/Utils.hpp."); + ColliderBit_error().raise(LOCAL_INFO, "Could not find jet recombination scheme in list available. Please add the missing option to the FJRecomScheme_map function in ColliderBit/src/Utils.cpp."); } return result; } From 51bfc93535bda3c5a983d06dfcb573b001699d8d Mon Sep 17 00:00:00 2001 From: Anders Kvellestad Date: Mon, 4 Dec 2023 16:28:54 +0100 Subject: [PATCH 27/27] Fixed small inconsistency and updated some comments. --- ColliderBit/include/gambit/ColliderBit/getPy8Collider.hpp | 2 +- ColliderBit/src/getHepMCEvent.cpp | 4 ++-- ColliderBit/src/getLHEvent.cpp | 2 +- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/ColliderBit/include/gambit/ColliderBit/getPy8Collider.hpp b/ColliderBit/include/gambit/ColliderBit/getPy8Collider.hpp index 978c414f59..04ddc385a6 100644 --- a/ColliderBit/include/gambit/ColliderBit/getPy8Collider.hpp +++ b/ColliderBit/include/gambit/ColliderBit/getPy8Collider.hpp @@ -126,7 +126,7 @@ namespace Gambit } result.jetcollection_taus = colOptions.getValue("jet_collection_taus"); - // Throw an error if the jetcollection_taus setting is not given and not using the antikt_R04 collection + // Throw an error if the "jet_collection_taus" setting does not match an entry in "jet_collections". if (std::find(jetcollection_names.begin(), jetcollection_names.end(), result.jetcollection_taus) == jetcollection_names.end()) { ColliderBit_error().raise(LOCAL_INFO,"Please provide the jet_collection_taus setting for jet collections."); diff --git a/ColliderBit/src/getHepMCEvent.cpp b/ColliderBit/src/getHepMCEvent.cpp index 21d0af1790..67a92c0429 100644 --- a/ColliderBit/src/getHepMCEvent.cpp +++ b/ColliderBit/src/getHepMCEvent.cpp @@ -214,8 +214,8 @@ namespace Gambit all_jet_collection_settings.push_back({key, algorithm, R, recombination_scheme, strategy}); } - jetcollection_taus = runOptions.getValueOrDef("antikt_R04", "jet_collection_taus"); - // Throw an error if the jetcollection_taus setting is not given and not using the antikt_R04 collection + jetcollection_taus = runOptions.getValue("jet_collection_taus"); + // Throw an error if the "jet_collection_taus" setting does not match an entry in "jet_collections". if (std::find(jetcollection_names.begin(), jetcollection_names.end(), jetcollection_taus) == jetcollection_names.end()) { ColliderBit_error().raise(LOCAL_INFO,"Please provide the jet_collection_taus setting for jet collections."); diff --git a/ColliderBit/src/getLHEvent.cpp b/ColliderBit/src/getLHEvent.cpp index 7c8f778129..b62c87e362 100644 --- a/ColliderBit/src/getLHEvent.cpp +++ b/ColliderBit/src/getLHEvent.cpp @@ -85,7 +85,7 @@ namespace Gambit // Note that jetcollection_taus is not used here as get_HEPUtils_event(...) has a much simpler jet definition than in Py8Conversions.hpp jetcollection_taus = colOptions.getValue("jet_collection_taus"); - // Throw an error if the jetcollection_taus setting is not given and not using the antikt_R04 collection + // Throw an error if the "jet_collection_taus" setting does not match an entry in "jet_collections". if (std::find(jetcollection_names.begin(), jetcollection_names.end(), jetcollection_taus) == jetcollection_names.end()) { ColliderBit_error().raise(LOCAL_INFO,"Please provide the jet_collection_taus setting for jet collections.");