Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-07-15 08:12:31 

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

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!