Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:09:44 

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