Factorize cluster shape parameter calculation

src/algorithms/calorimetry/CalorimeterClusterRecoCoG.cc

@@ -20,9 +20,6 @@
 #include <fmt/core.h>
 #include <podio/ObjectID.h>
 #include <podio/RelationRange.h>
-#include <Eigen/Core>
-#include <Eigen/Eigenvalues>
-#include <Eigen/Householder> // IWYU pragma: keep
 #include <cctype>
 #include <complex>
 #include <cstddef>
@@ -178,115 +175,6 @@ std::optional<edm4eic::MutableCluster> CalorimeterClusterRecoCoG::reconstruct(co
       debug("Bound cluster position to contributing hits due to {}", (overflow ? "overflow" : "underflow"));
     }
   }
-
-  // Additional convenience variables
-
-  //_______________________________________
-  // Calculate cluster profile:
-  //    radius,
-  //    dispersion (energy weighted radius),
-  //    theta-phi cluster widths (2D)
-  //    x-y-z cluster widths (3D)
-  float radius = 0, dispersion = 0, w_sum = 0;
-
-  Eigen::Matrix2f sum2_2D = Eigen::Matrix2f::Zero();
-  Eigen::Matrix3f sum2_3D = Eigen::Matrix3f::Zero();
-  Eigen::Vector2f sum1_2D = Eigen::Vector2f::Zero();
-  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.
-  edm4hep::Vector3f axis;
-
-  if (cl.getNhits() > 1) {
-    for (const auto& hit : pcl.getHits()) {
-      float w = weightFunc(hit.getEnergy(), totalE, logWeightBase, 0);
-
-      // theta, phi
-      Eigen::Vector2f pos2D( edm4hep::utils::anglePolar( hit.getPosition() ), edm4hep::utils::angleAzimuthal( hit.getPosition() ) );
-      // x, y, z
-      Eigen::Vector3f pos3D( hit.getPosition().x, hit.getPosition().y, hit.getPosition().z );
-
-      const auto delta = cl.getPosition() - hit.getPosition();
-      radius          += delta * delta;
-      dispersion      += delta * delta * w;
-
-      // Weighted Sum x*x, x*y, x*z, y*y, etc.
-      sum2_2D += w * pos2D * pos2D.transpose();
-      sum2_3D += w * pos3D * pos3D.transpose();
-
-      // Weighted Sum x, y, z
-      sum1_2D += w * pos2D;
-      sum1_3D += w * pos3D;
-
-      w_sum += w;
-    }
-
-    radius     = sqrt((1. / (cl.getNhits() - 1.)) * radius);
-    if( w_sum > 0 ) {
-      dispersion = sqrt( dispersion / w_sum );
-
-      // normalize matrices
-      sum2_2D /= w_sum;
-      sum2_3D /= w_sum;
-      sum1_2D /= w_sum;
-      sum1_3D /= w_sum;
-
-      // 2D and 3D covariance matrices
-      Eigen::Matrix2f cov2 = sum2_2D - sum1_2D * sum1_2D.transpose();
-      Eigen::Matrix3f cov3 = sum2_3D - sum1_3D * sum1_3D.transpose();
-
-      // 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, 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_eigen = eigenvectors.col(std::distance(
-            eigenValues_3D.begin(),
-            max_eigenvalue_it
-        ));
-      axis = {
-        axis_eigen(0,0).real(),
-        axis_eigen(1,0).real(),
-        axis_eigen(2,0).real(),
-      };
-    }
-  }
-
-  cl.addToShapeParameters( radius );
-  cl.addToShapeParameters( dispersion );
-  cl.addToShapeParameters( eigenValues_2D[0].real() ); // 2D theta-phi cluster width 1
-  cl.addToShapeParameters( eigenValues_2D[1].real() ); // 2D theta-phi cluster width 2
-  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
-
-
-  double dot_product = cl.getPosition() * axis;
-  if (dot_product < 0) {
-    axis = -1 * axis;
-  }
-
-  if (m_cfg.longitudinalShowerInfoAvailable) {
-    cl.setIntrinsicTheta(edm4hep::utils::anglePolar(axis));
-    cl.setIntrinsicPhi(edm4hep::utils::angleAzimuthal(axis));
-    // TODO intrinsicDirectionError
-  } else {
-    cl.setIntrinsicTheta(NAN);
-    cl.setIntrinsicPhi(NAN);
-  }
-
   return std::move(cl);
 }
 

src/algorithms/calorimetry/CalorimeterClusterRecoCoGConfig.h

@@ -23,8 +23,6 @@ namespace eicrecon {
         std::vector<double> logWeightBaseCoeffs{};
         double logWeightBase_Eref = 50 * dd4hep::MeV;
 
-        bool longitudinalShowerInfoAvailable = false;
-
         // Constrain the cluster position eta to be within
         // the eta of the contributing hits. This is useful to avoid edge effects
         // for endcaps.

src/algorithms/calorimetry/ClusterShapeCalculator.cc

@@ -0,0 +1,194 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2025 Chao Peng, Dhevan Gangadharan, Sebouh Paul, Derek Anderson
+
+// algorithm definition
+#include "ClusterShapeCalculator.h"
+
+#include <edm4hep/Vector3f.h>
+#include <edm4hep/utils/vector_utils.h>
+#include <Eigen/Core>
+#include <Eigen/Eigenvalues>
+#include <Eigen/Householder>
+#include <cmath>
+
+
+
+namespace eicrecon {
+
+  // --------------------------------------------------------------------------
+  //! Initialize algorithm
+  // --------------------------------------------------------------------------
+  void ClusterShapeCalculator::init() {
+
+    //... nothing to do ...//
+
+  }  // end 'init(dd4hep::Detector*)'
+
+
+
+  // --------------------------------------------------------------------------
+  //! Process inputs
+  // --------------------------------------------------------------------------
+  /*! Primary algorithm call: algorithm ingests a collection of clusters
+   *  and computes their cluster shape parameters.  Clusters are copied
+   *  onto output with computed shape parameters.  If associations are
+   *  provided, they are copied to the output.
+   *
+   *  Parameters calculated:
+   *    - radius,
+   *    - dispersion (energy weighted radius),
+   *    - theta-phi cluster widths (2D)
+   *    - x-y-z cluster widths (3D)
+   */
+  void ClusterShapeCalculator::process(
+    const ClusterShapeCalculator::Input& input,
+    const ClusterShapeCalculator::Output& output
+  ) const {
+
+    // grab inputs/outputs
+    const auto [in_clusters, in_associations] = input;
+    auto [out_clusters, out_associations] = output;
+
+    // exit if no clusters in collection
+    if (in_clusters->size() == 0) {
+      debug("No clusters in input collection.");
+      return;
+    }
+
+    // loop over input clusters
+    for (const auto& in_clust : *in_clusters) {
+
+      // copy input cluster
+      edm4eic::MutableCluster out_clust = in_clust.clone();
+
+      // ----------------------------------------------------------------------
+      // do shape parameter calculation
+      // ----------------------------------------------------------------------
+      {
+
+        // create addresses for quantities we'll need later
+        float radius = 0, dispersion = 0, w_sum = 0;
+
+        // set up matrices/vectors
+        Eigen::Matrix2f sum2_2D = Eigen::Matrix2f::Zero();
+        Eigen::Matrix3f sum2_3D = Eigen::Matrix3f::Zero();
+        Eigen::Vector2f sum1_2D = Eigen::Vector2f::Zero();
+        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.
+        edm4hep::Vector3f axis;
+        if (out_clust.getNhits() > 1) {
+          for (std::size_t iHit = 0; const auto& hit : out_clust.getHits()) {
+
+            // get weight of hit
+            const float w = out_clust.getHitContributions()[iHit] / hit.getEnergy();
+
+            // theta, phi
+            Eigen::Vector2f pos2D( edm4hep::utils::anglePolar( hit.getPosition() ), edm4hep::utils::angleAzimuthal( hit.getPosition() ) );
+            // x, y, z
+            Eigen::Vector3f pos3D( hit.getPosition().x, hit.getPosition().y, hit.getPosition().z );
+
+            const auto delta = out_clust.getPosition() - hit.getPosition();
+            radius     += delta * delta;
+            dispersion += delta * delta * w;
+
+            // Weighted Sum x*x, x*y, x*z, y*y, etc.
+            sum2_2D += w * pos2D * pos2D.transpose();
+            sum2_3D += w * pos3D * pos3D.transpose();
+
+            // Weighted Sum x, y, z
+            sum1_2D += w * pos2D;
+            sum1_3D += w * pos3D;
+
+            w_sum += w;
+            ++iHit;
+          }  // end hit loop
+
+          radius = sqrt((1. / (out_clust.getNhits() - 1.)) * radius);
+          if( w_sum > 0 ) {
+            dispersion = sqrt( dispersion / w_sum );
+
+            // normalize matrices
+            sum2_2D /= w_sum;
+            sum2_3D /= w_sum;
+            sum1_2D /= w_sum;
+            sum1_3D /= w_sum;
+
+            // 2D and 3D covariance matrices
+            Eigen::Matrix2f cov2 = sum2_2D - sum1_2D * sum1_2D.transpose();
+            Eigen::Matrix3f cov3 = sum2_3D - sum1_3D * sum1_3D.transpose();
+
+            // Solve for eigenvalues.  Corresponds to out_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, 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_eigen = eigenvectors.col(std::distance(
+                  eigenValues_3D.begin(),
+                  max_eigenvalue_it
+              ));
+            axis = {
+              axis_eigen(0,0).real(),
+              axis_eigen(1,0).real(),
+              axis_eigen(2,0).real(),
+            };
+          }  // end if weight sum is nonzero
+        }  // end if n hits > 1
+
+        // set shape parameters
+        out_clust.addToShapeParameters( radius );
+        out_clust.addToShapeParameters( dispersion );
+        out_clust.addToShapeParameters( eigenValues_2D[0].real() ); // 2D theta-phi out_cluster width 1
+        out_clust.addToShapeParameters( eigenValues_2D[1].real() ); // 2D theta-phi out_cluster width 2
+        out_clust.addToShapeParameters( eigenValues_3D[0].real() ); // 3D x-y-z out_cluster width 1
+        out_clust.addToShapeParameters( eigenValues_3D[1].real() ); // 3D x-y-z out_cluster width 2
+        out_clust.addToShapeParameters( eigenValues_3D[2].real() ); // 3D x-y-z out_cluster width 3
+
+        // check axis orientation
+        double dot_product = out_clust.getPosition() * axis;
+        if (dot_product < 0) {
+          axis = -1 * axis;
+        }
+
+        // set intrinsic theta/phi
+        if (m_cfg.longitudinalShowerInfoAvailable) {
+          out_clust.setIntrinsicTheta( edm4hep::utils::anglePolar(axis) );
+          out_clust.setIntrinsicPhi( edm4hep::utils::angleAzimuthal(axis) );
+          // TODO intrinsicDirectionError
+        } else {
+          out_clust.setIntrinsicTheta(NAN);
+          out_clust.setIntrinsicPhi(NAN);
+        }
+      }  // end shape parameter calculation
+
+      // ----------------------------------------------------------------------
+      // if provided, copy associations
+      // ----------------------------------------------------------------------
+      for (auto in_assoc : *in_associations) {
+        if (in_assoc.getRec() == in_clust) {
+          auto mc_par = in_assoc.getSim();
+          auto out_assoc = out_associations->create();
+          out_assoc.setRecID( out_clust.getObjectID().index );
+          out_assoc.setSimID( mc_par.getObjectID().index );
+          out_assoc.setRec( out_clust );
+          out_assoc.setSim( mc_par );
+          out_assoc.setWeight( in_assoc.getWeight() );
+        }
+      }  // end input association loop
+    }  // end input cluster loop
+  }  // end 'process(Input&, Output&)'
+
+}  // end eicrecon namespace

src/algorithms/calorimetry/ClusterShapeCalculator.h

@@ -0,0 +1,65 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2025 Chao Peng, Dhevan Gangadharan, Sebouh Paul, Derek Anderson
+
+#pragma once
+
+#include "algorithms/interfaces/WithPodConfig.h"
+#include "ClusterShapeCalculatorConfig.h"
+
+#include <algorithms/algorithm.h>
+#include <edm4eic/ClusterCollection.h>
+#include <edm4eic/MCRecoClusterParticleAssociationCollection.h>
+#include <optional>
+#include <string>
+#include <string_view>
+
+
+
+namespace eicrecon {
+
+  // --------------------------------------------------------------------------
+  //! Algorithm input/output
+  // --------------------------------------------------------------------------
+  using ClusterShapeCalculatorAlgorithm = algorithms::Algorithm<
+    algorithms::Input<
+      edm4eic::ClusterCollection,
+      std::optional<edm4eic::MCRecoClusterParticleAssociationCollection>
+    >,
+    algorithms::Output<
+      edm4eic::ClusterCollection,
+      std::optional<edm4eic::MCRecoClusterParticleAssociationCollection>
+    >
+  >;
+
+
+
+  // --------------------------------------------------------------------------
+  //! Calculate cluster shapes for provided clusters
+  // --------------------------------------------------------------------------
+  /*! An algorithm which takes a collection of clusters,
+   *  computes their cluster shape parameters, and saves
+   *  outputs the same clusters with computed parameters.
+   */
+  class ClusterShapeCalculator
+    : public ClusterShapeCalculatorAlgorithm
+    , public WithPodConfig<ClusterShapeCalculatorConfig>
+  {
+
+    public:
+
+      // ctor
+      ClusterShapeCalculator(std::string_view name) :
+        ClusterShapeCalculatorAlgorithm {
+          name,
+          {"inputClusters", "inputMCClusterAssociations"},
+          {"outputClusters", "outputMCClusterAssociations"},
+          "Computes cluster shape parameters"
+        } {}
+
+      // public methods
+      void init() final;
+      void process(const Input&, const Output&) const final;
+
+  };
+
+}  // namespace eicrecon