Convert ParticleIDs also from EDM4hep to LCIO

k4EDM4hep2LcioConv/include/k4EDM4hep2LcioConv/k4EDM4hep2LcioConv.h

@@ -85,6 +85,7 @@ namespace EDM4hep2LCIOConv {
     ObjectMapT<lcio::VertexImpl*, edm4hep::Vertex> vertices {};
     ObjectMapT<lcio::ReconstructedParticleImpl*, edm4hep::ReconstructedParticle> recoParticles {};
     ObjectMapT<lcio::MCParticleImpl*, edm4hep::MCParticle> mcParticles {};
+    ObjectMapT<lcio::ParticleIDImpl*, edm4hep::ParticleID> particleIDs {};
   };
 
   /**
@@ -321,6 +322,15 @@ namespace EDM4hep2LCIOConv {
     return convertMCParticles(mcparticle_coll, mc_particles_vec).release();
   }
 
+  /**
+   * Convert EDM4hep ParticleIDs to LCIO. NOTE: Since ParticleIDs cannot live in
+   * their own collections in LCIO this simply populates the pidMap that maps
+   * LCIO to EDM4hep particleIDs. **This just converts the data it is crucial to
+   * also resolve the relations afterwards!**
+   */
+  template<typename PidMapT>
+  void convertParticleIDs(const edm4hep::ParticleIDCollection* const edmCollection, PidMapT& pidMap, const int algoId);
+
   /**
    * Convert EDM4hep EventHeader to LCIO. This will directly populate the
    * corresponding information in the passed LCEvent. The input collection needs
@@ -389,6 +399,9 @@ namespace EDM4hep2LCIOConv {
   template<typename ClusterMapT, typename CaloHitMapT>
   void resolveRelationsClusters(ClusterMapT& clustersMap, const CaloHitMapT& caloHitMap);
 
+  template<typename PidMapT, typename RecoParticleMapT>
+  void resolveRelationsParticleIDs(PidMapT& pidMap, RecoParticleMapT& recoMap);
+
   /**
    * Resolve all relations in all converted objects that are held in the map.
    * Dispatch to the correpsonding implementation for all the types that have

k4EDM4hep2LcioConv/include/k4EDM4hep2LcioConv/k4EDM4hep2LcioConv.ipp

@@ -10,7 +10,6 @@ namespace EDM4hep2LCIOConv {
     TrackMapT& trackMap)
   {
     auto tracks = std::make_unique<lcio::LCCollectionVec>(lcio::LCIO::TRACK);
-
     // Loop over EDM4hep tracks converting them to lcio tracks.
     for (const auto& edm_tr : (*edmCollection)) {
       if (edm_tr.isAvailable()) {
@@ -520,6 +519,22 @@ namespace EDM4hep2LCIOConv {
     return mcparticles;
   }
 
+  template<typename PidMapT>
+  void convertParticleIDs(const edm4hep::ParticleIDCollection* const edmCollection, PidMapT& pidMap, const int algoId)
+  {
+    for (const auto& edmPid : (*edmCollection)) {
+      auto [lcioPid, _] = k4EDM4hep2LcioConv::detail::mapInsert(new lcio::ParticleIDImpl(), edmPid, pidMap).first;
+
+      lcioPid->setType(edmPid.getType());
+      lcioPid->setPDG(edmPid.getPDG());
+      lcioPid->setLikelihood(edmPid.getLikelihood());
+      lcioPid->setAlgorithmType(algoId);
+      for (const auto& param : edmPid.getParameters()) {
+        lcioPid->addParameter(param);
+      }
+    }
+  }
+
   template<typename MCParticleMapT, typename MCParticleLookupMapT>
   void resolveRelationsMCParticles(MCParticleMapT& mcParticlesMap, const MCParticleLookupMapT& lookupMap)
   {
@@ -708,6 +723,21 @@ namespace EDM4hep2LCIOConv {
     }
   }
 
+  template<typename PidMapT, typename RecoParticleMapT>
+  void resolveRelationsParticleIDs(PidMapT& pidMap, RecoParticleMapT& recoMap)
+  {
+    for (auto& [lcioPid, edmPid] : pidMap) {
+      const auto edmReco = edmPid.getParticle();
+      const auto lcioReco = k4EDM4hep2LcioConv::detail::mapLookupFrom(edmReco, recoMap);
+      if (lcioReco) {
+        lcioReco.value()->addParticleID(lcioPid);
+      }
+      else {
+        std::cerr << "Cannot find a reconstructed particle to attach a ParticleID to" << std::endl;
+      }
+    }
+  }
+
   template<typename ObjectMappingT>
   void resolveRelations(ObjectMappingT& collection_pairs)
   {
@@ -729,6 +759,7 @@ namespace EDM4hep2LCIOConv {
       lookup_pairs.vertices,
       lookup_pairs.clusters,
       lookup_pairs.tracks);
+    resolveRelationsParticleIDs(lookup_pairs.particleIDs, update_pairs.recoParticles);
     resolveRelationsVertices(update_pairs.vertices, lookup_pairs.recoParticles);
     resolveRelationsSimCaloHit(update_pairs.simCaloHits, lookup_pairs.mcParticles);
     resolveRelationsSimTrackerHits(update_pairs.simTrackerHits, lookup_pairs.mcParticles);

k4EDM4hep2LcioConv/src/k4EDM4hep2LcioConv.cpp

@@ -1,6 +1,6 @@
 #include "k4EDM4hep2LcioConv/k4EDM4hep2LcioConv.h"
 #include "edm4hep/Constants.h"
-#include "EVENT/MCParticle.h"
+#include "UTIL/PIDHandler.h"
 
 namespace EDM4hep2LCIOConv {
 
@@ -29,11 +29,22 @@ namespace EDM4hep2LCIOConv {
     return std::find(coll->begin(), coll->end(), collection_name) != coll->end();
   }
 
+  std::tuple<std::string, std::string> getPidAlgoName(const std::string& collectionName)
+  {
+    const auto pidPos = collectionName.find("_PID_");
+    return {collectionName.substr(pidPos + 5), collectionName.substr(0, pidPos)};
+  }
+
   std::unique_ptr<lcio::LCEventImpl> convertEvent(const podio::Frame& edmEvent, const podio::Frame& metadata)
   {
     auto lcioEvent = std::make_unique<lcio::LCEventImpl>();
     auto objectMappings = CollectionsPairVectors {};
 
+    // We simply collect all pairs of PID algorithm names and the RecoParticle
+    // collection to which they belong (assuming the naming convention
+    // introduced in the LCIO to EDM4hep conversion)
+    std::vector<std::tuple<std::string, std::string>> pidAlgoNames;
+
     const auto& collections = edmEvent.getAvailableCollections();
     for (const auto& name : collections) {
       const auto edmCollection = edmEvent.get(name);
@@ -92,11 +103,12 @@ namespace EDM4hep2LCIOConv {
       else if (auto coll = dynamic_cast<const edm4hep::EventHeaderCollection*>(edmCollection)) {
         convertEventHeader(coll, lcioEvent.get());
       }
-      else if (
-        dynamic_cast<const edm4hep::CaloHitContributionCollection*>(edmCollection) ||
-        dynamic_cast<const edm4hep::ParticleIDCollection*>(edmCollection)) {
-        // "converted" as part of FillMissingCollectoins at the end or as part
-        // of the reconstructed particle
+      else if (auto coll = dynamic_cast<const edm4hep::ParticleIDCollection*>(edmCollection)) {
+        convertParticleIDs(coll, objectMappings.particleIDs, pidAlgoNames.size());
+        pidAlgoNames.emplace_back(getPidAlgoName(name));
+      }
+      else if (dynamic_cast<const edm4hep::CaloHitContributionCollection*>(edmCollection)) {
+        // "converted" during relation resolving later
         continue;
       }
       else {
@@ -110,6 +122,12 @@ namespace EDM4hep2LCIOConv {
       }
     }
 
+    for (const auto& [algoName, recoCollName] : pidAlgoNames) {
+      std::cout << "Adding PID algo " << algoName << " to reco collection " << recoCollName << '\n';
+      UTIL::PIDHandler pidHandler(lcioEvent->getCollection(recoCollName));
+      std::cout << "assigned algo id " << pidHandler.addAlgorithm(algoName, {}) << '\n';
+    }
+
     resolveRelations(objectMappings);
 
     return lcioEvent;

tests/src/CompareEDM4hepEDM4hep.cc

@@ -334,7 +334,17 @@ bool compare(const edm4hep::ParticleIDCollection& origColl, const edm4hep::Parti
     const auto origPid = origColl[i];
     const auto pid = roundtripColl[i];
 
-    REQUIRE_SAME(origPid.getAlgorithmType(), pid.getAlgorithmType(), "algorithm type in ParticleID " << i);
+    // This might not be preserved in roundtripping
+    // REQUIRE_SAME(origPid.getAlgorithmType(), pid.getAlgorithmType(), "algorithm type in ParticleID " << i);
+
+    REQUIRE_SAME(origPid.getType(), pid.getType(), "type in ParticleID " << i);
+
+    const auto origParams = origPid.getParameters();
+    const auto params = pid.getParameters();
+    REQUIRE_SAME(origParams.size(), params.size(), "parameter sizes in ParticleID " << i);
+    for (size_t iP = 0; iP < params.size(); ++iP) {
+      REQUIRE_SAME(origParams[iP], params[iP], "parameter " << iP << " in ParticleID " << i);
+    }
 
     REQUIRE_SAME(
       origPid.getParticle().getObjectID(), pid.getParticle().getObjectID(), "related particle in ParticleID " << i);

tests/src/EDM4hep2LCIOUtilities.cc

@@ -364,13 +364,16 @@ std::vector<edm4hep::ParticleIDCollection> createParticleIDs(
   collections.reserve(recoIdcs.size());
 
   int algoId = 0;
+  float param = 0;
   for (const auto& idcs : recoIdcs) {
     algoId++;
     auto& coll = collections.emplace_back();
     for (const auto idx : idcs) {
       auto pid = coll.create();
       pid.setAlgorithmType(algoId);
+      pid.setType(idx);
       pid.setParticle(recoParticles[idx]);
+      pid.addToParameters(param++);
     }
   }