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 <TInterpreter.h>
0011 #include <TInterpreterValue.h>
0012 #include <edm4hep/Vector2f.h>
0013 #include <edm4hep/Vector3f.h>
0014 #include <fmt/format.h>
0015 #include <algorithm>
0016 #include <gsl/pointers>
0017 #include <map>
0018 #include <memory>
0019 #include <set>
0020 #include <sstream>
0021 #include <stdexcept>
0022 #include <string>
0023 #include <tuple>
0024 #include <utility>
0025 #include <vector>
0026
0027 #include "CalorimeterIslandCluster.h"
0028 #include "algorithms/calorimetry/CalorimeterIslandClusterConfig.h"
0029
0030 using namespace edm4eic;
0031
0032 namespace eicrecon {
0033
0034 unsigned int CalorimeterIslandCluster::function_id = 0;
0035
0036 static double Phi_mpi_pi(double phi) {
0037   return std::remainder(phi, 2 * M_PI);
0038 }
0039
0040 static edm4hep::Vector2f localDistXY(const CaloHit &h1, const CaloHit &h2) {
0041   const auto delta =h1.getLocal() - h2.getLocal();
0042   return {delta.x, delta.y};
0043 }
0044 static edm4hep::Vector2f localDistXZ(const CaloHit &h1, const CaloHit &h2) {
0045   const auto delta = h1.getLocal() - h2.getLocal();
0046   return {delta.x, delta.z};
0047 }
0048 static edm4hep::Vector2f localDistYZ(const CaloHit &h1, const CaloHit &h2) {
0049   const auto delta = h1.getLocal() - h2.getLocal();
0050   return {delta.y, delta.z};
0051 }
0052 static edm4hep::Vector2f dimScaledLocalDistXY(const CaloHit &h1, const CaloHit &h2) {
0053   const auto delta = h1.getLocal() - h2.getLocal();
0054
0055   const auto dimsum = h1.getDimension() + h2.getDimension();
0056
0057   return {2 * delta.x / dimsum.x, 2 * delta.y / dimsum.y};
0058 }
0059 static edm4hep::Vector2f globalDistRPhi(const CaloHit &h1, const CaloHit &h2) {
0060   using vector_type = decltype(edm4hep::Vector2f::a);
0061   return {
0062     static_cast<vector_type>(
0063       edm4hep::utils::magnitude(h1.getPosition()) - edm4hep::utils::magnitude(h2.getPosition())
0064     ),
0065     static_cast<vector_type>(
0066       Phi_mpi_pi(edm4hep::utils::angleAzimuthal(h1.getPosition()) - edm4hep::utils::angleAzimuthal(h2.getPosition()))
0067     )
0068   };
0069 }
0070 static edm4hep::Vector2f globalDistEtaPhi(const CaloHit &h1, const CaloHit &h2) {
0071   using vector_type = decltype(edm4hep::Vector2f::a);
0072   return {
0073     static_cast<vector_type>(
0074       edm4hep::utils::eta(h1.getPosition()) - edm4hep::utils::eta(h2.getPosition())
0075     ),
0076     static_cast<vector_type>(
0077       Phi_mpi_pi(edm4hep::utils::angleAzimuthal(h1.getPosition()) - edm4hep::utils::angleAzimuthal(h2.getPosition()))
0078     )
0079   };
0080 }
0081
0082 //------------------------
0083 // AlgorithmInit
0084 //------------------------
0085 void CalorimeterIslandCluster::init() {
0086
0087     static std::map<std::string,
0088                 std::tuple<std::function<edm4hep::Vector2f(const CaloHit&, const CaloHit&)>, std::vector<double>>>
0089     distMethods{
0090         {"localDistXY", {localDistXY, {dd4hep::mm, dd4hep::mm}}},        {"localDistXZ", {localDistXZ, {dd4hep::mm, dd4hep::mm}}},
0091         {"localDistYZ", {localDistYZ, {dd4hep::mm, dd4hep::mm}}},        {"dimScaledLocalDistXY", {dimScaledLocalDistXY, {1., 1.}}},
0092         {"globalDistRPhi", {globalDistRPhi, {dd4hep::mm, dd4hep::rad}}}, {"globalDistEtaPhi", {globalDistEtaPhi, {1., dd4hep::rad}}}
0093     };
0094
0095
0096     // set coordinate system
0097     auto set_dist_method = [this](std::pair<std::string, std::vector<double>> uprop) {
0098       if (uprop.second.size() == 0) {
0099         return false;
0100       }
0101       auto& [method, units] = distMethods[uprop.first];
0102       if (uprop.second.size() != units.size()) {
0103         warning("Expect {} values from {}, received {}. ignored it.", units.size(), uprop.first,  uprop.second.size());
0104         return false;
0105       } else {
0106         for (size_t i = 0; i < units.size(); ++i) {
0107           neighbourDist[i] = uprop.second[i] / units[i];
0108         }
0109         hitsDist = method;
0110         info("Clustering uses {} with distances <= [{}]", uprop.first, fmt::join(neighbourDist, ","));
0111       }
0112       return true;
0113     };
0114
0115     std::map<std::string, std::vector<double>> uprops{
0116             {"localDistXY", m_cfg.localDistXY},
0117             {"localDistXZ", m_cfg.localDistXZ},
0118             {"localDistYZ", m_cfg.localDistYZ},
0119             {"globalDistRPhi", m_cfg.globalDistRPhi},
0120             {"globalDistEtaPhi", m_cfg.globalDistEtaPhi},
0121             // default one should be the last one
0122             {"dimScaledLocalDistXY", m_cfg.dimScaledLocalDistXY}
0123     };
0124
0125     bool method_found = false;
0126
0127     // Adjacency matrix methods
0128     if (!m_cfg.adjacencyMatrix.empty()) {
0129       // sanity checks
0130       if (m_cfg.readout.empty()) {
0131         error("readoutClass is not provided, it is needed to know the fields in readout ids");
0132       }
0133       m_idSpec = m_detector->readout(m_cfg.readout).idSpec();
0134
0135       std::string func_name = fmt::format("_CalorimeterIslandCluster_{}", function_id++);
0136       std::ostringstream sstr;
0137       sstr << "bool " << func_name << "(double params[]){";
0138       unsigned int param_ix = 0;
0139       for(const auto &p : m_idSpec.fields()) {
0140         const std::string &name = p.first;
0141         const dd4hep::IDDescriptor::Field* field = p.second;
0142         sstr << "double " << name << "_1 = params[" << (param_ix++) << "];";
0143         sstr << "double " << name << "_2 = params[" << (param_ix++) << "];";
0144       }
0145       sstr << "return " << m_cfg.adjacencyMatrix << ";";
0146       sstr << "}";
0147       debug("Compiling {}", sstr.str());
0148
0149       TInterpreter *interp = TInterpreter::Instance();
0150       interp->ProcessLine(sstr.str().c_str());
0151       std::unique_ptr<TInterpreterValue> func_val { gInterpreter->MakeInterpreterValue() };
0152       interp->Evaluate(func_name.c_str(), *func_val);
0153       typedef bool (*func_t)(double params[]);
0154       func_t func = ((func_t)(func_val->GetAsPointer()));
0155
0156       is_neighbour = [this, func, param_ix](const CaloHit &h1, const CaloHit &h2) {
0157         std::vector<double> params;
0158         params.reserve(param_ix);
0159         for(const auto &p : m_idSpec.fields()) {
0160           const std::string &name = p.first;
0161           const dd4hep::IDDescriptor::Field* field = p.second;
0162           params.push_back(field->value(h1.getCellID()));
0163           params.push_back(field->value(h2.getCellID()));
0164           trace("{}_1 = {}", name, field->value(h1.getCellID()));
0165           trace("{}_2 = {}", name, field->value(h2.getCellID()));
0166         }
0167         return func(params.data());
0168       };
0169       method_found = true;
0170     }
0171
0172     // Coordinate distance methods
0173     if (not method_found) {
0174       for (auto& uprop : uprops) {
0175         if (set_dist_method(uprop)) {
0176           method_found = true;
0177
0178           is_neighbour = [this](const CaloHit &h1, const CaloHit &h2) {
0179             // in the same sector
0180             if (h1.getSector() == h2.getSector()) {
0181               auto dist = hitsDist(h1, h2);
0182               return (fabs(dist.a) <= neighbourDist[0]) && (fabs(dist.b) <= neighbourDist[1]);
0183               // different sector, local coordinates do not work, using global coordinates
0184             } else {
0185               // sector may have rotation (barrel), so z is included
0186               // (EDM4hep units are mm, so convert sectorDist to mm)
0187               return (edm4hep::utils::magnitude(h1.getPosition() - h2.getPosition()) <= m_cfg.sectorDist / dd4hep::mm);
0188             }
0189           };
0190
0191           info("Using clustering method: {}", uprop.first);
0192           break;
0193         }
0194       }
0195     }
0196
0197     if (not method_found) {
0198       throw std::runtime_error("Cannot determine the clustering coordinates");
0199     }
0200
0201     if (m_cfg.splitCluster) {
0202       auto transverseEnergyProfileMetric_it = std::find_if(distMethods.begin(), distMethods.end(), [&](auto &p) { return m_cfg.transverseEnergyProfileMetric == p.first; });
0203       if (transverseEnergyProfileMetric_it == distMethods.end()) {
0204           throw std::runtime_error(fmt::format("Unsupported value \"{}\" for \"transverseEnergyProfileMetric\"", m_cfg.transverseEnergyProfileMetric));
0205       }
0206       transverseEnergyProfileMetric = std::get<0>(transverseEnergyProfileMetric_it->second);
0207       std::vector<double> &units = std::get<1>(transverseEnergyProfileMetric_it->second);
0208       for (auto unit : units) {
0209         if (unit != units[0]) {
0210           throw std::runtime_error(fmt::format("Metric {} has incompatible dimension units", m_cfg.transverseEnergyProfileMetric));
0211         }
0212       }
0213       transverseEnergyProfileScaleUnits = units[0];
0214     }
0215
0216     return;
0217 }
0218
0219
0220 void CalorimeterIslandCluster::process(
0221       const CalorimeterIslandCluster::Input& input,
0222       const CalorimeterIslandCluster::Output& output) const {
0223
0224     const auto [hits] = input;
0225     auto [proto_clusters] = output;
0226
0227     // group neighboring hits
0228     std::vector<std::set<std::size_t>> groups;
0229
0230     std::vector<bool> visits(hits->size(), false);
0231     for (size_t i = 0; i < hits->size(); ++i) {
0232
0233       {
0234         const auto& hit = (*hits)[i];
0235         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);
0236       }
0237       // already in a group
0238       if (visits[i]) {
0239         continue;
0240       }
0241       groups.emplace_back();
0242       // create a new group, and group all the neighboring hits
0243       bfs_group(*hits, groups.back(), i, visits);
0244     }
0245
0246     for (auto& group : groups) {
0247       if (group.empty()) {
0248         continue;
0249       }
0250       auto maxima = find_maxima(*hits, group, !m_cfg.splitCluster);
0251       split_group(*hits, group, maxima, proto_clusters);
0252
0253       debug("hits in a group: {}, local maxima: {}", group.size(), maxima.size());
0254     }
0255 }
0256
0257 } // namespace eicrecon