feat: Add a fastjet plugin

Core/include/Acts/EventData/TrackProxy.hpp

@@ -320,6 +320,22 @@ class TrackProxy {
   /// @return the global momentum vector
   Vector3 momentum() const { return absoluteMomentum() * direction(); }
 
+  /// Get the four-momentum vector: (px, py, pz, e)
+  /// @return the four-momentum vector
+  Vector4 fourMomentum() const {
+    Vector4 p4 = Vector4::Zero();
+
+    Vector3 p3 = momentum();
+    p4[eMom0] = p3[eMom0];
+    p4[eMom1] = p3[eMom1];
+    p4[eMom2] = p3[eMom2];
+
+    float m = particleHypothesis().mass();
+    p4[eEnergy] = std::sqrt(m * m + p3.squaredNorm());
+
+    return p4;
+  }
+
   /// Return the number of track states associated to this track
   /// @note This is calculated by iterating over the track states which is
   ///       somewhat expensive. Consider caching this value if you need It

Plugins/FastJet/CMakeLists.txt

@@ -1,17 +1,18 @@
-add_library(ActsPluginFastJet INTERFACE)
+add_library(ActsPluginFastJet SHARED src/TrackJets.cpp)
 
 target_include_directories(
     ActsPluginFastJet
-    INTERFACE
+    PUBLIC
         $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include>
         $<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>
 )
-target_link_libraries(ActsPluginFastJet INTERFACE FastJet)
+
+target_link_libraries(ActsPluginFastJet PUBLIC ActsCore FastJet)
 
 install(
     TARGETS ActsPluginFastJet
     EXPORT ActsPluginFastJetTargets
     LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR}
 )
-# Add once fastjet plugin has headers
-# install(DIRECTORY include/Acts DESTINATION ${CMAKE_INSTALL_INCLUDEDIR})
+
+install(DIRECTORY include/Acts DESTINATION ${CMAKE_INSTALL_INCLUDEDIR})

Plugins/FastJet/include/Acts/Plugins/FastJet/TrackJets.hpp

@@ -0,0 +1,88 @@
+// This file is part of the ACTS project.
+//
+// Copyright (C) 2016 CERN for the benefit of the ACTS project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at https://mozilla.org/MPL/2.0/.
+
+#pragma once
+
+#include "Acts/Definitions/Common.hpp"
+#include "Acts/Definitions/Units.hpp"
+
+#include <optional>
+#include <vector>
+
+#include <fastjet/JetDefinition.hh>
+#include <fastjet/PseudoJet.hh>
+
+namespace Acts::FastJet {
+
+/// Default jet definition: Anti-kt with a radius of 0.4
+const fastjet::JetDefinition DefaultJetDefinition =
+    fastjet::JetDefinition(fastjet::antikt_algorithm, 0.4);
+
+template <typename TrackContainer>
+class InputTracks {
+ public:
+  /// Constructor; saves a reference to the container
+  /// @param tracks the TrackContainer
+  explicit InputTracks(TrackContainer& tracks) : m_tracks{tracks} {}
+
+  /// Get vector for 4-momenta from the track collection
+  /// @return vector of fastjet::PseudoJet, one per track
+  std::vector<fastjet::PseudoJet> fourMomenta() const;
+
+  /// Get the tracks making up a track-jet
+  ///
+  /// @param jet the jet from which to get the constituent tracks
+  /// @param coreR optional radius inside which to get the tracks
+  ///
+  /// @return a vector of TrackProxy
+  std::vector<typename TrackContainer::TrackProxy> tracksInJet(
+      const fastjet::PseudoJet& jet, std::optional<float> coreR = {});
+
+ private:
+  TrackContainer& m_tracks;
+};
+
+class TrackJetSequence {
+ public:
+  /// Factory function to create a sequence of track jets
+  ///
+  /// @param tracks the input tracks
+  /// @jetDef the jet definition to use, defaults to "DefaultJetDefinition"
+  static TrackJetSequence create(
+      std::vector<fastjet::PseudoJet>& fourMomenta,
+      fastjet::JetDefinition jetDef = DefaultJetDefinition);
+
+  static TrackJetSequence create(
+      std::vector<fastjet::PseudoJet>&& fourMomenta,
+      fastjet::JetDefinition jetDef = DefaultJetDefinition) {
+    return create(fourMomenta, std::move(jetDef));
+  }
+
+  /// Get all the track jets passing the pT & eta cuts
+  ///
+  /// @param ptMin the minimum jet pT in GeV
+  /// @param etaMax the maximum jet absolute eta
+  ///
+  /// @return a vector of fastjet::PseudoJet objects
+  std::vector<fastjet::PseudoJet> jets(float ptMin = 20 *
+                                                     Acts::UnitConstants::GeV,
+                                       float etaMax = 2.5);
+
+ private:
+  /// Main constructor. Users should call "TrackJetSequence::create" instead
+  ///
+  /// @param clusterSeq the fastjet::ClusterSequence object
+  explicit TrackJetSequence(const fastjet::ClusterSequence& clusterSeq)
+      : m_clusterSeq{clusterSeq} {}
+
+  fastjet::ClusterSequence m_clusterSeq{};
+};
+
+}  // namespace Acts::FastJet
+
+#include "TrackJets.ipp"

Plugins/FastJet/include/Acts/Plugins/FastJet/TrackJets.ipp

@@ -0,0 +1,49 @@
+// This file is part of the ACTS project.
+//
+// Copyright (C) 2016 CERN for the benefit of the ACTS project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at https://mozilla.org/MPL/2.0/.
+
+#pragma once
+
+#include "Acts/Definitions/Algebra.hpp"
+
+namespace Acts::FastJet {
+
+template <typename TrackContainer>
+std::vector<fastjet::PseudoJet> InputTracks<TrackContainer>::fourMomenta()
+    const {
+  std::vector<fastjet::PseudoJet> inputs;
+  for (std::size_t i = 0; i < m_tracks.size(); i++) {
+    Acts::Vector4 p = m_tracks.getTrack(i).fourMomentum();
+    inputs.emplace_back(p[Acts::eMom0], p[Acts::eMom1], p[Acts::eMom2],
+                        p[Acts::eEnergy]);
+    inputs.back().set_user_index(i);
+  }
+  return inputs;
+}
+
+template <typename TrackContainer>
+std::vector<typename TrackContainer::TrackProxy>
+InputTracks<TrackContainer>::tracksInJet(const fastjet::PseudoJet& jet,
+                                         std::optional<float> coreR) {
+  fastjet::Selector sel = fastjet::SelectorIdentity();
+  if (coreR.has_value()) {
+    if (*coreR < 0) {
+      throw std::invalid_argument("coreR must be positive!");
+    }
+    sel = fastjet::SelectorCircle(*coreR);
+    sel.set_reference(jet);
+  }
+
+  std::vector<typename TrackContainer::TrackProxy> tracks;
+  for (fastjet::PseudoJet& cst : sel(jet.constituents())) {
+    tracks.push_back(m_tracks.getTrack(cst.user_index()));
+  }
+
+  return tracks;
+}
+
+}  // namespace Acts::FastJet

Plugins/FastJet/src/TrackJets.cpp

@@ -0,0 +1,25 @@
+// This file is part of the ACTS project.
+//
+// Copyright (C) 2016 CERN for the benefit of the ACTS project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at https://mozilla.org/MPL/2.0/.
+
+#include "Acts/Plugins/FastJet/TrackJets.hpp"
+
+namespace Acts::FastJet {
+
+TrackJetSequence TrackJetSequence::create(
+    std::vector<fastjet::PseudoJet>& inputs, fastjet::JetDefinition jetDef) {
+  fastjet::ClusterSequence cs(inputs, jetDef);
+  return TrackJetSequence(std::move(cs));
+}
+
+std::vector<fastjet::PseudoJet> TrackJetSequence::jets(float ptMin,
+                                                       float etaMax) {
+  fastjet::Selector sel_eta = fastjet::SelectorAbsEtaMax(etaMax);
+  return sel_eta(m_clusterSeq.inclusive_jets(ptMin));
+}
+
+}  // namespace Acts::FastJet

Tests/UnitTests/Plugins/CMakeLists.txt

@@ -2,6 +2,7 @@ add_subdirectory_if(ActSVG ACTS_BUILD_PLUGIN_ACTSVG)
 add_subdirectory_if(Detray ACTS_BUILD_PLUGIN_TRACCC)
 add_subdirectory_if(DD4hep ACTS_BUILD_PLUGIN_DD4HEP)
 add_subdirectory_if(ExaTrkX ACTS_BUILD_PLUGIN_EXATRKX)
+add_subdirectory_if(FastJet ACTS_BUILD_PLUGIN_FASTJET)
 add_subdirectory_if(Geant4 ACTS_BUILD_PLUGIN_GEANT4)
 add_subdirectory_if(GeoModel ACTS_BUILD_PLUGIN_GEOMODEL)
 add_subdirectory_if(Json ACTS_BUILD_PLUGIN_JSON)

Tests/UnitTests/Plugins/FastJet/CMakeLists.txt

@@ -0,0 +1,2 @@
+set(unittest_extra_libraries ActsPluginFastJet)
+add_unittest(TrackJetsTests TrackJetsTests.cpp)

Tests/UnitTests/Plugins/FastJet/TrackJetsTests.cpp

@@ -0,0 +1,171 @@
+// This file is part of the ACTS project.
+//
+// Copyright (C) 2016 CERN for the benefit of the ACTS project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at https://mozilla.org/MPL/2.0/.
+
+#include <boost/test/unit_test.hpp>
+
+#include <Acts/Plugins/FastJet/TrackJets.hpp>
+
+class Track {
+ public:
+  static constexpr float mass = 139.57061 * Acts::UnitConstants::MeV;
+
+  Track(float pt, float eta, float phi) {
+    Acts::Vector3 p3 = Acts::Vector3::Zero();
+    p3[0] = pt * std::cos(phi);
+    p3[1] = pt * std::sin(phi);
+    p3[2] = pt * std::sinh(eta);
+    float e = std::sqrt(mass * mass + p3.squaredNorm());
+    m_fourMom[0] = p3[0];
+    m_fourMom[1] = p3[1];
+    m_fourMom[2] = p3[2];
+    m_fourMom[3] = e;
+  }
+
+  Acts::Vector4 fourMomentum() const { return m_fourMom; }
+
+ private:
+  Acts::Vector4 m_fourMom{};
+};
+
+bool operator==(Track const& lhs, Track const& rhs) {
+  return lhs.fourMomentum() == rhs.fourMomentum();
+}
+
+class TrackContainer {
+ public:
+  using TrackProxy = Track;
+
+  TrackContainer() = default;
+  void insert(Track track) { m_vec.push_back(std::move(track)); }
+  std::size_t size() { return m_vec.size(); }
+
+  using ConstTrackProxy = const Track&;
+  ConstTrackProxy getTrack(std::size_t i) {
+    if (i < size()) {
+      return m_vec[i];
+    }
+    throw std::runtime_error("Too few tracks");
+  }
+
+ private:
+  std::vector<Track> m_vec{};
+};
+
+BOOST_AUTO_TEST_CASE(SingleTrack) {
+  TrackContainer tracks;
+  tracks.insert(Track(100, 0, 0));
+
+  Acts::FastJet::InputTracks inputTracks(tracks);
+
+  Acts::FastJet::TrackJetSequence jetSeq =
+      Acts::FastJet::TrackJetSequence::create(inputTracks.fourMomenta());
+  std::vector<fastjet::PseudoJet> jets = jetSeq.jets();
+
+  BOOST_CHECK_EQUAL(jets.size(), 1);
+  BOOST_CHECK_EQUAL(jets[0].constituents().size(), 1);
+  BOOST_CHECK_EQUAL(jets[0].constituents()[0].user_index(), 0);
+  BOOST_CHECK_CLOSE(jets[0].pt(), 100, 1e-3);
+  BOOST_CHECK_CLOSE(jets[0].eta(), 0, 1e-3);
+  BOOST_CHECK_CLOSE(jets[0].phi(), 0, 1e-3);
+  BOOST_CHECK_CLOSE(jets[0].m(), Track::mass, 1);
+}
+
+BOOST_AUTO_TEST_CASE(TwoTracksTwoJets) {
+  TrackContainer tracks;
+  tracks.insert(Track(100, 0, 0.0));
+  tracks.insert(Track(100, 0, std::numbers::pi));
+
+  Acts::FastJet::InputTracks inputTracks(tracks);
+
+  Acts::FastJet::TrackJetSequence jetSeq =
+      Acts::FastJet::TrackJetSequence::create(inputTracks.fourMomenta());
+  std::vector<fastjet::PseudoJet> jets = jetSeq.jets();
+
+  BOOST_CHECK_EQUAL(jets.size(), 2);
+
+  std::vector<Track> trks_0 = inputTracks.tracksInJet(jets[0]);
+  BOOST_CHECK_EQUAL(trks_0.size(), 1);
+  BOOST_CHECK(trks_0[0] == tracks.getTrack(0) ||
+              trks_0[0] == tracks.getTrack(1));
+
+  std::vector<Track> trks_1 = inputTracks.tracksInJet(jets[1]);
+  BOOST_CHECK_EQUAL(trks_1.size(), 1);
+  BOOST_CHECK(trks_1[0] == tracks.getTrack(0) ||
+              trks_1[0] == tracks.getTrack(1));
+  BOOST_CHECK(trks_0[0] != trks_1[0]);
+}
+
+BOOST_AUTO_TEST_CASE(TwoTracksOneJet) {
+  TrackContainer tracks;
+  tracks.insert(Track(100, 0, 0.0));
+  tracks.insert(Track(100, 0, 0.2));
+
+  Acts::FastJet::InputTracks inputTracks(tracks);
+
+  Acts::FastJet::TrackJetSequence jetSeq =
+      Acts::FastJet::TrackJetSequence::create(inputTracks.fourMomenta());
+  std::vector<fastjet::PseudoJet> jets = jetSeq.jets();
+
+  BOOST_CHECK_EQUAL(jets.size(), 1);
+
+  std::vector<Track> trks_0 = inputTracks.tracksInJet(jets[0]);
+  BOOST_CHECK_EQUAL(trks_0.size(), 2);
+  BOOST_CHECK(trks_0[0] == tracks.getTrack(0) ||
+              trks_0[0] == tracks.getTrack(1));
+  BOOST_CHECK(trks_0[1] == tracks.getTrack(0) ||
+              trks_0[1] == tracks.getTrack(1));
+  BOOST_CHECK(trks_0[0] != trks_0[1]);
+}
+
+BOOST_AUTO_TEST_CASE(TracksInJetCore) {
+  TrackContainer tracks;
+  tracks.insert(Track(100, 0, 0));
+  tracks.insert(Track(10, 0.05, 0));
+  tracks.insert(Track(10, -0.05, 0));
+  tracks.insert(Track(10, 0.2, 0));
+  tracks.insert(Track(10, -0.2, 0));
+
+  Acts::FastJet::InputTracks inputTracks(tracks);
+
+  Acts::FastJet::TrackJetSequence jetSeq =
+      Acts::FastJet::TrackJetSequence::create(inputTracks.fourMomenta());
+  std::vector<fastjet::PseudoJet> jets = jetSeq.jets();
+
+  BOOST_REQUIRE_EQUAL(jets.size(), 1);
+
+  std::vector<Track> trks = inputTracks.tracksInJet(jets[0], 0.1);
+  BOOST_CHECK_EQUAL(trks.size(), 3);
+
+  BOOST_CHECK(std::find(trks.begin(), trks.end(), tracks.getTrack(0)) !=
+              trks.end());
+  BOOST_CHECK(std::find(trks.begin(), trks.end(), tracks.getTrack(1)) !=
+              trks.end());
+  BOOST_CHECK(std::find(trks.begin(), trks.end(), tracks.getTrack(2)) !=
+              trks.end());
+  BOOST_CHECK(std::find(trks.begin(), trks.end(), tracks.getTrack(3)) ==
+              trks.end());
+  BOOST_CHECK(std::find(trks.begin(), trks.end(), tracks.getTrack(4)) ==
+              trks.end());
+}
+
+BOOST_AUTO_TEST_CASE(EmptyTrackContainer) {
+  Acts::FastJet::TrackJetSequence jetSeq =
+      Acts::FastJet::TrackJetSequence::create(
+          std::vector<fastjet::PseudoJet>());
+  BOOST_CHECK_EQUAL(jetSeq.jets().size(), 0);
+}
+
+BOOST_AUTO_TEST_CASE(InvalidCoreRadius) {
+  TrackContainer tracks;
+  tracks.insert(Track(100, 0, 0));
+  Acts::FastJet::InputTracks inputTracks(tracks);
+  Acts::FastJet::TrackJetSequence jetSeq =
+      Acts::FastJet::TrackJetSequence::create(inputTracks.fourMomenta());
+  BOOST_CHECK_THROW(inputTracks.tracksInJet(jetSeq.jets()[0], -1.0),
+                    std::invalid_argument);
+}