Rich cluster and ring reconstruction

- add fuzzy k clustering algorithms
- rename the clustering algorithm

JugReco/CMakeLists.txt

@@ -38,6 +38,8 @@ gaudi_add_module(JugRecoPlugins
   src/components/CrystalEndcapsReco.cpp 
   src/components/ClusterRecoCoG.cpp 
   src/components/PhotoMultiplierReco.cpp 
+  src/components/PhotoRingClusters.cpp 
+  src/components/FuzzyKClusters.cpp 
   LINK_LIBRARIES GaudiAlgLib GaudiKernel JugBase ROOT NPDet::DD4podIO EICD::eicd DDRec Acts)
 
 target_compile_options(JugRecoPlugins PRIVATE -Wno-suggest-override)

JugReco/src/components/FuzzyKClusters.cpp

+/*  Fuzzy K Clustering Algorithms
+ *
+ *  Author: Chao Peng (ANL)
+ *  Date: 10/14/2020
+ *
+ */
+
+#include "FuzzyKClusters.h"
+#include <exception>
+#include <iostream>
+#include <cmath>
+
+
+using namespace fkc;
+using namespace Eigen;
+
+
+// =================================================================================================
+//  KMeans Algorithm
+// =================================================================================================
+
+KMeans::KMeans()
+: n_iters(0), variance(0.)
+{
+    // place holder
+}
+
+KMeans::~KMeans()
+{
+    // place holder
+}
+
+MatrixXd KMeans::Fit(const MatrixXd &data, int k, double q, double epsilon, int max_iters)
+{
+    auto res = Initialize(data, k, q);
+
+    for (n_iters = 0; n_iters < max_iters; ++n_iters) {
+        auto old_mems = mems;
+        Distances(res, data);
+        Memberships(q);
+        FormClusters(res, data, q);
+
+        if ((old_mems - mems).cwiseAbs().maxCoeff() < epsilon) {
+            break;
+        }
+    }
+
+    return res;
+}
+
+// initialize and guess the clusters
+MatrixXd KMeans::Initialize(const MatrixXd &data, int k, double q)
+{
+    // resize matrices
+    dists.resize(k, data.rows());
+
+    // guess the cluster centers
+    mems = MatrixXd::Random(k, data.rows());
+    for (int j = 0; j < mems.cols(); ++j) {
+        auto csum = mems.col(j).sum();
+        for (int i = 0; i < mems.rows(); ++i) {
+            mems(i, j) = mems(i, j)/csum;
+        }
+    }
+
+    MatrixXd clusters(k, data.cols());
+    FormClusters(clusters, data, q);
+    return clusters;
+}
+
+// distance matrix (num_clusters, num_data)
+void KMeans::Distances(const MatrixXd &centroids, const MatrixXd &data)
+{
+    for (int i = 0; i < centroids.rows(); ++i) {
+        for (int j = 0; j < data.rows(); ++j) {
+            dists(i, j) = (centroids.row(i) - data.row(j)).cwiseAbs2().sum();
+        }
+    }
+}
+
+// membership matrix (num_clusters, num_data)
+void KMeans::Memberships(double q)
+{
+    // coeffcient-wise operation
+    auto d = dists.array().pow(-1.0/(q - 1.0)).matrix();
+
+    for (int j = 0; j < d.cols(); ++j) {
+        auto dsum = d.col(j).sum();
+        for (int i = 0; i < d.rows(); ++i) {
+            mems(i, j) = d(i, j)/dsum;
+        }
+    }
+}
+
+// rebuild clusters
+void KMeans::FormClusters(MatrixXd &clusters, const MatrixXd &data, double q)
+{
+    auto weights = mems.array().pow(q).matrix();
+    for (int i = 0; i < clusters.rows(); ++i) {
+        clusters.row(i) *= 0;
+        for (int j = 0; j < data.rows(); ++j) {
+            clusters.row(i) += data.row(j)*weights(i, j);
+        }
+        clusters.row(i) /= weights.row(i).sum();
+    }
+}
+
+
+// =================================================================================================
+//  KRings Algorithm, extended from KMeans
+//  Reference:
+//      [1] Y. H. Man and I. Gath,
+//          "Detection and separation of ring-shaped clusters using fuzzy clustering,"
+//          in IEEE Transactions on Pattern Analysis and Machine Intelligence,
+//          vol. 16, no. 8, pp. 855-861, Aug. 1994, doi: 10.1109/34.308484.
+// =================================================================================================
+
+KRings::KRings()
+: KMeans()
+{
+    // place holder
+}
+
+KRings::~KRings()
+{
+    // place holder
+}
+
+MatrixXd KRings::Fit(const MatrixXd &data, int k, double q, double epsilon, int max_iters)
+{
+    auto res = Initialize(data, k, q);
+
+    for (n_iters = 0; n_iters < max_iters; ++n_iters) {
+        auto old_mems = mems;
+        Distances(res, data);
+        Memberships(q);
+        FormRadii(res, q);
+        FormClusters(res, data, q);
+
+        if ((old_mems - mems).cwiseAbs().maxCoeff() < epsilon) {
+            break;
+        }
+    }
+
+    return res;
+}
+
+// initialize and guess the clusters
+MatrixXd KRings::Initialize(const MatrixXd &data, int k, double q)
+{
+    MatrixXd clusters(k, data.cols() + 1);
+    auto centers = clusters.leftCols(data.cols());
+
+    // call KMeans to help initialization
+    KMeans fkm;
+    centers = fkm.Fit(data, k, q, 1e-4, 5);
+
+    dists.resize(k, data.rows());
+    dists_euc = fkm.GetDistances().cwiseSqrt();
+    mems = fkm.GetMemberships();
+    FormRadii(clusters, q);
+    return clusters;
+}
+
+// distance matrix (num_clusters, num_data)
+void KRings::Distances(const MatrixXd &centroids, const MatrixXd &data)
+{
+    auto const centers = centroids.leftCols(centroids.cols() - 1);
+    auto const radii = centroids.rightCols(1);
+
+    for (int i = 0; i < centroids.rows(); ++i) {
+        for (int j = 0; j < data.rows(); ++j) {
+            dists_euc(i, j) = std::sqrt((centers.row(i) - data.row(j)).cwiseAbs2().sum());
+            dists(i, j) = std::pow(dists_euc(i, j) - radii(i, 0), 2);
+        }
+    }
+}
+
+// rebuild clusters radii
+void KRings::FormRadii(MatrixXd &clusters, double q)
+{
+    auto radii = clusters.rightCols(1);
+    auto weights = mems.array().pow(q).matrix();
+
+    for (int i = 0; i < weights.rows(); ++i) {
+        radii(i, 0) = 0;
+        for (int j = 0; j < weights.cols(); ++j) {
+            radii(i, 0) += weights(i, j)*dists_euc(i, j);
+        }
+        radii(i, 0) /= weights.row(i).sum();
+    }
+}
+
+// rebuild clusters centers
+void KRings::FormClusters(MatrixXd &clusters, const MatrixXd &data, double q)
+{
+    auto centers = clusters.leftCols(data.cols());
+    const auto &radii = clusters.rightCols(1);
+    auto weights = mems.array().pow(q).matrix();
+
+    for (int i = 0; i < weights.rows(); ++i) {
+        MatrixXd icenter = centers.row(i);
+        centers.row(i) *= 0;
+        for (int j = 0; j < weights.cols(); ++j) {
+            double scale = radii(i, 0)/dists_euc(i, j);
+            centers.row(i) += weights(i, j)*(data.row(j) - (data.row(j) - icenter)*scale);
+        }
+        centers.row(i) /= weights.row(i).sum();
+    }
+}
+
+
+// =================================================================================================
+//  KEllipses Algorithm, extended from KRings
+//  Reference:
+//      [1] I. Gath and D. Hoory, Pattern Recognition Letters 16 (1995) 727-741,
+//          https://doi.org/10.1016/0167-8655(95)00030-K.
+// =================================================================================================
+// @TODO
+

JugReco/src/components/FuzzyKClusters.h

+#pragma once
+
+#include <Eigen/Dense>
+
+namespace fkc {
+
+class KMeans
+{
+public:
+    KMeans();
+    virtual ~KMeans();
+
+    virtual Eigen::MatrixXd Fit(const Eigen::MatrixXd &data, int k,
+                                double q = 2.0, double epsilon = 1e-4, int max_iters = 1000);
+
+    int NIters() const { return n_iters; }
+    double Variance() const { return variance; }
+    const Eigen::MatrixXd &GetDistances() const { return dists; }
+    const Eigen::MatrixXd &GetMemberships() const { return mems; }
+
+    Eigen::MatrixXd &GetDistances() { return dists; }
+    Eigen::MatrixXd &GetMemberships() { return mems; }
+
+protected:
+    virtual Eigen::MatrixXd Initialize(const Eigen::MatrixXd &data, int k, double q);
+    virtual void Distances(const Eigen::MatrixXd &centroids, const Eigen::MatrixXd &data);
+    virtual void Memberships(double q);
+    virtual void FormClusters(Eigen::MatrixXd &clusters, const Eigen::MatrixXd &data, double q);
+
+protected:
+    int n_iters;
+    double variance;
+    Eigen::MatrixXd dists, mems;
+};
+
+class KRings : public KMeans
+{
+public:
+    KRings();
+    ~KRings();
+
+    virtual Eigen::MatrixXd Fit(const Eigen::MatrixXd &data, int k,
+                                double q = 2.0, double epsilon = 1e-4, int max_iters = 1000);
+
+protected:
+    virtual Eigen::MatrixXd Initialize(const Eigen::MatrixXd &data, int k, double q);
+    virtual void Distances(const Eigen::MatrixXd &centroids, const Eigen::MatrixXd &data);
+    virtual void FormClusters(Eigen::MatrixXd &clusters, const Eigen::MatrixXd &data, double q);
+    virtual void FormRadii(Eigen::MatrixXd &clusters, double g);
+
+protected:
+    Eigen::MatrixXd dists_euc;
+};
+
+} // namespace fkc
+

JugReco/src/components/PhotoMultiplierReco.cpp

 /*  General PhotoMultiplier Reconstruction
  *
- *  Apply the given quantum efficiency for photon detection
- *  Converts the number of s to signal amplitude
+ *  Estimate the number of photo-electrons and convert timeStamp to time
+ *  Collect cell information
  *
  *  Author: Chao Peng (ANL)
  *  Date: 10/03/2020
@@ -40,7 +40,7 @@ public:
     DataHandle<eic::PMTHitCollection>
         m_outputHitCollection{"outputHitCollection", Gaudi::DataHandle::Writer, this};
     Gaudi::Property<double> m_timeStep{this, "timeStep", 0.0625*ns};
-    Gaudi::Property<double> m_minNpe{this, "minNpe", 0.5};
+    Gaudi::Property<double> m_minNpe{this, "minNpe", 0.0};
     Gaudi::Property<double> m_speMean{this, "speMean", 80.0};
     Gaudi::Property<double> m_pedMean{this, "pedMean", 200.0};
     /// Pointer to the geometry service

JugReco/src/components/PhotoRingClusters.cpp

+/*  Clustering Algorithm for Ring Imaging Cherenkov (RICH) events
+ *
+ *  Author: Chao Peng (ANL)
+ *  Date: 10/04/2020
+ *
+ */
+
+#include <algorithm>
+
+#include "Gaudi/Property.h"
+#include "GaudiAlg/GaudiAlgorithm.h"
+#include "GaudiKernel/ToolHandle.h"
+#include "GaudiAlg/Transformer.h"
+#include "GaudiAlg/GaudiTool.h"
+#include "GaudiKernel/RndmGenerators.h"
+#include "GaudiKernel/PhysicalConstants.h"
+
+#include "DDRec/CellIDPositionConverter.h"
+#include "DDRec/SurfaceManager.h"
+#include "DDRec/Surface.h"
+
+// FCCSW
+#include "JugBase/DataHandle.h"
+#include "JugBase/IGeoSvc.h"
+
+// Event Model related classes
+#include "eicd/PMTHitCollection.h"
+#include "eicd/RIChClusterCollection.h"
+#include "FuzzyKClusters.h"
+
+using namespace Gaudi::Units;
+using namespace Eigen;
+
+
+namespace Jug::Reco {
+class PhotoRingClusters : public GaudiAlgorithm
+{
+public:
+    DataHandle<eic::PMTHitCollection>
+        m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this};
+    DataHandle<eic::RIChClusterCollection>
+        m_outputClusterCollection{"outputClusterCollection", Gaudi::DataHandle::Writer, this};
+    Gaudi::Property<double> m_minNpe{this, "minNpe", 0.0};
+    /// Pointer to the geometry service
+    SmartIF<IGeoSvc> m_geoSvc;
+
+    // ill-formed: using GaudiAlgorithm::GaudiAlgorithm;
+    PhotoRingClusters(const std::string& name, ISvcLocator* svcLoc)
+        : GaudiAlgorithm(name, svcLoc)
+    {
+        declareProperty("inputHitCollection",   m_inputHitCollection,   "");
+        declareProperty("outputClusterCollection",  m_outputClusterCollection,  "");
+    }
+
+    StatusCode initialize() override
+    {
+        if (GaudiAlgorithm::initialize().isFailure()) {
+            return StatusCode::FAILURE;
+        }
+        m_geoSvc = service("GeoSvc");
+        if (!m_geoSvc) {
+            error() << "Unable to locate Geometry Service. "
+                    << "Make sure you have GeoSvc and SimSvc in the right order in the configuration." << endmsg;
+            return StatusCode::FAILURE;
+        }
+        return StatusCode::SUCCESS;
+    }
+
+    StatusCode execute() override
+    {
+        // input collections
+	    const auto &rawhits = *m_inputHitCollection.get();
+        // Create output collections
+        auto &clusters = *m_outputClusterCollection.createAndPut();
+
+        // algorithm
+        auto alg = fkc::KRings();
+
+        MatrixXd data(rawhits.size(), 2);
+        for (int i = 0; i < data.rows(); ++i) {
+            data.row(i) << rawhits[i].local_x(), rawhits[i].local_y();
+        }
+
+        return StatusCode::SUCCESS;
+    }
+};
+
+DECLARE_COMPONENT(PhotoRingClusters)
+
+} // namespace Jug::Reco
+
+