TrackFinding/src/GraphBasedSeedingAlgorithm.cpp

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 "ActsExamples/TrackFinding/GraphBasedSeedingAlgorithm.hpp"
0010
0011 #include "Acts/EventData/SpacePointColumns.hpp"
0012 #include "Acts/EventData/SpacePointContainer2.hpp"
0013 #include "Acts/Geometry/GeometryContext.hpp"
0014 #include "Acts/Geometry/GeometryIdentifier.hpp"
0015 #include "Acts/Seeding2/GbtsGeometry.hpp"
0016 #include "Acts/Seeding2/GbtsLayerConnection.hpp"
0017 #include "Acts/Seeding2/GbtsTrackingFilter.hpp"
0018 #include "ActsExamples/EventData/IndexSourceLink.hpp"
0019
0020 #include <cmath>
0021 #include <fstream>
0022 #include <iostream>
0023 #include <map>
0024 #include <numbers>
0025 #include <sstream>
0026 #include <vector>
0027
0028 namespace ActsExamples {
0029
0030 GraphBasedSeedingAlgorithm::GraphBasedSeedingAlgorithm(
0031     const Config &cfg, std::unique_ptr<const Acts::Logger> logger)
0032     : IAlgorithm("GraphBasedSeedingAlgorithm", std::move(logger)), m_cfg(cfg) {
0033   // initialise the space point, seed and cluster handles
0034   m_inputSpacePoints.initialize(m_cfg.inputSpacePoints);
0035   m_outputSeeds.initialize(m_cfg.outputSeeds);
0036   m_inputClusters.initialize(m_cfg.inputClusters);
0037
0038   // create the TrigInDetSiLayers (Logical Layers),
0039   // as well as a map that tracks there index in m_layerGeometry
0040   auto layerGeometry =
0041       layerNumbering(Acts::GeometryContext::dangerouslyDefaultConstruct());
0042
0043   // create the connection objects
0044   Acts::Experimental::GbtsLayerConnectionMap layerConnectionMap(
0045       m_cfg.seedFinderConfig.connectorInputFile,
0046       m_cfg.seedFinderConfig.lrtMode);
0047
0048   // option that allows for adding custom eta binning (default is at 0.2)
0049   if (m_cfg.seedFinderConfig.etaBinWidthOverride != 0.0f) {
0050     layerConnectionMap.etaBinWidth = m_cfg.seedFinderConfig.etaBinWidthOverride;
0051   }
0052
0053   // initialise the object that holds all the geometry information needed for
0054   // the algorithm
0055   auto geometry = std::make_shared<Acts::Experimental::GbtsGeometry>(
0056       layerGeometry, layerConnectionMap);
0057
0058   m_finder = Acts::Experimental::GraphBasedTrackSeeder(
0059       Acts::Experimental::GraphBasedTrackSeeder::DerivedConfig(
0060           m_cfg.seedFinderConfig),
0061       geometry, this->logger().cloneWithSuffix("GbtsFinder"));
0062
0063   m_filter = Acts::Experimental::GbtsTrackingFilter(m_cfg.trackingFilterConfig,
0064                                                     geometry);
0065
0066   // parse the mapping file and turn into map
0067   m_actsGbtsMap = makeActsGbtsMap();
0068
0069   printConfig();
0070 }
0071
0072 ProcessCode GraphBasedSeedingAlgorithm::execute(
0073     const AlgorithmContext &ctx) const {
0074   // initialise input space points from handle and define new container
0075   const SpacePointContainer &spacePoints = m_inputSpacePoints(ctx);
0076
0077   // take space points, add variables needed for GBTS and add them to new
0078   // container due to how space point container works, we need to keep the
0079   // container and the external columns we added alive this is done by using a
0080   // tuple of the core container and the two extra columns
0081   auto coreSpacePoints = makeSpContainer(spacePoints, m_actsGbtsMap);
0082
0083   // used to reserve size of nodes 2D vector in core
0084   std::uint32_t maxLayers = m_layerIdMap.size();
0085
0086   // ROI file:Defines what region in detector we are interested in, currently
0087   // set to entire detector
0088   Acts::Experimental::GbtsRoiDescriptor internalRoi(
0089       0, -4.5, 4.5, 0, -std::numbers::pi, std::numbers::pi, 0, -150., 150.);
0090
0091   // create the seeds
0092   Acts::SeedContainer2 seeds =
0093       m_finder->createSeeds(coreSpacePoints, internalRoi, maxLayers, *m_filter);
0094
0095   // update seed space point indices to original space point container
0096   for (auto seed : seeds) {
0097     for (auto &spIndex : seed.spacePointIndices()) {
0098       spIndex = coreSpacePoints.at(spIndex).copyFromIndex();
0099     }
0100   }
0101
0102   m_outputSeeds(ctx, std::move(seeds));
0103
0104   return ProcessCode::SUCCESS;
0105 }
0106
0107 std::map<GraphBasedSeedingAlgorithm::ActsIDs,
0108          GraphBasedSeedingAlgorithm::GbtsIDs>
0109 GraphBasedSeedingAlgorithm::makeActsGbtsMap() const {
0110   std::map<ActsIDs, GbtsIDs> actsToGbtsMap;
0111
0112   // prepare the acts to gbts mapping file
0113   // 0 in this file refers to no Gbts ID
0114   std::ifstream data(m_cfg.layerMappingFile);
0115   std::string line;
0116   // row = physical module, column = ACTS ID components
0117   std::vector<std::vector<std::string>> parsedCsv;
0118   while (std::getline(data, line)) {
0119     std::stringstream lineStream(line);
0120     std::string cell;
0121     std::vector<std::string> parsedRow;
0122     while (std::getline(lineStream, cell, ',')) {
0123       parsedRow.push_back(cell);
0124     }
0125
0126     parsedCsv.push_back(parsedRow);
0127   }
0128
0129   // file in format ACTS_vol,ACTS_lay,ACTS_mod,gbtsId
0130   for (auto i : parsedCsv) {
0131     std::uint32_t actsVol = stoi(i[0]);
0132     std::uint32_t actsLay = stoi(i[1]);
0133     std::uint32_t actsMod = stoi(i[2]);
0134     std::uint32_t gbts = stoi(i[5]);
0135     std::uint32_t etaMod = stoi(i[6]);
0136     std::uint32_t actsJoint = actsVol * 100 + actsLay;
0137     ActsIDs actsId{static_cast<std::uint64_t>(actsJoint),
0138                    static_cast<std::uint64_t>(actsMod)};
0139     GbtsIDs gbtsId{static_cast<std::uint32_t>(gbts),
0140                    static_cast<std::uint32_t>(etaMod), 0};
0141     actsToGbtsMap.insert({{actsId}, {gbtsId}});
0142   }
0143
0144   return actsToGbtsMap;
0145 }
0146
0147 Acts::SpacePointContainer2 GraphBasedSeedingAlgorithm::makeSpContainer(
0148     const SpacePointContainer &spacePoints,
0149     std::map<ActsIDs, GbtsIDs> map) const {
0150   Acts::SpacePointContainer2 coreSpacePoints(
0151       Acts::SpacePointColumns::X | Acts::SpacePointColumns::Y |
0152       Acts::SpacePointColumns::Z | Acts::SpacePointColumns::R |
0153       Acts::SpacePointColumns::Phi | Acts::SpacePointColumns::CopyFromIndex);
0154
0155   // add new column for layer ID and clusterwidth
0156   auto layerColomn = coreSpacePoints.createColumn<std::uint32_t>("layerId");
0157   auto clusterWidthColomn = coreSpacePoints.createColumn<float>("clusterWidth");
0158   auto localPositionColomn =
0159       coreSpacePoints.createColumn<float>("localPositionY");
0160   coreSpacePoints.reserve(spacePoints.size());
0161
0162   // for loop filling space point container and assigning layer ID's using the
0163   // map, also assigning cluster width and local position (currently false input
0164   // as not in examples but will be added in future)
0165   for (const auto &spacePoint : spacePoints) {
0166     // Gbts space point vector
0167     // loop over space points, call on map
0168     const auto &sourceLink = spacePoint.sourceLinks();
0169
0170     // warning if source link empty
0171     if (sourceLink.empty()) {
0172       // warning in officaial acts format
0173       ACTS_WARNING("warning source link vector is empty");
0174       continue;
0175     }
0176
0177     const auto &indexSourceLink = sourceLink.front().get<IndexSourceLink>();
0178
0179     std::uint32_t actsVolId = indexSourceLink.geometryId().volume();
0180     std::uint32_t actsLayId = indexSourceLink.geometryId().layer();
0181     std::uint32_t actsModId = indexSourceLink.geometryId().sensitive();
0182
0183     // dont want strips or HGTD
0184     if (actsVolId == 2 || actsVolId == 22 || actsVolId == 23 ||
0185         actsVolId == 24) {
0186       continue;
0187     }
0188
0189     // Search for vol, lay and module=0, if doesn't esist (end) then search
0190     // for full thing vol*100+lay as first number in pair then 0 or mod id
0191     auto actsJointId = actsVolId * 100 + actsLayId;
0192
0193     // here the key needs to be pair of(vol*100+lay, 0)
0194     ActsIDs key{static_cast<std::uint64_t>(actsJointId), 0};
0195     auto find = map.find(key);
0196
0197     // if end then make new key of (vol*100+lay, modid)
0198     if (find == map.end()) {
0199       key = ActsIDs{static_cast<std::uint64_t>(actsJointId),
0200                     static_cast<std::uint64_t>(actsModId)};  // mod ID
0201       find = map.find(key);
0202     }
0203
0204     // warning if key not in map
0205     if (find == map.end()) {
0206       ACTS_WARNING("Key not found in Gbts map for volume id: "
0207                    << actsVolId << " and layer id: " << actsLayId);
0208       continue;
0209     }
0210
0211     // now should be pixel with Gbts ID:
0212     // new map the item is a pair so want first from it
0213     std::uint32_t gbtsId = std::get<0>(find->second);
0214
0215     if (gbtsId == 0) {
0216       ACTS_WARNING("No assigned Gbts ID for key for volume id: "
0217                    << actsVolId << " and layer id: " << actsLayId);
0218     }
0219
0220     // access IDs from map
0221
0222     auto newSp = coreSpacePoints.createSpacePoint();
0223
0224     newSp.x() = spacePoint.x();
0225     newSp.y() = spacePoint.y();
0226     newSp.z() = spacePoint.z();
0227     newSp.r() = spacePoint.r();
0228     newSp.phi() = std::atan2(spacePoint.y(), spacePoint.x());
0229     newSp.copyFromIndex() = spacePoint.index();
0230     newSp.extra(layerColomn) = std::get<2>(find->second);
0231     // false input as this is not available in examples
0232     newSp.extra(clusterWidthColomn) = 0;
0233     newSp.extra(localPositionColomn) = 0;
0234   }
0235
0236   ACTS_VERBOSE("space point collection successfully assigned layerId's");
0237
0238   return coreSpacePoints;
0239 }
0240
0241 std::vector<Acts::Experimental::GbtsLayerDescription>
0242 GraphBasedSeedingAlgorithm::layerNumbering(const Acts::GeometryContext &gctx) {
0243   std::vector<Acts::Experimental::GbtsLayerDescription> inputVector;
0244   std::vector<std::size_t> countVector;
0245
0246   m_cfg.trackingGeometry->visitSurfaces([this, &inputVector, &countVector,
0247                                          &gctx](const Acts::Surface *surface) {
0248     Acts::GeometryIdentifier geoId = surface->geometryId();
0249     auto actsVolId = geoId.volume();
0250     auto actsLayId = geoId.layer();
0251     auto mod_id = geoId.sensitive();
0252     auto bounds_vect = surface->bounds().values();
0253     auto center = surface->center(gctx);
0254
0255     // make bounds global
0256     Acts::Vector3 globalFakeMom(1, 1, 1);
0257     Acts::Vector2 min_bound_local =
0258         Acts::Vector2(bounds_vect[0], bounds_vect[1]);
0259     Acts::Vector2 max_bound_local =
0260         Acts::Vector2(bounds_vect[2], bounds_vect[3]);
0261     Acts::Vector3 min_bound_global =
0262         surface->localToGlobal(gctx, min_bound_local, globalFakeMom);
0263     Acts::Vector3 max_bound_global =
0264         surface->localToGlobal(gctx, max_bound_local, globalFakeMom);
0265
0266     // checking that not wrong way round
0267     if (min_bound_global(0) > max_bound_global(0)) {
0268       min_bound_global.swap(max_bound_global);
0269     }
0270
0271     float rc = 0.0;
0272     float minBound = 100000.0;
0273     float maxBound = -100000.0;
0274
0275     // convert to Gbts ID
0276     auto actsJointId = actsVolId * 100 + actsLayId;
0277     // here the key needs to be pair of(vol*100+lay, 0)
0278     auto key = ActsIDs{actsJointId, 0};
0279     auto find = m_actsGbtsMap.find(key);
0280     // initialise first to avoid FLTUND later
0281     std::uint32_t gbtsId = 0;
0282     // new map, item is pair want first
0283     gbtsId = std::get<0>(find->second);
0284     // if end then make new key of (vol*100+lay, modid)
0285     if (find == m_actsGbtsMap.end()) {
0286       key = ActsIDs{actsJointId, mod_id};  // mod ID
0287       find = m_actsGbtsMap.find(key);
0288       gbtsId = std::get<0>(find->second);
0289     }
0290
0291     Acts::Experimental::GbtsLayerType barrelEc =
0292         Acts::Experimental::GbtsLayerType::Barrel;  // a variable that says if
0293                                                     // barrrel, 0 = barrel
0294     std::uint32_t etaMod = std::get<1>(find->second);
0295
0296     // assign barrelEc depending on Gbts_layer
0297     if (79 < gbtsId && gbtsId < 85) {  // 80s, barrel
0298       barrelEc = Acts::Experimental::GbtsLayerType::Barrel;
0299     } else if (89 < gbtsId && gbtsId < 99) {  // 90s positive
0300       barrelEc = Acts::Experimental::GbtsLayerType::Endcap;
0301     } else {  // 70s negative
0302       barrelEc = Acts::Experimental::GbtsLayerType::Endcap;
0303     }
0304
0305     if (barrelEc == Acts::Experimental::GbtsLayerType::Barrel) {
0306       rc = std::sqrt(center(0) * center(0) +
0307                      center(1) * center(1));  // barrel center in r
0308       // bounds of z
0309       if (min_bound_global(2) < minBound) {
0310         minBound = min_bound_global(2);
0311       }
0312       if (max_bound_global(2) > maxBound) {
0313         maxBound = max_bound_global(2);
0314       }
0315     } else if (barrelEc == Acts::Experimental::GbtsLayerType::Endcap) {
0316       rc = center(2);  // not barrel center in Z
0317       // bounds of r
0318       float min = std::sqrt(min_bound_global(0) * min_bound_global(0) +
0319                             min_bound_global(1) * min_bound_global(1));
0320       float max = std::sqrt(max_bound_global(0) * max_bound_global(0) +
0321                             max_bound_global(1) * max_bound_global(1));
0322       if (min < minBound) {
0323         minBound = min;
0324       }
0325       if (max > maxBound) {
0326         maxBound = max;
0327       }
0328     } else {
0329       throw std::runtime_error(
0330           "Invalid barrel/endcap assignment for GbtsLayer");
0331     }
0332
0333     std::int32_t combinedId = gbtsId * 1000 + etaMod;
0334
0335     const auto currentIndex =
0336         find_if(inputVector.begin(), inputVector.end(),
0337                 [combinedId](auto n) { return n.id == combinedId; });
0338     if (currentIndex != inputVector.end()) {  // not end so does exist
0339       std::size_t index = std::distance(inputVector.begin(), currentIndex);
0340       inputVector[index].refCoord += rc;
0341       inputVector[index].minBound =
0342           std::min(inputVector[index].minBound, minBound);
0343       inputVector[index].maxBound =
0344           std::max(inputVector[index].maxBound, maxBound);
0345       countVector[index] += 1;  // increase count at the index
0346
0347     } else {  // end so doesn't exists
0348       // make new if one with Gbts ID doesn't exist:
0349       Acts::Experimental::GbtsLayerDescription newGbtsId(
0350           combinedId, barrelEc, rc, minBound, maxBound);
0351       inputVector.push_back(newGbtsId);
0352       // so the element exists and not divinding by 0
0353       countVector.push_back(1);
0354
0355       // tracking the index of each GbtsLayerDescription as there added
0356
0357       // so layer ID refers to actual index and not size of vector
0358       std::uint32_t layerId = countVector.size() - 1;
0359       std::get<2>(find->second) = layerId;
0360       m_layerIdMap.insert({combinedId, layerId});
0361     }
0362     // look up for every combined ID to see if it has a layer
0363     if (auto findLayer = m_layerIdMap.find(combinedId);
0364         findLayer == m_layerIdMap.end()) {
0365       ACTS_WARNING("No assigned Layer ID for combined ID: " << combinedId);
0366     } else {
0367       std::get<2>(find->second) = findLayer->second;
0368     }
0369
0370     // add to file each time,
0371     // print to csv for each module, no repeats so dont need to make
0372     // map for averaging
0373     if (m_cfg.fillModuleCsv) {
0374       std::fstream fout;
0375       fout.open("ACTS_modules.csv", std::ios::out | std::ios::app);
0376       fout << actsVolId << ", "  // vol
0377            << actsLayId << ", "  // lay
0378            << mod_id << ", "     // module
0379            << gbtsId << ","      // Gbts id
0380            << etaMod << ","      // etaMod
0381            << center(2) << ", "  // z
0382            << std::sqrt(center(0) * center(0) + center(1) * center(1))  // r
0383            << "\n";
0384     }
0385   });
0386
0387   for (std::size_t i = 0; i < inputVector.size(); i++) {
0388     inputVector[i].refCoord = inputVector[i].refCoord / countVector[i];
0389   }
0390
0391   return inputVector;
0392 }
0393
0394 void GraphBasedSeedingAlgorithm::printConfig() const {
0395   ACTS_DEBUG("===== GraphBasedTrackSeeder =====");
0396   const auto &cfg1 = m_cfg.seedFinderConfig;
0397   ACTS_DEBUG("BeamSpotCorrection: " << cfg1.beamSpotCorrection);
0398   ACTS_DEBUG("connectorInputFile: " << cfg1.connectorInputFile);
0399   ACTS_DEBUG("lutInputFile: " << cfg1.lutInputFile);
0400   ACTS_DEBUG("lrtMode: " << cfg1.lrtMode);
0401   ACTS_DEBUG("useMl: " << cfg1.useMl);
0402   ACTS_DEBUG("matchBeforeCreate: " << cfg1.matchBeforeCreate);
0403   ACTS_DEBUG("useOldTunings: " << cfg1.useOldTunings);
0404   ACTS_DEBUG("tauRatioCut: " << cfg1.tauRatioCut);
0405   ACTS_DEBUG("tauRatioPrecut: " << cfg1.tauRatioPrecut);
0406   ACTS_DEBUG("etaBinWidthOverride: " << cfg1.etaBinWidthOverride);
0407   ACTS_DEBUG("nMaxPhiSlice: " << cfg1.nMaxPhiSlice);
0408   ACTS_DEBUG("minPt: " << cfg1.minPt);
0409   ACTS_DEBUG("ptCoeff: " << cfg1.ptCoeff);
0410   ACTS_DEBUG("useEtaBinning: " << cfg1.useEtaBinning);
0411   ACTS_DEBUG("doubletFilterRZ: " << cfg1.doubletFilterRZ);
0412   ACTS_DEBUG("nMaxEdges: " << cfg1.nMaxEdges);
0413   ACTS_DEBUG("minDeltaRadius: " << cfg1.minDeltaRadius);
0414   ACTS_DEBUG("edgeMaskMinEta: " << cfg1.edgeMaskMinEta);
0415   ACTS_DEBUG("hitShareThreshold: " << cfg1.hitShareThreshold);
0416   ACTS_DEBUG("maxEndcapClusterWidth: " << cfg1.maxEndcapClusterWidth);
0417   ACTS_DEBUG("===== GbtsTrackFilter =====");
0418   const auto &cfg2 = m_cfg.trackingFilterConfig;
0419   ACTS_DEBUG("sigmaMS: " << cfg2.sigmaMS);
0420   ACTS_DEBUG("radLen: " << cfg2.radLen);
0421   ACTS_DEBUG("sigmaX: " << cfg2.sigmaX);
0422   ACTS_DEBUG("sigmaY: " << cfg2.sigmaY);
0423   ACTS_DEBUG("weightX: " << cfg2.weightX);
0424   ACTS_DEBUG("weightY: " << cfg2.weightY);
0425   ACTS_DEBUG("maxDChi2X: " << cfg2.maxDChi2X);
0426   ACTS_DEBUG("maxDChi2Y: " << cfg2.maxDChi2Y);
0427   ACTS_DEBUG("addHit: " << cfg2.addHit);
0428   ACTS_DEBUG("maxCurvature: " << cfg2.maxCurvature);
0429   ACTS_DEBUG("maxZ0: " << cfg2.maxZ0);
0430   ACTS_DEBUG("================================");
0431 }
0432
0433 }  // namespace ActsExamples