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