src/Utilities/StripModulePairing.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/Utilities/StripModulePairing.hpp"
0010
0011 #include "Acts/Geometry/GeometryContext.hpp"
0012 #include "ActsExamples/EventData/GeometryContainers.hpp"
0013
0014 ActsExamples::StripModulePairMap ActsExamples::pairStripModules(
0015     const Acts::TrackingGeometry& trackingGeometry,
0016     const std::vector<Acts::GeometryIdentifier>& stripGeometrySelection,
0017     const Acts::Logger& logger) {
0018   StripModulePairMap result;
0019
0020   ACTS_INFO("Strip space point geometry selection:");
0021   for (const auto& geoId : stripGeometrySelection) {
0022     ACTS_INFO("  " << geoId);
0023   }
0024
0025   // We need to use a default geometry context here to access the center
0026   // coordinates of modules.
0027   const auto gctx = Acts::GeometryContext::dangerouslyDefaultConstruct();
0028
0029   // Build strip partner map, i.e., which modules are stereo partners
0030   // As a heuristic we assume that the stereo partners are the modules
0031   // which have the shortest mutual distance
0032   std::vector<const Acts::Surface*> allSensitivesVector;
0033   trackingGeometry.visitSurfaces(
0034       [&](const auto surface) { allSensitivesVector.push_back(surface); },
0035       true);
0036   std::ranges::sort(allSensitivesVector, detail::CompareGeometryId{},
0037                     detail::GeometryIdGetter{});
0038   GeometryIdMultiset<const Acts::Surface*> allSensitives(
0039       allSensitivesVector.begin(), allSensitivesVector.end());
0040
0041   for (auto selector : stripGeometrySelection) {
0042     // Apply volume/layer range
0043     auto rangeLayer =
0044         selectLowestNonZeroGeometryObject(allSensitives, selector);
0045
0046     // Apply selector on extra if extra != 0
0047     auto range = rangeLayer | std::views::filter([&](auto srf) {
0048                    return srf->geometryId().extra() != 0
0049                               ? srf->geometryId().extra() == selector.extra()
0050                               : true;
0051                  });
0052
0053     const auto sizeBefore = result.size();
0054     const std::size_t nSurfaces = std::distance(range.begin(), range.end());
0055
0056     if (nSurfaces < 2) {
0057       ACTS_WARNING("Only " << nSurfaces << " surfaces for selector " << selector
0058                            << ", skip");
0059       continue;
0060     }
0061     ACTS_DEBUG("Found " << nSurfaces << " surfaces for selector " << selector);
0062
0063     // Very dumb all-to-all search
0064     for (auto mod1 : range) {
0065       if (result.contains(mod1->geometryId())) {
0066         continue;
0067       }
0068
0069       const Acts::Surface* partner = nullptr;
0070       double minDist = std::numeric_limits<double>::max();
0071
0072       for (auto mod2 : range) {
0073         if (mod1 == mod2) {
0074           continue;
0075         }
0076         auto c1 = mod1->center(gctx);
0077         auto c2 = mod2->center(gctx);
0078         if (minDist > (c1 - c2).norm()) {
0079           minDist = (c1 - c2).norm();
0080           partner = mod2;
0081         }
0082       }
0083
0084       ACTS_VERBOSE("Found stereo pair: " << mod1->geometryId() << " <-> "
0085                                          << partner->geometryId());
0086       ACTS_VERBOSE("- " << mod1->center(gctx).transpose() << " <-> "
0087                         << partner->center(gctx).transpose());
0088       const auto [it1, success1] =
0089           result.insert({mod1->geometryId(), partner->geometryId()});
0090       const auto [it2, success2] =
0091           result.insert({partner->geometryId(), mod1->geometryId()});
0092       if (!success1 || !success2) {
0093         throw std::runtime_error("error inserting in map");
0094       }
0095     }
0096
0097     const std::size_t sizeAfter = result.size();
0098     const std::size_t missing = nSurfaces - (sizeAfter - sizeBefore);
0099     if (missing > 0) {
0100       ACTS_WARNING("Did not find a stereo partner for " << missing
0101                                                         << " surfaces");
0102     }
0103   }
0104
0105   return result;
0106 }