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