Core/Seeding/SeedFinderTest.cpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2024 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 "Acts/Definitions/Algebra.hpp"
0010 #include "Acts/Definitions/Units.hpp"
0011 #include "Acts/Geometry/Extent.hpp"
0012 #include "Acts/Seeding/BinnedGroup.hpp"
0013 #include "Acts/Seeding/Seed.hpp"
0014 #include "Acts/Seeding/SeedFilter.hpp"
0015 #include "Acts/Seeding/SeedFilterConfig.hpp"
0016 #include "Acts/Seeding/SeedFinder.hpp"
0017 #include "Acts/Seeding/SeedFinderConfig.hpp"
0018 #include "Acts/Seeding/SpacePointGrid.hpp"
0019 #include "Acts/Utilities/Grid.hpp"
0020 #include "Acts/Utilities/GridBinFinder.hpp"
0021
0022 #include <algorithm>
0023 #include <chrono>
0024 #include <cmath>
0025 #include <cstdlib>
0026 #include <fstream>
0027 #include <iostream>
0028 #include <iterator>
0029 #include <memory>
0030 #include <string>
0031 #include <tuple>
0032 #include <utility>
0033 #include <vector>
0034
0035 #include "ATLASCuts.hpp"
0036 #include "SpacePoint.hpp"
0037
0038 using namespace Acts::UnitLiterals;
0039
0040 std::vector<const SpacePoint*> readFile(const std::string& filename) {
0041   std::string line;
0042   int layer = 0;
0043   std::vector<const SpacePoint*> readSP;
0044
0045   std::ifstream spFile(filename);
0046   if (spFile.is_open()) {
0047     while (!spFile.eof()) {
0048       std::getline(spFile, line);
0049       std::stringstream ss(line);
0050       std::string linetype;
0051       ss >> linetype;
0052       if (linetype == "lxyz") {
0053         float x = 0, y = 0, z = 0, varianceR = 0, varianceZ = 0;
0054         std::optional<float> t, varianceT;
0055         ss >> layer >> x >> y >> z >> varianceR >> varianceZ;
0056         const float r = std::hypot(x, y);
0057
0058         float cov = varianceZ * varianceZ * .08333;
0059         if (cov < varianceR) {
0060           cov = varianceR;
0061         }
0062
0063         if (std::abs(z) > 450.) {
0064           varianceZ = 9. * cov;
0065           varianceR = .06;
0066         } else {
0067           varianceR = 9. * cov;
0068           varianceZ = .06;
0069         }
0070
0071         SpacePoint* sp = new SpacePoint{
0072             x, y, z, r, layer, varianceR, varianceZ, t, varianceT};
0073         //     if(r < 200.){
0074         //       sp->setClusterList(1,0);
0075         //     }
0076         readSP.push_back(sp);
0077       }
0078     }
0079   }
0080   return readSP;
0081 }
0082
0083 int main(int argc, char** argv) {
0084   std::string file{"sp.txt"};
0085   bool help(false);
0086   bool quiet(false);
0087
0088   int opt = -1;
0089   while ((opt = getopt(argc, argv, "hf:q")) != -1) {
0090     switch (opt) {
0091       case 'f':
0092         file = optarg;
0093         break;
0094       case 'q':
0095         quiet = true;
0096         break;
0097       case 'h':
0098         help = true;
0099         [[fallthrough]];
0100       default: /* '?' */
0101         std::cerr << "Usage: " << argv[0] << " [-hq] [-f FILENAME]\n";
0102         if (help) {
0103           std::cout << "      -h : this help" << std::endl;
0104           std::cout
0105               << "      -f FILE : read spacepoints from FILE. Default is \""
0106               << file << "\"" << std::endl;
0107           std::cout << "      -q : don't print out all found seeds"
0108                     << std::endl;
0109         }
0110
0111         exit(EXIT_FAILURE);
0112     }
0113   }
0114
0115   std::ifstream f(file);
0116   if (!f.good()) {
0117     std::cerr << "input file \"" << file << "\" does not exist\n";
0118     exit(EXIT_FAILURE);
0119   }
0120
0121   auto start_read = std::chrono::system_clock::now();
0122   std::vector<const SpacePoint*> spVec = readFile(file);
0123   auto end_read = std::chrono::system_clock::now();
0124   std::chrono::duration<double> elapsed_read = end_read - start_read;
0125
0126   std::cout << "read " << spVec.size() << " SP from file " << file << " in "
0127             << elapsed_read.count() << "s" << std::endl;
0128
0129   Acts::SeedFinderConfig<SpacePoint> config;
0130   // silicon detector max
0131   config.rMax = 160._mm;
0132   config.deltaRMin = 5._mm;
0133   config.deltaRMax = 160._mm;
0134   config.deltaRMinTopSP = config.deltaRMin;
0135   config.deltaRMinBottomSP = config.deltaRMin;
0136   config.deltaRMaxTopSP = config.deltaRMax;
0137   config.deltaRMaxBottomSP = config.deltaRMax;
0138   config.collisionRegionMin = -250._mm;
0139   config.collisionRegionMax = 250._mm;
0140   config.zMin = -2800._mm;
0141   config.zMax = 2800._mm;
0142   config.maxSeedsPerSpM = 5;
0143   // 2.7 eta
0144   config.cotThetaMax = 7.40627;
0145   config.sigmaScattering = 1.00000;
0146
0147   config.minPt = 500._MeV;
0148
0149   config.impactMax = 10._mm;
0150
0151   config.useVariableMiddleSPRange = false;
0152
0153   Acts::SeedFinderOptions options;
0154   options.beamPos = {-.5_mm, -.5_mm};
0155   options.bFieldInZ = 2_T;
0156
0157   int numPhiNeighbors = 1;
0158
0159   // extent used to store r range for middle spacepoint
0160   Acts::Extent rRangeSPExtent;
0161
0162   config.useVariableMiddleSPRange = false;
0163   const Acts::Range1D<float> rMiddleSPRange;
0164
0165   std::vector<std::pair<int, int>> zBinNeighborsTop;
0166   std::vector<std::pair<int, int>> zBinNeighborsBottom;
0167
0168   auto bottomBinFinder = std::make_unique<Acts::GridBinFinder<2ul>>(
0169       Acts::GridBinFinder<2ul>(numPhiNeighbors, zBinNeighborsBottom));
0170   auto topBinFinder = std::make_unique<Acts::GridBinFinder<2ul>>(
0171       Acts::GridBinFinder<2ul>(numPhiNeighbors, zBinNeighborsTop));
0172   Acts::SeedFilterConfig sfconf;
0173   Acts::ATLASCuts<SpacePoint> atlasCuts = Acts::ATLASCuts<SpacePoint>();
0174   config.seedFilter = std::make_unique<Acts::SeedFilter<SpacePoint>>(
0175       Acts::SeedFilter<SpacePoint>(sfconf, &atlasCuts));
0176   Acts::SeedFinder<SpacePoint, Acts::CylindricalSpacePointGrid<SpacePoint>>
0177       a;  // test creation of unconfigured finder
0178   a = Acts::SeedFinder<SpacePoint, Acts::CylindricalSpacePointGrid<SpacePoint>>(
0179       config);
0180
0181   // covariance tool, sets covariances per spacepoint as required
0182   auto ct = [=](const SpacePoint& sp, float, float, float) {
0183     Acts::Vector3 position(sp.x(), sp.y(), sp.z());
0184     Acts::Vector2 covariance(sp.varianceR, sp.varianceZ);
0185     return std::make_tuple(position, covariance, sp.t());
0186   };
0187
0188   // setup spacepoint grid config
0189   Acts::CylindricalSpacePointGridConfig gridConf;
0190   gridConf.minPt = config.minPt;
0191   gridConf.rMax = config.rMax;
0192   gridConf.zMax = config.zMax;
0193   gridConf.zMin = config.zMin;
0194   gridConf.deltaRMax = config.deltaRMax;
0195   gridConf.cotThetaMax = config.cotThetaMax;
0196   // setup spacepoint grid options
0197   Acts::CylindricalSpacePointGridOptions gridOpts;
0198   gridOpts.bFieldInZ = options.bFieldInZ;
0199   // create grid with bin sizes according to the configured geometry
0200
0201   Acts::CylindricalSpacePointGrid<SpacePoint> grid =
0202       Acts::CylindricalSpacePointGridCreator::createGrid<SpacePoint>(gridConf,
0203                                                                      gridOpts);
0204   Acts::CylindricalSpacePointGridCreator::fillGrid(
0205       config, options, grid, spVec.begin(), spVec.end(), ct, rRangeSPExtent);
0206   auto spGroup = Acts::CylindricalBinnedGroup<SpacePoint>(
0207       std::move(grid), *bottomBinFinder, *topBinFinder);
0208
0209   std::vector<std::vector<Acts::Seed<SpacePoint>>> seedVector;
0210   decltype(a)::SeedingState state;
0211   auto start = std::chrono::system_clock::now();
0212   for (auto [bottom, middle, top] : spGroup) {
0213     auto& v = seedVector.emplace_back();
0214     a.createSeedsForGroup(options, state, spGroup.grid(), std::back_inserter(v),
0215                           bottom, middle, top, rMiddleSPRange);
0216   }
0217   auto end = std::chrono::system_clock::now();
0218   std::chrono::duration<double> elapsed_seconds = end - start;
0219   std::cout << "time to create seeds: " << elapsed_seconds.count() << std::endl;
0220   std::cout << "Number of regions: " << seedVector.size() << std::endl;
0221   int numSeeds = 0;
0222   for (auto& outVec : seedVector) {
0223     numSeeds += outVec.size();
0224   }
0225   std::cout << "Number of seeds generated: " << numSeeds << std::endl;
0226   if (!quiet) {
0227     for (auto& regionVec : seedVector) {
0228       for (std::size_t i = 0; i < regionVec.size(); i++) {
0229         const Acts::Seed<SpacePoint>* seed = &regionVec[i];
0230         const SpacePoint* sp = seed->sp()[0];
0231         std::cout << " (" << sp->x() << ", " << sp->y() << ", " << sp->z()
0232                   << ") ";
0233         sp = seed->sp()[1];
0234         std::cout << sp->layer << " (" << sp->x() << ", " << sp->y() << ", "
0235                   << sp->z() << ") ";
0236         sp = seed->sp()[2];
0237         std::cout << sp->layer << " (" << sp->x() << ", " << sp->y() << ", "
0238                   << sp->z() << ") ";
0239         std::cout << std::endl;
0240       }
0241     }
0242   }
0243   return 0;
0244 }