Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-07-26 08:18:59

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 #include "Acts/Seeding/GbtsTrainingTool.hpp"
0010 
0011 #include <cmath>
0012 #include <fstream>
0013 #include <iostream>
0014 #include <stdexcept>
0015 #include <string>
0016 
0017 namespace Acts::Experimental {
0018 
0019 GbtsLayerConnectionTool::LayerDescription::LayerDescription(
0020     float minR_, float maxR_, float minZ_, float maxZ_, std::int32_t gbtsId_)
0021     : minR(minR_), maxR(maxR_), minZ(minZ_), maxZ(maxZ_), gbtsId(gbtsId_) {}
0022 
0023 GbtsLayerConnectionTool::GbtsLayerConnectionTool(
0024     const Config& config, std::unique_ptr<const Logger> logger)
0025     : m_cfg(config), m_logger(std::move(logger)) {
0026   if (m_cfg.detectorGeometry.empty()) {
0027     throw std::runtime_error("File does not exist or could not be opened");
0028   }
0029 
0030   // reserves
0031   const std::uint32_t pairReserve = m_cfg.detectorGeometry.size();
0032   m_layerPairs.reserve(pairReserve * (pairReserve - 1));
0033 
0034   // create map linked pairs of GBTS ids with number of transitions between
0035   // layers
0036   for (std::uint32_t i = 0; i < m_cfg.detectorGeometry.size(); i++) {
0037     for (std::uint32_t j = 0; j < m_cfg.detectorGeometry.size(); j++) {
0038       if (i == j) {
0039         continue;
0040       }
0041 
0042       LayerIdPair pair;
0043       pair.first = m_cfg.detectorGeometry[i].gbtsId;
0044       pair.second = m_cfg.detectorGeometry[j].gbtsId;
0045 
0046       // key = GBTS ids of pair, value = number of transitions
0047       m_layerPairs.emplace(pair, 0);
0048     }
0049   }
0050 }
0051 
0052 void GbtsLayerConnectionTool::addTrack(
0053     const std::vector<HitCoordinates>& track) {
0054   if (track.size() < 2) {
0055     ACTS_WARNING("Track only has one measurement, skipping");
0056     return;
0057   }
0058 
0059   // container for gbts IDs of the track
0060   std::vector<std::optional<std::int32_t>> layerGbtsIds{};
0061   layerGbtsIds.reserve(track.size());
0062 
0063   // find GBTS ids for all measurements in a track
0064   for (const auto& measurement : track) {
0065     const auto gbtsId = findGbtsIdByCoord(measurement);
0066     if (!gbtsId) {
0067       ACTS_WARNING("No Gbts Layer for coordinates with r: "
0068                    << measurement.r << " and z: " << measurement.z);
0069     }
0070     layerGbtsIds.emplace_back(gbtsId);
0071   }
0072 
0073   // update map with track layer transitions
0074   for (std::uint32_t id = 0; id + 1 < layerGbtsIds.size(); id++) {
0075     const auto& index1 = layerGbtsIds[id];
0076     const auto& index2 = layerGbtsIds[id + 1];
0077 
0078     // skip nonexistent layers ids
0079     if (!index1 || !index2) {
0080       continue;
0081     }
0082 
0083     if (index1.value() == index2.value()) {
0084       ACTS_WARNING("Track transitions between same layer, skipping");
0085 
0086       continue;
0087     }
0088 
0089     m_layerPairs[{index1.value(), index2.value()}] += 1;
0090   }
0091 
0092   m_totalTracks++;
0093 }
0094 
0095 GbtsLayerConnectionTool::LayerIdPairs
0096 GbtsLayerConnectionTool::createConnectionTable(
0097     const std::string& outputFileLocation) const {
0098   if (m_totalTracks == 0) {
0099     throw std::runtime_error(
0100         "Warning: no tracks were added when creating connection table");
0101   }
0102 
0103   // define output text file
0104   std::ofstream outputFile(outputFileLocation);
0105 
0106   // obtain total incoming transitions for each src layer (used as denominator
0107   // of probability)
0108   std::vector<std::uint32_t> srcTotals;
0109   srcTotals.resize(m_cfg.detectorGeometry.size(), 0);
0110 
0111   for (const auto& [layerPair, nTransitions] : m_layerPairs) {
0112     std::uint32_t layerIndex = getIndexByGbtsId(layerPair.first);
0113     srcTotals[layerIndex] += nTransitions;
0114   }
0115 
0116   // find transitions that pass probability cut and add to temp container
0117   LayerIdPairs tempPairs;
0118   for (const auto& [layerPair, nTransitions] : m_layerPairs) {
0119     const std::uint32_t srcIndex = getIndexByGbtsId(layerPair.first);
0120 
0121     float probability{};
0122 
0123     if (srcTotals[srcIndex] == 0) {
0124       // avoid NAN error with 0/0
0125       probability = 0;
0126     } else {
0127       probability = static_cast<float>(nTransitions) / srcTotals[srcIndex];
0128     }
0129 
0130     const bool passCut = (m_cfg.probThreshold == -1)
0131                              ? (probability != 0)
0132                              : (probability >= m_cfg.probThreshold);
0133 
0134     if (passCut) {
0135       tempPairs.emplace(layerPair.first, layerPair.second);
0136     }
0137   }
0138 
0139   // if symmetrizing connection table, add transitions that mirror found ones
0140   if (m_cfg.doSymmetrization) {
0141     for (const auto& layerPair : tempPairs) {
0142       // find swapped ids
0143       const auto srcSwappedId = oppositeSideLayer(layerPair.first);
0144       const auto dstSwappedId = oppositeSideLayer(layerPair.second);
0145 
0146       if (!srcSwappedId || !dstSwappedId) {
0147         ACTS_WARNING("Cannot find oppisite side layer, skipping");
0148         continue;
0149       }
0150       // search set to see if swapped pair has already been added
0151       const bool notAdded =
0152           (tempPairs.count({srcSwappedId.value(), dstSwappedId.value()}) == 0);
0153 
0154       // if not already added, add to output file
0155       if (notAdded) {
0156         tempPairs.emplace(srcSwappedId.value(), dstSwappedId.value());
0157       }
0158     }
0159   }
0160 
0161   return tempPairs;
0162 }
0163 
0164 std::optional<std::int32_t> GbtsLayerConnectionTool::findGbtsIdByCoord(
0165     const HitCoordinates& hit) const {
0166   for (const auto& layer : m_cfg.detectorGeometry) {
0167     const float zMin = layer.minZ - m_cfg.zMinTol;
0168     const float zMax = layer.maxZ + m_cfg.zMaxTol;
0169     const float rMin = layer.minR - m_cfg.rMinTol;
0170     const float rMax = layer.maxR + m_cfg.rMaxTol;
0171 
0172     if (zMin <= hit.z && hit.z <= zMax) {
0173       if (rMin <= hit.r && hit.r <= rMax) {
0174         return layer.gbtsId;
0175       }
0176     }
0177   }
0178 
0179   return std::nullopt;
0180 }
0181 
0182 std::uint32_t GbtsLayerConnectionTool::getIndexByGbtsId(
0183     std::int32_t gbtsId) const {
0184   for (std::uint32_t idx = 0; idx < m_cfg.detectorGeometry.size(); idx++) {
0185     if (gbtsId == m_cfg.detectorGeometry[idx].gbtsId) {
0186       return idx;
0187     }
0188   }
0189 
0190   throw std::runtime_error("index not found for GBTS ID");
0191 }
0192 
0193 std::optional<std::int32_t> GbtsLayerConnectionTool::oppositeSideLayer(
0194     std::int32_t layerId) const {
0195   const std::uint32_t layerIndex = getIndexByGbtsId(layerId);
0196 
0197   const auto& layer = m_cfg.detectorGeometry[layerIndex];
0198 
0199   bool switchedMinZ{};
0200   bool switchedMaxZ{};
0201 
0202   bool sameMaxR{};
0203   bool sameMinR{};
0204 
0205   for (const auto& switchedLayer : m_cfg.detectorGeometry) {
0206     switchedMinZ = (switchedLayer.minZ == -layer.maxZ);
0207     switchedMaxZ = (switchedLayer.maxZ == -layer.minZ);
0208 
0209     sameMaxR = (switchedLayer.maxR == layer.maxR);
0210     sameMinR = (switchedLayer.minR == layer.minR);
0211 
0212     if (switchedMinZ && switchedMaxZ && sameMaxR && sameMinR) {
0213       return switchedLayer.gbtsId;
0214     }
0215   }
0216 
0217   return std::nullopt;
0218 }
0219 
0220 }  // namespace Acts::Experimental