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/Geometry/GeometryIdentifier.hpp"
0012 #include "Acts/Seeding/SeedFilter.hpp"
0013 #include "ActsExamples/EventData/IndexSourceLink.hpp"
0014 #include "ActsExamples/EventData/Measurement.hpp"
0015 #include "ActsExamples/EventData/ProtoTrack.hpp"
0016 #include "ActsExamples/EventData/SimSeed.hpp"
0017 #include "ActsExamples/Framework/WhiteBoard.hpp"
0018
0019 #include <fstream>
0020 #include <iostream>
0021 #include <map>
0022 #include <numbers>
0023 #include <random>
0024 #include <sstream>
0025 #include <vector>
0026
0027 template class Acts::GbtsLayer<ActsExamples::SimSpacePoint>;
0028 template class Acts::GbtsGeometry<ActsExamples::SimSpacePoint>;
0029 template class Acts::GbtsNode<ActsExamples::SimSpacePoint>;
0030 template class Acts::GbtsEtaBin<ActsExamples::SimSpacePoint>;
0031 template struct Acts::GbtsSP<ActsExamples::SimSpacePoint>;
0032 template class Acts::GbtsDataStorage<ActsExamples::SimSpacePoint>;
0033 template class Acts::GbtsEdge<ActsExamples::SimSpacePoint>;
0034
0035 // constructor:
0036 ActsExamples::GbtsSeedingAlgorithm::GbtsSeedingAlgorithm(
0037     ActsExamples::GbtsSeedingAlgorithm::Config cfg, Acts::Logging::Level lvl)
0038     : ActsExamples::IAlgorithm("SeedingAlgorithm", lvl), m_cfg(std::move(cfg)) {
0039   // fill config struct
0040   m_cfg.layerMappingFile = m_cfg.layerMappingFile;
0041
0042   m_cfg.seedFilterConfig = m_cfg.seedFilterConfig.toInternalUnits();
0043
0044   m_cfg.seedFinderConfig =
0045       m_cfg.seedFinderConfig.toInternalUnits().calculateDerivedQuantities();
0046
0047   m_cfg.seedFinderOptions =
0048       m_cfg.seedFinderOptions.toInternalUnits().calculateDerivedQuantities(
0049           m_cfg.seedFinderConfig);
0050
0051   for (const auto &spName : m_cfg.inputSpacePoints) {
0052     if (spName.empty()) {
0053       throw std::invalid_argument("Invalid space point input collection");
0054     }
0055
0056     auto &handle = m_inputSpacePoints.emplace_back(
0057         std::make_unique<ReadDataHandle<SimSpacePointContainer>>(
0058             this,
0059             "InputSpacePoints#" + std::to_string(m_inputSpacePoints.size())));
0060     handle->initialize(spName);
0061   }
0062
0063   m_outputSeeds.initialize(m_cfg.outputSeeds);
0064
0065   m_inputClusters.initialize(m_cfg.inputClusters);
0066
0067   // map
0068   m_cfg.ActsGbtsMap = makeActsGbtsMap();
0069   // input trig vector
0070   m_cfg.seedFinderConfig.m_layerGeometry = LayerNumbering();
0071
0072   std::ifstream input_ifstream(
0073       m_cfg.seedFinderConfig.connector_input_file.c_str(), std::ifstream::in);
0074
0075   // connector
0076   std::unique_ptr<Acts::GbtsConnector> inputConnector =
0077       std::make_unique<Acts::GbtsConnector>(input_ifstream);
0078
0079   m_gbtsGeo = std::make_unique<Acts::GbtsGeometry<SimSpacePoint>>(
0080       m_cfg.seedFinderConfig.m_layerGeometry, inputConnector);
0081
0082 }  // this is not Gbts config type because it is a member of the algs config,
0083    // which is of type Gbts cofig
0084
0085 // execute:
0086 ActsExamples::ProcessCode ActsExamples::GbtsSeedingAlgorithm::execute(
0087     const AlgorithmContext &ctx) const {
0088   std::vector<Acts::GbtsSP<SimSpacePoint>> GbtsSpacePoints =
0089       MakeGbtsSpacePoints(ctx, m_cfg.ActsGbtsMap);
0090
0091   for (auto sp : GbtsSpacePoints) {
0092     ACTS_DEBUG("Gbts space points:  Gbts_id: "
0093                << sp.gbtsID << " z: " << sp.SP->z() << " r: " << sp.SP->r()
0094                << " ACTS volume:  "
0095                << sp.SP->sourceLinks()
0096                       .front()
0097                       .get<IndexSourceLink>()
0098                       .geometryId()
0099                       .volume());
0100   }
0101
0102   // this is now calling on a core algorithm
0103   Acts::SeedFinderGbts<SimSpacePoint> finder(
0104       m_cfg.seedFinderConfig, *m_gbtsGeo,
0105       logger().cloneWithSuffix("GbtdFinder"));
0106
0107   // output of function needed for seed
0108
0109   finder.loadSpacePoints(GbtsSpacePoints);
0110
0111   // trigGbts file :
0112   Acts::RoiDescriptor internalRoi(0, -4.5, 4.5, 0, -std::numbers::pi,
0113                                   std::numbers::pi, 0, -150., 150.);
0114   // ROI file:
0115   //  Acts::RoiDescriptor internalRoi(0, -5, 5, 0, -std::numbers::pi,
0116   //  std::numbers::pi, 0, -225., 225.);
0117
0118   // new version returns seeds
0119   SimSeedContainer seeds = finder.createSeeds(internalRoi, *m_gbtsGeo);
0120
0121   m_outputSeeds(ctx, std::move(seeds));
0122
0123   return ActsExamples::ProcessCode::SUCCESS;
0124 }
0125
0126 std::map<std::pair<int, int>, std::pair<int, int>>
0127 ActsExamples::GbtsSeedingAlgorithm::makeActsGbtsMap() const {
0128   std::map<std::pair<int, int>, std::pair<int, int>> ActsGbts;
0129   std::ifstream data(
0130       m_cfg.layerMappingFile);  // 0 in this file refers to no Gbts ID
0131   std::string line;
0132   std::vector<std::vector<std::string>> parsedCsv;
0133   while (std::getline(data, line)) {
0134     std::stringstream lineStream(line);
0135     std::string cell;
0136     std::vector<std::string> parsedRow;
0137     while (std::getline(lineStream, cell, ',')) {
0138       parsedRow.push_back(cell);
0139     }
0140
0141     parsedCsv.push_back(parsedRow);
0142   }
0143   // file in format ACTS_vol,ACTS_lay,ACTS_mod,Gbts_id
0144   for (auto i : parsedCsv) {
0145     int ACTS_vol = stoi(i[0]);
0146     int ACTS_lay = stoi(i[1]);
0147     int ACTS_mod = stoi(i[2]);
0148     int Gbts = stoi(i[5]);
0149     int eta_mod = stoi(i[6]);
0150     int ACTS_joint = ACTS_vol * 100 + ACTS_lay;
0151     ActsGbts.insert({{ACTS_joint, ACTS_mod}, {Gbts, eta_mod}});
0152   }
0153
0154   return ActsGbts;
0155 }
0156
0157 std::vector<Acts::GbtsSP<ActsExamples::SimSpacePoint>>
0158 ActsExamples::GbtsSeedingAlgorithm::MakeGbtsSpacePoints(
0159     const AlgorithmContext &ctx,
0160     std::map<std::pair<int, int>, std::pair<int, int>> map) const {
0161   // create space point vectors
0162   std::vector<Acts::GbtsSP<ActsExamples::SimSpacePoint>> gbtsSpacePoints;
0163   gbtsSpacePoints.reserve(
0164       m_inputSpacePoints.size());  // not sure if this is enough
0165
0166   // for loop filling space
0167   for (const auto &isp : m_inputSpacePoints) {
0168     for (const auto &spacePoint : (*isp)(ctx)) {
0169       // Gbts space point vector
0170       // loop over space points, call on map
0171       const auto &sourceLink = spacePoint.sourceLinks();
0172       const auto &indexSourceLink = sourceLink.front().get<IndexSourceLink>();
0173
0174       // warning if source link empty
0175       if (sourceLink.empty()) {
0176         // warning in officaial acts format
0177         ACTS_WARNING("warning source link vector is empty");
0178         continue;
0179       }
0180       int ACTS_vol_id = indexSourceLink.geometryId().volume();
0181       int ACTS_lay_id = indexSourceLink.geometryId().layer();
0182       int ACTS_mod_id = indexSourceLink.geometryId().sensitive();
0183
0184       // dont want strips or HGTD
0185       if (ACTS_vol_id == 2 || ACTS_vol_id == 22 || ACTS_vol_id == 23 ||
0186           ACTS_vol_id == 24) {
0187         continue;
0188       }
0189
0190       // Search for vol, lay and module=0, if doesn't esist (end) then search
0191       // for full thing vol*100+lay as first number in pair then 0 or mod id
0192       auto ACTS_joint_id = ACTS_vol_id * 100 + ACTS_lay_id;
0193       auto key = std::make_pair(
0194           ACTS_joint_id,
0195           0);  // here the key needs to be pair of(vol*100+lay, 0)
0196       auto Find = map.find(key);
0197
0198       if (Find ==
0199           map.end()) {  // if end then make new key of (vol*100+lay, modid)
0200         key = std::make_pair(ACTS_joint_id, ACTS_mod_id);  // 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                      << ACTS_vol_id << " and layer id: " << ACTS_lay_id);
0208         continue;
0209       }
0210
0211       // now should be pixel with Gbts ID:
0212       int Gbts_id =
0213           Find->second
0214               .first;  // new map the item is a pair so want first from it
0215
0216       if (Gbts_id == 0) {
0217         ACTS_WARNING("No assigned Gbts ID for key for volume id: "
0218                      << ACTS_vol_id << " and layer id: " << ACTS_lay_id);
0219       }
0220
0221       // access IDs from map
0222       int eta_mod = Find->second.second;
0223       int combined_id = Gbts_id * 1000 + eta_mod;
0224
0225       // fill Gbts vector with current sapce point and ID
0226       gbtsSpacePoints.emplace_back(&spacePoint, Gbts_id, combined_id);
0227     }
0228   }
0229   ACTS_VERBOSE("Space points successfully assigned Gbts ID");
0230
0231   return gbtsSpacePoints;
0232 }
0233
0234 std::vector<Acts::TrigInDetSiLayer>
0235 ActsExamples::GbtsSeedingAlgorithm::LayerNumbering() const {
0236   std::vector<Acts::TrigInDetSiLayer> input_vector;
0237   std::vector<std::size_t> count_vector;
0238
0239   m_cfg.trackingGeometry->visitSurfaces([this, &input_vector, &count_vector](
0240                                             const Acts::Surface *surface) {
0241     Acts::GeometryIdentifier geoid = surface->geometryId();
0242     auto ACTS_vol_id = geoid.volume();
0243     auto ACTS_lay_id = geoid.layer();
0244     auto mod_id = geoid.sensitive();
0245     auto bounds_vect = surface->bounds().values();
0246     auto center = surface->center(Acts::GeometryContext());
0247
0248     // make bounds global
0249     Acts::Vector3 globalFakeMom(1, 1, 1);
0250     Acts::Vector2 min_bound_local =
0251         Acts::Vector2(bounds_vect[0], bounds_vect[1]);
0252     Acts::Vector2 max_bound_local =
0253         Acts::Vector2(bounds_vect[2], bounds_vect[3]);
0254     Acts::Vector3 min_bound_global = surface->localToGlobal(
0255         Acts::GeometryContext(), min_bound_local, globalFakeMom);
0256     Acts::Vector3 max_bound_global = surface->localToGlobal(
0257         Acts::GeometryContext(), max_bound_local, globalFakeMom);
0258
0259     if (min_bound_global(0) >
0260         max_bound_global(0)) {  // checking that not wrong way round
0261       min_bound_global.swap(max_bound_global);
0262     }
0263
0264     float rc = 0.0;
0265     float minBound = 100000.0;
0266     float maxBound = -100000.0;
0267
0268     // convert to Gbts ID
0269     auto ACTS_joint_id = ACTS_vol_id * 100 + ACTS_lay_id;
0270     auto key =
0271         std::make_pair(ACTS_joint_id,
0272                        0);  // here the key needs to be pair of(vol*100+lay, 0)
0273     auto Find = m_cfg.ActsGbtsMap.find(key);
0274     int Gbts_id = 0;               // initialise first to avoid FLTUND later
0275     Gbts_id = Find->second.first;  // new map, item is pair want first
0276     if (Find ==
0277         m_cfg.ActsGbtsMap
0278             .end()) {  // if end then make new key of (vol*100+lay, modid)
0279       key = std::make_pair(ACTS_joint_id, mod_id);  // mod ID
0280       Find = m_cfg.ActsGbtsMap.find(key);
0281       Gbts_id = Find->second.first;
0282     }
0283
0284     short barrel_ec = 0;  // a variable that says if barrrel, 0 = barrel
0285     int eta_mod = Find->second.second;
0286
0287     // assign barrel_ec depending on Gbts_layer
0288     if (79 < Gbts_id && Gbts_id < 85) {  // 80s, barrel
0289       barrel_ec = 0;
0290     } else if (89 < Gbts_id && Gbts_id < 99) {  // 90s positive
0291       barrel_ec = 2;
0292     } else {  // 70s negative
0293       barrel_ec = -2;
0294     }
0295
0296     if (barrel_ec == 0) {
0297       rc = sqrt(center(0) * center(0) +
0298                 center(1) * center(1));  // barrel center in r
0299       // bounds of z
0300       if (min_bound_global(2) < minBound) {
0301         minBound = min_bound_global(2);
0302       }
0303       if (max_bound_global(2) > maxBound) {
0304         maxBound = max_bound_global(2);
0305       }
0306     } else {
0307       rc = center(2);  // not barrel center in Z
0308       // bounds of r
0309       float min = sqrt(min_bound_global(0) * min_bound_global(0) +
0310                        min_bound_global(1) * min_bound_global(1));
0311       float max = sqrt(max_bound_global(0) * max_bound_global(0) +
0312                        max_bound_global(1) * max_bound_global(1));
0313       if (min < minBound) {
0314         minBound = min;
0315       }
0316       if (max > maxBound) {
0317         maxBound = max;
0318       }
0319     }
0320
0321     int combined_id = Gbts_id * 1000 + eta_mod;
0322     auto current_index =
0323         find_if(input_vector.begin(), input_vector.end(),
0324                 [combined_id](auto n) { return n.m_subdet == combined_id; });
0325     if (current_index != input_vector.end()) {  // not end so does exist
0326       std::size_t index = std::distance(input_vector.begin(), current_index);
0327       input_vector[index].m_refCoord += rc;
0328       input_vector[index].m_minBound += minBound;
0329       input_vector[index].m_maxBound += maxBound;
0330       count_vector[index] += 1;  // increase count at the index
0331
0332     } else {  // end so doesn't exists
0333       // make new if one with Gbts ID doesn't exist:
0334       Acts::TrigInDetSiLayer new_Gbts_ID(combined_id, barrel_ec, rc, minBound,
0335                                          maxBound);
0336       input_vector.push_back(new_Gbts_ID);
0337       count_vector.push_back(
0338           1);  // so the element exists and not divinding by 0
0339     }
0340
0341     if (m_cfg.fill_module_csv) {
0342       std::fstream fout;
0343       fout.open("ACTS_modules.csv",
0344                 std::ios::out | std::ios::app);  // add to file each time
0345       // print to csv for each module, no repeats so dont need to make
0346       // map for averaging
0347       fout << ACTS_vol_id << ", "                                  // vol
0348            << ACTS_lay_id << ", "                                  // lay
0349            << mod_id << ", "                                       // module
0350            << Gbts_id << ","                                       // Gbts id
0351            << eta_mod << ","                                       // eta_mod
0352            << center(2) << ", "                                    // z
0353            << sqrt(center(0) * center(0) + center(1) * center(1))  // r
0354            << "\n";
0355     }
0356   });
0357
0358   for (std::size_t i = 0; i < input_vector.size(); i++) {
0359     input_vector[i].m_refCoord = input_vector[i].m_refCoord / count_vector[i];
0360   }
0361
0362   return input_vector;
0363 }