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