Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-18 09:11:35 

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