Digitization/src/DigitizationConfigurator.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/Digitization/DigitizationConfigurator.hpp"
0010
0011 #include "Acts/Surfaces/AnnulusBounds.hpp"
0012 #include "Acts/Surfaces/DiscTrapezoidBounds.hpp"
0013 #include "Acts/Surfaces/RadialBounds.hpp"
0014 #include "Acts/Surfaces/RectangleBounds.hpp"
0015 #include "Acts/Surfaces/Surface.hpp"
0016 #include "Acts/Surfaces/SurfaceBounds.hpp"
0017 #include "Acts/Surfaces/TrapezoidBounds.hpp"
0018 #include "Acts/Utilities/IAxis.hpp"
0019 #include "Acts/Utilities/IMultiAxis.hpp"
0020 #include "Acts/Utilities/Zip.hpp"
0021 #include "ActsExamples/Digitization/SmearingConfig.hpp"
0022 #include "ActsFatras/Digitization/UncorrelatedHitSmearer.hpp"
0023
0024 #include <algorithm>
0025 #include <cmath>
0026 #include <memory>
0027 #include <stdexcept>
0028 #include <vector>
0029
0030 namespace ActsExamples {
0031
0032 namespace {
0033
0034 /// @note This does not really compare if the configs are equal, therefore
0035 /// it is no operator==. The contained std::function types cannot really
0036 /// be checked for equality.
0037 bool segmentationMaybeEqual(const std::shared_ptr<const Acts::IMultiAxis> &a,
0038                             const std::shared_ptr<const Acts::IMultiAxis> &b) {
0039   if (a == b) {
0040     return true;
0041   }
0042   if (a == nullptr || b == nullptr) {
0043     return false;
0044   }
0045   // The axis equality does not cover the axis directions
0046   return *a == *b &&
0047          std::ranges::equal(*a, *b, [](const auto &aa, const auto &ba) {
0048            return aa.getDirection() == ba.getDirection();
0049          });
0050 }
0051
0052 bool digiConfigMaybeEqual(DigiComponentsConfig &a, DigiComponentsConfig &b) {
0053   // Check smearing config
0054   for (const auto &[as, bs] :
0055        Acts::zip(a.smearingDigiConfig.params, b.smearingDigiConfig.params)) {
0056     if (as.index != bs.index) {
0057       return false;
0058     }
0059   }
0060   if (a.smearingDigiConfig.maxRetries != b.smearingDigiConfig.maxRetries) {
0061     return false;
0062   }
0063   // Check geometric config
0064   const auto &ag = a.geometricDigiConfig;
0065   const auto &bg = b.geometricDigiConfig;
0066   return (ag.indices == bg.indices &&
0067           segmentationMaybeEqual(ag.segmentation, bg.segmentation) &&
0068           ag.thickness == bg.thickness && ag.threshold == bg.threshold &&
0069           ag.digital == bg.digital);
0070 }
0071
0072 }  // namespace
0073
0074 void DigitizationConfigurator::operator()(const Acts::Surface *surface) {
0075   if (!surface->isSensitive()) {
0076     return;
0077   }
0078
0079   Acts::GeometryIdentifier geoId = surface->geometryId();
0080   const auto dInputConfig = inputDigiComponents.find(geoId);
0081   if (dInputConfig == inputDigiComponents.end()) {
0082     return;
0083   }
0084
0085   // The output config, copy over the smearing part
0086   DigiComponentsConfig dOutputConfig;
0087   dOutputConfig.smearingDigiConfig = dInputConfig->smearingDigiConfig;
0088
0089   if (!dInputConfig->geometricDigiConfig.indices.empty()) {
0090     // Copy over what can be done
0091     dOutputConfig.geometricDigiConfig.indices =
0092         dInputConfig->geometricDigiConfig.indices;
0093     dOutputConfig.geometricDigiConfig.thickness =
0094         dInputConfig->geometricDigiConfig.thickness;
0095     dOutputConfig.geometricDigiConfig.chargeSmearer =
0096         dInputConfig->geometricDigiConfig.chargeSmearer;
0097     dOutputConfig.geometricDigiConfig.digital =
0098         dInputConfig->geometricDigiConfig.digital;
0099
0100     const Acts::SurfaceBounds &sBounds = surface->bounds();
0101     const std::vector<double> boundValues = sBounds.values();
0102
0103     const auto &inputSegmentation =
0104         dInputConfig->geometricDigiConfig.segmentation;
0105     if (inputSegmentation == nullptr) {
0106       throw std::invalid_argument(
0107           "DigitizationConfigurator: geometric digitization configured "
0108           "without a segmentation");
0109     }
0110     std::vector<std::unique_ptr<const Acts::IAxis>> outputAxes;
0111
0112     switch (sBounds.type()) {
0113       // The module is a rectangle module
0114       case Acts::SurfaceBounds::eRectangle: {
0115         if (inputSegmentation->getAxis(0).getDirection() ==
0116             Acts::AxisDirection::AxisX) {
0117           const double minX = boundValues[Acts::RectangleBounds::eMinX];
0118           const double maxX = boundValues[Acts::RectangleBounds::eMaxX];
0119
0120           const unsigned int nBins = inputSegmentation->getAxis(0).getNBins();
0121
0122           outputAxes.push_back(Acts::IAxis::createEquidistant(
0123               Acts::AxisBoundaryType::Bound, minX, maxX, nBins,
0124               Acts::AxisDirection::AxisX));
0125         }
0126         if (inputSegmentation->getAxis(0).getDirection() ==
0127                 Acts::AxisDirection::AxisY ||
0128             inputSegmentation->getNAxes() == 2) {
0129           const unsigned int accessBin =
0130               inputSegmentation->getNAxes() == 2 ? 1 : 0;
0131
0132           const double minY = boundValues[Acts::RectangleBounds::eMinY];
0133           const double maxY = boundValues[Acts::RectangleBounds::eMaxY];
0134
0135           const unsigned int nBins =
0136               inputSegmentation->getAxis(accessBin).getNBins();
0137
0138           outputAxes.push_back(Acts::IAxis::createEquidistant(
0139               Acts::AxisBoundaryType::Bound, minY, maxY, nBins,
0140               Acts::AxisDirection::AxisY));
0141         }
0142       } break;
0143
0144       // The module is a trapezoid module
0145       case Acts::SurfaceBounds::eTrapezoid: {
0146         if (inputSegmentation->getAxis(0).getDirection() ==
0147             Acts::AxisDirection::AxisX) {
0148           const double maxX =
0149               std::max(boundValues[Acts::TrapezoidBounds::eHalfLengthXnegY],
0150                        boundValues[Acts::TrapezoidBounds::eHalfLengthXposY]);
0151
0152           const unsigned int nBins = inputSegmentation->getAxis(0).getNBins();
0153
0154           outputAxes.push_back(Acts::IAxis::createEquidistant(
0155               Acts::AxisBoundaryType::Bound, -maxX, maxX, nBins,
0156               Acts::AxisDirection::AxisX));
0157         }
0158         if (inputSegmentation->getAxis(0).getDirection() ==
0159                 Acts::AxisDirection::AxisY ||
0160             inputSegmentation->getNAxes() == 2) {
0161           const unsigned int accessBin =
0162               inputSegmentation->getNAxes() == 2 ? 1 : 0;
0163
0164           const double maxY = boundValues[Acts::TrapezoidBounds::eHalfLengthY];
0165
0166           const unsigned int nBins =
0167               inputSegmentation->getAxis(accessBin).getNBins();
0168
0169           outputAxes.push_back(Acts::IAxis::createEquidistant(
0170               Acts::AxisBoundaryType::Bound, -maxY, maxY, nBins,
0171               Acts::AxisDirection::AxisY));
0172         }
0173       } break;
0174
0175       // The module is an annulus module
0176       case Acts::SurfaceBounds::eAnnulus: {
0177         if (inputSegmentation->getAxis(0).getDirection() ==
0178             Acts::AxisDirection::AxisR) {
0179           const double minR = boundValues[Acts::AnnulusBounds::eMinR];
0180           const double maxR = boundValues[Acts::AnnulusBounds::eMaxR];
0181
0182           const unsigned int nBins = inputSegmentation->getAxis(0).getNBins();
0183
0184           outputAxes.push_back(Acts::IAxis::createEquidistant(
0185               Acts::AxisBoundaryType::Bound, minR, maxR, nBins,
0186               Acts::AxisDirection::AxisR));
0187         }
0188         if (inputSegmentation->getAxis(0).getDirection() ==
0189                 Acts::AxisDirection::AxisPhi ||
0190             inputSegmentation->getNAxes() == 2) {
0191           const double averagePhi =
0192               boundValues[Acts::AnnulusBounds::eAveragePhi];
0193           const double minPhi =
0194               averagePhi - boundValues[Acts::AnnulusBounds::eMinPhiRel];
0195           const double maxPhi =
0196               averagePhi + boundValues[Acts::AnnulusBounds::eMaxPhiRel];
0197
0198           const unsigned int nBins = inputSegmentation->getAxis(0).getNBins();
0199
0200           outputAxes.push_back(Acts::IAxis::createEquidistant(
0201               Acts::AxisBoundaryType::Bound, minPhi, maxPhi, nBins,
0202               Acts::AxisDirection::AxisPhi));
0203         }
0204       } break;
0205
0206       // The module is a Disc Trapezoid
0207       case Acts::SurfaceBounds::eDiscTrapezoid: {
0208         const double minR = boundValues[Acts::DiscTrapezoidBounds::eMinR];
0209         const double maxR = boundValues[Acts::DiscTrapezoidBounds::eMaxR];
0210
0211         if (inputSegmentation->getAxis(0).getDirection() ==
0212             Acts::AxisDirection::AxisR) {
0213           const unsigned int nBins = inputSegmentation->getAxis(0).getNBins();
0214
0215           outputAxes.push_back(Acts::IAxis::createEquidistant(
0216               Acts::AxisBoundaryType::Bound, minR, maxR, nBins,
0217               Acts::AxisDirection::AxisR));
0218         }
0219         if (inputSegmentation->getAxis(0).getDirection() ==
0220                 Acts::AxisDirection::AxisPhi ||
0221             inputSegmentation->getNAxes() == 2) {
0222           const unsigned int accessBin =
0223               inputSegmentation->getNAxes() == 2 ? 1 : 0;
0224
0225           const double hxMinR =
0226               boundValues[Acts::DiscTrapezoidBounds::eHalfLengthXminR];
0227           const double hxMaxR =
0228               boundValues[Acts::DiscTrapezoidBounds::eHalfLengthXmaxR];
0229
0230           const double averagePhi =
0231               boundValues[Acts::DiscTrapezoidBounds::eAveragePhi];
0232           const double alphaMinR = std::atan2(minR, hxMinR);
0233           const double alphaMaxR = std::atan2(maxR, hxMaxR);
0234           const double alpha = std::max(alphaMinR, alphaMaxR);
0235
0236           const unsigned int nBins =
0237               inputSegmentation->getAxis(accessBin).getNBins();
0238
0239           outputAxes.push_back(Acts::IAxis::createEquidistant(
0240               Acts::AxisBoundaryType::Bound, averagePhi - alpha,
0241               averagePhi + alpha, nBins, Acts::AxisDirection::AxisPhi));
0242         }
0243       } break;
0244
0245       case Acts::SurfaceBounds::eDisc: {
0246         if (inputSegmentation->getAxis(0).getDirection() ==
0247             Acts::AxisDirection::AxisR) {
0248           const double minR = boundValues[Acts::RadialBounds::eMinR];
0249           const double maxR = boundValues[Acts::RadialBounds::eMaxR];
0250
0251           const unsigned int nBins = inputSegmentation->getAxis(0).getNBins();
0252
0253           outputAxes.push_back(Acts::IAxis::createEquidistant(
0254               Acts::AxisBoundaryType::Bound, minR, maxR, nBins,
0255               Acts::AxisDirection::AxisR));
0256         }
0257         if (inputSegmentation->getAxis(0).getDirection() ==
0258                 Acts::AxisDirection::AxisPhi ||
0259             inputSegmentation->getNAxes() == 2) {
0260           const unsigned int accessBin =
0261               inputSegmentation->getNAxes() == 2 ? 1 : 0;
0262
0263           const double averagePhi =
0264               boundValues[Acts::RadialBounds::eAveragePhi];
0265           const double halfPhiSector =
0266               boundValues[Acts::RadialBounds::eHalfPhiSector];
0267           const double minPhi = averagePhi - halfPhiSector;
0268           const double maxPhi = averagePhi + halfPhiSector;
0269
0270           const unsigned int nBins =
0271               inputSegmentation->getAxis(accessBin).getNBins();
0272
0273           outputAxes.push_back(Acts::IAxis::createEquidistant(
0274               Acts::AxisBoundaryType::Bound, minPhi, maxPhi, nBins,
0275               Acts::AxisDirection::AxisPhi));
0276         }
0277       } break;
0278
0279       default:
0280         break;
0281     }
0282
0283     // Set the adapted segmentation class
0284     if (outputAxes.size() == 1) {
0285       dOutputConfig.geometricDigiConfig.segmentation =
0286           Acts::IMultiAxis::create(*outputAxes[0]);
0287     } else if (outputAxes.size() == 2) {
0288       dOutputConfig.geometricDigiConfig.segmentation =
0289           Acts::IMultiAxis::create(*outputAxes[0], *outputAxes[1]);
0290     }
0291   }
0292
0293   // Compactify the output map where possible
0294   if (compactify) {
0295     // Check for a representing volume configuration, insert if not
0296     // present
0297     Acts::GeometryIdentifier volGeoId =
0298         Acts::GeometryIdentifier().withVolume(geoId.volume());
0299
0300     auto volRep = volumeLayerComponents.find(volGeoId);
0301     if (volRep != volumeLayerComponents.end() &&
0302         digiConfigMaybeEqual(dOutputConfig, volRep->second)) {
0303       // return if the volume representation already covers this one
0304       return;
0305     }
0306     volumeLayerComponents[volGeoId] = dOutputConfig;
0307     outputDigiComponents.push_back({volGeoId, dOutputConfig});
0308
0309     // Check for a representing layer configuration, insert if not present
0310     Acts::GeometryIdentifier volLayGeoId =
0311         Acts::GeometryIdentifier(volGeoId).withLayer(geoId.layer());
0312     const auto volLayRep = volumeLayerComponents.find(volLayGeoId);
0313     if (volLayRep != volumeLayerComponents.end() &&
0314         digiConfigMaybeEqual(dOutputConfig, volLayRep->second)) {
0315       return;
0316     }
0317     volumeLayerComponents[volLayGeoId] = dOutputConfig;
0318     outputDigiComponents.push_back({volLayGeoId, dOutputConfig});
0319   }
0320
0321   // Insert into the output list
0322   outputDigiComponents.push_back({geoId, dOutputConfig});
0323 }
0324
0325 }  // namespace ActsExamples