added PrimaryVertices factory (subCollection of CentralTrackVertices)

src/algorithms/reco/PrimaryVertices.cc

@@ -0,0 +1,102 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2024 Daniel Brandenburg
+
+#include <Math/GenVector/LorentzVector.h>
+#include <Math/GenVector/PxPyPzM4D.h>
+#include <Math/Vector4Dfwd.h>
+#include <edm4eic/VertexCollection.h>
+#include <edm4eic/unit_system.h>
+#include <edm4hep/Vector3f.h>
+#include <fmt/core.h>
+#include <podio/ObjectID.h>
+#include <iterator>
+#include <map>
+#include <utility>
+
+#include "algorithms/reco/PrimaryVertices.h"
+#include "algorithms/reco/PrimaryVerticesConfig.h"
+
+namespace eicrecon {
+
+  /**
+   * @brief Initialize the PrimaryVertices Algorithm
+   *
+   * @param logger
+   */
+  void PrimaryVertices::init(std::shared_ptr<spdlog::logger>& logger) {
+    m_log = logger;
+  }
+
+  /**
+   * @brief Produce a list of primary vertex candidates
+   *
+   * @param rcvtx  - input collection of all vertex candidates
+   * @return std::unique_ptr<edm4eic::VertexCollection>
+   */
+  std::unique_ptr<edm4eic::VertexCollection> PrimaryVertices::execute(
+                const edm4eic::VertexCollection *rcvtx
+        ){
+
+    // this map will store intermediate results
+    // so that we can sort them before filling output
+    // collection
+    std::map<int, edm4eic::Vertex> primaryVertexMap;
+
+    // our output collection of primary vertex
+    // ordered by N_trk = associatedParticle array size
+    auto out_primary_vertices =  std::make_unique<edm4eic::VertexCollection>();
+    out_primary_vertices->setSubsetCollection();
+
+    m_log->trace( "We have {} candidate vertices",
+        rcvtx->size()
+      );
+
+    for ( const auto& vtx: *rcvtx ) {
+
+      const auto &v = vtx.getPosition();
+
+      // some basic vertex selection
+      if ( sqrt( v.x*v.x + v.y*v.y ) / edm4eic::unit::mm > m_cfg.maxVr ||
+           fabs( v.z ) / edm4eic::unit::mm > m_cfg.maxVz )
+           continue;
+
+      if ( vtx.getChi2() > m_cfg.maxChi2 ) continue;
+
+      //
+      int N_trk = vtx.getAssociatedParticles().size();
+      m_log->trace( "\t N_trk = {}", N_trk );
+      primaryVertexMap[ N_trk ] = vtx;
+    } // vertex loop
+
+    // map sorts in descending order by default
+    // sort by descending
+    bool first = true;
+    // for (auto kv : primaryVertexMap) {
+    for (auto kv = primaryVertexMap.rbegin(); kv != primaryVertexMap.rend(); ++kv) {
+
+      int N_trk = kv->first;
+      // Do not save primary candidates that
+      // are not within range
+      if ( N_trk > m_cfg.maxNtrk
+        || N_trk < m_cfg.minNtrk ){
+        continue;
+      }
+
+      // For logging and development
+      // report the highest N_trk candidate chosen
+      if ( first ){
+        m_log->trace( "Max N_trk Candidate:" );
+        m_log->trace( "\t N_trk = {}", N_trk );
+        m_log->trace( "\t Primary vertex has xyz=( {}, {}, {} )", kv->second.getPosition().x, kv->second.getPosition().y, kv->second.getPosition().z );
+        first = false;
+      }
+      out_primary_vertices->push_back( kv->second );
+    } // reverse loop on primaryVertexMap
+
+
+    // Return primary vertex ranked
+    // in order from largest N_trk to smallest
+    return out_primary_vertices;
+  }
+
+} // namespace eicrecon

src/algorithms/reco/PrimaryVertices.h

@@ -0,0 +1,29 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2024 Daniel Brandenburg
+
+#pragma once
+
+#include <edm4eic/VertexCollection.h>
+#include <spdlog/logger.h>
+#include <memory>
+
+#include "algorithms/reco/PrimaryVerticesConfig.h"
+#include "algorithms/interfaces/WithPodConfig.h"
+
+
+namespace eicrecon {
+
+    class PrimaryVertices : public WithPodConfig<PrimaryVerticesConfig>{
+
+    public:
+
+        void init(std::shared_ptr<spdlog::logger>& logger);
+        std::unique_ptr<edm4eic::VertexCollection> execute(
+                const edm4eic::VertexCollection *rcvtx
+        );
+
+    private:
+        std::shared_ptr<spdlog::logger> m_log;
+
+    };
+} // namespace eicrecon

src/algorithms/reco/PrimaryVerticesConfig.h

@@ -0,0 +1,23 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2024 Daniel Brandenburg
+
+#pragma once
+
+#include <string>
+#include <DD4hep/DD4hepUnits.h>
+
+namespace eicrecon {
+
+  struct PrimaryVerticesConfig {
+
+    // For now these are wide open
+    // In the future the cut should depend
+    // on the generator settings
+    float maxVr = 50.0; // mm
+    float maxVz = 500.0; // mm
+    float maxChi2 = 10000.0; //
+    int   minNtrk = 1;  // >=
+    int   maxNtrk = 1000000; // <=
+  };
+
+}  // end eicrecon namespace

src/global/reco/PrimaryVertices_factory.h

@@ -0,0 +1,47 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2024 Daniel Brandenburg
+
+#pragma once
+
+#include <JANA/JEvent.h>
+#include <memory>
+#include <string>
+#include <utility>
+#include <vector>
+
+#include "algorithms/reco/PrimaryVertices.h"
+#include "algorithms/reco/PrimaryVerticesConfig.h"
+#include "extensions/jana/JOmniFactory.h"
+
+namespace eicrecon {
+
+class PrimaryVertices_factory :
+        public JOmniFactory<PrimaryVertices_factory, PrimaryVerticesConfig> {
+
+public:
+    using AlgoT = eicrecon::PrimaryVertices;
+private:
+    std::unique_ptr<AlgoT> m_algo;
+
+    PodioInput<edm4eic::Vertex> m_rc_vertices_input {this};
+
+    // Declare outputs
+    PodioOutput<edm4eic::Vertex> m_out_primary_vertices {this};
+
+public:
+    void Configure() {
+        m_algo = std::make_unique<AlgoT>();
+        m_algo->init(logger());
+        m_algo->applyConfig(config());
+    }
+
+    void ChangeRun(int64_t run_number) {
+    }
+
+    void Process(int64_t run_number, uint64_t event_number) {
+        auto output = m_algo->execute(m_rc_vertices_input());
+        m_out_primary_vertices() = std::move(output);
+    }
+};
+
+} // eicrecon

src/global/reco/reco.cc

@@ -47,6 +47,7 @@
 #include "global/reco/ReconstructedElectrons_factory.h"
 #include "global/reco/ScatteredElectronsEMinusPz_factory.h"
 #include "global/reco/ScatteredElectronsTruth_factory.h"
+#include "global/reco/PrimaryVertices_factory.h"
 
 extern "C" {
 void InitPlugin(JApplication *app) {
@@ -384,6 +385,19 @@ void InitPlugin(JApplication *app) {
             app
     ));
 
+    app->Add(new JOmniFactoryGeneratorT<PrimaryVertices_factory>(
+        "PrimaryVertices",
+        {
+          "CentralTrackVertices"
+        },
+        {
+          "PrimaryVertices"
+        },
+        {
+        },
+        app
+    ));
+
 
 }
 } // extern "C"

src/services/io/podio/JEventProcessorPODIO.cc

@@ -221,6 +221,7 @@ JEventProcessorPODIO::JEventProcessorPODIO() {
             "ReconstructedElectrons",
             "ScatteredElectronsTruth",
             "ScatteredElectronsEMinusPz",
+            "PrimaryVertices",
 #if EDM4EIC_VERSION_MAJOR >= 6
             "HadronicFinalState",
 #endif