TrackFinding/src/GbtsSeedingAlgorithm.cpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2020 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 http://mozilla.org/MPL/2.0/.
0008
0009 #include "ActsExamples/TrackFinding/GbtsSeedingAlgorithm.hpp"
0010
0011 #include "Acts/Geometry/GeometryIdentifier.hpp"
0012 #include "Acts/Seeding/Seed.hpp"
0013 #include "Acts/Seeding/SeedFilter.hpp"
0014 #include "ActsExamples/EventData/IndexSourceLink.hpp"
0015 #include "ActsExamples/EventData/Measurement.hpp"
0016 #include "ActsExamples/EventData/ProtoTrack.hpp"
0017 #include "ActsExamples/EventData/SimSeed.hpp"
0018 #include "ActsExamples/Framework/WhiteBoard.hpp"
0019
0020 #include <fstream>
0021 #include <iostream>
0022 #include <map>
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_inputSourceLinks.initialize(m_cfg.inputSourceLinks);
0066
0067   m_inputClusters.initialize(m_cfg.inputClusters);
0068
0069   m_cfg.seedFinderConfig.seedFilter =
0070       std::make_unique<Acts::SeedFilter<SimSpacePoint>>(
0071           Acts::SeedFilter<SimSpacePoint>(m_cfg.seedFilterConfig));
0072
0073   // map
0074   m_cfg.ActsGbtsMap = makeActsGbtsMap();
0075   // input trig vector
0076   m_cfg.seedFinderConfig.m_layerGeometry = LayerNumbering();
0077
0078   std::ifstream input_ifstream(
0079       m_cfg.seedFinderConfig.connector_input_file.c_str(), std::ifstream::in);
0080
0081   // connector
0082   std::unique_ptr<Acts::GbtsConnector> inputConnector =
0083       std::make_unique<Acts::GbtsConnector>(input_ifstream);
0084
0085   m_gbtsGeo = std::make_unique<Acts::GbtsGeometry<SimSpacePoint>>(
0086       m_cfg.seedFinderConfig.m_layerGeometry, inputConnector);
0087
0088 }  // this is not Gbts config type because it is a member of the algs config,
0089    // which is of type Gbts cofig
0090
0091 // execute:
0092 ActsExamples::ProcessCode ActsExamples::GbtsSeedingAlgorithm::execute(
0093     const AlgorithmContext &ctx) const {
0094   std::vector<Acts::GbtsSP<SimSpacePoint>> GbtsSpacePoints =
0095       MakeGbtsSpacePoints(ctx, m_cfg.ActsGbtsMap);
0096
0097   // cluster width
0098   //  const ClusterContainer* clusters = &m_inputClusters(ctx) ;
0099
0100   // for (const auto& sp : GbtsSpacePoints){
0101   //   const auto& sl = sp.SP->sourceLinks().front().get<IndexSourceLink>() ;
0102   //   const auto& cluster = clusters->at(sl.index()) ;
0103   //   std::cout << "testing 0: " << cluster.sizeLoc0 << " 1: " <<
0104   //   cluster.sizeLoc1 << std::endl ;
0105
0106   // }
0107
0108   for (auto sp : GbtsSpacePoints) {
0109     ACTS_DEBUG("Gbts space points: "
0110                << " Gbts_id: " << sp.gbtsID << " z: " << sp.SP->z()
0111                << " r: " << sp.SP->r() << " ACTS volume:  "
0112                << sp.SP->sourceLinks()
0113                       .front()
0114                       .get<IndexSourceLink>()
0115                       .geometryId()
0116                       .volume()
0117                << "\n");
0118   }
0119
0120   // this is now calling on a core algorithm
0121   Acts::SeedFinderGbts<SimSpacePoint> finder(m_cfg.seedFinderConfig,
0122                                              *m_gbtsGeo);
0123
0124   // need this function as create_coords is needed for seeds
0125   std::function<std::pair<Acts::Vector3, Acts::Vector2>(
0126       const SimSpacePoint *sp)>
0127       create_coordinates = [](const SimSpacePoint *sp) {
0128         Acts::Vector3 position(sp->x(), sp->y(), sp->z());
0129         Acts::Vector2 variance(sp->varianceR(), sp->varianceZ());
0130         return std::make_pair(position, variance);
0131       };
0132   // output of function needed for seed
0133
0134   finder.loadSpacePoints(GbtsSpacePoints);
0135
0136   // trigGbts file :
0137   Acts::RoiDescriptor internalRoi(0, -4.5, 4.5, 0, -M_PI, M_PI, 0, -150.0,
0138                                   150.0);
0139   // ROI file:
0140   //  Acts::RoiDescriptor internalRoi(0, -5, 5, 0, -M_PI, M_PI, 0, -225.0,
0141   //                                  225.0);
0142
0143   // new version returns seeds
0144   SimSeedContainer seeds = finder.createSeeds(internalRoi, *m_gbtsGeo);
0145
0146   m_outputSeeds(ctx, std::move(seeds));
0147
0148   return ActsExamples::ProcessCode::SUCCESS;
0149 }
0150
0151 std::map<std::pair<int, int>, std::pair<int, int>>
0152 ActsExamples::GbtsSeedingAlgorithm::makeActsGbtsMap() const {
0153   std::map<std::pair<int, int>, std::pair<int, int>> ActsGbts;
0154   std::ifstream data(
0155       m_cfg.layerMappingFile);  // 0 in this file refers to no Gbts ID
0156   std::string line;
0157   std::vector<std::vector<std::string>> parsedCsv;
0158   while (std::getline(data, line)) {
0159     std::stringstream lineStream(line);
0160     std::string cell;
0161     std::vector<std::string> parsedRow;
0162     while (std::getline(lineStream, cell, ',')) {
0163       parsedRow.push_back(cell);
0164     }
0165
0166     parsedCsv.push_back(parsedRow);
0167   }
0168   // file in format ACTS_vol,ACTS_lay,ACTS_mod,Gbts_id
0169   for (auto i : parsedCsv) {
0170     int ACTS_vol = stoi(i[0]);
0171     int ACTS_lay = stoi(i[1]);
0172     int ACTS_mod = stoi(i[2]);
0173     int Gbts = stoi(i[5]);
0174     int eta_mod = stoi(i[6]);
0175     int ACTS_joint = ACTS_vol * 100 + ACTS_lay;
0176     ActsGbts.insert({{ACTS_joint, ACTS_mod}, {Gbts, eta_mod}});
0177   }
0178
0179   return ActsGbts;
0180 }
0181
0182 std::vector<Acts::GbtsSP<ActsExamples::SimSpacePoint>>
0183 ActsExamples::GbtsSeedingAlgorithm::MakeGbtsSpacePoints(
0184     const AlgorithmContext &ctx,
0185     std::map<std::pair<int, int>, std::pair<int, int>> map) const {
0186   // create space point vectors
0187   std::vector<const ActsExamples::SimSpacePoint *> spacePoints;
0188   std::vector<Acts::GbtsSP<ActsExamples::SimSpacePoint>> gbtsSpacePoints;
0189   gbtsSpacePoints.reserve(
0190       m_inputSpacePoints.size());  // not sure if this is enough
0191
0192   // for loop filling space
0193   for (const auto &isp : m_inputSpacePoints) {
0194     for (const auto &spacePoint : (*isp)(ctx)) {
0195       // fill original space point vector
0196       spacePoints.push_back(&spacePoint);
0197
0198       // Gbts space point vector
0199       // loop over space points, call on map
0200       const auto &source_link = spacePoint.sourceLinks();
0201       const auto &index_source_link =
0202           source_link.front().get<IndexSourceLink>();
0203
0204       // warning if source link empty
0205       if (source_link.empty()) {
0206         // warning in officaial acts format
0207         ACTS_WARNING("warning source link vector is empty");
0208         continue;
0209       }
0210       int ACTS_vol_id = index_source_link.geometryId().volume();
0211       int ACTS_lay_id = index_source_link.geometryId().layer();
0212       int ACTS_mod_id = index_source_link.geometryId().sensitive();
0213
0214       // dont want strips or HGTD
0215       if (ACTS_vol_id == 2 || ACTS_vol_id == 22 || ACTS_vol_id == 23 ||
0216           ACTS_vol_id == 24) {
0217         continue;
0218       }
0219
0220       // Search for vol, lay and module=0, if doesn't esist (end) then search
0221       // for full thing vol*100+lay as first number in pair then 0 or mod id
0222       auto ACTS_joint_id = ACTS_vol_id * 100 + ACTS_lay_id;
0223       auto key = std::make_pair(
0224           ACTS_joint_id,
0225           0);  // here the key needs to be pair of(vol*100+lay, 0)
0226       auto Find = map.find(key);
0227
0228       if (Find ==
0229           map.end()) {  // if end then make new key of (vol*100+lay, modid)
0230         key = std::make_pair(ACTS_joint_id, ACTS_mod_id);  // mod ID
0231         Find = map.find(key);
0232       }
0233
0234       // warning if key not in map
0235       if (Find == map.end()) {
0236         ACTS_WARNING("Key not found in Gbts map for volume id: "
0237                      << ACTS_vol_id << " and layer id: " << ACTS_lay_id);
0238         continue;
0239       }
0240
0241       // now should be pixel with Gbts ID:
0242       int Gbts_id =
0243           Find->second
0244               .first;  // new map the item is a pair so want first from it
0245
0246       if (Gbts_id == 0) {
0247         ACTS_WARNING("No assigned Gbts ID for key for volume id: "
0248                      << ACTS_vol_id << " and layer id: " << ACTS_lay_id);
0249       }
0250
0251       // access IDs from map
0252       int eta_mod = Find->second.second;
0253       int combined_id = Gbts_id * 1000 + eta_mod;
0254
0255       // fill Gbts vector with current sapce point and ID
0256       gbtsSpacePoints.emplace_back(&spacePoint, Gbts_id, combined_id);
0257     }
0258   }
0259   ACTS_VERBOSE("Space points successfully assigned Gbts ID");
0260
0261   return gbtsSpacePoints;
0262 }
0263
0264 std::vector<Acts::TrigInDetSiLayer>
0265 ActsExamples::GbtsSeedingAlgorithm::LayerNumbering() const {
0266   std::vector<Acts::TrigInDetSiLayer> input_vector;
0267   std::vector<std::size_t> count_vector;
0268
0269   m_cfg.trackingGeometry->visitSurfaces([this, &input_vector, &count_vector](
0270                                             const Acts::Surface *surface) {
0271     Acts::GeometryIdentifier geoid = surface->geometryId();
0272     auto ACTS_vol_id = geoid.volume();
0273     auto ACTS_lay_id = geoid.layer();
0274     auto mod_id = geoid.sensitive();
0275     auto bounds_vect = surface->bounds().values();
0276     auto center = surface->center(Acts::GeometryContext());
0277
0278     // make bounds global
0279     Acts::Vector3 globalFakeMom(1, 1, 1);
0280     Acts::Vector2 min_bound_local =
0281         Acts::Vector2(bounds_vect[0], bounds_vect[1]);
0282     Acts::Vector2 max_bound_local =
0283         Acts::Vector2(bounds_vect[2], bounds_vect[3]);
0284     Acts::Vector3 min_bound_global = surface->localToGlobal(
0285         Acts::GeometryContext(), min_bound_local, globalFakeMom);
0286     Acts::Vector3 max_bound_global = surface->localToGlobal(
0287         Acts::GeometryContext(), max_bound_local, globalFakeMom);
0288
0289     if (min_bound_global(0) >
0290         max_bound_global(0)) {  // checking that not wrong way round
0291       min_bound_global.swap(max_bound_global);
0292     }
0293
0294     float rc = 0.0;
0295     float minBound = 100000.0;
0296     float maxBound = -100000.0;
0297
0298     // convert to Gbts ID
0299     auto ACTS_joint_id = ACTS_vol_id * 100 + ACTS_lay_id;
0300     auto key =
0301         std::make_pair(ACTS_joint_id,
0302                        0);  // here the key needs to be pair of(vol*100+lay, 0)
0303     auto Find = m_cfg.ActsGbtsMap.find(key);
0304     int Gbts_id = 0;               // initialise first to avoid FLTUND later
0305     Gbts_id = Find->second.first;  // new map, item is pair want first
0306     if (Find ==
0307         m_cfg.ActsGbtsMap
0308             .end()) {  // if end then make new key of (vol*100+lay, modid)
0309       key = std::make_pair(ACTS_joint_id, mod_id);  // mod ID
0310       Find = m_cfg.ActsGbtsMap.find(key);
0311       Gbts_id = Find->second.first;
0312     }
0313
0314     short barrel_ec = 0;  // a variable that says if barrrel, 0 = barrel
0315     int eta_mod = Find->second.second;
0316
0317     // assign barrel_ec depending on Gbts_layer
0318     if (79 < Gbts_id && Gbts_id < 85) {  // 80s, barrel
0319       barrel_ec = 0;
0320     } else if (89 < Gbts_id && Gbts_id < 99) {  // 90s positive
0321       barrel_ec = 2;
0322     } else {  // 70s negative
0323       barrel_ec = -2;
0324     }
0325
0326     if (barrel_ec == 0) {
0327       rc = sqrt(center(0) * center(0) +
0328                 center(1) * center(1));  // barrel center in r
0329       // bounds of z
0330       if (min_bound_global(2) < minBound) {
0331         minBound = min_bound_global(2);
0332       }
0333       if (max_bound_global(2) > maxBound) {
0334         maxBound = max_bound_global(2);
0335       }
0336     } else {
0337       rc = center(2);  // not barrel center in Z
0338       // bounds of r
0339       float min = sqrt(min_bound_global(0) * min_bound_global(0) +
0340                        min_bound_global(1) * min_bound_global(1));
0341       float max = sqrt(max_bound_global(0) * max_bound_global(0) +
0342                        max_bound_global(1) * max_bound_global(1));
0343       if (min < minBound) {
0344         minBound = min;
0345       }
0346       if (max > maxBound) {
0347         maxBound = max;
0348       }
0349     }
0350
0351     int combined_id = Gbts_id * 1000 + eta_mod;
0352     auto current_index =
0353         find_if(input_vector.begin(), input_vector.end(),
0354                 [combined_id](auto n) { return n.m_subdet == combined_id; });
0355     if (current_index != input_vector.end()) {  // not end so does exist
0356       std::size_t index = std::distance(input_vector.begin(), current_index);
0357       input_vector[index].m_refCoord += rc;
0358       input_vector[index].m_minBound += minBound;
0359       input_vector[index].m_maxBound += maxBound;
0360       count_vector[index] += 1;  // increase count at the index
0361
0362     } else {  // end so doesn't exists
0363       // make new if one with Gbts ID doesn't exist:
0364       Acts::TrigInDetSiLayer new_Gbts_ID(combined_id, barrel_ec, rc, minBound,
0365                                          maxBound);
0366       input_vector.push_back(new_Gbts_ID);
0367       count_vector.push_back(
0368           1);  // so the element exists and not divinding by 0
0369     }
0370
0371     if (m_cfg.fill_module_csv) {
0372       std::fstream fout;
0373       fout.open("ACTS_modules.csv",
0374                 std::ios::out | std::ios::app);  // add to file each time
0375       // print to csv for each module, no repeats so dont need to make
0376       // map for averaging
0377       fout << ACTS_vol_id << ", "                                  // vol
0378            << ACTS_lay_id << ", "                                  // lay
0379            << mod_id << ", "                                       // module
0380            << Gbts_id << ","                                       // Gbts id
0381            << eta_mod << ","                                       // eta_mod
0382            << center(2) << ", "                                    // z
0383            << sqrt(center(0) * center(0) + center(1) * center(1))  // r
0384            << "\n";
0385     }
0386   });
0387
0388   for (long unsigned int i = 0; i < input_vector.size(); i++) {
0389     input_vector[i].m_refCoord = input_vector[i].m_refCoord / count_vector[i];
0390   }
0391
0392   return input_vector;
0393 }