algorithms/calorimetry/ImagingTopoCluster.h

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 Chao Peng, Sylvester Joosten, Whitney Armstrong, Wouter Deconinck
0003
0004 /*
0005  *  Topological Cell Clustering Algorithm for Imaging Calorimetry
0006  *  1. group all the adjacent pixels
0007  *
0008  *  Author: Chao Peng (ANL), 06/02/2021
0009  *  Original reference: https://arxiv.org/pdf/1603.02934.pdf
0010  *
0011  *  Modifications:
0012  *
0013  *  Wouter Deconinck (Manitoba), 08/24/2024
0014  *  - converted hit storage model from std::vector to std::set sorted on layer
0015  *    where only hits remaining to be assigned to a group are in the set
0016  *  - erase hits that are too low in energy to be part of a cluster
0017  *  - converted group storage model from std::set to std::list to allow adding
0018  *    hits while keeping iterators valid
0019  *
0020  */
0021 #pragma once
0022
0023 #include <algorithms/algorithm.h>
0024 // Event Model related classes
0025 #include <edm4eic/CalorimeterHitCollection.h>
0026 #include <edm4eic/ProtoClusterCollection.h>
0027 #include <array>
0028 #include <cstddef>
0029 #include <list>
0030 #include <set>
0031 #include <string>
0032 #include <string_view>
0033
0034 #include "ImagingTopoClusterConfig.h"
0035 #include "algorithms/interfaces/WithPodConfig.h"
0036
0037 namespace eicrecon {
0038
0039 using ImagingTopoClusterAlgorithm =
0040     algorithms::Algorithm<algorithms::Input<edm4eic::CalorimeterHitCollection>,
0041                           algorithms::Output<edm4eic::ProtoClusterCollection>>;
0042
0043 class ImagingTopoCluster : public ImagingTopoClusterAlgorithm,
0044                            public WithPodConfig<ImagingTopoClusterConfig> {
0045
0046 public:
0047   ImagingTopoCluster(std::string_view name)
0048       : ImagingTopoClusterAlgorithm{
0049             name,
0050             {"inputHitCollection"},
0051             {"outputProtoClusterCollection"},
0052             "Topological cell clustering algorithm for imaging calorimetry."} {}
0053
0054 private:
0055   // unitless counterparts of the input parameters
0056   std::array<double, 2> sameLayerDistXY{0, 0};
0057   std::array<double, 2> diffLayerDistXY{0, 0};
0058   std::array<double, 3> sameLayerDistXYZ{0, 0, 0};
0059   std::array<double, 3> diffLayerDistXYZ{0, 0, 0};
0060   std::array<double, 2> sameLayerDistEtaPhi{0, 0};
0061   std::array<double, 2> diffLayerDistEtaPhi{0, 0};
0062   std::array<double, 2> sameLayerDistTZ{0, 0};
0063   std::array<double, 2> diffLayerDistTZ{0, 0};
0064   double sectorDist{0};
0065   double minClusterHitEdep{0};
0066   double minClusterCenterEdep{0};
0067   double minClusterEdep{0};
0068
0069 public:
0070   void init();
0071   void process(const Input& input, const Output& output) const final;
0072
0073 private:
0074   // helper function to group hits
0075   bool is_neighbour(const edm4eic::CalorimeterHit& h1, const edm4eic::CalorimeterHit& h2) const;
0076
0077   // grouping function with Breadth-First Search
0078   // note: template to allow Compare only known in local scope of caller
0079   template <typename Compare>
0080   void bfs_group(const edm4eic::CalorimeterHitCollection& hits,
0081                  std::set<std::size_t, Compare>& indices, std::list<std::size_t>& group,
0082                  const std::size_t idx) const {
0083
0084     // loop over group as it grows, until the end is stable and we reach it
0085     for (auto idx1 = group.begin(); idx1 != group.end(); ++idx1) {
0086       // check neighbours (note comments on loop over set above)
0087       for (auto idx2 = indices.begin(); idx2 != indices.end();
0088            indices.empty() ? idx2 = indices.end() : idx2) {
0089
0090         // skip idx1 and original idx
0091         // (we cannot erase idx since it would invalidate iterator in calling scope)
0092         if (*idx2 == *idx1 || *idx2 == idx) {
0093           idx2++;
0094           continue;
0095         }
0096
0097         // skip rest of list of hits when we're past relevant layers
0098         //if (hits[*idx2].getLayer() - hits[*idx1].getLayer() > m_cfg.neighbourLayersRange) {
0099         //  break;
0100         //}
0101
0102         // not energetic enough for cluster hit
0103         if (hits[*idx2].getEnergy() < m_cfg.minClusterHitEdep) {
0104           idx2 = indices.erase(idx2);
0105           continue;
0106         }
0107
0108         if (is_neighbour(hits[*idx1], hits[*idx2])) {
0109           group.push_back(*idx2);
0110           idx2 = indices.erase(idx2); // takes role of idx2++
0111         } else {
0112           idx2++;
0113         }
0114       }
0115     }
0116   }
0117 };
0118
0119 } // namespace eicrecon