Acts/Seeding/GbtsDataStorage.hpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #pragma once
0010
0011 // TODO: update to C++17 style
0012 #include "Acts/Definitions/Algebra.hpp"
0013 #include "Acts/Seeding/GbtsGeometry.hpp"
0014
0015 #include <algorithm>
0016 #include <iostream>
0017 #include <map>
0018 #include <numbers>
0019 #include <vector>
0020
0021 namespace Acts {
0022
0023 constexpr std::size_t MAX_SEG_PER_NODE = 1000;  // was 30
0024 constexpr std::size_t N_SEG_CONNS = 6;          // was 6
0025
0026 // new sp struct
0027 template <typename space_point_t>
0028 struct GbtsSP {
0029   const space_point_t *SP;  // want inside to have pointer
0030   int gbtsID;
0031   int combined_ID;
0032   GbtsSP(const space_point_t *sp, int id, int combined_id)
0033       : SP(sp), gbtsID(id), combined_ID{combined_id} {
0034     if (SP->sourceLinks().size() == 1) {  // pixels have 1 SL
0035       m_isPixel = true;
0036     } else {
0037       m_isPixel = false;
0038     }
0039     m_phi = std::atan(SP->x() / SP->y());
0040   };
0041   bool isPixel() const { return m_isPixel; }
0042   bool isSCT() const { return !m_isPixel; }
0043   float phi() const { return m_phi; }
0044   bool m_isPixel;
0045   float m_phi;
0046 };
0047
0048 template <typename space_point_t>
0049 class GbtsNode {
0050  public:
0051   GbtsNode(const GbtsSP<space_point_t> &spGbts, float minT = -100.0,
0052            float maxT = 100.0)
0053       : m_spGbts(spGbts), m_minCutOnTau(minT), m_maxCutOnTau(maxT) {}
0054
0055   inline void addIn(int i) {
0056     if (m_in.size() < MAX_SEG_PER_NODE) {
0057       m_in.push_back(i);
0058     }
0059   }
0060
0061   inline void addOut(int i) {
0062     if (m_out.size() < MAX_SEG_PER_NODE) {
0063       m_out.push_back(i);
0064     }
0065   }
0066
0067   inline bool isConnector() const {
0068     if (m_in.empty() || m_out.empty()) {
0069       return false;
0070     }
0071     return true;
0072   }
0073
0074   inline bool isFull() const {
0075     if (m_in.size() == MAX_SEG_PER_NODE && m_out.size() == MAX_SEG_PER_NODE) {
0076       return true;
0077     } else {
0078       return false;
0079     }
0080   }
0081
0082   const GbtsSP<space_point_t> &m_spGbts;
0083
0084   std::vector<unsigned int> m_in;  // indices of the edges in the edge storage
0085   std::vector<unsigned int> m_out;
0086   float m_minCutOnTau, m_maxCutOnTau;
0087 };
0088
0089 template <typename space_point_t>
0090 class GbtsEtaBin {
0091  public:
0092   GbtsEtaBin() { m_vn.clear(); }
0093
0094   void sortByPhi() {
0095     std::ranges::sort(m_vn, [](const auto &n1, const auto &n2) {
0096       return (n1->m_spGbts.phi() < n2->m_spGbts.phi());
0097     });
0098   }
0099
0100   bool empty() const { return m_vn.empty(); }
0101
0102   void generatePhiIndexing(float dphi) {
0103     for (unsigned int nIdx = 0; nIdx < m_vn.size(); nIdx++) {
0104       GbtsNode<space_point_t> &pN = *m_vn.at(nIdx);
0105       // float phi = pN->m_sp.phi();
0106       // float phi = (std::atan(pN->m_sp.x() / pN->m_sp.y()));
0107       float phi = pN.m_spGbts.phi();
0108       if (phi <= std::numbers::pi_v<float> - dphi) {
0109         continue;
0110       }
0111
0112       m_vPhiNodes.push_back(std::pair<float, unsigned int>(
0113           phi - static_cast<float>(2. * std::numbers::pi), nIdx));
0114     }
0115
0116     for (unsigned int nIdx = 0; nIdx < m_vn.size(); nIdx++) {
0117       GbtsNode<space_point_t> &pN = *m_vn.at(nIdx);
0118       float phi = pN.m_spGbts.phi();
0119       m_vPhiNodes.push_back(std::pair<float, unsigned int>(phi, nIdx));
0120     }
0121
0122     for (unsigned int nIdx = 0; nIdx < m_vn.size(); nIdx++) {
0123       GbtsNode<space_point_t> &pN = *m_vn.at(nIdx);
0124       float phi = pN.m_spGbts.phi();
0125       if (phi >= -std::numbers::pi_v<float> + dphi) {
0126         break;
0127       }
0128       m_vPhiNodes.push_back(std::pair<float, unsigned int>(
0129           phi + static_cast<float>(2. * std::numbers::pi), nIdx));
0130     }
0131   }
0132
0133   std::vector<std::unique_ptr<GbtsNode<space_point_t>>> m_vn;
0134   std::vector<std::pair<float, unsigned int>> m_vPhiNodes;
0135 };
0136
0137 template <typename space_point_t>
0138 class GbtsDataStorage {
0139  public:
0140   GbtsDataStorage(const GbtsGeometry<space_point_t> &g) : m_geo(g) {
0141     m_etaBins.reserve(g.num_bins());
0142     for (int k = 0; k < g.num_bins(); k++) {
0143       m_etaBins.emplace_back(GbtsEtaBin<space_point_t>());
0144     }
0145   }
0146
0147   int addSpacePoint(const GbtsSP<space_point_t> &sp, bool useClusterWidth) {
0148     const GbtsLayer<space_point_t> *pL =
0149         m_geo.getGbtsLayerByKey(sp.combined_ID);
0150
0151     if (pL == nullptr) {
0152       return -1;
0153     }
0154
0155     int binIndex = pL->getEtaBin(sp.SP->z(), sp.SP->r());
0156
0157     if (binIndex == -1) {
0158       return -2;
0159     }
0160
0161     bool isBarrel = (pL->m_layer.m_type == 0);
0162
0163     if (isBarrel) {
0164       float min_tau = -100.0;
0165       float max_tau = 100.0;
0166       // can't do this bit yet as dont have cluster width
0167       if (useClusterWidth) {
0168         float cluster_width = 1;  // temporary while cluster width not available
0169         min_tau = 6.7 * (cluster_width - 0.2);
0170         max_tau =
0171             1.6 + 0.15 / (cluster_width + 0.2) + 6.1 * (cluster_width - 0.2);
0172       }
0173
0174       m_etaBins.at(binIndex).m_vn.push_back(
0175           std::make_unique<GbtsNode<space_point_t>>(
0176               sp, min_tau, max_tau));  // adding ftf member to nodes
0177     } else {
0178       if (useClusterWidth) {
0179         float cluster_width = 1;  // temporary while cluster width not available
0180         if (cluster_width > 0.2) {
0181           return -3;
0182         }
0183       }
0184       m_etaBins.at(binIndex).m_vn.push_back(
0185           std::make_unique<GbtsNode<space_point_t>>(sp));
0186     }
0187
0188     return 0;
0189   }
0190
0191   // for safety to prevent passing as copy
0192   GbtsDataStorage(const GbtsDataStorage &) = delete;
0193   GbtsDataStorage &operator=(const GbtsDataStorage &) = delete;
0194
0195   unsigned int numberOfNodes() const {
0196     unsigned int n = 0;
0197
0198     for (auto &b : m_etaBins) {
0199       n += b.m_vn.size();
0200     }
0201     return n;
0202   }
0203
0204   void getConnectingNodes(std::vector<const GbtsNode<space_point_t> *> &vn) {
0205     vn.clear();
0206     vn.reserve(numberOfNodes());
0207     for (const auto &b : m_etaBins) {
0208       for (const auto &n : b.m_vn) {
0209         if (n->m_in.empty()) {
0210           continue;
0211         }
0212         if (n->m_out.empty()) {
0213           continue;
0214         }
0215         vn.push_back(n.get());
0216       }
0217     }
0218   }
0219
0220   void sortByPhi() {
0221     for (auto &b : m_etaBins) {
0222       b.sortByPhi();
0223     }
0224   }
0225
0226   void generatePhiIndexing(float dphi) {
0227     for (auto &b : m_etaBins) {
0228       b.generatePhiIndexing(dphi);
0229     }
0230   }
0231
0232   const GbtsEtaBin<space_point_t> &getEtaBin(int idx) const {
0233     if (idx >= static_cast<int>(m_etaBins.size())) {
0234       idx = idx - 1;
0235     }
0236     return m_etaBins.at(idx);
0237   }
0238
0239  protected:
0240   const GbtsGeometry<space_point_t> &m_geo;
0241
0242   std::vector<GbtsEtaBin<space_point_t>> m_etaBins;
0243 };
0244
0245 template <typename space_point_t>
0246 class GbtsEdge {
0247  public:
0248   struct CompareLevel {
0249    public:
0250     bool operator()(const GbtsEdge *pS1, const GbtsEdge *pS2) {
0251       return pS1->m_level > pS2->m_level;
0252     }
0253   };
0254
0255   GbtsEdge(GbtsNode<space_point_t> *n1, GbtsNode<space_point_t> *n2, float p1,
0256            float p2, float p3, float p4)
0257       : m_n1(n1), m_n2(n2), m_level(1), m_next(1) {
0258     m_p[0] = p1;
0259     m_p[1] = p2;
0260     m_p[2] = p3;
0261     m_p[3] = p4;
0262   }
0263
0264   GbtsEdge() : m_n1(nullptr), m_n2(nullptr), m_level(-1), m_next(-1) {}
0265
0266   // GbtsEdge(const GbtsEdge<space_point_t> &e)
0267   //     : m_n1(e.m_n1), m_n2(e.m_n2) {}
0268
0269   // inline void initialize(GbtsNode<space_point_t> *n1,
0270   //                        GbtsNode<space_point_t> *n2) {
0271   //   m_n1 = n1;
0272   //   m_n2 = n2;
0273   //   m_level = 1;
0274   //   m_next = 1;
0275   //   m_nNei = 0;
0276   // }
0277
0278   GbtsNode<space_point_t> *m_n1{nullptr};
0279   GbtsNode<space_point_t> *m_n2{nullptr};
0280
0281   signed char m_level{}, m_next{};
0282
0283   unsigned char m_nNei{0};
0284   float m_p[4]{};
0285
0286   unsigned int m_vNei[N_SEG_CONNS]{};  // global indices of the connected edges
0287 };
0288
0289 }  // namespace Acts