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, 2> sameLayerDistEtaPhi{0, 0};
0059   std::array<double, 2> diffLayerDistEtaPhi{0, 0};
0060   std::array<double, 2> sameLayerDistTZ{0, 0};
0061   std::array<double, 2> diffLayerDistTZ{0, 0};
0062   double sectorDist{0};
0063   double minClusterHitEdep{0};
0064   double minClusterCenterEdep{0};
0065   double minClusterEdep{0};
0066
0067 public:
0068   void init();
0069   void process(const Input& input, const Output& output) const final;
0070
0071 private:
0072   // helper function to group hits
0073   bool is_neighbour(const edm4eic::CalorimeterHit& h1, const edm4eic::CalorimeterHit& h2) const;
0074
0075   // grouping function with Breadth-First Search
0076   // note: template to allow Compare only known in local scope of caller
0077   template <typename Compare>
0078   void bfs_group(const edm4eic::CalorimeterHitCollection& hits,
0079                  std::set<std::size_t, Compare>& indices, std::list<std::size_t>& group,
0080                  const std::size_t idx) const {
0081
0082     // loop over group as it grows, until the end is stable and we reach it
0083     for (auto idx1 = group.begin(); idx1 != group.end(); ++idx1) {
0084       // check neighbours (note comments on loop over set above)
0085       for (auto idx2 = indices.begin(); idx2 != indices.end();
0086            indices.empty() ? idx2 = indices.end() : idx2) {
0087
0088         // skip idx1 and original idx
0089         // (we cannot erase idx since it would invalidate iterator in calling scope)
0090         if (*idx2 == *idx1 || *idx2 == idx) {
0091           idx2++;
0092           continue;
0093         }
0094
0095         // skip rest of list of hits when we're past relevant layers
0096         //if (hits[*idx2].getLayer() - hits[*idx1].getLayer() > m_cfg.neighbourLayersRange) {
0097         //  break;
0098         //}
0099
0100         // not energetic enough for cluster hit
0101         if (hits[*idx2].getEnergy() < m_cfg.minClusterHitEdep) {
0102           idx2 = indices.erase(idx2);
0103           continue;
0104         }
0105
0106         if (is_neighbour(hits[*idx1], hits[*idx2])) {
0107           group.push_back(*idx2);
0108           idx2 = indices.erase(idx2); // takes role of idx2++
0109         } else {
0110           idx2++;
0111         }
0112       }
0113     }
0114   }
0115 };
0116
0117 } // namespace eicrecon