Back to home page

EIC code displayed by LXR

 
 

    


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 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2021 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 http://mozilla.org/MPL/2.0/.
0008 
0009 #pragma once
0010 
0011 // TODO: update to C++17 style
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;  // combined ID
0024   int m_type;    // 0: barrel, +/-n : endcap
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) {  // barrel
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 {  // endcap
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;  // increasing them slightly to avoid range_check
0069                        // exceptions
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) {  // barrel
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 {  // endcap
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) {  // barrel
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 {  // endcap
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) {  // barrel
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 {  // endcap
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);  // index in the global storage
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) {  // barrel -> barrel
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) {  // barrel -> endcap
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;  // eta-bin indices
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     // calculating bin tables in the connector...
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;  // n2 : the new connectors
0288         unsigned int dst = (*cIt)->m_dst;  // n1
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++) {    // loop over bins in Layer 1
0307           for (int b2 = 0; b2 < nSrcBins; b2++) {  // loop over bins in Layer 2
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   // for safety to prevent passing as copy
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;  // this should be combined ID
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 }  // namespace Acts