Digitization/src/DigitizationAlgorithm.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/DigitizationAlgorithm.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/TrackParametrization.hpp"
0013 #include "Acts/Geometry/GeometryIdentifier.hpp"
0014 #include "Acts/Utilities/BinUtility.hpp"
0015 #include "ActsExamples/Digitization/ModuleClusters.hpp"
0016 #include "ActsExamples/EventData/GeometryContainers.hpp"
0017 #include "ActsExamples/EventData/Index.hpp"
0018 #include "ActsExamples/Framework/AlgorithmContext.hpp"
0019
0020 #include <algorithm>
0021 #include <array>
0022 #include <limits>
0023 #include <ostream>
0024 #include <stdexcept>
0025 #include <string>
0026 #include <utility>
0027
0028 namespace ActsExamples {
0029
0030 DigitizationAlgorithm::DigitizationAlgorithm(Config config,
0031                                              Acts::Logging::Level level)
0032     : IAlgorithm("DigitizationAlgorithm", level), m_cfg(std::move(config)) {
0033   if (m_cfg.inputSimHits.empty()) {
0034     throw std::invalid_argument("Missing simulated hits input collection");
0035   }
0036   if (m_cfg.surfaceByIdentifier.empty()) {
0037     throw std::invalid_argument("Missing Surface-GeometryID association map");
0038   }
0039   if (!m_cfg.randomNumbers) {
0040     throw std::invalid_argument("Missing random numbers tool");
0041   }
0042   if (m_cfg.digitizationConfigs.empty()) {
0043     throw std::invalid_argument("Missing digitization configuration");
0044   }
0045
0046   if (m_cfg.doClusterization) {
0047     if (m_cfg.outputMeasurements.empty()) {
0048       throw std::invalid_argument("Missing measurements output collection");
0049     }
0050     if (m_cfg.outputClusters.empty()) {
0051       throw std::invalid_argument("Missing cluster output collection");
0052     }
0053     if (m_cfg.outputMeasurementParticlesMap.empty()) {
0054       throw std::invalid_argument(
0055           "Missing hit-to-particles map output collection");
0056     }
0057     if (m_cfg.outputMeasurementSimHitsMap.empty()) {
0058       throw std::invalid_argument(
0059           "Missing hit-to-simulated-hits map output collection");
0060     }
0061     if (m_cfg.outputParticleMeasurementsMap.empty()) {
0062       throw std::invalid_argument(
0063           "Missing particle-to-measurements map output collection");
0064     }
0065     if (m_cfg.outputSimHitMeasurementsMap.empty()) {
0066       throw std::invalid_argument(
0067           "Missing particle-to-simulated-hits map output collection");
0068     }
0069
0070     m_outputMeasurements.initialize(m_cfg.outputMeasurements);
0071     m_outputClusters.initialize(m_cfg.outputClusters);
0072     m_outputMeasurementParticlesMap.initialize(
0073         m_cfg.outputMeasurementParticlesMap);
0074     m_outputMeasurementSimHitsMap.initialize(m_cfg.outputMeasurementSimHitsMap);
0075     m_outputParticleMeasurementsMap.initialize(
0076         m_cfg.outputParticleMeasurementsMap);
0077     m_outputSimHitMeasurementsMap.initialize(m_cfg.outputSimHitMeasurementsMap);
0078   }
0079
0080   if (m_cfg.doOutputCells) {
0081     if (m_cfg.outputCells.empty()) {
0082       throw std::invalid_argument("Missing cell output collection");
0083     }
0084
0085     m_outputCells.initialize(m_cfg.outputCells);
0086   }
0087
0088   m_inputHits.initialize(m_cfg.inputSimHits);
0089
0090   // Create the digitizers from the configuration
0091   std::vector<std::pair<Acts::GeometryIdentifier, Digitizer>> digitizerInput;
0092
0093   for (std::size_t i = 0; i < m_cfg.digitizationConfigs.size(); ++i) {
0094     GeometricConfig geoCfg;
0095     Acts::GeometryIdentifier geoId = m_cfg.digitizationConfigs.idAt(i);
0096
0097     const auto& digiCfg = m_cfg.digitizationConfigs.valueAt(i);
0098     geoCfg = digiCfg.geometricDigiConfig;
0099     // Copy so we can sort in-place
0100     SmearingConfig smCfg = digiCfg.smearingDigiConfig;
0101
0102     std::vector<Acts::BoundIndices> indices;
0103     for (auto& gcf : smCfg.params) {
0104       indices.push_back(gcf.index);
0105     }
0106     indices.insert(indices.begin(), geoCfg.indices.begin(),
0107                    geoCfg.indices.end());
0108
0109     // Make sure the configured input parameter indices are sorted and unique
0110     std::ranges::sort(indices);
0111
0112     auto dup = std::adjacent_find(indices.begin(), indices.end());
0113     if (dup != indices.end()) {
0114       std::invalid_argument(
0115           "Digitization configuration contains duplicate parameter indices");
0116     }
0117
0118     switch (smCfg.params.size()) {
0119       case 0u:
0120         digitizerInput.emplace_back(geoId, makeDigitizer<0u>(digiCfg));
0121         break;
0122       case 1u:
0123         digitizerInput.emplace_back(geoId, makeDigitizer<1u>(digiCfg));
0124         break;
0125       case 2u:
0126         digitizerInput.emplace_back(geoId, makeDigitizer<2u>(digiCfg));
0127         break;
0128       case 3u:
0129         digitizerInput.emplace_back(geoId, makeDigitizer<3u>(digiCfg));
0130         break;
0131       case 4u:
0132         digitizerInput.emplace_back(geoId, makeDigitizer<4u>(digiCfg));
0133         break;
0134       default:
0135         throw std::invalid_argument("Unsupported smearer size");
0136     }
0137   }
0138
0139   m_digitizers = Acts::GeometryHierarchyMap<Digitizer>(digitizerInput);
0140 }
0141
0142 ProcessCode DigitizationAlgorithm::execute(const AlgorithmContext& ctx) const {
0143   // Retrieve input
0144   const auto& simHits = m_inputHits(ctx);
0145   ACTS_DEBUG("Loaded " << simHits.size() << " sim hits");
0146
0147   // Prepare output containers
0148   // need list here for stable addresses
0149   MeasurementContainer measurements;
0150   ClusterContainer clusters;
0151
0152   IndexMultimap<SimBarcode> measurementParticlesMap;
0153   IndexMultimap<Index> measurementSimHitsMap;
0154   measurements.reserve(simHits.size());
0155   measurementParticlesMap.reserve(simHits.size());
0156   measurementSimHitsMap.reserve(simHits.size());
0157
0158   // Setup random number generator
0159   auto rng = m_cfg.randomNumbers->spawnGenerator(ctx);
0160
0161   // Some statistics
0162   std::size_t skippedHits = 0;
0163
0164   // Some algorithms do the clusterization themselves such as the traccc chain.
0165   // Thus we need to store the cell data from the simulation.
0166   CellsMap cellsMap;
0167
0168   ACTS_DEBUG("Starting loop over modules ...");
0169   for (const auto& simHitsGroup : groupByModule(simHits)) {
0170     // Manual pair unpacking instead of using
0171     //   auto [moduleGeoId, moduleSimHits] : ...
0172     // otherwise clang on macos complains that it is unable to capture the local
0173     // binding in the lambda used for visiting the smearer below.
0174     Acts::GeometryIdentifier moduleGeoId = simHitsGroup.first;
0175     const auto& moduleSimHits = simHitsGroup.second;
0176
0177     auto surfaceItr = m_cfg.surfaceByIdentifier.find(moduleGeoId);
0178
0179     if (surfaceItr == m_cfg.surfaceByIdentifier.end()) {
0180       // this is either an invalid geometry id or a misconfigured smearer
0181       // setup; both cases can not be handled and should be fatal.
0182       ACTS_ERROR("Could not find surface " << moduleGeoId
0183                                            << " for configured smearer");
0184       return ProcessCode::ABORT;
0185     }
0186
0187     const Acts::Surface* surfacePtr = surfaceItr->second;
0188
0189     auto digitizerItr = m_digitizers.find(moduleGeoId);
0190     if (digitizerItr == m_digitizers.end()) {
0191       ACTS_VERBOSE("No digitizer present for module " << moduleGeoId);
0192       continue;
0193     } else {
0194       ACTS_VERBOSE("Digitizer found for module " << moduleGeoId);
0195     }
0196
0197     // Run the digitizer. Iterate over the hits for this surface inside the
0198     // visitor so we do not need to lookup the variant object per-hit.
0199     std::visit(
0200         [&](const auto& digitizer) {
0201           ModuleClusters moduleClusters(
0202               digitizer.geometric.segmentation, digitizer.geometric.indices,
0203               m_cfg.doMerge, m_cfg.mergeNsigma, m_cfg.mergeCommonCorner);
0204
0205           for (auto h = moduleSimHits.begin(); h != moduleSimHits.end(); ++h) {
0206             const auto& simHit = *h;
0207             const auto simHitIdx = simHits.index_of(h);
0208
0209             DigitizedParameters dParameters;
0210
0211             if (simHit.depositedEnergy() < m_cfg.minEnergyDeposit) {
0212               ACTS_VERBOSE("Skip hit because energy deposit to small");
0213               continue;
0214             }
0215
0216             // Geometric part - 0, 1, 2 local parameters are possible
0217             if (!digitizer.geometric.indices.empty()) {
0218               ACTS_VERBOSE("Configured to geometric digitize "
0219                            << digitizer.geometric.indices.size()
0220                            << " parameters.");
0221               const auto& cfg = digitizer.geometric;
0222               Acts::Vector3 driftDir = cfg.drift(simHit.position(), rng);
0223               auto channelsRes = m_channelizer.channelize(
0224                   simHit, *surfacePtr, ctx.geoContext, driftDir,
0225                   cfg.segmentation, cfg.thickness);
0226               if (!channelsRes.ok() || channelsRes->empty()) {
0227                 ACTS_DEBUG(
0228                     "Geometric channelization did not work, skipping this "
0229                     "hit.");
0230                 continue;
0231               }
0232               ACTS_VERBOSE("Activated " << channelsRes->size()
0233                                         << " channels for this hit.");
0234               dParameters =
0235                   localParameters(digitizer.geometric, *channelsRes, rng);
0236             }
0237
0238             // Smearing part - (optionally) rest
0239             if (!digitizer.smearing.indices.empty()) {
0240               ACTS_VERBOSE("Configured to smear "
0241                            << digitizer.smearing.indices.size()
0242                            << " parameters.");
0243               auto res =
0244                   digitizer.smearing(rng, simHit, *surfacePtr, ctx.geoContext);
0245               if (!res.ok()) {
0246                 ++skippedHits;
0247                 ACTS_DEBUG("Problem in hit smearing, skip hit ("
0248                            << res.error().message() << ")");
0249                 continue;
0250               }
0251               const auto& [par, cov] = res.value();
0252               for (Eigen::Index ip = 0; ip < par.rows(); ++ip) {
0253                 dParameters.indices.push_back(digitizer.smearing.indices[ip]);
0254                 dParameters.values.push_back(par[ip]);
0255                 dParameters.variances.push_back(cov(ip, ip));
0256               }
0257             }
0258
0259             // Check on success - threshold could have eliminated all channels
0260             if (dParameters.values.empty()) {
0261               ACTS_VERBOSE(
0262                   "Parameter digitization did not yield a measurement.");
0263               continue;
0264             }
0265
0266             moduleClusters.add(std::move(dParameters), simHitIdx);
0267           }
0268
0269           auto digitizeParametersResult = moduleClusters.digitizedParameters();
0270
0271           // Store the cell data into a map.
0272           if (m_cfg.doOutputCells) {
0273             std::vector<Cluster::Cell> cells;
0274             for (const auto& [dParameters, simHitsIdxs] :
0275                  digitizeParametersResult) {
0276               for (const auto& cell : dParameters.cluster.channels) {
0277                 cells.push_back(cell);
0278               }
0279             }
0280             cellsMap.insert({moduleGeoId, std::move(cells)});
0281           }
0282
0283           if (m_cfg.doClusterization) {
0284             for (auto& [dParameters, simHitsIdxs] : digitizeParametersResult) {
0285               auto measurement =
0286                   createMeasurement(measurements, moduleGeoId, dParameters);
0287
0288               dParameters.cluster.globalPosition = measurementGlobalPosition(
0289                   dParameters, *surfacePtr, ctx.geoContext);
0290               clusters.emplace_back(std::move(dParameters.cluster));
0291
0292               for (auto [i, simHitIdx] : Acts::enumerate(simHitsIdxs)) {
0293                 measurementParticlesMap.emplace_hint(
0294                     measurementParticlesMap.end(), measurement.index(),
0295                     simHits.nth(simHitIdx)->particleId());
0296                 measurementSimHitsMap.emplace_hint(measurementSimHitsMap.end(),
0297                                                    measurement.index(),
0298                                                    simHitIdx);
0299               }
0300             }
0301           }
0302         },
0303         *digitizerItr);
0304   }
0305
0306   if (skippedHits > 0) {
0307     ACTS_WARNING(
0308         skippedHits
0309         << " skipped in Digitization. Enable DEBUG mode to see more details.");
0310   }
0311
0312   if (m_cfg.doClusterization) {
0313     ACTS_DEBUG("Created " << measurements.size() << " measurements, "
0314                           << clusters.size() << " clusters" << " from "
0315                           << simHits.size() << " sim hits.");
0316
0317     m_outputMeasurements(ctx, std::move(measurements));
0318     m_outputClusters(ctx, std::move(clusters));
0319
0320     // invert them before they are moved
0321     m_outputParticleMeasurementsMap(
0322         ctx, invertIndexMultimap(measurementParticlesMap));
0323     m_outputSimHitMeasurementsMap(ctx,
0324                                   invertIndexMultimap(measurementSimHitsMap));
0325
0326     m_outputMeasurementParticlesMap(ctx, std::move(measurementParticlesMap));
0327     m_outputMeasurementSimHitsMap(ctx, std::move(measurementSimHitsMap));
0328   }
0329
0330   if (m_cfg.doOutputCells) {
0331     m_outputCells(ctx, std::move(cellsMap));
0332   }
0333
0334   return ProcessCode::SUCCESS;
0335 }
0336
0337 DigitizedParameters DigitizationAlgorithm::localParameters(
0338     const GeometricConfig& geoCfg,
0339     const std::vector<ActsFatras::Segmentizer::ChannelSegment>& channels,
0340     RandomEngine& rng) const {
0341   DigitizedParameters dParameters;
0342
0343   const auto& binningData = geoCfg.segmentation.binningData();
0344
0345   // For digital readout, the weight needs to be split in x and y
0346   std::array<double, 2u> pos = {0., 0.};
0347   std::array<double, 2u> totalWeight = {0., 0.};
0348   std::array<std::size_t, 2u> bmin = {std::numeric_limits<std::size_t>::max(),
0349                                       std::numeric_limits<std::size_t>::max()};
0350   std::array<std::size_t, 2u> bmax = {0, 0};
0351
0352   // The component digital store
0353   std::array<std::set<std::size_t>, 2u> componentChannels;
0354
0355   // Combine the channels
0356   for (const auto& ch : channels) {
0357     auto bin = ch.bin;
0358     double charge = geoCfg.charge(ch.activation, rng);
0359     // Loop and check
0360     if (charge > geoCfg.threshold) {
0361       double weight = geoCfg.digital ? 1. : charge;
0362       for (std::size_t ib = 0; ib < 2; ++ib) {
0363         if (geoCfg.digital && geoCfg.componentDigital) {
0364           // only fill component of this row/column if not yet filled
0365           if (!componentChannels[ib].contains(bin[ib])) {
0366             totalWeight[ib] += weight;
0367             pos[ib] += weight * binningData[ib].center(bin[ib]);
0368             componentChannels[ib].insert(bin[ib]);
0369           }
0370         } else {
0371           totalWeight[ib] += weight;
0372           pos[ib] += weight * binningData[ib].center(bin[ib]);
0373         }
0374         // min max channels
0375         bmin[ib] = std::min(bmin[ib], static_cast<std::size_t>(bin[ib]));
0376         bmax[ib] = std::max(bmax[ib], static_cast<std::size_t>(bin[ib]));
0377       }
0378       // Create a copy of the channel, as activation may change
0379       auto chdig = ch;
0380       chdig.bin = ch.bin;
0381       chdig.activation = charge;
0382       dParameters.cluster.channels.push_back(chdig);
0383     }
0384   }
0385   if (totalWeight[0] > 0. && totalWeight[1] > 0.) {
0386     pos[0] /= totalWeight[0];
0387     pos[1] /= totalWeight[1];
0388     dParameters.indices = geoCfg.indices;
0389     for (auto idx : dParameters.indices) {
0390       dParameters.values.push_back(pos[idx]);
0391     }
0392     std::size_t size0 = (bmax[0] - bmin[0] + 1);
0393     std::size_t size1 = (bmax[1] - bmin[1] + 1);
0394
0395     dParameters.variances = geoCfg.variances({size0, size1}, bmin);
0396     dParameters.cluster.sizeLoc0 = size0;
0397     dParameters.cluster.sizeLoc1 = size1;
0398   }
0399
0400   return dParameters;
0401 }
0402
0403 }  // namespace ActsExamples