Warning, file /include/Acts/Seeding/GbtsGeometry.hpp was not indexed
or was modified since last indexation (in which case cross-reference links may be missing, inaccurate or erroneous).
0001
0002
0003
0004
0005
0006
0007
0008
0009 #pragma once
0010
0011
0012 #include "Acts/TrackFinding/GbtsConnector.hpp"
0013
0014 #include <algorithm>
0015 #include <cmath>
0016 #include <map>
0017 #include <memory>
0018 #include <vector>
0019
0020 namespace Acts {
0021 class TrigInDetSiLayer {
0022 public:
0023 int m_subdet;
0024 int m_type;
0025 float m_refCoord;
0026 float m_minBound, m_maxBound;
0027
0028 TrigInDetSiLayer(int subdet, short int type, float center, float min,
0029 float max)
0030 : m_subdet(subdet),
0031 m_type(type),
0032 m_refCoord(center),
0033 m_minBound(min),
0034 m_maxBound(max) {}
0035 };
0036
0037 template <typename space_point_t>
0038 class GbtsLayer {
0039 public:
0040 GbtsLayer(const TrigInDetSiLayer &ls, float ew, int bin0)
0041 : m_layer(ls), m_etaBinWidth(ew) {
0042 if (m_layer.m_type == 0) {
0043 m_r1 = m_layer.m_refCoord;
0044 m_r2 = m_layer.m_refCoord;
0045 m_z1 = m_layer.m_minBound;
0046 m_z2 = m_layer.m_maxBound;
0047 } else {
0048 m_r1 = m_layer.m_minBound;
0049 m_r2 = m_layer.m_maxBound;
0050 m_z1 = m_layer.m_refCoord;
0051 m_z2 = m_layer.m_refCoord;
0052 }
0053
0054 float t1 = m_z1 / m_r1;
0055 float eta1 = -std::log(sqrt(1 + t1 * t1) - t1);
0056
0057 float t2 = m_z2 / m_r2;
0058 float eta2 = -std::log(sqrt(1 + t2 * t2) - t2);
0059
0060 m_minEta = eta1;
0061 m_maxEta = eta2;
0062
0063 if (m_maxEta < m_minEta) {
0064 m_minEta = eta2;
0065 m_maxEta = eta1;
0066 }
0067
0068 m_maxEta += 1e-6;
0069
0070 m_minEta -= 1e-6;
0071
0072 float deltaEta = m_maxEta - m_minEta;
0073
0074 int binCounter = bin0;
0075
0076 if (deltaEta < m_etaBinWidth) {
0077 m_nBins = 1;
0078 m_bins.push_back(binCounter++);
0079 m_etaBin = deltaEta;
0080 if (m_layer.m_type == 0) {
0081 m_minRadius.push_back(m_layer.m_refCoord - 2.0);
0082 m_maxRadius.push_back(m_layer.m_refCoord + 2.0);
0083 m_minBinCoord.push_back(m_layer.m_minBound);
0084 m_maxBinCoord.push_back(m_layer.m_maxBound);
0085 } else {
0086 m_minRadius.push_back(m_layer.m_minBound - 2.0);
0087 m_maxRadius.push_back(m_layer.m_maxBound + 2.0);
0088 m_minBinCoord.push_back(m_layer.m_minBound);
0089 m_maxBinCoord.push_back(m_layer.m_maxBound);
0090 }
0091 } else {
0092 float nB = static_cast<int>(deltaEta / m_etaBinWidth);
0093 m_nBins = nB;
0094 if (deltaEta - m_etaBinWidth * nB > 0.5 * m_etaBinWidth) {
0095 m_nBins++;
0096 }
0097 m_etaBin = deltaEta / m_nBins;
0098
0099 if (m_nBins == 1) {
0100 m_bins.push_back(binCounter++);
0101 if (m_layer.m_type == 0) {
0102 m_minRadius.push_back(m_layer.m_refCoord - 2.0);
0103 m_maxRadius.push_back(m_layer.m_refCoord + 2.0);
0104 m_minBinCoord.push_back(m_layer.m_minBound);
0105 m_maxBinCoord.push_back(m_layer.m_maxBound);
0106 } else {
0107 m_minRadius.push_back(m_layer.m_minBound - 2.0);
0108 m_maxRadius.push_back(m_layer.m_maxBound + 2.0);
0109 m_minBinCoord.push_back(m_layer.m_minBound);
0110 m_maxBinCoord.push_back(m_layer.m_maxBound);
0111 }
0112 } else {
0113 float eta = m_minEta + 0.5 * m_etaBin;
0114
0115 for (int i = 1; i <= m_nBins; i++) {
0116 m_bins.push_back(binCounter++);
0117
0118 float e1 = eta - 0.5 * m_etaBin;
0119 float e2 = eta + 0.5 * m_etaBin;
0120
0121 if (m_layer.m_type == 0) {
0122 m_minRadius.push_back(m_layer.m_refCoord - 2.0);
0123 m_maxRadius.push_back(m_layer.m_refCoord + 2.0);
0124 float z1 = m_layer.m_refCoord * std::sinh(e1);
0125 m_minBinCoord.push_back(z1);
0126 float z2 = m_layer.m_refCoord * std::sinh(e2);
0127 m_maxBinCoord.push_back(z2);
0128 } else {
0129 float r = m_layer.m_refCoord / std::sinh(e1);
0130 m_minBinCoord.push_back(r);
0131 m_minRadius.push_back(r - 2.0);
0132 r = m_layer.m_refCoord / std::sinh(e2);
0133 m_maxBinCoord.push_back(r);
0134 m_maxRadius.push_back(r + 2.0);
0135 }
0136
0137 eta += m_etaBin;
0138 }
0139 }
0140 }
0141 }
0142
0143 int getEtaBin(float zh, float rh) const {
0144 if (m_bins.size() == 1) {
0145 return m_bins.at(0);
0146 }
0147 float t1 = zh / rh;
0148 float eta = -std::log(std::sqrt(1 + t1 * t1) - t1);
0149
0150 int idx = static_cast<int>((eta - m_minEta) / m_etaBin);
0151
0152 if (idx < 0) {
0153 idx = 0;
0154 }
0155 if (idx >= static_cast<int>(m_bins.size())) {
0156 idx = static_cast<int>(m_bins.size()) - 1;
0157 }
0158 return m_bins.at(idx);
0159 }
0160
0161 float getMinBinRadius(int idx) const {
0162 if (idx >= static_cast<int>(m_minRadius.size())) {
0163 idx = idx - 1;
0164 }
0165 if (idx < 0) {
0166 idx = 0;
0167 }
0168 return m_minRadius.at(idx);
0169 }
0170
0171 float getMaxBinRadius(int idx) const {
0172 if (idx >= static_cast<int>(m_maxRadius.size())) {
0173 idx = idx - 1;
0174 }
0175 if (idx < 0) {
0176 idx = 0;
0177 }
0178 return m_maxRadius.at(idx);
0179 }
0180
0181 int num_bins() const { return m_bins.size(); }
0182
0183 bool verifyBin(const GbtsLayer<space_point_t> *pL, int b1, int b2,
0184 float min_z0, float max_z0) const {
0185 float z1min = m_minBinCoord.at(b1);
0186 float z1max = m_maxBinCoord.at(b1);
0187 float r1 = m_layer.m_refCoord;
0188
0189 if (m_layer.m_type == 0 && pL->m_layer.m_type == 0) {
0190
0191 const float tol = 5.0;
0192
0193 float min_b2 = pL->m_minBinCoord.at(b2);
0194 float max_b2 = pL->m_maxBinCoord.at(b2);
0195
0196 float r2 = pL->m_layer.m_refCoord;
0197
0198 float A = r2 / (r2 - r1);
0199 float B = r1 / (r2 - r1);
0200
0201 float z0_min = z1min * A - max_b2 * B;
0202 float z0_max = z1max * A - min_b2 * B;
0203
0204 if (z0_max < min_z0 - tol || z0_min > max_z0 + tol) {
0205 return false;
0206 }
0207 return true;
0208 }
0209
0210 if (m_layer.m_type == 0 && pL->m_layer.m_type != 0) {
0211
0212 const float tol = 10.0;
0213
0214 float z2 = pL->m_layer.m_refCoord;
0215 float r2max = pL->m_maxBinCoord.at(b2);
0216 float r2min = pL->m_minBinCoord.at(b2);
0217
0218 if (r2max <= r1) {
0219 return false;
0220 }
0221 if (r2min <= r1) {
0222 r2min = r1 + 1e-3;
0223 }
0224
0225 float z0_max = 0.0;
0226 float z0_min = 0.0;
0227
0228 if (z2 > 0) {
0229 z0_max = (z1max * r2max - z2 * r1) / (r2max - r1);
0230 z0_min = (z1min * r2min - z2 * r1) / (r2min - r1);
0231 } else {
0232 z0_max = (z1max * r2min - z2 * r1) / (r2min - r1);
0233 z0_min = (z1min * r2max - z2 * r1) / (r2max - r1);
0234 }
0235
0236 if (z0_max < min_z0 - tol || z0_min > max_z0 + tol) {
0237 return false;
0238 }
0239 return true;
0240 }
0241
0242 return true;
0243 }
0244
0245 const TrigInDetSiLayer &m_layer;
0246 std::vector<int> m_bins;
0247 std::vector<float> m_minRadius;
0248 std::vector<float> m_maxRadius;
0249 std::vector<float> m_minBinCoord;
0250 std::vector<float> m_maxBinCoord;
0251
0252 float m_minEta{}, m_maxEta{};
0253
0254 protected:
0255 float m_etaBinWidth{}, m_phiBinWidth{};
0256
0257 float m_r1{}, m_z1{}, m_r2{}, m_z2{};
0258 float m_nBins{};
0259 float m_etaBin{};
0260 };
0261
0262 template <typename space_point_t>
0263 class GbtsGeometry {
0264 public:
0265 GbtsGeometry(const std::vector<TrigInDetSiLayer> &layers,
0266 std::unique_ptr<Acts::GbtsConnector> &conn)
0267
0268 : m_connector(std::move(conn)) {
0269 const float min_z0 = -168.0;
0270 const float max_z0 = 168.0;
0271
0272 m_etaBinWidth = m_connector->m_etaBin;
0273 for (const auto &layer : layers) {
0274 const GbtsLayer<space_point_t> *pL = addNewLayer(layer, m_nEtaBins);
0275 m_nEtaBins += pL->num_bins();
0276 }
0277
0278
0279
0280 for (std::map<int, std::vector<GbtsConnection *>>::const_iterator it =
0281 m_connector->m_connMap.begin();
0282 it != m_connector->m_connMap.end(); ++it) {
0283 const std::vector<GbtsConnection *> &vConn = (*it).second;
0284
0285 for (std::vector<GbtsConnection *>::const_iterator cIt = vConn.begin();
0286 cIt != vConn.end(); ++cIt) {
0287 unsigned int src = (*cIt)->m_src;
0288 unsigned int dst = (*cIt)->m_dst;
0289
0290 const GbtsLayer<space_point_t> *pL1 = getGbtsLayerByKey(dst);
0291 const GbtsLayer<space_point_t> *pL2 = getGbtsLayerByKey(src);
0292
0293 if (pL1 == nullptr) {
0294 std::cout << " skipping invalid dst layer " << dst << std::endl;
0295 continue;
0296 }
0297 if (pL2 == nullptr) {
0298 std::cout << " skipping invalid src layer " << src << std::endl;
0299 continue;
0300 }
0301 int nSrcBins = pL2->m_bins.size();
0302 int nDstBins = pL1->m_bins.size();
0303
0304 (*cIt)->m_binTable.resize(nSrcBins * nDstBins, 0);
0305
0306 for (int b1 = 0; b1 < nDstBins; b1++) {
0307 for (int b2 = 0; b2 < nSrcBins; b2++) {
0308 if (!pL1->verifyBin(pL2, b1, b2, min_z0, max_z0)) {
0309 continue;
0310 }
0311 int address = b1 + b2 * nDstBins;
0312 (*cIt)->m_binTable.at(address) = 1;
0313 }
0314 }
0315 }
0316 }
0317 }
0318
0319 GbtsGeometry() = default;
0320
0321
0322 GbtsGeometry(const GbtsGeometry &) = delete;
0323 GbtsGeometry &operator=(const GbtsGeometry &) = delete;
0324
0325 ~GbtsGeometry() {
0326 for (typename std::vector<GbtsLayer<space_point_t> *>::iterator it =
0327 m_layArray.begin();
0328 it != m_layArray.end(); ++it) {
0329 delete (*it);
0330 }
0331
0332 m_layMap.clear();
0333 m_layArray.clear();
0334 }
0335
0336 const GbtsLayer<space_point_t> *getGbtsLayerByKey(unsigned int key) const {
0337 typename std::map<unsigned int, GbtsLayer<space_point_t> *>::const_iterator
0338 it = m_layMap.find(key);
0339 if (it == m_layMap.end()) {
0340 return nullptr;
0341 }
0342
0343 return (*it).second;
0344 }
0345
0346 const GbtsLayer<space_point_t> *getGbtsLayerByIndex(int idx) const {
0347 return m_layArray.at(idx);
0348 }
0349
0350 int num_bins() const { return m_nEtaBins; }
0351
0352 Acts::GbtsConnector *connector() const { return m_connector.get(); }
0353
0354 protected:
0355 const GbtsLayer<space_point_t> *addNewLayer(const TrigInDetSiLayer &l,
0356 int bin0) {
0357 unsigned int layerKey = l.m_subdet;
0358 float ew = m_etaBinWidth;
0359
0360 GbtsLayer<space_point_t> *pHL = new GbtsLayer<space_point_t>(l, ew, bin0);
0361
0362 m_layMap.insert(
0363 std::pair<unsigned int, GbtsLayer<space_point_t> *>(layerKey, pHL));
0364 m_layArray.push_back(pHL);
0365 return pHL;
0366 }
0367
0368 float m_etaBinWidth{};
0369
0370 std::map<unsigned int, GbtsLayer<space_point_t> *> m_layMap;
0371 std::vector<GbtsLayer<space_point_t> *> m_layArray;
0372
0373 int m_nEtaBins{0};
0374
0375 std::unique_ptr<Acts::GbtsConnector> m_connector;
0376 };
0377
0378 }