HepUtils update: Multiple Jet Collections and Cluster sequence access

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, 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())
+      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;
     }

ColliderBit/include/gambit/ColliderBit/colliders/BaseCollider.hpp

@@ -21,6 +21,8 @@
 #include <string>
 #include <vector>
 
+#include "gambit/ColliderBit/Utils.hpp"
+
 namespace Gambit
 {
   namespace ColliderBit
@@ -33,7 +35,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 +69,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<jet_collection_settings> all_jet_collection_settings;
+
+        /// Key for jet collection used in adding taus
+        std::string jetcollection_taus;
 
     };
 

ColliderBit/include/gambit/ColliderBit/colliders/Pythia8/Py8EventConversions.hpp

@@ -14,12 +14,14 @@
 ///  \author Pat Scott
 ///  \author Martin White
 ///  \author Are Raklev  June 2021
+///  \author Chris Chang Nov 2023
 ///
 ///  *********************************************
 
 #pragma once
 
 #include "gambit/ColliderBit/colliders/EventConversionUtils.hpp"
+#include "gambit/ColliderBit/colliders/BaseCollider.hpp"
 
 //#define COLLIDERBIT_DEBUG
 
@@ -35,7 +37,7 @@ namespace Gambit
     /// Convert a hadron-level EventT into an unsmeared HEPUtils::Event
     /// @todo Overlap between jets and prompt containers: need some isolation in MET calculation
     template<typename EventT>
-    void convertParticleEvent(const EventT& pevt, HEPUtils::Event& result, double antiktR, double jet_pt_min)
+    void convertParticleEvent(const EventT& pevt, HEPUtils::Event& result, std::vector<jet_collection_settings> all_jet_collection_settings, str jetcollection_taus, double jet_pt_min)
     {
       result.clear();
 
@@ -113,7 +115,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 +133,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 +147,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,91 +195,103 @@ namespace Gambit
       }
 
       /// Jet finding
-      /// @todo Choose jet algorithm via detector _settings? Run several algs?
-      const FJNS::JetDefinition jet_def(FJNS::antikt_algorithm, antiktR);
-      FJNS::ClusterSequence cseq(jetparticles, jet_def);
-      std::vector<FJNS::PseudoJet> 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);
+
+        /// 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<const FJNS::ClusterSequence> CSeqBasePtr = result.emplace_clusterseq(jetparticles, jet_def, jetcollection.key);
+        /// Get the resulting pseudojets
+        std::vector<FJNS::PseudoJet> pjets = sorted_by_pt(CSeqBasePtr->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
-        result.add_jet(new HEPUtils::Jet(HEPUtils::mk_p4(pj), isB, isC, isW, isZ, ish));
+          // Add jet to collection including tags and PseudoJet
+          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
@@ -299,11 +318,12 @@ 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;
         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
@@ -312,7 +332,7 @@ namespace Gambit
 
     /// Convert a partonic (no hadrons) EventT into an unsmeared HEPUtils::Event
     template<typename EventT>
-    void convertPartonEvent(const EventT& pevt, HEPUtils::Event& result, double antiktR, double jet_pt_min)
+    void convertPartonEvent(const EventT& pevt, HEPUtils::Event& result, std::vector<jet_collection_settings> all_jet_collection_settings, str jetcollection_taus, double jet_pt_min)
     {
       result.clear();
 
@@ -389,40 +409,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<FJNS::PseudoJet> 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);
+        std::shared_ptr<const FJNS::ClusterSequence> CSeqBasePtr = result.emplace_clusterseq(jetparticles, jet_def, jetcollection.key);
+        std::vector<FJNS::PseudoJet> pjets = sorted_by_pt(CSeqBasePtr->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);
         }
       }