TrackFinding/src/GbtsSeedingAlgorithm.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/GbtsSeedingAlgorithm.hpp"
0010
0011 #include "Acts/EventData/SeedContainer2.hpp"
0012 #include "Acts/EventData/SpacePointContainer2.hpp"
0013 #include "Acts/Geometry/GeometryIdentifier.hpp"
0014 #include "ActsExamples/EventData/IndexSourceLink.hpp"
0015 #include "ActsExamples/EventData/SimSeed.hpp"
0016
0017 #include <cmath>
0018 #include <fstream>
0019 #include <iostream>
0020 #include <map>
0021 #include <numbers>
0022 #include <sstream>
0023 #include <stdexcept>
0024 #include <vector>
0025
0026 ActsExamples::GbtsSeedingAlgorithm::GbtsSeedingAlgorithm(
0027     ActsExamples::GbtsSeedingAlgorithm::Config cfg, Acts::Logging::Level lvl)
0028     : ActsExamples::IAlgorithm("SeedingAlgorithm", lvl), m_cfg(std::move(cfg)) {
0029   // initialise the spacepoint, seed and cluster handles
0030   m_inputSpacePoints.initialize(
0031       m_cfg.inputSpacePoints);  // TO DO: change bindings so it only gives a
0032                                 // string instead of a vector
0033   m_outputSeeds.initialize(m_cfg.outputSeeds);
0034   m_inputClusters.initialize(m_cfg.inputClusters);
0035
0036   // parse the mapping file and turn into map
0037   m_cfg.ActsGbtsMap = makeActsGbtsMap();
0038
0039   // create the TrigInDetSiLayers (Logical Layers),
0040   // as well as a map that tracks there index in m_layerGeometry
0041   m_layerGeometry = LayerNumbering();
0042
0043   // parse connection file
0044   std::ifstream input_ifstream(
0045       m_cfg.seedFinderConfig.connectorInputFile.c_str(), std::ifstream::in);
0046
0047   if (input_ifstream.peek() == std::ifstream::traits_type::eof()) {
0048     ACTS_WARNING("Cannot find layer connections file ");
0049     throw std::runtime_error("connection file not found");
0050
0051   }
0052
0053   // create the connection objects
0054   else {
0055     m_connector = std::make_unique<Acts::Experimental::GbtsConnector>(
0056         input_ifstream, m_cfg.seedFinderConfig.LRTmode);
0057
0058     // option that allows for adding custom eta binning (default is at 0.2)
0059     if (m_cfg.seedFinderConfig.etaBinOverride != 0.0f) {
0060       m_connector->m_etaBin = m_cfg.seedFinderConfig.etaBinOverride;
0061     }
0062   }
0063
0064   // initialise the object that holds all the geometry information needed for
0065   // the algorithm
0066   m_gbtsGeo = std::make_unique<Acts::Experimental::GbtsGeometry>(
0067       m_layerGeometry, m_connector);
0068
0069   // manually convert min Pt as no conversion available in ACTS Examples
0070   // (currently inputs as 0.9 GeV but need 900 MeV)
0071   m_cfg.seedFinderConfig.minPt = m_cfg.seedFinderConfig.minPt * 1000;
0072   printSeedFinderGbtsConfig(m_cfg.seedFinderConfig);
0073 }
0074
0075 // execute:
0076 ActsExamples::ProcessCode ActsExamples::GbtsSeedingAlgorithm::execute(
0077     const AlgorithmContext &ctx) const {
0078   // take spacepoints, add variables needed for GBTS and add them to new
0079   // container due to how spacepoint container works, we need to keep the
0080   // container and the external columns we added alive this is done by using a
0081   // tuple of the core container and the two extra columns
0082   auto SpContainerComponents = MakeSpContainer(ctx, m_cfg.ActsGbtsMap);
0083
0084   // this is now calling on a core algorithm
0085   Acts::Experimental::SeedFinderGbts finder(
0086       m_cfg.seedFinderConfig, m_gbtsGeo.get(), &m_layerGeometry,
0087       logger().cloneWithSuffix("GbtdFinder"));
0088
0089   // used to reserve size of nodes 2D vector in core
0090   int max_layers = m_LayeridMap.size();
0091
0092   // ROI file:Defines what region in detector we are interested in, currently
0093   // set to entire detector
0094   //  Acts::Experimental::RoiDescriptor internalRoi(0, -5, 5, 0,
0095   //  -std::numbers::pi, std::numbers::pi, 0, -225., 225.);
0096   Acts::Experimental::RoiDescriptor internalRoi(
0097       0, -4.5, 4.5, 0, -std::numbers::pi, std::numbers::pi, 0, -150., 150.);
0098
0099   // create the seeds
0100
0101   Acts::SeedContainer2 seeds =
0102       finder.CreateSeeds(internalRoi, SpContainerComponents, max_layers);
0103
0104   // move seeds to simseedcontainer to be used down stream taking fist middle
0105   // and last sps currently as simseeds need to be hard types so only 3
0106   // spacepoint can be added but in future we should be able to have any length
0107   // seed
0108   SimSeedContainer seedContainerForStorage;
0109   seedContainerForStorage.reserve(seeds.size());
0110   for (const auto &seed : seeds) {
0111     auto sps = seed.spacePointIndices();
0112     unsigned int indices = sps.size() - 1;
0113     std::size_t mid = static_cast<std::size_t>(std::round(indices / 2.0));
0114     seedContainerForStorage.emplace_back(
0115         *std::get<0>(SpContainerComponents)
0116              .at(sps[0])  // first spacepoint
0117              .sourceLinks()[0]
0118              .get<const SimSpacePoint *>(),
0119         *std::get<0>(SpContainerComponents)
0120              .at(sps[mid])  // middle spacepoint
0121              .sourceLinks()[0]
0122              .get<const SimSpacePoint *>(),
0123         *std::get<0>(SpContainerComponents)
0124              .at(sps[indices])  // last spacepoint
0125              .sourceLinks()[0]
0126              .get<const SimSpacePoint *>());
0127
0128     seedContainerForStorage.back().setVertexZ(seed.vertexZ());
0129     seedContainerForStorage.back().setQuality(seed.quality());
0130   }
0131
0132   m_outputSeeds(ctx, std::move(seedContainerForStorage));
0133
0134   return ActsExamples::ProcessCode::SUCCESS;
0135 }
0136
0137 std::map<std::pair<int, int>, std::tuple<int, int, int>>
0138 ActsExamples::GbtsSeedingAlgorithm::makeActsGbtsMap() const {
0139   std::map<std::pair<int, int>, std::tuple<int, int, int>> ActsGbts;
0140
0141   // prepare the acts to gbts mapping file
0142   std::ifstream data(
0143       m_cfg.layerMappingFile);  // 0 in this file refers to no Gbts ID
0144   std::string line;
0145   std::vector<std::vector<std::string>>
0146       parsedCsv;  // row = physical module, column = ACTS ID components
0147   while (std::getline(data, line)) {
0148     std::stringstream lineStream(line);
0149     std::string cell;
0150     std::vector<std::string> parsedRow;
0151     while (std::getline(lineStream, cell, ',')) {
0152       parsedRow.push_back(cell);
0153     }
0154
0155     parsedCsv.push_back(parsedRow);
0156   }
0157
0158   // file in format ACTS_vol,ACTS_lay,ACTS_mod,Gbts_id
0159   for (auto i : parsedCsv) {
0160     int ACTS_vol = stoi(i[0]);
0161     int ACTS_lay = stoi(i[1]);
0162     int ACTS_mod = stoi(i[2]);
0163     int Gbts = stoi(i[5]);
0164     int eta_mod = stoi(i[6]);
0165     int ACTS_joint = ACTS_vol * 100 + ACTS_lay;
0166     ActsGbts.insert({{ACTS_joint, ACTS_mod}, {Gbts, eta_mod, 0}});
0167   }
0168
0169   return ActsGbts;
0170 }
0171
0172 Acts::Experimental::SPContainerComponentsType
0173 ActsExamples::GbtsSeedingAlgorithm::MakeSpContainer(
0174     const AlgorithmContext &ctx,
0175     std::map<std::pair<int, int>, std::tuple<int, int, int>> map) const {
0176   // new seeding container test
0177   // initialise input spacepoints from handle and define new container
0178   const SimSpacePointContainer &spacePoints = m_inputSpacePoints(ctx);
0179   Acts::SpacePointContainer2 coreSpacePoints(
0180
0181       Acts::SpacePointColumns::SourceLinks | Acts::SpacePointColumns::X |
0182       Acts::SpacePointColumns::Y | Acts::SpacePointColumns::Z |
0183       Acts::SpacePointColumns::R | Acts::SpacePointColumns::Phi);
0184
0185   // add new column for layer ID and clusterwidth
0186   auto LayerColoumn = coreSpacePoints.createColumn<int>("LayerID");
0187   auto ClusterWidthColoumn =
0188       coreSpacePoints.createColumn<float>("Cluster_Width");
0189   auto LocalPositionColoumn =
0190       coreSpacePoints.createColumn<float>("LocalPositionY");
0191   coreSpacePoints.reserve(spacePoints.size());
0192
0193   // for loop filling space
0194
0195   for (const auto &spacePoint : spacePoints) {
0196     // Gbts space point vector
0197     // loop over space points, call on map
0198     const auto &sourceLink = spacePoint.sourceLinks();
0199
0200     // warning if source link empty
0201     if (sourceLink.empty()) {
0202       // warning in officaial acts format
0203       ACTS_WARNING("warning source link vector is empty");
0204       continue;
0205     }
0206
0207     const auto &indexSourceLink = sourceLink.front().get<IndexSourceLink>();
0208
0209     int ACTS_vol_id = indexSourceLink.geometryId().volume();
0210     int ACTS_lay_id = indexSourceLink.geometryId().layer();
0211     int ACTS_mod_id = indexSourceLink.geometryId().sensitive();
0212
0213     // dont want strips or HGTD
0214     if (ACTS_vol_id == 2 || ACTS_vol_id == 22 || ACTS_vol_id == 23 ||
0215         ACTS_vol_id == 24) {
0216       continue;
0217     }
0218
0219     // Search for vol, lay and module=0, if doesn't esist (end) then search
0220     // for full thing vol*100+lay as first number in pair then 0 or mod id
0221     auto ACTS_joint_id = ACTS_vol_id * 100 + ACTS_lay_id;
0222     auto key =
0223         std::make_pair(ACTS_joint_id,
0224                        0);  // here the key needs to be pair of(vol*100+lay, 0)
0225     auto Find = map.find(key);
0226
0227     if (Find ==
0228         map.end()) {  // if end then make new key of (vol*100+lay, modid)
0229       key = std::make_pair(ACTS_joint_id, ACTS_mod_id);  // mod ID
0230       Find = map.find(key);
0231     }
0232
0233     // warning if key not in map
0234     if (Find == map.end()) {
0235       ACTS_WARNING("Key not found in Gbts map for volume id: "
0236                    << ACTS_vol_id << " and layer id: " << ACTS_lay_id);
0237       continue;
0238     }
0239
0240     // now should be pixel with Gbts ID:
0241     int Gbts_id = std::get<0>(
0242         Find->second);  // new map the item is a pair so want first from it
0243
0244     if (Gbts_id == 0) {
0245       ACTS_WARNING("No assigned Gbts ID for key for volume id: "
0246                    << ACTS_vol_id << " and layer id: " << ACTS_lay_id);
0247     }
0248
0249     // access IDs from map
0250
0251     auto newSp = coreSpacePoints.createSpacePoint();
0252
0253     newSp.assignSourceLinks(
0254         std::array<Acts::SourceLink, 1>{Acts::SourceLink(&spacePoint)});
0255     newSp.x() = spacePoint.x();
0256     newSp.y() = spacePoint.y();
0257     newSp.z() = spacePoint.z();
0258     newSp.r() = spacePoint.r();
0259     newSp.phi() = std::atan2(spacePoint.y(), spacePoint.x());
0260     newSp.extra(LayerColoumn) = std::get<2>(Find->second);
0261     // false input as this is not available in examples
0262     newSp.extra(ClusterWidthColoumn) = 0;
0263     newSp.extra(LocalPositionColoumn) = 0;
0264   }
0265
0266   ACTS_VERBOSE("Space point collection successfully assigned LayerID's");
0267
0268   return std::make_tuple(std::move(coreSpacePoints), LayerColoumn.asConst(),
0269                          ClusterWidthColoumn.asConst(),
0270                          LocalPositionColoumn.asConst());
0271 }
0272
0273 std::vector<Acts::Experimental::TrigInDetSiLayer>
0274 ActsExamples::GbtsSeedingAlgorithm::LayerNumbering() const {
0275   std::vector<Acts::Experimental::TrigInDetSiLayer> input_vector{};
0276   std::vector<std::size_t> count_vector{};
0277
0278   m_cfg.trackingGeometry->visitSurfaces([this, &input_vector, &count_vector](
0279                                             const Acts::Surface *surface) {
0280     Acts::GeometryIdentifier geoid = surface->geometryId();
0281     auto ACTS_vol_id = geoid.volume();
0282     auto ACTS_lay_id = geoid.layer();
0283     auto mod_id = geoid.sensitive();
0284     auto bounds_vect = surface->bounds().values();
0285     auto center = surface->center(Acts::GeometryContext());
0286
0287     // make bounds global
0288     Acts::Vector3 globalFakeMom(1, 1, 1);
0289     Acts::Vector2 min_bound_local =
0290         Acts::Vector2(bounds_vect[0], bounds_vect[1]);
0291     Acts::Vector2 max_bound_local =
0292         Acts::Vector2(bounds_vect[2], bounds_vect[3]);
0293     Acts::Vector3 min_bound_global = surface->localToGlobal(
0294         Acts::GeometryContext(), min_bound_local, globalFakeMom);
0295     Acts::Vector3 max_bound_global = surface->localToGlobal(
0296         Acts::GeometryContext(), max_bound_local, globalFakeMom);
0297
0298     if (min_bound_global(0) >
0299         max_bound_global(0)) {  // checking that not wrong way round
0300       min_bound_global.swap(max_bound_global);
0301     }
0302
0303     float rc = 0.0;
0304     float minBound = 100000.0;
0305     float maxBound = -100000.0;
0306
0307     // convert to Gbts ID
0308     auto ACTS_joint_id = ACTS_vol_id * 100 + ACTS_lay_id;
0309     auto key =
0310         std::make_pair(ACTS_joint_id,
0311                        0);  // here the key needs to be pair of(vol*100+lay, 0)
0312     auto Find = m_cfg.ActsGbtsMap.find(key);
0313     int Gbts_id = 0;  // initialise first to avoid FLTUND later
0314     Gbts_id = std::get<0>(Find->second);  // new map, item is pair want first
0315     if (Find ==
0316         m_cfg.ActsGbtsMap
0317             .end()) {  // if end then make new key of (vol*100+lay, modid)
0318       key = std::make_pair(ACTS_joint_id, mod_id);  // mod ID
0319       Find = m_cfg.ActsGbtsMap.find(key);
0320       Gbts_id = std::get<0>(Find->second);
0321     }
0322
0323     short barrel_ec = 0;  // a variable that says if barrrel, 0 = barrel
0324     int eta_mod = std::get<1>(Find->second);
0325
0326     // assign barrel_ec depending on Gbts_layer
0327     if (79 < Gbts_id && Gbts_id < 85) {  // 80s, barrel
0328       barrel_ec = 0;
0329     } else if (89 < Gbts_id && Gbts_id < 99) {  // 90s positive
0330       barrel_ec = 2;
0331     } else {  // 70s negative
0332       barrel_ec = -2;
0333     }
0334
0335     if (barrel_ec == 0) {
0336       rc = sqrt(center(0) * center(0) +
0337                 center(1) * center(1));  // barrel center in r
0338       // bounds of z
0339       if (min_bound_global(2) < minBound) {
0340         minBound = min_bound_global(2);
0341       }
0342       if (max_bound_global(2) > maxBound) {
0343         maxBound = max_bound_global(2);
0344       }
0345     } else {
0346       rc = center(2);  // not barrel center in Z
0347       // bounds of r
0348       float min = sqrt(min_bound_global(0) * min_bound_global(0) +
0349                        min_bound_global(1) * min_bound_global(1));
0350       float max = sqrt(max_bound_global(0) * max_bound_global(0) +
0351                        max_bound_global(1) * max_bound_global(1));
0352       if (min < minBound) {
0353         minBound = min;
0354       }
0355       if (max > maxBound) {
0356         maxBound = max;
0357       }
0358     }
0359
0360     int combined_id = Gbts_id * 1000 + eta_mod;
0361
0362     auto current_index =
0363         find_if(input_vector.begin(), input_vector.end(),
0364                 [combined_id](auto n) { return n.m_subdet == combined_id; });
0365     if (current_index != input_vector.end()) {  // not end so does exist
0366       std::size_t index = std::distance(input_vector.begin(), current_index);
0367       input_vector[index].m_refCoord += rc;
0368       input_vector[index].m_minBound += minBound;
0369       input_vector[index].m_maxBound += maxBound;
0370       count_vector[index] += 1;  // increase count at the index
0371
0372     } else {  // end so doesn't exists
0373       // make new if one with Gbts ID doesn't exist:
0374       Acts::Experimental::TrigInDetSiLayer new_Gbts_ID(combined_id, barrel_ec,
0375                                                        rc, minBound, maxBound);
0376       input_vector.push_back(new_Gbts_ID);
0377       count_vector.push_back(
0378           1);  // so the element exists and not divinding by 0
0379
0380       // tracking the index of each TrigInDetSiLayer as there added to the
0381       // vector
0382       int LayerID =
0383           count_vector.size() -
0384           1;  // so layer ID refers to actual index and not size of vector
0385       std::get<2>(Find->second) = LayerID;
0386       m_LayeridMap.insert({combined_id, LayerID});
0387     }
0388     // look up for every combined ID to see if it has a layer
0389     auto FindLayer = m_LayeridMap.find(combined_id);
0390     if (FindLayer == m_LayeridMap.end()) {
0391       ACTS_WARNING("No assigned Layer ID for combined ID: " << combined_id);
0392     } else {
0393       std::get<2>(Find->second) = FindLayer->second;
0394     }
0395
0396     if (m_cfg.fill_module_csv) {
0397       std::fstream fout;
0398       fout.open("ACTS_modules.csv",
0399                 std::ios::out | std::ios::app);  // add to file each time
0400       // print to csv for each module, no repeats so dont need to make
0401       // map for averaging
0402       fout << ACTS_vol_id << ", "                                  // vol
0403            << ACTS_lay_id << ", "                                  // lay
0404            << mod_id << ", "                                       // module
0405            << Gbts_id << ","                                       // Gbts id
0406            << eta_mod << ","                                       // eta_mod
0407            << center(2) << ", "                                    // z
0408            << sqrt(center(0) * center(0) + center(1) * center(1))  // r
0409            << "\n";
0410     }
0411   });
0412
0413   for (std::size_t i = 0; i < input_vector.size(); i++) {
0414     input_vector[i].m_refCoord = input_vector[i].m_refCoord / count_vector[i];
0415   }
0416
0417   return input_vector;
0418 }
0419
0420 void ActsExamples::GbtsSeedingAlgorithm::printSeedFinderGbtsConfig(
0421     const Acts::Experimental::SeedFinderGbtsConfig &cfg) {
0422   ACTS_DEBUG("===== SeedFinderGbtsConfig =====");
0423
0424   ACTS_DEBUG("BeamSpotCorrection: " << cfg.BeamSpotCorrection
0425                                     << " (default: false)");
0426   ACTS_DEBUG("connectorInputFile: " << cfg.connectorInputFile
0427                                     << " (default: empty string)");
0428   ACTS_DEBUG("lutInputFile: " << cfg.lutInputFile
0429                               << " (default: empty string)");
0430   ACTS_DEBUG("LRTmode: " << cfg.LRTmode << " (default: false)");
0431   ACTS_DEBUG("useML: " << cfg.useML << " (default: false)");
0432   ACTS_DEBUG("matchBeforeCreate: " << cfg.matchBeforeCreate
0433                                    << " (default: false)");
0434   ACTS_DEBUG("useOldTunings: " << cfg.useOldTunings << " (default: false)");
0435   ACTS_DEBUG("tau_ratio_cut: " << cfg.tau_ratio_cut << " (default: 0.007)");
0436   ACTS_DEBUG("tau_ratio_precut: " << cfg.tau_ratio_precut
0437                                   << " (default: 0.009f)");
0438   ACTS_DEBUG("etaBinOverride: " << cfg.etaBinOverride << " (default: 0.0)");
0439   ACTS_DEBUG("nMaxPhiSlice: " << cfg.nMaxPhiSlice << " (default: 53)");
0440   ACTS_DEBUG("minPt: " << cfg.minPt
0441                        << " (default: 1.0 * Acts::UnitConstants::GeV)");
0442   ACTS_DEBUG("phiSliceWidth: " << cfg.phiSliceWidth << " (default: derived)");
0443   ACTS_DEBUG("ptCoeff: " << cfg.ptCoeff
0444                          << " (default: 0.29997 * 1.9972 / 2.0)");
0445   ACTS_DEBUG("useEtaBinning: " << cfg.useEtaBinning << " (default: true)");
0446   ACTS_DEBUG("doubletFilterRZ: " << cfg.doubletFilterRZ << " (default: true)");
0447   ACTS_DEBUG("nMaxEdges: " << cfg.nMaxEdges << " (default: 2000000)");
0448   ACTS_DEBUG("minDeltaRadius: " << cfg.minDeltaRadius << " (default: 2.0)");
0449   ACTS_DEBUG("sigmaMS: " << cfg.sigmaMS << " (default: 0.016)");
0450   ACTS_DEBUG("radLen: " << cfg.radLen << " (default: 0.025)");
0451   ACTS_DEBUG("sigma_x: " << cfg.sigma_x << " (default: 0.08)");
0452   ACTS_DEBUG("sigma_y: " << cfg.sigma_y << " (default: 0.25)");
0453   ACTS_DEBUG("weight_x: " << cfg.weight_x << " (default: 0.5)");
0454   ACTS_DEBUG("weight_y: " << cfg.weight_y << " (default: 0.5)");
0455   ACTS_DEBUG("maxDChi2_x: " << cfg.maxDChi2_x << " (default: 5.0)");
0456   ACTS_DEBUG("maxDChi2_y: " << cfg.maxDChi2_y << " (default: 6.0)");
0457   ACTS_DEBUG("add_hit: " << cfg.add_hit << " (default: 14.0)");
0458   ACTS_DEBUG("max_curvature: " << cfg.max_curvature << " (default: 1e-3f)");
0459   ACTS_DEBUG("max_z0: " << cfg.max_z0 << " (default: 170.0)");
0460   ACTS_DEBUG("edge_mask_min_eta: " << cfg.edge_mask_min_eta
0461                                    << " (default: 1.5)");
0462   ACTS_DEBUG("hit_share_threshold: " << cfg.hit_share_threshold
0463                                      << " (default: 0.49)");
0464   ACTS_DEBUG("max_endcap_clusterwidth: " << cfg.max_endcap_clusterwidth
0465                                          << " (default: 0.35)");
0466
0467   ACTS_DEBUG("================================");
0468 }