PodioSource: Adding remove, merge and createPseudoJets analyzers (#427)

* PodioSource: Adding remove, merge and createPseudoJets analyzers

* Fixing formatting

* Using for loop instead of std::find

* Explicitelly turning on WITH_PODIO_DATASOURCE

* Adding EDM4HEP::edm4hep

* Building FCCAnalyses as shared library

* Formatting

.github/workflows/test.yml

@@ -28,7 +28,7 @@ jobs:
            mkdir -p build install; \
            source ${{ matrix.STACK }}; \
            cd build; \
-           cmake -DCMAKE_INSTALL_PREFIX=../install -DCMAKE_CXX_STANDARD=17 -DCMAKE_CXX_FLAGS=" -fdiagnostics-color=always " -DWITH_DD4HEP=ON -DWITH_ACTS=OFF -DWITH_ONNX=ON -G Ninja ..;'
+           cmake -DCMAKE_INSTALL_PREFIX=../install -DCMAKE_CXX_FLAGS=" -fdiagnostics-color=always " -DWITH_DD4HEP=ON -DWITH_ACTS=OFF -DWITH_ONNX=ON -DWITH_PODIO_DATASOURCE=ON -G Ninja ..;'
     - name: Compile
       run: |
         docker exec CI_container /bin/bash -c 'cd ./Package; \
@@ -52,4 +52,4 @@ jobs:
         docker exec CI_container /bin/bash -c 'cd ./Package; \
           source ${{ matrix.STACK }}; \
           source ./setup.sh; \
-          fccanalysis run examples/FCCee/higgs/mH-recoil/mumu/analysis_stage1.py --output myoutput.root --test'
+          fccanalysis run examples/FCCee/higgs/mH-recoil/mumu/analysis_stage1.py --output myoutput.root --test'

analyzers/dataframe/CMakeLists.txt

@@ -33,6 +33,8 @@ if(NOT WITH_ACTS)
 endif()
 
 if(NOT WITH_PODIO_DATASOURCE)
+  list(FILTER headers EXCLUDE REGEX "JetClusteringUtilsSource.h")
+  list(FILTER sources EXCLUDE REGEX "JetClusteringUtilsSource.cc")
   list(FILTER headers EXCLUDE REGEX "LinkSource.h")
   list(FILTER headers EXCLUDE REGEX "ReconstructedParticleSource.h")
   list(FILTER sources EXCLUDE REGEX "ReconstructedParticleSource.cc")
@@ -46,7 +48,7 @@ message(STATUS "FCCAnalyses sources:\n   ${sources}")
 message(STATUS "CMAKE_CURRENT_SOURCE_DIR  ${CMAKE_CURRENT_SOURCE_DIR}")
 message(STATUS "CMAKE_INSTALL_INCLUDEDIR  ${CMAKE_INSTALL_INCLUDEDIR}")
 
-add_library(FCCAnalyses SHARED ${sources} ${headers})
+add_library(FCCAnalyses SHARED ${sources})
 target_include_directories(FCCAnalyses PUBLIC
                            $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}>
                            $<BUILD_INTERFACE:${CMAKE_SOURCE_DIR}/addons>
@@ -65,12 +67,12 @@ target_link_libraries(FCCAnalyses PUBLIC
                       ROOT::MathCore
                       ROOT::ROOTVecOps
                       ROOT::ROOTDataFrame
-                      EDM4HEP::edm4hep
-                      EDM4HEP::edm4hepDict
-                      EDM4HEP::utils
                       podio::podio
                       podio::podioRootIO
                       podio::podioDataSource
+                      EDM4HEP::edm4hep
+                      EDM4HEP::edm4hepDict
+                      EDM4HEP::utils
                       ${ADDONS_LIBRARIES}
                       ${DELPHES_LIBRARY}
                       gfortran # todo: why necessary?
@@ -110,4 +112,4 @@ if (${ROOT_VERSION} GREATER 6)
     install(FILES
           "${PROJECT_BINARY_DIR}/analyzers/dataframe/libFCCAnalyses_rdict.pcm"
                 DESTINATION "${INSTALL_LIB_DIR}" COMPONENT dev)
-    endif()
+endif()

analyzers/dataframe/FCCAnalyses/JetClusteringUtilsSource.h

@@ -0,0 +1,32 @@
+#ifndef ANALYZERS_SOURCE_JET_CLUSTERING_UTILS_H
+#define ANALYZERS_SOURCE_JET_CLUSTERING_UTILS_H
+
+// EDM4hep
+#include "edm4hep/ReconstructedParticleCollection.h"
+
+// FastJet
+#include "FastJet/JetClustering.h"
+#include "fastjet/JetDefinition.hh"
+
+/**
+ * @brief  Jet clustering tools and utilities.
+ *
+ * This namespace contains a set of functions and utilities to perform jet
+ * clustering using FastJet clustering algorithms with pseudoJets created from
+ * EDM4hep collections of particles.
+ *
+ * Most of the analyzers work with FastJet pseudoJets.
+ */
+namespace FCCAnalyses ::PodioSource ::JetClustering {
+/**
+ * @brief  Create FastJet pseudoJets for later usage by the jet clustering
+ *         algorithm(s).
+ *
+ * @param[in] inColl  Input collection of the reconstructed particles.
+ * @return  Vector of pseudoJets.
+ */
+std::vector<fastjet::PseudoJet>
+createPseudoJets(const edm4hep::ReconstructedParticleCollection &inColl);
+} // namespace FCCAnalyses::PodioSource::JetClustering
+
+#endif /* ANALYZERS_SOURCE_JET_CLUSTERING_UTILS_H */

analyzers/dataframe/FCCAnalyses/ReconstructedParticleSource.h

@@ -190,6 +190,54 @@ getCharge(const edm4hep::ReconstructedParticleCollection &inColl);
 edm4hep::ReconstructedParticleCollection
 sortByPt(const edm4hep::ReconstructedParticleCollection &inColl);
 
+/**
+ * @brief Remove a particle from a collection.
+ *
+ * If the matching parameter is false them only the particle with the same ID
+ * will be removed from the collection.
+ * If the matching parameter is true then ID is ignored and any particle
+ * matching with the one provided will be removed.
+ *
+ * @param[in] inColl  Collection of input particles.
+ * @param[in] inPartToBeRemoved  Input particle to be removed.
+ * @param[in] matching  Use matching instead of IDs/
+ * @return  Particles remaining.
+ */
+edm4hep::ReconstructedParticleCollection
+remove(const edm4hep::ReconstructedParticleCollection &inColl,
+       const edm4hep::ReconstructedParticle &inPartToBeRemoved,
+       const bool matching = false);
+
+/**
+ * @brief Remove multiple particles from a collection.
+ *
+ * If the matching parameter is false them only the particles with the same IDs
+ * will be removed from the collection.
+ * If the matching parameter is true then ID is ignored and any particle
+ * matching with the particles provided in the remove collection will be
+ * removed.
+ *
+ * @param[in] inColl  Collection of input particles.
+ * @param[in] inPartsToBeRemoved  Collection of input particles to be removed.
+ * @param[in] matching  Use matching instead of IDs/
+ * @return  Particles remaining.
+ */
+edm4hep::ReconstructedParticleCollection
+remove(const edm4hep::ReconstructedParticleCollection &inColl,
+       const edm4hep::ReconstructedParticleCollection &inPartsToBeRemoved,
+       const bool matching = false);
+
+/**
+ * @brief Merge two particle collections together.
+ *
+ * @param[in] inColl1  First collection of input particles.
+ * @param[in] inColl2  Second collection of input particles.
+ * @return  Particles in merged collection.
+ */
+edm4hep::ReconstructedParticleCollection
+merge(const edm4hep::ReconstructedParticleCollection &inColl1,
+      const edm4hep::ReconstructedParticleCollection &inColl2);
+
 /**
  * \brief Build two particle resonances from an arbitrary list of input
  *        reconstructed particles.

analyzers/dataframe/src/JetClusteringUtilsSource.cc

@@ -0,0 +1,18 @@
+#include "FCCAnalyses/JetClusteringUtilsSource.h"
+
+namespace FCCAnalyses ::PodioSource ::JetClustering {
+// ----------------------------------------------------------------------------
+std::vector<fastjet::PseudoJet>
+createPseudoJets(const edm4hep::ReconstructedParticleCollection &inColl) {
+  std::vector<fastjet::PseudoJet> result;
+  unsigned index = 0;
+  for (const auto &particle : inColl) {
+    result.emplace_back(particle.getMomentum().x, particle.getMomentum().y,
+                        particle.getMomentum().z, particle.getEnergy());
+    result.back().set_user_index(index);
+    ++index;
+  }
+  return result;
+}
+
+} // namespace FCCAnalyses::PodioSource::JetClustering

analyzers/dataframe/src/ReconstructedParticleSource.cc

@@ -213,6 +213,107 @@ sortByPt(const edm4hep::ReconstructedParticleCollection &inColl) {
   return outColl;
 }
 
+// ----------------------------------------------------------------------------
+edm4hep::ReconstructedParticleCollection
+remove(const edm4hep::ReconstructedParticleCollection &inColl,
+       const edm4hep::ReconstructedParticle &inPartToBeRemoved,
+       const bool matching) {
+  edm4hep::ReconstructedParticleCollection inPartsToBeRemoved;
+  inPartsToBeRemoved.setSubsetCollection();
+  inPartsToBeRemoved.push_back(inPartToBeRemoved);
+
+  return remove(inColl, inPartsToBeRemoved, matching);
+}
+
+// ----------------------------------------------------------------------------
+edm4hep::ReconstructedParticleCollection
+remove(const edm4hep::ReconstructedParticleCollection &inColl,
+       const edm4hep::ReconstructedParticleCollection &inPartsToBeRemoved,
+       const bool matching) {
+  edm4hep::ReconstructedParticleCollection outColl;
+  outColl.setSubsetCollection();
+
+  if (!matching) {
+    for (const auto &particle : inColl) {
+      // TODO: Use std::find for this
+      // if (std::find(inPartsToBeRemoved.begin(), inPartsToBeRemoved.end(),
+      //               particle) != inPartsToBeRemoved.end()) {
+      bool removePart = false;
+      for (const auto &partToBeRemoved : inPartsToBeRemoved) {
+        if (particle == partToBeRemoved) {
+          removePart = true;
+        }
+      }
+      if (removePart) {
+        continue;
+      }
+      outColl.push_back(particle);
+    }
+  } else {
+    float epsilon = 1e-6;
+    for (const auto &particle1 : inColl) {
+      bool removePart = false;
+      for (const auto &particle2 : inPartsToBeRemoved) {
+        float massDiff = 0.;
+        if (particle1.getMass() > 0.) {
+          massDiff = std::fabs(particle1.getMass() - particle2.getMass()) /
+                     particle1.getMass();
+        }
+        float pXDiff = 0.;
+        if (particle1.getMomentum().x != 0.) {
+          pXDiff = std::fabs(
+              (particle1.getMomentum().x - particle2.getMomentum().x) /
+              particle1.getMomentum().x);
+        }
+        float pYDiff = 0.;
+        if (particle1.getMomentum().y != 0.) {
+          pYDiff = std::fabs(
+              (particle1.getMomentum().y - particle2.getMomentum().y) /
+              particle1.getMomentum().y);
+        }
+        float pZDiff = 0.;
+        if (particle1.getMomentum().z != 0.) {
+          pZDiff = std::fabs(
+              (particle1.getMomentum().z - particle2.getMomentum().z) /
+              particle1.getMomentum().z);
+        }
+        float chargeDiff =
+            std::fabs(particle1.getCharge() - particle2.getCharge());
+        int32_t pdgDiff = std::abs(particle1.getPDG() - particle1.getPDG());
+
+        if (massDiff < epsilon && pXDiff < epsilon && pYDiff < epsilon &&
+            pZDiff < epsilon && chargeDiff < epsilon && pdgDiff < 1) {
+          removePart = true;
+        }
+      }
+      if (removePart) {
+        continue;
+      }
+      outColl.push_back(particle1);
+    }
+  }
+
+  return outColl;
+}
+
+// ----------------------------------------------------------------------------
+edm4hep::ReconstructedParticleCollection
+merge(const edm4hep::ReconstructedParticleCollection &inColl1,
+      const edm4hep::ReconstructedParticleCollection &inColl2) {
+  edm4hep::ReconstructedParticleCollection outColl;
+  outColl.setSubsetCollection();
+
+  for (const auto &particle : inColl1) {
+    outColl.push_back(particle);
+  }
+
+  for (const auto &particle : inColl2) {
+    outColl.push_back(particle);
+  }
+
+  return outColl;
+}
+
 // --------------------------------------------------------------------------
 resonanceBuilder::resonanceBuilder(float resonanceMass)
     : m_resonanceMass(resonanceMass) {

tests/unittest/CMakeLists.txt

@@ -4,19 +4,28 @@ find_catch_instance()
 # list of labels that we want to ignore
 set(filter_tests "")
 
-add_executable(unittest unittest.cpp
-                        myutils.cpp
-                        algorithms.cpp
-                        ReconstructedParticle.cpp
+list(APPEND unittest_src unittest.cpp
+                         myutils.cpp
+                         algorithms.cpp
+                         ReconstructedParticle.cpp
 )
-target_link_libraries(unittest PUBLIC FCCAnalyses
-                                      gfortran
-                               PRIVATE Catch2::Catch2WithMain)
+
+if(WITH_PODIO_DATASOURCE)
+  list(APPEND unittest_src JetClusteringUtilsSource.cpp
+                           ReconstructedParticleSource.cpp
+  )
+endif()
+
+add_executable(unittest ${unittest_src})
+add_dependencies(unittest FCCAnalyses)
+target_link_libraries(unittest PRIVATE Catch2::Catch2WithMain
+                               PUBLIC FCCAnalyses)
 target_include_directories(unittest PUBLIC ${VDT_INCLUDE_DIR})
 target_compile_definitions(unittest PUBLIC "-DTEST_INPUT_DATA_DIR=${TEST_INPUT_DATA_DIR}")
 include(Catch)
 catch_discover_tests(unittest
     WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}
     TEST_PREFIX "UT_" # make it possible to filter easily with -R ^UT
     TEST_SPEC ${filter_tests} # discover only tests that are known to not fail
+    TEST_DL_PATHS $<TARGET_FILE_DIR:FCCAnalyses> # Path to the current FCCAnalyses shared library
 )

tests/unittest/JetClusteringUtilsSource.cpp

@@ -0,0 +1,35 @@
+#include "FCCAnalyses/JetClusteringUtilsSource.h"
+
+// Catch2
+#include "catch2/catch_test_macros.hpp"
+#include <catch2/catch_approx.hpp>
+
+TEST_CASE("create-pseudojets", "[JetClusteringUtilsSource]") {
+  edm4hep::ReconstructedParticleCollection pColl;
+  edm4hep::MutableReconstructedParticle p1;
+  p1.setPDG(11);
+  p1.setMomentum({20., 30., 40.});
+  p1.setMass(5.486e-4);
+  p1.setEnergy(53.85164807);
+  pColl.push_back(p1);
+  edm4hep::MutableReconstructedParticle p2;
+  p2.setPDG(13);
+  p2.setMomentum({10., 20., 30.});
+  p2.setMass(1.05658e-1);
+  p2.setEnergy(37.416723);
+  pColl.push_back(p2);
+  auto res = FCCAnalyses::PodioSource::JetClustering::createPseudoJets(pColl);
+  REQUIRE(res.size() == 2);
+  REQUIRE(res[0].px() == Catch::Approx(p1.getMomentum().x));
+  REQUIRE(res[0].py() == Catch::Approx(p1.getMomentum().y));
+  REQUIRE(res[0].pz() == Catch::Approx(p1.getMomentum().z));
+  REQUIRE(res[0].e() == Catch::Approx(p1.getEnergy()));
+  REQUIRE(res[0].m2() ==
+          Catch::Approx(p1.getMass() * p1.getMass()).margin(1.e-3));
+  REQUIRE(res[1].px() == Catch::Approx(p2.getMomentum().x));
+  REQUIRE(res[1].py() == Catch::Approx(p2.getMomentum().y));
+  REQUIRE(res[1].pz() == Catch::Approx(p2.getMomentum().z));
+  REQUIRE(res[1].e() == Catch::Approx(p2.getEnergy()));
+  REQUIRE(res[1].m2() ==
+          Catch::Approx(p2.getMass() * p2.getMass()).margin(1.e-3));
+}

tests/unittest/ReconstructedParticle.cpp

@@ -1,8 +1,9 @@
+// FCCAnalyses
 #include "FCCAnalyses/ReconstructedParticle.h"
 
 // Catch2
+#include "catch2/catch_approx.hpp"
 #include "catch2/catch_test_macros.hpp"
-#include <catch2/catch_approx.hpp>
 
 // EDM4hep
 #include "edm4hep/EDM4hepVersion.h"