algorithms/calorimetry/CalorimeterIslandCluster.cc

0001 // Copyright (C) 2022, 2023 Chao Peng, Wouter Deconinck, Sylvester Joosten, Dmitry Kalinkin, David Lawrence
0002 // SPDX-License-Identifier: LGPL-3.0-or-later
0003
0004 // References:
0005 //   https://cds.cern.ch/record/687345/files/note01_034.pdf
0006 //   https://www.jlab.org/primex/weekly_meetings/primexII/slides_2012_01_20/island_algorithm.pdf
0007
0008 #include <DD4hep/Readout.h>
0009 #include <Evaluator/DD4hepUnits.h>
0010 #include <algorithms/service.h>
0011 #include <edm4hep/Vector2f.h>
0012 #include <edm4hep/Vector3f.h>
0013 #include <fmt/format.h>
0014 #include <algorithm>
0015 #include <gsl/pointers>
0016 #include <map>
0017 #include <set>
0018 #include <stdexcept>
0019 #include <string>
0020 #include <tuple>
0021 #include <unordered_map>
0022 #include <utility>
0023 #include <vector>
0024
0025 #include "CalorimeterIslandCluster.h"
0026 #include "algorithms/calorimetry/CalorimeterIslandClusterConfig.h"
0027 #include "services/evaluator/EvaluatorSvc.h"
0028
0029 using namespace edm4eic;
0030
0031 namespace eicrecon {
0032
0033 static double Phi_mpi_pi(double phi) {
0034   return std::remainder(phi, 2 * M_PI);
0035 }
0036
0037 static edm4hep::Vector2f localDistXY(const CaloHit &h1, const CaloHit &h2) {
0038   const auto delta =h1.getLocal() - h2.getLocal();
0039   return {delta.x, delta.y};
0040 }
0041 static edm4hep::Vector2f localDistXZ(const CaloHit &h1, const CaloHit &h2) {
0042   const auto delta = h1.getLocal() - h2.getLocal();
0043   return {delta.x, delta.z};
0044 }
0045 static edm4hep::Vector2f localDistYZ(const CaloHit &h1, const CaloHit &h2) {
0046   const auto delta = h1.getLocal() - h2.getLocal();
0047   return {delta.y, delta.z};
0048 }
0049 static edm4hep::Vector2f dimScaledLocalDistXY(const CaloHit &h1, const CaloHit &h2) {
0050   const auto delta = h1.getLocal() - h2.getLocal();
0051
0052   const auto dimsum = h1.getDimension() + h2.getDimension();
0053
0054   return {2 * delta.x / dimsum.x, 2 * delta.y / dimsum.y};
0055 }
0056 static edm4hep::Vector2f globalDistRPhi(const CaloHit &h1, const CaloHit &h2) {
0057   using vector_type = decltype(edm4hep::Vector2f::a);
0058   return {
0059     static_cast<vector_type>(
0060       edm4hep::utils::magnitude(h1.getPosition()) - edm4hep::utils::magnitude(h2.getPosition())
0061     ),
0062     static_cast<vector_type>(
0063       Phi_mpi_pi(edm4hep::utils::angleAzimuthal(h1.getPosition()) - edm4hep::utils::angleAzimuthal(h2.getPosition()))
0064     )
0065   };
0066 }
0067 static edm4hep::Vector2f globalDistEtaPhi(const CaloHit &h1, const CaloHit &h2) {
0068   using vector_type = decltype(edm4hep::Vector2f::a);
0069   return {
0070     static_cast<vector_type>(
0071       edm4hep::utils::eta(h1.getPosition()) - edm4hep::utils::eta(h2.getPosition())
0072     ),
0073     static_cast<vector_type>(
0074       Phi_mpi_pi(edm4hep::utils::angleAzimuthal(h1.getPosition()) - edm4hep::utils::angleAzimuthal(h2.getPosition()))
0075     )
0076   };
0077 }
0078
0079 //------------------------
0080 // AlgorithmInit
0081 //------------------------
0082 void CalorimeterIslandCluster::init() {
0083
0084     static std::map<std::string,
0085                 std::tuple<std::function<edm4hep::Vector2f(const CaloHit&, const CaloHit&)>, std::vector<double>>>
0086     distMethods{
0087         {"localDistXY", {localDistXY, {dd4hep::mm, dd4hep::mm}}},        {"localDistXZ", {localDistXZ, {dd4hep::mm, dd4hep::mm}}},
0088         {"localDistYZ", {localDistYZ, {dd4hep::mm, dd4hep::mm}}},        {"dimScaledLocalDistXY", {dimScaledLocalDistXY, {1., 1.}}},
0089         {"globalDistRPhi", {globalDistRPhi, {dd4hep::mm, dd4hep::rad}}}, {"globalDistEtaPhi", {globalDistEtaPhi, {1., dd4hep::rad}}}
0090     };
0091
0092
0093     // set coordinate system
0094     auto set_dist_method = [this](std::pair<std::string, std::vector<double>> uprop) {
0095       if (uprop.second.size() == 0) {
0096         return false;
0097       }
0098       auto& [method, units] = distMethods[uprop.first];
0099       if (uprop.second.size() != units.size()) {
0100         warning("Expect {} values from {}, received {}. ignored it.", units.size(), uprop.first,  uprop.second.size());
0101         return false;
0102       } else {
0103         for (size_t i = 0; i < units.size(); ++i) {
0104           neighbourDist[i] = uprop.second[i] / units[i];
0105         }
0106         hitsDist = method;
0107         info("Clustering uses {} with distances <= [{}]", uprop.first, fmt::join(neighbourDist, ","));
0108       }
0109       return true;
0110     };
0111
0112     std::vector<std::pair<std::string, std::vector<double>>> uprops{
0113             {"localDistXY", m_cfg.localDistXY},
0114             {"localDistXZ", m_cfg.localDistXZ},
0115             {"localDistYZ", m_cfg.localDistYZ},
0116             {"globalDistRPhi", m_cfg.globalDistRPhi},
0117             {"globalDistEtaPhi", m_cfg.globalDistEtaPhi},
0118             // default one should be the last one
0119             {"dimScaledLocalDistXY", m_cfg.dimScaledLocalDistXY}
0120     };
0121
0122     auto& serviceSvc = algorithms::ServiceSvc::instance();
0123
0124     std::function hit_pair_to_map = [this](const edm4eic::CalorimeterHit &h1, const edm4eic::CalorimeterHit &h2) {
0125       std::unordered_map<std::string, double> params;
0126       for(const auto &p : m_idSpec.fields()) {
0127         const std::string &name = p.first;
0128         const dd4hep::IDDescriptor::Field* field = p.second;
0129         params.emplace(name + "_1", field->value(h1.getCellID()));
0130         params.emplace(name + "_2", field->value(h2.getCellID()));
0131         trace("{}_1 = {}", name, field->value(h1.getCellID()));
0132         trace("{}_2 = {}", name, field->value(h2.getCellID()));
0133       }
0134       return params;
0135     };
0136
0137     if (m_cfg.readout.empty()) {
0138       if ((!m_cfg.adjacencyMatrix.empty()) || (!m_cfg.peakNeighbourhoodMatrix.empty())) {
0139         throw std::runtime_error("'readout' is not provided, it is needed to know the fields in readout ids");
0140       }
0141     } else {
0142       m_idSpec = m_detector->readout(m_cfg.readout).idSpec();
0143     }
0144
0145     bool method_found = false;
0146
0147     // Adjacency matrix methods
0148     if (!m_cfg.adjacencyMatrix.empty()) {
0149       is_neighbour = serviceSvc.service<EvaluatorSvc>("EvaluatorSvc")->compile(m_cfg.adjacencyMatrix, hit_pair_to_map);
0150       method_found = true;
0151     }
0152
0153     // Coordinate distance methods
0154     if (not method_found) {
0155       for (auto& uprop : uprops) {
0156         if (set_dist_method(uprop)) {
0157           method_found = true;
0158
0159           is_neighbour = [this](const CaloHit &h1, const CaloHit &h2) {
0160             // in the same sector
0161             if (h1.getSector() == h2.getSector()) {
0162               auto dist = hitsDist(h1, h2);
0163               return (fabs(dist.a) <= neighbourDist[0]) && (fabs(dist.b) <= neighbourDist[1]);
0164               // different sector, local coordinates do not work, using global coordinates
0165             } else {
0166               // sector may have rotation (barrel), so z is included
0167               // (EDM4hep units are mm, so convert sectorDist to mm)
0168               return (edm4hep::utils::magnitude(h1.getPosition() - h2.getPosition()) <= m_cfg.sectorDist / dd4hep::mm);
0169             }
0170           };
0171
0172           break;
0173         }
0174       }
0175     }
0176
0177     if (not method_found) {
0178       throw std::runtime_error("Cannot determine the clustering coordinates");
0179     }
0180
0181     if (m_cfg.splitCluster) {
0182       if (!m_cfg.peakNeighbourhoodMatrix.empty()) {
0183         is_maximum_neighbourhood = serviceSvc.service<EvaluatorSvc>("EvaluatorSvc")->compile(m_cfg.peakNeighbourhoodMatrix, hit_pair_to_map);
0184       } else {
0185         is_maximum_neighbourhood = is_neighbour;
0186       }
0187
0188       auto transverseEnergyProfileMetric_it = std::find_if(distMethods.begin(), distMethods.end(), [&](auto &p) { return m_cfg.transverseEnergyProfileMetric == p.first; });
0189       if (transverseEnergyProfileMetric_it == distMethods.end()) {
0190           throw std::runtime_error(fmt::format("Unsupported value \"{}\" for \"transverseEnergyProfileMetric\"", m_cfg.transverseEnergyProfileMetric));
0191       }
0192       transverseEnergyProfileMetric = std::get<0>(transverseEnergyProfileMetric_it->second);
0193       std::vector<double> &units = std::get<1>(transverseEnergyProfileMetric_it->second);
0194       for (auto unit : units) {
0195         if (unit != units[0]) {
0196           throw std::runtime_error(fmt::format("Metric {} has incompatible dimension units", m_cfg.transverseEnergyProfileMetric));
0197         }
0198       }
0199       transverseEnergyProfileScaleUnits = units[0];
0200     }
0201
0202     return;
0203 }
0204
0205
0206 void CalorimeterIslandCluster::process(
0207       const CalorimeterIslandCluster::Input& input,
0208       const CalorimeterIslandCluster::Output& output) const {
0209
0210     const auto [hits] = input;
0211     auto [proto_clusters] = output;
0212
0213     // group neighboring hits
0214     std::vector<std::set<std::size_t>> groups;
0215
0216     std::vector<bool> visits(hits->size(), false);
0217     for (size_t i = 0; i < hits->size(); ++i) {
0218
0219       {
0220         const auto& hit = (*hits)[i];
0221         debug("hit {:d}: energy = {:.4f} MeV, local = ({:.4f}, {:.4f}) mm, global=({:.4f}, {:.4f}, {:.4f}) mm", i, hit.getEnergy() * 1000., hit.getLocal().x, hit.getLocal().y, hit.getPosition().x,  hit.getPosition().y, hit.getPosition().z);
0222       }
0223       // already in a group
0224       if (visits[i]) {
0225         continue;
0226       }
0227       groups.emplace_back();
0228       // create a new group, and group all the neighboring hits
0229       bfs_group(*hits, groups.back(), i, visits);
0230     }
0231
0232     for (auto& group : groups) {
0233       if (group.empty()) {
0234         continue;
0235       }
0236       auto maxima = find_maxima(*hits, group, !m_cfg.splitCluster);
0237       split_group(*hits, group, maxima, proto_clusters);
0238
0239       debug("hits in a group: {}, local maxima: {}", group.size(), maxima.size());
0240     }
0241 }
0242
0243 } // namespace eicrecon