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