Upstream FastJet algorithm from MuonCollider fork

CMakeLists.txt

@@ -26,6 +26,7 @@ set(${PROJECT_NAME}_VERSION_MINOR 4)
 set(${PROJECT_NAME}_VERSION_PATCH 0)
 
 set(${PROJECT_NAME}_VERSION "${${PROJECT_NAME}_VERSION_MAJOR}.${${PROJECT_NAME}_VERSION_MINOR}.${${PROJECT_NAME}_VERSION_PATCH}")
+set(CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/cmake ${CMAKE_MODULE_PATH})
 
 option(BUILD_TRACKING "Build conformal tracking" ON)
 
@@ -37,6 +38,8 @@ find_package(k4FWCore 1.4 REQUIRED)
 find_package(ROOT REQUIRED COMPONENTS RIO Tree)
 find_package(k4SimGeant4 REQUIRED)
 
+find_package(FastJet REQUIRED)
+
 include(cmake/Key4hepConfig.cmake)
 
 include(GNUInstallDirs)

README.md

@@ -1,3 +1,14 @@
 # k4Reco
 
 Gaudi algorithms for reconstruction using EDM4hep natively.
+
+## Example steering files
+
+Example steering files that show how to configure and run the algorithms live in the
+`options` directory of the corresponding subfolder of [`k4Reco`](k4Reco), e.g.
+[`k4Reco/Tracking/options`](k4Reco/Tracking/options). The ones that are only used to
+run and validate the algorithms in the tests are kept next to the tests instead, e.g.
+[`test/FastJet`](test/FastJet) for the `FastJetAlg`.
+
+In either case these files are examples and the configurations that we use for testing,
+they are not blessed production configurations.

cmake/FindFastJet.cmake

@@ -0,0 +1,29 @@
+# - Locate FastJet library
+# Defines:
+#
+#  FASTJET_FOUND
+#  FASTJET_INCLUDE_DIR
+#  FASTJET_INCLUDE_DIRS (not cached)
+#  FASTJET_LIBRARY
+#  FASTJET_LIBRARIES (not cached)
+#  FASTJET_LIBRARY_DIRS (not cached)
+
+find_path(FASTJET_INCLUDE_DIR fastjet/version.hh
+          HINTS $ENV{FASTJET_ROOT_DIR}/include ${FASTJET_ROOT_DIR}/include)
+
+find_library(FASTJET_LIBRARY NAMES fastjet
+             HINTS $ENV{FASTJET_ROOT_DIR}/lib ${FASTJET_ROOT_DIR}/lib)
+
+find_library(FASTJETPLUGINS_LIBRARY NAMES fastjetplugins
+             HINTS $ENV{FASTJET_ROOT_DIR}/lib ${FASTJET_ROOT_DIR}/lib)
+
+# handle the QUIETLY and REQUIRED arguments and set FASTJET_FOUND to TRUE if
+# all listed variables are TRUE
+INCLUDE(FindPackageHandleStandardArgs)
+FIND_PACKAGE_HANDLE_STANDARD_ARGS(FastJet DEFAULT_MSG FASTJET_INCLUDE_DIR FASTJET_LIBRARY)
+
+mark_as_advanced(FASTJET_FOUND FASTJET_INCLUDE_DIR FASTJET_LIBRARY)
+
+set(FASTJET_INCLUDE_DIRS ${FASTJET_INCLUDE_DIR})
+set(FASTJET_LIBRARIES ${FASTJET_LIBRARY} ${FASTJETPLUGINS_LIBRARY})
+get_filename_component(FASTJET_LIBRARY_DIRS ${FASTJET_LIBRARY} PATH)

k4Reco/CMakeLists.txt

@@ -19,6 +19,7 @@ limitations under the License.
 
 set(sources
   DDPlanarDigi/components/DDPlanarDigi.cpp
+  FastJet/src/FastJetAlg.cpp
   GaudiLumiCalClusterer/components/GaudiLumiCalClusterer.cpp
   GaudiLumiCalClusterer/src/LumiCalClusterer.cpp
   GaudiLumiCalClusterer/src/LumiCalClusterer_getCalHits.cpp
@@ -44,6 +45,7 @@ set(libs
         ROOT::Core
         ROOT::MathCore
         ROOT::Tree
+        ${FASTJET_LIBRARIES}
 )
 
 if(BUILD_TRACKING)
@@ -123,6 +125,11 @@ target_include_directories(k4RecoPlugins PUBLIC
   $<BUILD_INTERFACE:${CMAKE_CURRENT_LIST_DIR}/Tracking/include>
   $<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>)
 
+# Only the FastJetAlg sources use these, nothing links against a Gaudi module
+target_include_directories(k4RecoPlugins PRIVATE
+  $<BUILD_INTERFACE:${CMAKE_CURRENT_LIST_DIR}/FastJet/include>)
+target_include_directories(k4RecoPlugins SYSTEM PRIVATE ${FASTJET_INCLUDE_DIRS})
+
 if(BUILD_TRACKING)
   target_include_directories(k4RecoPlugins
     PUBLIC
@@ -136,7 +143,12 @@ if(BUILD_TRACKING)
   )
 endif()
 
-install(TARGETS GaudiTrkUtils k4RecoPlugins
+set(install_targets k4RecoPlugins)
+if(BUILD_TRACKING)
+  list(APPEND install_targets GaudiTrkUtils)
+endif()
+
+install(TARGETS ${install_targets}
   EXPORT k4RecoTargets
   RUNTIME DESTINATION "${CMAKE_INSTALL_BINDIR}" COMPONENT bin
   LIBRARY DESTINATION "${CMAKE_INSTALL_LIBDIR}" COMPONENT shlib

k4Reco/FastJet/include/EClusterMode.h

@@ -0,0 +1,55 @@
+/*
+ * Copyright (c) 2020-2024 Key4hep-Project.
+ *
+ * This file is part of Key4hep.
+ * See https://key4hep.github.io/key4hep-doc/ for further info.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+#ifndef K4RECO_FASTJET_ECLUSTERMODE_H
+#define K4RECO_FASTJET_ECLUSTERMODE_H
+
+#include <iostream>
+
+// The enum, name and value of the enum for the Cluster Mode
+namespace k4Reco::FastJet {
+enum EClusterMode {
+  NONE = 0,
+  FJ_exclusive_yCut = 1,     // exclusive clustering mode implemented in FastJet
+  FJ_exclusive_nJets = 2,    // exclusive clustering mode implemented in FastJet
+  FJ_inclusive = 4,          // inclusive "-"
+  OWN_inclusiveIteration = 8 // use FJ inclusive Clustering, but iterate until we have the desired number of jets
+};
+std::ostream& operator<<(std::ostream& out, EClusterMode& m) {
+  switch (m) {
+  case OWN_inclusiveIteration:
+    out << "InclusiveIterativeNJets";
+    break;
+  case FJ_inclusive:
+    out << "Inclusive";
+    break;
+  case FJ_exclusive_nJets:
+    out << "ExclusiveNJets";
+    break;
+  case FJ_exclusive_yCut:
+    out << "ExclusiveYCut";
+    break;
+  default:
+    out << "unknown";
+    break;
+  }
+  return out;
+}
+
+} // namespace k4Reco::FastJet
+#endif

k4Reco/FastJet/include/FastJetAlg.hxx

@@ -0,0 +1,218 @@
+/*
+ * Copyright (c) 2020-2024 Key4hep-Project.
+ *
+ * This file is part of Key4hep.
+ * See https://key4hep.github.io/key4hep-doc/ for further info.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+/*
+ * FastJetAlg.h
+ *
+ *  Created on: 25.05.2010
+ *      Author: Lars Weuste (MPP Munich) - weuste@mpp.mpg.de
+ *			iterative inclusive algorithm based on design by Marco Battaglia (CERN) -
+ * Marco.Battaglia@cern.ch Converted to Gaudi on: 25.08.2025 Conversion: Samuel Ferraro - samuel.rowles.ferraro@cern.ch
+ */
+
+#ifndef FASTJETALG_H_
+#define FASTJETALG_H_
+
+#include "EClusterMode.h"
+
+#include <edm4hep/ReconstructedParticleCollection.h>
+#include <k4FWCore/Transformer.h>
+
+// FastJet
+#include <fastjet/CDFJetCluPlugin.hh>
+#include <fastjet/CDFMidPointPlugin.hh>
+#include <fastjet/ClusterSequence.hh>
+#include <fastjet/EECambridgePlugin.hh>
+#include <fastjet/JadePlugin.hh>
+#include <fastjet/JetDefinition.hh>
+#include <fastjet/NestedDefsPlugin.hh>
+#include <fastjet/PseudoJet.hh>
+#include <fastjet/SISConePlugin.hh>
+#include <fastjet/SISConeSphericalPlugin.hh>
+
+#include <map>
+#include <stdexcept>
+#include <string>
+#include <vector>
+
+namespace k4Reco::FastJet {
+constexpr static int ITERATIVE_INCLUSIVE_MAX_ITERATIONS = 20;
+
+static const std::map<std::string, fastjet::JetAlgorithm> NAME_TO_ALGORITHM_MAP = {
+    {"kt_algorithm", fastjet::kt_algorithm},
+    {"cambridge_algorithm", fastjet::cambridge_algorithm},
+    {"antikt_algorithm", fastjet::antikt_algorithm},
+    {"genkt_algorithm", fastjet::genkt_algorithm},
+    {"cambridge_for_passive_algorithm", fastjet::cambridge_for_passive_algorithm},
+    {"genkt_for_passive_algorithm", fastjet::genkt_for_passive_algorithm},
+    {"ee_kt_algorithm", fastjet::ee_kt_algorithm},
+    {"ee_genkt_algorithm", fastjet::ee_genkt_algorithm},
+};
+static const std::map<std::string, fastjet::RecombinationScheme> NAME_TO_RECO_SCHEME_MAP = {
+    {"E_scheme", fastjet::E_scheme},         {"pt_scheme", fastjet::pt_scheme},   {"pt2_scheme", fastjet::pt2_scheme},
+    {"Et_scheme", fastjet::Et_scheme},       {"Et2_scheme", fastjet::Et2_scheme}, {"BIpt_scheme", fastjet::BIpt_scheme},
+    {"BIpt2_scheme", fastjet::BIpt2_scheme},
+};
+static const std::map<std::string, int> NAME_TO_NR_PARAMS_MAP = {
+    {"kt_algorithm", 1},
+    {"cambridge_algorithm", 1},
+    {"antikt_algorithm", 1},
+    {"genkt_algorithm", 2},
+    {"cambridge_for_passive_algorithm", 1},
+    {"genkt_for_passive_algorithm", 1},
+    {"ee_kt_algorithm", 0},
+    {"ee_genkt_algorithm", 2},
+    {"SISConePlugin", 2},
+    {"SISConeSphericalPlugin", 2},
+};
+static const std::map<std::string, int> NAME_TO_CLUSTER_MODE_MAP = {
+    {"kt_algorithm", FJ_inclusive | FJ_exclusive_nJets | FJ_exclusive_yCut | OWN_inclusiveIteration},
+    {"cambridge_algorithm", FJ_inclusive | FJ_exclusive_nJets | FJ_exclusive_yCut | OWN_inclusiveIteration},
+    {"antikt_algorithm", FJ_inclusive | OWN_inclusiveIteration},
+    {"genkt_algorithm", FJ_inclusive | OWN_inclusiveIteration | FJ_exclusive_nJets | FJ_exclusive_yCut},
+    {"cambridge_for_passive_algorithm", FJ_inclusive | OWN_inclusiveIteration | FJ_exclusive_nJets | FJ_exclusive_yCut},
+    {"genkt_for_passive_algorithm", FJ_inclusive | OWN_inclusiveIteration},
+    {"ee_kt_algorithm", FJ_exclusive_nJets | FJ_exclusive_yCut},
+    {"ee_genkt_algorithm", FJ_inclusive | FJ_exclusive_nJets | FJ_exclusive_yCut},
+    {"SISConePlugin", FJ_inclusive | OWN_inclusiveIteration},
+    {"SISConeSphericalPlugin", FJ_inclusive | OWN_inclusiveIteration},
+};
+
+/// The parameters that are expected for a given clustering mode, together with
+/// the EClusterMode value they map to
+struct ClusterModeParams {
+  EClusterMode mode;                   ///< The cluster mode corresponding to this name
+  std::vector<std::string> paramNames; ///< The names of the expected clusteringParams (in order)
+};
+static const std::map<std::string, ClusterModeParams> NAME_TO_CLUSTER_MODE_PARAMS_MAP = {
+    {"Inclusive", {FJ_inclusive, {"minPt"}}},
+    {"InclusiveIterativeNJets", {OWN_inclusiveIteration, {"nrJets", "minE"}}},
+    {"ExclusiveNJets", {FJ_exclusive_nJets, {"nrJets"}}},
+    {"ExclusiveYCut", {FJ_exclusive_yCut, {"yCut"}}},
+};
+
+class JetDefinitionFactory {
+  using creatorFunc = std::function<std::unique_ptr<fastjet::JetDefinition>(
+      fastjet::JetAlgorithm, const std::vector<float>&, fastjet::RecombinationScheme, fastjet::Strategy)>;
+  std::map<std::string, creatorFunc> registry;
+
+public:
+  JetDefinitionFactory();
+  std::unique_ptr<fastjet::JetDefinition> create(const std::string& type, fastjet::JetAlgorithm m_jetAlgoType,
+                                                 const std::vector<float>& params,
+                                                 fastjet::RecombinationScheme m_jetRecoScheme,
+                                                 fastjet::Strategy m_strategy);
+
+private:
+  static std::unique_ptr<fastjet::JetDefinition> useZeroParams(fastjet::JetAlgorithm m_jetAlgoType,
+                                                               const std::vector<float>&,
+                                                               fastjet::RecombinationScheme m_jetRecoScheme,
+                                                               fastjet::Strategy m_strategy) {
+    return std::make_unique<fastjet::JetDefinition>(m_jetAlgoType, m_jetRecoScheme, m_strategy);
+  }
+  static std::unique_ptr<fastjet::JetDefinition> useOneParams(fastjet::JetAlgorithm m_jetAlgoType,
+                                                              const std::vector<float>& params,
+                                                              fastjet::RecombinationScheme m_jetRecoScheme,
+                                                              fastjet::Strategy m_strategy) {
+    return std::make_unique<fastjet::JetDefinition>(m_jetAlgoType, params.at(0), m_jetRecoScheme, m_strategy);
+  }
+  static std::unique_ptr<fastjet::JetDefinition> useTwoParams(fastjet::JetAlgorithm m_jetAlgoType,
+                                                              const std::vector<float>& params,
+                                                              fastjet::RecombinationScheme m_jetRecoScheme,
+                                                              fastjet::Strategy m_strategy) {
+    return std::make_unique<fastjet::JetDefinition>(m_jetAlgoType, params.at(0), params.at(1), m_jetRecoScheme,
+                                                    m_strategy);
+  }
+};
+
+using PseudoJetList = std::vector<fastjet::PseudoJet>;
+} // namespace k4Reco::FastJet
+
+struct FastJetAlg : k4FWCore::MultiTransformer<
+                        std::tuple<edm4hep::ReconstructedParticleCollection, edm4hep::ReconstructedParticleCollection>(
+                            const edm4hep::ReconstructedParticleCollection&)> {
+public:
+  FastJetAlg(const std::string& name, ISvcLocator* svcLoc);
+
+  /// Validates the configured algorithm, clustering mode and parameters and builds the jet definition
+  StatusCode initialize();
+
+  /// Clusters the input particles into jets, returning the jets and their constituents
+  std::tuple<edm4hep::ReconstructedParticleCollection, edm4hep::ReconstructedParticleCollection>
+  operator()(const edm4hep::ReconstructedParticleCollection& inputCollection) const;
+
+private:
+  Gaudi::Property<std::string> m_jetAlgoName{this, "algorithm", "kt_algorithm",
+                                             "Name of the algorithm to use for making jets. E.g. kt_algorithm. Full "
+                                             "list of algorithms can be seen in the FastJet code."};
+  Gaudi::Property<std::vector<float>> m_jetAlgoParams{this,
+                                                      "algorithmParameters",
+                                                      {0.7},
+                                                      "Parameters required by each specific algorithm. The amount of "
+                                                      "parameters and interpretation varies per algorithm."};
+  Gaudi::Property<std::string> m_clusterModeName{this, "clusteringMode", "Inclusive",
+                                                 "Clustering mode for the algorithm. Not all clustering modes are "
+                                                 "available for all algorithms. See Fast Jet code for availabilities."};
+  Gaudi::Property<std::vector<float>> m_clusterModeParams{
+      this,
+      "clusteringParams",
+      {0.0},
+      "Parameters for the clusteringMode. \nclusterMode \"Inclusive\" takes <minPt>\nclusterMode "
+      "\"InclusiveIterativeNJets\" takes <nrJets> <minE>\n"
+      "clusterMode \"ExclusiveNJets\" takes <nrJets>\nclusterMode \"ExclusiveYCut\" takes <yCut>. Note: not all modes "
+      "are available for all algorithms."
+      "Some parameters are input as floats, despite conversion to integers."};
+
+  Gaudi::Property<std::string> m_jetRecoSchemeName{
+      this, "recombinationScheme", std::string("E_scheme"),
+      "The recombination scheme used when merging 2 particles. Usually there is no need to use anything else than "
+      "4-Vector addition: E_scheme."};
+
+  // jet algorithm
+  std::unique_ptr<fastjet::JetDefinition> m_jetAlgo;
+  fastjet::JetAlgorithm m_jetAlgoType;
+
+  // clustering mode
+  k4Reco::FastJet::EClusterMode m_clusterMode;
+
+  // jet reco scheme
+  fastjet::RecombinationScheme m_jetRecoScheme;
+
+  // jet strategy
+  std::string m_strategyName;
+  fastjet::Strategy m_strategy;
+
+  // parameters
+  unsigned m_requestedNumberOfJets;
+  double m_yCut;
+  double m_minPt;
+  double m_minE;
+
+private:
+  fastjet::JetAlgorithm getAlgoType() const;
+  bool validateParams();
+  bool validateClusterModes() const;
+  bool validateClusterModeParams() const;
+  std::unique_ptr<k4Reco::FastJet::JetDefinitionFactory> theJetDefinitionFactory =
+      std::make_unique<k4Reco::FastJet::JetDefinitionFactory>();
+};
+
+DECLARE_COMPONENT(FastJetAlg)
+
+#endif /* FASTJETALG_H_ */