Use JANA's new first-class PODIO integration (#578)

### Briefly, what does this PR introduce? This PR allows eicrecon to use PODIO associations correctly, resolving a major source of confusion, entropy, and bugs. As of v2.1.0, the JANA2 framework has first-class support for PODIO. This means that users can directly read and write PODIO collections. Even if users continue to use the classic JANA vector-of-pointers style, under the hood the PODIO objects are always registered to a collection and that collection is registered to a frame. JANA understands and abides by PODIO's memory ownership semantics, including supporting subset collections. This PR also converts three problematic factories into multifactories, which support multiple outputs. These factories used direct object ID rewriting in order to set up valid associations, and/or had particularly thorny PODIO ownership issues. They are now models for how to write factories in the future. These factories are: MatchClusters_factory, ParticlesWithTruthPID_factory, and TrackingResult_factory. ### What kind of change does this PR introduce? - [ ] Bug fix (issue #__) - [ ] New feature (issue #__) - [ ] Documentation update - [ ] Other: __ ### Please check if this PR fulfills the following: - [ ] Tests for the changes have been added - [ ] Documentation has been added / updated - [x] Changes have been communicated to collaborators ### Does this PR introduce breaking changes? What changes might users need to make to their code? - The way that users should write future factories has changed. Users should look to TrackingResult_factory as an example of how to write factories with multiple outputs. - Output files are now written using Podio's Frame writer, which is incompatible with Podio's older EventStore writer. ### Does this PR change default behavior? - JEventProcessorPODIO used to maintain a list of factories that have thrown an exception. After it crashed once, it would never be called again. Now, it will be called again on every event. The reason for this change in behavior is that PODIO v0.16.3 has a bug where collection IDs depend on insertion order, which causes the PODIO frame writer to crash. Excepting factories still produce an empty collection (thus preserving the insertion order), which the JEventProcessorPODIO then persists. --------- Co-authored-by: Nathan Brei <[email protected]>
eic · Apr 14, 2023 · e14e81b · e14e81b
1 parent 11e728d
commit e14e81b
Show file tree

Hide file tree

Showing 151 changed files with 1,449 additions and 3,419 deletions.
diff --git a/.gitignore b/.gitignore
@@ -2,6 +2,7 @@
 
 # Autogenerated file
 src/services/io/podio/datamodel_glue.h
+src/services/io/podio/datamodel_includes.h
 
 # Indirectory builds
 include/*

diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -87,6 +87,9 @@ find_package(JANA REQUIRED)
 message(STATUS "${CMAKE_PROJECT_NAME}: JANA2 CMake   : ${JANA_DIR}")
 message(STATUS "${CMAKE_PROJECT_NAME}: JANA2 includes: ${JANA_INCLUDE_DIR}")
 message(STATUS "${CMAKE_PROJECT_NAME}: JANA2 library : ${JANA_LIBRARY}")
+add_compile_definitions(HAVE_PODIO)
+# TODO: NWB: Do this correctly in JANA via target_compile_definitions for v2.1.1
+# TODO: NWB: Maybe we want these to be prefixed, e.g. JANA2_HAVE_PODIO
 
 # PODIO, EDM4HEP, EDM4EIC event models
 find_package(podio REQUIRED)

diff --git a/src/algorithms/reco/MatchClusters.cc b/src/algorithms/reco/MatchClusters.cc
@@ -25,15 +25,15 @@
 #include "edm4eic/TrackParametersCollection.h"
 #include "edm4eic/vector_utils.h"
 
-#include "ParticlesWithAssociation.h"
+
 
 namespace eicrecon {
 
     void MatchClusters::init(std::shared_ptr<spdlog::logger> logger) {
         m_log = logger;
     }
 
-    ParticlesWithAssociation *MatchClusters::execute(
+    std::tuple<edm4eic::ReconstructedParticleCollection*, edm4eic::MCRecoParticleAssociationCollection*> MatchClusters::execute(
             std::vector<const edm4hep::MCParticle *> mcparticles,
             std::vector<edm4eic::ReconstructedParticle *> inparts,
             std::vector<edm4eic::MCRecoParticleAssociation *> inpartsassoc,
@@ -42,17 +42,12 @@ namespace eicrecon {
 
         m_log->debug("Processing cluster info for new event");
 
-
-        // Resulting reconstructed particles
-        std::vector<edm4eic::ReconstructedParticle *> outparts;
-
-        // Resulting associations
-        std::vector<edm4eic::MCRecoParticleAssociation *> outpartsassoc;
-
+        // Resulting reconstructed particles and associations
+        edm4eic::ReconstructedParticleCollection* outparts = new edm4eic::ReconstructedParticleCollection();
+        edm4eic::MCRecoParticleAssociationCollection* outpartsassoc = new edm4eic::MCRecoParticleAssociationCollection();
 
         m_log->debug("Step 0/2: Getting indexed list of clusters...");
 
-
         // get an indexed map of all clusters
         auto clusterMap = indexedClusters(cluster_collections, cluster_assoc_collections);
 
@@ -64,8 +59,7 @@ namespace eicrecon {
             m_log->debug(" --> Processing charged particle {}, PDG {}, energy {}", inpart->getObjectID().index,
                          inpart->getPDG(), inpart->getEnergy());
 
-            auto outpart = new edm4eic::ReconstructedParticle(*inpart);
-            outparts.push_back(outpart);
+            auto outpart = outparts->create();
 
             int mcID = -1;
 
@@ -91,13 +85,12 @@ namespace eicrecon {
             }
 
             // create truth associations
-            auto assoc = edm4eic::MutableMCRecoParticleAssociation();
-            assoc.setRecID(outpart->getObjectID().index);
+            auto assoc = outpartsassoc->create();
+            assoc.setRecID(outpart.getObjectID().index);
             assoc.setSimID(mcID);
             assoc.setWeight(1.0);
-            assoc.setRec(*outpart);
-            //assoc.setSim(mcparticles[mcID]);
-            outpartsassoc.push_back(new edm4eic::MCRecoParticleAssociation(assoc));
+            assoc.setRec(outpart);
+            assoc.setRec(outpart);
         }
 
         // 2. Now loop over all remaining clusters and add neutrals. Also add in Hcal energy
@@ -125,7 +118,7 @@ namespace eicrecon {
             m_log->debug(" --> Reconstructed neutral particle with PDG: {}, energy: {}", outpart.getPDG(),
                          outpart.getEnergy());
 
-            outparts.push_back(new edm4eic::ReconstructedParticle(outpart));
+            outparts->push_back(outpart);
 
             // Create truth associations
             auto assoc = edm4eic::MutableMCRecoParticleAssociation();
@@ -135,7 +128,8 @@ namespace eicrecon {
             assoc.setRec(outpart);
             //assoc.setSim(mcparticles[mcID]);
         }
-        return new ParticlesWithAssociation(std::move(outparts), std::move(outpartsassoc));
+
+        return {outparts, outpartsassoc};
     }
 
 

diff --git a/src/algorithms/reco/MatchClusters.h b/src/algorithms/reco/MatchClusters.h
@@ -24,7 +24,6 @@
 #include <edm4eic/TrackParametersCollection.h>
 #include <edm4eic/vector_utils.h>
 
-#include "ParticlesWithAssociation.h"
 
 namespace eicrecon {
 
@@ -33,27 +32,27 @@ namespace eicrecon {
 
     public:
 
+        using MatchingResults = std::tuple<edm4eic::ReconstructedParticleCollection*, edm4eic::MCRecoParticleAssociationCollection*>;
+
         void init(std::shared_ptr<spdlog::logger> logger);
 
-        ParticlesWithAssociation *execute(
-                std::vector<const edm4hep::MCParticle *> mcparticles,
-                std::vector<edm4eic::ReconstructedParticle *> inparts,            // TODO fix const
-                std::vector<edm4eic::MCRecoParticleAssociation *> inpartsassoc,   // TODO fix const
-                const std::vector<std::vector<const edm4eic::Cluster*>> &cluster_collections,
-                const std::vector<std::vector<const edm4eic::MCRecoClusterParticleAssociation*>> &cluster_assoc_collections);
+        MatchingResults execute(
+            std::vector<const edm4hep::MCParticle *> mcparticles,
+            std::vector<edm4eic::ReconstructedParticle *> inparts,            // TODO fix const
+            std::vector<edm4eic::MCRecoParticleAssociation *> inpartsassoc,   // TODO fix const
+            const std::vector<std::vector<const edm4eic::Cluster*>> &cluster_collections,
+            const std::vector<std::vector<const edm4eic::MCRecoClusterParticleAssociation*>> &cluster_assoc_collections);
 
     private:
 
         std::shared_ptr<spdlog::logger> m_log;
 
-
         // get a map of mcID --> cluster
         // input: cluster_collections --> list of handles to all cluster collections
         std::map<int, const edm4eic::Cluster*> indexedClusters(
                 const std::vector<std::vector<const edm4eic::Cluster*>> &cluster_collections,
                 const std::vector<std::vector<const edm4eic::MCRecoClusterParticleAssociation*>> &associations_collections);
 
-
         // reconstruct a neutral cluster
         // (for now assuming the vertex is at (0,0,0))
         edm4eic::ReconstructedParticle

diff --git a/src/algorithms/tracking/ParticlesFromTrackFit.cc b/src/algorithms/tracking/ParticlesFromTrackFit.cc
@@ -29,7 +29,7 @@ void eicrecon::Reco::ParticlesFromTrackFit::init(std::shared_ptr<spdlog::logger>
     m_log = log;
 }
 
-ParticlesFromTrackFitResult* eicrecon::Reco::ParticlesFromTrackFit::execute(const std::vector<const eicrecon::TrackingResultTrajectory *> &trajectories) {
+ParticlesFromTrackFitResultNew eicrecon::Reco::ParticlesFromTrackFit::execute(const std::vector<const eicrecon::TrackingResultTrajectory *> &trajectories) {
 
     // create output collections
     auto rec_parts = std::make_unique<edm4eic::ReconstructedParticleCollection >();
@@ -137,5 +137,5 @@ ParticlesFromTrackFitResult* eicrecon::Reco::ParticlesFromTrackFit::execute(cons
         });
     }
 
-    return new ParticlesFromTrackFitResult(std::move(rec_parts), std::move(track_pars));
+    return std::make_pair(std::move(rec_parts), std::move(track_pars));
 }
diff --git a/src/algorithms/tracking/ParticlesFromTrackFit.h b/src/algorithms/tracking/ParticlesFromTrackFit.h
@@ -10,6 +10,9 @@
 #include <algorithms/tracking/ParticlesFromTrackFitResult.h>
 
 
+using ParticlesFromTrackFitResultNew = std::pair<std::unique_ptr<edm4eic::ReconstructedParticleCollection>,
+                                                 std::unique_ptr<edm4eic::TrackParametersCollection>>;
+
 namespace eicrecon::Reco {
 
     /** Extract the particles form fit trajectories.
@@ -24,7 +27,7 @@ namespace eicrecon::Reco {
     public:
         void init(std::shared_ptr<spdlog::logger> log);
 
-        ParticlesFromTrackFitResult *execute(const std::vector<const eicrecon::TrackingResultTrajectory *> &trajectories);
+        ParticlesFromTrackFitResultNew execute(const std::vector<const eicrecon::TrackingResultTrajectory *> &trajectories);
 
     };
 } // namespace Jug::Reco
diff --git a/src/algorithms/tracking/ParticlesWithTruthPID.cpp b/src/algorithms/tracking/ParticlesWithTruthPID.cpp
@@ -3,19 +3,8 @@
 
 #include "ParticlesWithTruthPID.h"
 
-// Event Model related classes
-#include <edm4hep/MutableMCParticle.h>
-#include <edm4eic/ReconstructedParticle.h>
-#include <edm4eic/MutableReconstructedParticle.h>
-#include <edm4eic/MutableMCRecoParticleAssociation.h>
 #include <edm4eic/vector_utils.h>
 
-// Getting access to PODIO indexes
-#include <extensions/podio_access/accessor.h>
-#include <edm4eic/ReconstructedParticleObj.h>
-
-using ReconstructedParticleObjPtr = edm4eic::ReconstructedParticleObj*;
-ALLOW_ACCESS(edm4eic::MutableReconstructedParticle, m_obj, ReconstructedParticleObjPtr);
 
 
 namespace eicrecon {
@@ -25,41 +14,40 @@ namespace eicrecon {
 
     }
 
-    ParticlesWithAssociation *ParticlesWithTruthPID::process(
-            const std::vector<const edm4hep::MCParticle *> &mc_particles,
-            const std::vector<const edm4eic::TrackParameters *> &track_params) {
+    ParticlesWithAssociationNew ParticlesWithTruthPID::process(
+            const edm4hep::MCParticleCollection* mc_particles,
+            const edm4eic::TrackParametersCollection* track_params) {
+
         // input collection
 
         /// Resulting reconstructed particles
-        std::vector<edm4eic::ReconstructedParticle *> reco_particles;
-
-        /// Resulting associations
-        std::vector<edm4eic::MCRecoParticleAssociation *> associations;
+        auto reco_particles = std::make_unique<edm4eic::ReconstructedParticleCollection>();
+        auto associations = std::make_unique<edm4eic::MCRecoParticleAssociationCollection>();
 
         const double sinPhiOver2Tolerance = sin(0.5 * m_cfg.phiTolerance);
         tracePhiToleranceOnce(sinPhiOver2Tolerance, m_cfg.phiTolerance);
 
-        std::vector<bool> mc_prt_is_consumed(mc_particles.size(), false);         // MCParticle is already consumed flag
+        std::vector<bool> mc_prt_is_consumed(mc_particles->size(), false);         // MCParticle is already consumed flag
 
-        for (const auto &trk: track_params) {
-            const auto mom = edm4eic::sphericalToVector(1.0 / std::abs(trk->getQOverP()), trk->getTheta(),
-                                                        trk->getPhi());
-            const auto charge_rec = trk->getCharge();
+        for (const auto &trk: *track_params) {
+            const auto mom = edm4eic::sphericalToVector(1.0 / std::abs(trk.getQOverP()), trk.getTheta(),
+                                                        trk.getPhi());
+            const auto charge_rec = trk.getCharge();
 
 
             m_log->debug("Match:  [id]   [mom]   [theta]  [phi]    [charge]  [PID]");
             m_log->debug(" Track : {:<4} {:<8.3f} {:<8.3f} {:<8.2f} {:<4}",
-                         trk->getObjectID().index, edm4eic::magnitude(mom), edm4eic::anglePolar(mom), edm4eic::angleAzimuthal(mom), charge_rec);
+                         trk.getObjectID().index, edm4eic::magnitude(mom), edm4eic::anglePolar(mom), edm4eic::angleAzimuthal(mom), charge_rec);
 
             // utility variables for matching
             int best_match = -1;
             double best_delta = std::numeric_limits<double>::max();
-            for (size_t ip = 0; ip < mc_particles.size(); ++ip) {
-                const auto &mc_part = mc_particles[ip];
-                const auto &p = mc_part->getMomentum();
+            for (size_t ip = 0; ip < mc_particles->size(); ++ip) {
+                const auto &mc_part = (*mc_particles)[ip];
+                const auto &p = mc_part.getMomentum();
 
-                m_log->trace("  MCParticle with id={:<4} mom={:<8.3f} charge={}", mc_part->getObjectID().index,
-                             edm4eic::magnitude(p), mc_part->getCharge());
+                m_log->trace("  MCParticle with id={:<4} mom={:<8.3f} charge={}", mc_part.getObjectID().index,
+                             edm4eic::magnitude(p), mc_part.getCharge());
 
                 // Check if used
                 if (mc_prt_is_consumed[ip]) {
@@ -68,19 +56,19 @@ namespace eicrecon {
                 }
 
                 // Check if non-primary
-                if (mc_part->getGeneratorStatus() > 1) {
+                if (mc_part.getGeneratorStatus() > 1) {
                     m_log->trace("    Ignoring. GeneratorStatus > 1 => Non-primary particle");
                     continue;
                 }
 
                 // Check if neutral
-                if (mc_part->getCharge() == 0) {
+                if (mc_part.getCharge() == 0) {
                     m_log->trace("    Ignoring. Neutral particle");
                     continue;
                 }
 
                 // Check opposite charge
-                if (mc_part->getCharge() * charge_rec < 0) {
+                if (mc_part.getCharge() * charge_rec < 0) {
                     m_log->trace("    Ignoring. Opposite charge particle");
                     continue;
                 }
@@ -123,22 +111,16 @@ namespace eicrecon {
             if (best_match >= 0) {
                 m_log->trace("Best match is found and is: {}", best_match);
                 mc_prt_is_consumed[best_match] = true;
-                const auto &best_mc_part = mc_particles[best_match];
-                best_pid = best_mc_part->getPDG();
+                const auto &best_mc_part = (*mc_particles)[best_match];
+                best_pid = best_mc_part.getPDG();
                 referencePoint = {
-                        static_cast<float>(best_mc_part->getVertex().x),
-                        static_cast<float>(best_mc_part->getVertex().y),
-                        static_cast<float>(best_mc_part->getVertex().z)}; // @TODO: not sure if vertex/reference point makes sense here
+                        static_cast<float>(best_mc_part.getVertex().x),
+                        static_cast<float>(best_mc_part.getVertex().y),
+                        static_cast<float>(best_mc_part.getVertex().z)}; // @TODO: not sure if vertex/reference point makes sense here
                 // time                 = mcpart.getTime();
-                mass = best_mc_part->getMass();
+                mass = best_mc_part.getMass();
             }
 
-            // Getting access to private PODIO member to set the index in the collection
-            // Evil, Evil, Evil, Evil is here TODO change asap
-            // Nothing is private in this world anymore!
-            auto obj = ACCESS(rec_part, m_obj);
-            obj->id.index = (int) reco_particles.size();     // Set the index in the collection
-
             rec_part.setType(static_cast<int16_t>(best_match >= 0 ? 0 : -1)); // @TODO: determine type codes
             rec_part.setEnergy((float) std::hypot(edm4eic::magnitude(mom), mass));
             rec_part.setMomentum(mom);
@@ -148,27 +130,32 @@ namespace eicrecon {
             rec_part.setGoodnessOfPID(1); // perfect PID
             rec_part.setPDG(best_pid);
             // rec_part.covMatrix()  // @TODO: covariance matrix on 4-momentum
+            // Add reconstructed particle to collection BEFORE doing association
+            reco_particles->push_back(rec_part);
+
             // Also write MC <--> truth particle association if match was found
             if (best_match >= 0) {
                 auto rec_assoc = edm4eic::MutableMCRecoParticleAssociation();
                 rec_assoc.setRecID(rec_part.getObjectID().index);
-                rec_assoc.setSimID(mc_particles[best_match]->getObjectID().index);
+                rec_assoc.setSimID((*mc_particles)[best_match].getObjectID().index);
                 rec_assoc.setWeight(1);
                 rec_assoc.setRec(rec_part);
-                auto sim = mc_particles[best_match]->clone();
+                auto sim = (*mc_particles)[best_match];
                 rec_assoc.setSim(sim);
-                associations.emplace_back(new edm4eic::MCRecoParticleAssociation(rec_assoc));
+
+                // Add association to collection
+                associations->push_back(rec_assoc);
 
                 if (m_log->level() <= spdlog::level::debug) {
 
-                    const auto &mcpart = mc_particles[best_match];
-                    const auto &p = mcpart->getMomentum();
+                    const auto &mcpart = (*mc_particles)[best_match];
+                    const auto &p = mcpart.getMomentum();
                     const auto p_mag = edm4eic::magnitude(p);
                     const auto p_phi = edm4eic::angleAzimuthal(p);
                     const auto p_theta = edm4eic::anglePolar(p);
                     m_log->debug(" MCPart: {:<4} {:<8.3f} {:<8.3f} {:<8.2f} {:<6}",
-                                 mcpart->getObjectID().index, p_mag, p_theta, p_phi, mcpart->getCharge(),
-                                 mcpart->getPDG());
+                                 mcpart.getObjectID().index, p_mag, p_theta, p_phi, mcpart.getCharge(),
+                                 mcpart.getPDG());
 
                     m_log->debug(" Assoc: id={} SimId={} RecId={}",
                                  rec_assoc.getObjectID().index, rec_assoc.getSimID(), rec_assoc.getSimID());
@@ -178,12 +165,10 @@ namespace eicrecon {
                 m_log->debug(" MCPart: Did not find a good match");
             }
 
-            // Add particle to the output vector
-            reco_particles.push_back(new edm4eic::ReconstructedParticle(rec_part));
         }
 
         // Assembling the results
-        return new ParticlesWithAssociation(std::move(reco_particles), std::move(associations));
+        return std::make_pair(std::move(reco_particles), std::move(associations));
     }
 
     void ParticlesWithTruthPID::tracePhiToleranceOnce(const double sinPhiOver2Tolerance, double phiTolerance) {