Merge branch 'main' into fix_CellID_FarClustering

eic · Jun 4, 2024 · 997ac48 · 997ac48
2 parents a626503 + 1413111
commit 997ac48
Show file tree

Hide file tree

Showing 66 changed files with 2,286 additions and 2,003 deletions.
diff --git a/.github/iwyu.imp b/.github/iwyu.imp
@@ -30,6 +30,7 @@
   { include: ['<JANA/JApplicationFwd.h>', private, '<JANA/JApplication.h>', public] },
 
   # the rest
+  { include: ['<bits/chrono.h>', private, '<chrono>', public] },
   { include: ['<bits/std_abs.h>', private, '<math.h>', public] },
   { include: ['<bits/utility.h>', private, '<tuple>', public] },
   # https://github.com/include-what-you-use/include-what-you-use/issues/166

diff --git a/.github/workflows/linux-eic-shell.yml b/.github/workflows/linux-eic-shell.yml
@@ -931,7 +931,7 @@ jobs:
           GITHUB_SHA=${{ github.event.pull_request.head.sha || github.sha }}
           GITHUB_PR=${{ github.event.pull_request.number }}
           EICRECON_VERSION="${{ github.event.pull_request.head.ref || github.ref_name }}"
-          PIPELINE_NAME_CONTAINER=${{ github.repository }}: ${{ github.event.pull_request.title || github.ref_name }}
+          PIPELINE_NAME=${{ github.repository }}: ${{ github.event.pull_request.title || github.ref_name }}
     - name: Post a GitHub CI status
       run: |
         gh api \

diff --git a/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.cc b/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.cc
@@ -10,6 +10,7 @@
 
 #include <Evaluator/DD4hepUnits.h>
 #include <boost/algorithm/string/join.hpp>
+#include <boost/iterator/iterator_facade.hpp>
 #include <boost/range/adaptor/map.hpp>
 #include <edm4eic/CalorimeterHitCollection.h>
 #include <edm4hep/CaloHitContributionCollection.h>
@@ -26,6 +27,7 @@
 #include <complex>
 #include <cstddef>
 #include <gsl/pointers>
+#include <iterator>
 #include <limits>
 #include <map>
 #include <optional>
@@ -246,12 +248,13 @@ std::optional<edm4eic::MutableCluster> CalorimeterClusterRecoCoG::reconstruct(co
   Eigen::Vector3f sum1_3D = Eigen::Vector3f::Zero();
   Eigen::Vector2cf eigenValues_2D = Eigen::Vector2cf::Zero();
   Eigen::Vector3cf eigenValues_3D = Eigen::Vector3cf::Zero();
-
+  //the axis is the direction of the eigenvalue corresponding to the largest eigenvalue.
+  double axis_x=0, axis_y=0, axis_z=0;
   if (cl.getNhits() > 1) {
 
     for (const auto& hit : pcl.getHits()) {
 
-      float w = weightFunc(hit.getEnergy(), cl.getEnergy(), m_cfg.logWeightBase, 0);
+      float w = weightFunc(hit.getEnergy(), totalE, logWeightBase, 0);
 
       // theta, phi
       Eigen::Vector2f pos2D( edm4hep::utils::anglePolar( hit.getPosition() ), edm4hep::utils::angleAzimuthal( hit.getPosition() ) );
@@ -273,8 +276,8 @@ std::optional<edm4eic::MutableCluster> CalorimeterClusterRecoCoG::reconstruct(co
       w_sum += w;
     }
 
+    radius     = sqrt((1. / (cl.getNhits() - 1.)) * radius);
     if( w_sum > 0 ) {
-      radius     = sqrt((1. / (cl.getNhits() - 1.)) * radius);
       dispersion = sqrt( dispersion / w_sum );
 
       // normalize matrices
@@ -289,11 +292,33 @@ std::optional<edm4eic::MutableCluster> CalorimeterClusterRecoCoG::reconstruct(co
 
       // Solve for eigenvalues.  Corresponds to cluster's 2nd moments (widths)
       Eigen::EigenSolver<Eigen::Matrix2f> es_2D(cov2, false); // set to true for eigenvector calculation
-      Eigen::EigenSolver<Eigen::Matrix3f> es_3D(cov3, false); // set to true for eigenvector calculation
+      Eigen::EigenSolver<Eigen::Matrix3f> es_3D(cov3, true); // set to true for eigenvector calculation
 
       // eigenvalues of symmetric real matrix are always real
       eigenValues_2D = es_2D.eigenvalues();
       eigenValues_3D = es_3D.eigenvalues();
+      //find the eigenvector corresponding to the largest eigenvalue
+      auto eigenvectors= es_3D.eigenvectors();
+      auto max_eigenvalue_it = std::max_element(
+        eigenValues_3D.begin(),
+        eigenValues_3D.end(),
+        [](auto a, auto b) {
+            return std::real(a) < std::real(b);
+        }
+      );
+      auto axis = eigenvectors.col(std::distance(
+            eigenValues_3D.begin(),
+            max_eigenvalue_it
+        ));
+      axis_x=axis(0,0).real();
+      axis_y=axis(1,0).real();
+      axis_z=axis(2,0).real();
+      double norm=sqrt(axis_x*axis_x+axis_y*axis_y+axis_z*axis_z);
+      if (norm!=0){
+        axis_x/=norm;
+        axis_y/=norm;
+        axis_z/=norm;
+      }
     }
   }
 
@@ -304,7 +329,10 @@ std::optional<edm4eic::MutableCluster> CalorimeterClusterRecoCoG::reconstruct(co
   cl.addToShapeParameters( eigenValues_3D[0].real() ); // 3D x-y-z cluster width 1
   cl.addToShapeParameters( eigenValues_3D[1].real() ); // 3D x-y-z cluster width 2
   cl.addToShapeParameters( eigenValues_3D[2].real() ); // 3D x-y-z cluster width 3
-
+  //last 3 shape parameters are the components of the axis direction
+  cl.addToShapeParameters( axis_x );
+  cl.addToShapeParameters( axis_y );
+  cl.addToShapeParameters( axis_z );
   return std::move(cl);
 }
 

diff --git a/src/algorithms/calorimetry/ImagingClusterReco.h b/src/algorithms/calorimetry/ImagingClusterReco.h
@@ -23,6 +23,7 @@
 // Event Model related classes
 #include <edm4hep/MCParticleCollection.h>
 #include <edm4hep/SimCalorimeterHitCollection.h>
+#include <edm4hep/utils/vector_utils.h>
 #include <edm4eic/CalorimeterHitCollection.h>
 #include <edm4eic/ClusterCollection.h>
 #include <edm4eic/MCRecoClusterParticleAssociationCollection.h>

diff --git a/src/algorithms/digi/SiliconTrackerDigi.cc b/src/algorithms/digi/SiliconTrackerDigi.cc
@@ -20,10 +20,7 @@
 
 namespace eicrecon {
 
-void SiliconTrackerDigi::init(std::shared_ptr<spdlog::logger>& logger) {
-    // set logger
-    m_log = logger;
-
+void SiliconTrackerDigi::init() {
     // Create random gauss function
     m_gauss = [&](){
         return m_random.Gaus(0, m_cfg.timeResolution);
@@ -50,20 +47,20 @@ void SiliconTrackerDigi::process(
         double result_time = sim_hit.getTime() + time_smearing;
         auto hit_time_stamp = (std::int32_t) (result_time * 1e3);
 
-        m_log->debug("--------------------");
-        m_log->debug("Hit cellID   = {}", sim_hit.getCellID());
-        m_log->debug("   position  = ({:.2f}, {:.2f}, {:.2f})", sim_hit.getPosition().x, sim_hit.getPosition().y, sim_hit.getPosition().z);
-        m_log->debug("   xy_radius = {:.2f}", std::hypot(sim_hit.getPosition().x, sim_hit.getPosition().y));
-        m_log->debug("   momentum  = ({:.2f}, {:.2f}, {:.2f})", sim_hit.getMomentum().x, sim_hit.getMomentum().y, sim_hit.getMomentum().z);
-        m_log->debug("   edep = {:.2f}", sim_hit.getEDep());
-        m_log->debug("   time = {:.4f}[ns]", sim_hit.getTime());
-        m_log->debug("   particle time = {}[ns]", sim_hit.getMCParticle().getTime());
-        m_log->debug("   time smearing: {:.4f}, resulting time = {:.4f} [ns]", time_smearing, result_time);
-        m_log->debug("   hit_time_stamp: {} [~ps]", hit_time_stamp);
+        debug("--------------------");
+        debug("Hit cellID   = {}", sim_hit.getCellID());
+        debug("   position  = ({:.2f}, {:.2f}, {:.2f})", sim_hit.getPosition().x, sim_hit.getPosition().y, sim_hit.getPosition().z);
+        debug("   xy_radius = {:.2f}", std::hypot(sim_hit.getPosition().x, sim_hit.getPosition().y));
+        debug("   momentum  = ({:.2f}, {:.2f}, {:.2f})", sim_hit.getMomentum().x, sim_hit.getMomentum().y, sim_hit.getMomentum().z);
+        debug("   edep = {:.2f}", sim_hit.getEDep());
+        debug("   time = {:.4f}[ns]", sim_hit.getTime());
+        debug("   particle time = {}[ns]", sim_hit.getMCParticle().getTime());
+        debug("   time smearing: {:.4f}, resulting time = {:.4f} [ns]", time_smearing, result_time);
+        debug("   hit_time_stamp: {} [~ps]", hit_time_stamp);
 
 
         if (sim_hit.getEDep() < m_cfg.threshold) {
-            m_log->debug("  edep is below threshold of {:.2f} [keV]", m_cfg.threshold / dd4hep::keV);
+            debug("  edep is below threshold of {:.2f} [keV]", m_cfg.threshold / dd4hep::keV);
             continue;
         }
 
@@ -77,7 +74,7 @@ void SiliconTrackerDigi::process(
         } else {
             // There is previous values in the cell
             auto& hit = cell_hit_map[sim_hit.getCellID()];
-            m_log->debug("  Hit already exists in cell ID={}, prev. hit time: {}", sim_hit.getCellID(), hit.getTimeStamp());
+            debug("  Hit already exists in cell ID={}, prev. hit time: {}", sim_hit.getCellID(), hit.getTimeStamp());
 
             // keep earliest time for hit
             auto time_stamp = hit.getTimeStamp();

diff --git a/src/algorithms/digi/SiliconTrackerDigi.h b/src/algorithms/digi/SiliconTrackerDigi.h
@@ -5,12 +5,10 @@
 
 #include <TRandomGen.h>
 #include <algorithms/algorithm.h>
-#include <edm4eic/RawTrackerHitCollection.h>
 #include <edm4eic/MCRecoTrackerHitAssociationCollection.h>
+#include <edm4eic/RawTrackerHitCollection.h>
 #include <edm4hep/SimTrackerHitCollection.h>
-#include <spdlog/logger.h>
 #include <functional>
-#include <memory>
 #include <string>
 #include <string_view>
 
@@ -19,45 +17,35 @@
 
 namespace eicrecon {
 
-  using SiliconTrackerDigiAlgorithm = algorithms::Algorithm<
-    algorithms::Input<
-      edm4hep::SimTrackerHitCollection
-    >,
-    algorithms::Output<
-      edm4eic::RawTrackerHitCollection,
-      edm4eic::MCRecoTrackerHitAssociationCollection
-    >
-  >;
-
-  class SiliconTrackerDigi
-  : public SiliconTrackerDigiAlgorithm,
-    public WithPodConfig<SiliconTrackerDigiConfig> {
-
-  public:
-    SiliconTrackerDigi(std::string_view name)
-      : SiliconTrackerDigiAlgorithm{name,
-                            {"inputHitCollection"},
-                            {"outputRawHitCollection","outputHitAssociations"},
-                            "Apply threshold, digitize within ADC range, "
-                            "convert time with smearing resolution."} {}
+using SiliconTrackerDigiAlgorithm =
+    algorithms::Algorithm<algorithms::Input<edm4hep::SimTrackerHitCollection>,
+                          algorithms::Output<edm4eic::RawTrackerHitCollection,
+                                             edm4eic::MCRecoTrackerHitAssociationCollection>>;
 
-    void init(std::shared_ptr<spdlog::logger>& logger);
-    void process(const Input&, const Output&) const final;
+class SiliconTrackerDigi : public SiliconTrackerDigiAlgorithm,
+                           public WithPodConfig<SiliconTrackerDigiConfig> {
 
-  private:
-    /** algorithm logger */
-    std::shared_ptr<spdlog::logger> m_log;
+public:
+  SiliconTrackerDigi(std::string_view name)
+      : SiliconTrackerDigiAlgorithm{name,
+                                    {"inputHitCollection"},
+                                    {"outputRawHitCollection", "outputHitAssociations"},
+                                    "Apply threshold, digitize within ADC range, "
+                                    "convert time with smearing resolution."} {}
 
-    /** Random number generation*/
-    TRandomMixMax m_random;
-    std::function<double()> m_gauss;
+  void init() final;
+  void process(const Input&, const Output&) const final;
 
-    // FIXME replace with standard random engine
-    //std::default_random_engine generator; // TODO: need something more appropriate here
-    //std::normal_distribution<double> m_normDist; // defaults to mean=0, sigma=1
+private:
+  /** Random number generation*/
+  TRandomMixMax m_random;
+  std::function<double()> m_gauss;
 
-    //algorithms::Generator m_rng = algorithms::RandomSvc::instance().generator();
+  // FIXME replace with standard random engine
+  // std::default_random_engine generator; // TODO: need something more appropriate here
+  // std::normal_distribution<double> m_normDist; // defaults to mean=0, sigma=1
 
-  };
+  // algorithms::Generator m_rng = algorithms::RandomSvc::instance().generator();
+};
 
-} // eicrecon
+} // namespace eicrecon
diff --git a/src/algorithms/fardetectors/CMakeLists.txt b/src/algorithms/fardetectors/CMakeLists.txt
@@ -11,4 +11,8 @@ plugin_glob_all(${PLUGIN_NAME})
 
 # Find dependencies
 plugin_add_dd4hep(${PLUGIN_NAME})
+plugin_add_eigen3(${PLUGIN_NAME})
 plugin_add_event_model(${PLUGIN_NAME})
+plugin_add_cern_root(${PLUGIN_NAME})
+
+plugin_link_libraries(${PLUGIN_NAME} ROOT::TMVA)
diff --git a/src/algorithms/fardetectors/FarDetectorLinearTracking.cc b/src/algorithms/fardetectors/FarDetectorLinearTracking.cc
@@ -25,9 +25,7 @@
 
 namespace eicrecon {
 
-    void FarDetectorLinearTracking::init(std::shared_ptr<spdlog::logger>& logger) {
-
-      m_log      = logger;
+    void FarDetectorLinearTracking::init() {
 
       // For changing how strongly each layer hit is in contributing to the fit
       m_layerWeights = Eigen::VectorXd::Constant(m_cfg.n_layer,1);
@@ -49,7 +47,7 @@ namespace eicrecon {
         // Check the number of input collections is correct
         int nCollections = inputhits.size();
         if(nCollections!=m_cfg.n_layer){
-          m_log->error("Wrong number of input collections passed to algorithm");
+          error("Wrong number of input collections passed to algorithm");
           return;
         }
 
@@ -58,7 +56,7 @@ namespace eicrecon {
         // TODO - Implement more sensible solution
         for(const auto& layerHits: inputhits){
           if((*layerHits).size()>m_cfg.layer_hits_max){
-            m_log->info("Too many hits in layer");
+            info("Too many hits in layer");
             return;
           }
         }
@@ -160,9 +158,9 @@ namespace eicrecon {
 
       double angle = std::acos(hitDiff.dot(m_optimumDirection));
 
-      m_log->debug("Vector: x={}, y={}, z={}",hitDiff.x(),hitDiff.y(),hitDiff.z());
-      m_log->debug("Optimum: x={}, y={}, z={}",m_optimumDirection.x(),m_optimumDirection.y(),m_optimumDirection.z());
-      m_log->debug("Angle: {}, Tolerance {}",angle,m_cfg.step_angle_tolerance);
+      debug("Vector: x={}, y={}, z={}",hitDiff.x(),hitDiff.y(),hitDiff.z());
+      debug("Optimum: x={}, y={}, z={}",m_optimumDirection.x(),m_optimumDirection.y(),m_optimumDirection.z());
+      debug("Angle: {}, Tolerance {}",angle,m_cfg.step_angle_tolerance);
 
       if(angle>m_cfg.step_angle_tolerance) return false;
 

diff --git a/src/algorithms/fardetectors/FarDetectorLinearTracking.h b/src/algorithms/fardetectors/FarDetectorLinearTracking.h
@@ -3,65 +3,54 @@
 
 #pragma once
 
+#include <Eigen/Core>
 #include <algorithms/algorithm.h>
 #include <algorithms/interfaces/WithPodConfig.h>
 #include <edm4eic/TrackSegmentCollection.h>
 #include <edm4hep/TrackerHitCollection.h>
 #include <gsl/pointers>
-#include <memory>
 #include <string>
 #include <string_view>
 #include <vector>
-#include <Eigen/Core>
-#include <spdlog/logger.h>
+
 #include "FarDetectorLinearTrackingConfig.h"
 
 namespace eicrecon {
 
-    using FarDetectorLinearTrackingAlgorithm = algorithms::Algorithm<
-        algorithms::Input<
-            std::vector<edm4hep::TrackerHitCollection>
-        >,
-        algorithms::Output<
-            edm4eic::TrackSegmentCollection
-        >
-    >;
-
-    class FarDetectorLinearTracking
-    : public FarDetectorLinearTrackingAlgorithm,
-      public WithPodConfig<FarDetectorLinearTrackingConfig>  {
-
-    public:
-        FarDetectorLinearTracking(std::string_view name)
-            : FarDetectorLinearTrackingAlgorithm{name,
-                {"inputHitCollections"},
-                {"outputTrackSegments"},
-                "Fit track segments from hits in the tracker layers"} {}
+using FarDetectorLinearTrackingAlgorithm =
+    algorithms::Algorithm<algorithms::Input<std::vector<edm4hep::TrackerHitCollection>>,
+                          algorithms::Output<edm4eic::TrackSegmentCollection>>;
 
-        /** One time initialization **/
-        void init(std::shared_ptr<spdlog::logger>& logger);
+class FarDetectorLinearTracking : public FarDetectorLinearTrackingAlgorithm,
+                                  public WithPodConfig<FarDetectorLinearTrackingConfig> {
 
-        /** Event by event processing **/
-        void process(const Input&, const Output&) const final;
+public:
+  FarDetectorLinearTracking(std::string_view name)
+      : FarDetectorLinearTrackingAlgorithm{name,
+                                           {"inputHitCollections"},
+                                           {"outputTrackSegments"},
+                                           "Fit track segments from hits in the tracker layers"} {}
 
-    private:
-        std::shared_ptr<spdlog::logger> m_log;
+  /** One time initialization **/
+  void init() final;
 
-        Eigen::VectorXd m_layerWeights;
+  /** Event by event processing **/
+  void process(const Input&, const Output&) const final;
 
-        Eigen::Vector3d m_optimumDirection;
+private:
+  Eigen::VectorXd m_layerWeights;
 
-        void buildMatrixRecursive(int level,
-                                 Eigen::MatrixXd* hitMatrix,
-                                 const std::vector<gsl::not_null<const edm4hep::TrackerHitCollection*>>& hits,
-                                 gsl::not_null<edm4eic::TrackSegmentCollection*> outputTracks) const;
+  Eigen::Vector3d m_optimumDirection;
 
-        void checkHitCombination(Eigen::MatrixXd* hitMatrix,
-                                gsl::not_null<edm4eic::TrackSegmentCollection*> outputTracks) const;
+  void
+  buildMatrixRecursive(int level, Eigen::MatrixXd* hitMatrix,
+                       const std::vector<gsl::not_null<const edm4hep::TrackerHitCollection*>>& hits,
+                       gsl::not_null<edm4eic::TrackSegmentCollection*> outputTracks) const;
 
-        bool checkHitPair(const Eigen::Vector3d& hit1,
-                          const Eigen::Vector3d& hit2) const;
+  void checkHitCombination(Eigen::MatrixXd* hitMatrix,
+                           gsl::not_null<edm4eic::TrackSegmentCollection*> outputTracks) const;
 
-    };
+  bool checkHitPair(const Eigen::Vector3d& hit1, const Eigen::Vector3d& hit2) const;
+};
 
-} // eicrecon
+} // namespace eicrecon