Digitization/src/DigitizationAlgorithm.cpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2021 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/DigitizationAlgorithm.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/TrackParametrization.hpp"
0013 #include "Acts/Geometry/GeometryIdentifier.hpp"
0014 #include "Acts/Geometry/TrackingGeometry.hpp"
0015 #include "Acts/Surfaces/Surface.hpp"
0016 #include "Acts/Utilities/BinUtility.hpp"
0017 #include "Acts/Utilities/Result.hpp"
0018 #include "ActsExamples/Digitization/ModuleClusters.hpp"
0019 #include "ActsExamples/EventData/GeometryContainers.hpp"
0020 #include "ActsExamples/EventData/Index.hpp"
0021 #include "ActsExamples/EventData/IndexSourceLink.hpp"
0022 #include "ActsExamples/EventData/SimHit.hpp"
0023 #include "ActsExamples/Framework/AlgorithmContext.hpp"
0024 #include "ActsExamples/Utilities/GroupBy.hpp"
0025 #include "ActsExamples/Utilities/Range.hpp"
0026 #include "ActsFatras/EventData/Barcode.hpp"
0027 #include "ActsFatras/EventData/Hit.hpp"
0028
0029 #include <algorithm>
0030 #include <array>
0031 #include <cmath>
0032 #include <cstdint>
0033 #include <ostream>
0034 #include <set>
0035 #include <stdexcept>
0036 #include <string>
0037 #include <utility>
0038
0039 ActsExamples::DigitizationAlgorithm::DigitizationAlgorithm(
0040     DigitizationConfig config, Acts::Logging::Level level)
0041     : ActsExamples::IAlgorithm("DigitizationAlgorithm", level),
0042       m_cfg(std::move(config)) {
0043   if (m_cfg.inputSimHits.empty()) {
0044     throw std::invalid_argument("Missing simulated hits input collection");
0045   }
0046   if (m_cfg.outputMeasurements.empty()) {
0047     throw std::invalid_argument("Missing measurements output collection");
0048   }
0049   if (m_cfg.outputSourceLinks.empty()) {
0050     throw std::invalid_argument("Missing source links output collection");
0051   }
0052   if (m_cfg.outputMeasurementParticlesMap.empty()) {
0053     throw std::invalid_argument(
0054         "Missing hit-to-particles map output collection");
0055   }
0056   if (m_cfg.outputMeasurementSimHitsMap.empty()) {
0057     throw std::invalid_argument(
0058         "Missing hit-to-simulated-hits map output collection");
0059   }
0060   if (!m_cfg.trackingGeometry) {
0061     throw std::invalid_argument("Missing tracking geometry");
0062   }
0063   if (!m_cfg.randomNumbers) {
0064     throw std::invalid_argument("Missing random numbers tool");
0065   }
0066
0067   if (m_cfg.digitizationConfigs.empty()) {
0068     throw std::invalid_argument("Missing digitization configuration");
0069   }
0070
0071   m_simContainerReadHandle.initialize(m_cfg.inputSimHits);
0072   m_sourceLinkWriteHandle.initialize(m_cfg.outputSourceLinks);
0073   m_measurementWriteHandle.initialize(m_cfg.outputMeasurements);
0074   m_clusterWriteHandle.initialize(m_cfg.outputClusters);
0075   m_measurementParticlesMapWriteHandle.initialize(
0076       m_cfg.outputMeasurementParticlesMap);
0077   m_measurementSimHitsMapWriteHandle.initialize(
0078       m_cfg.outputMeasurementSimHitsMap);
0079
0080   // Create the digitizers from the configuration
0081   std::vector<std::pair<Acts::GeometryIdentifier, Digitizer>> digitizerInput;
0082
0083   for (std::size_t i = 0; i < m_cfg.digitizationConfigs.size(); ++i) {
0084     GeometricConfig geoCfg;
0085     Acts::GeometryIdentifier geoId = m_cfg.digitizationConfigs.idAt(i);
0086
0087     const auto& digiCfg = m_cfg.digitizationConfigs.valueAt(i);
0088     geoCfg = digiCfg.geometricDigiConfig;
0089     // Copy so we can sort in-place
0090     SmearingConfig smCfg = digiCfg.smearingDigiConfig;
0091
0092     std::vector<Acts::BoundIndices> indices;
0093     for (auto& gcf : smCfg) {
0094       indices.push_back(gcf.index);
0095     }
0096     indices.insert(indices.begin(), geoCfg.indices.begin(),
0097                    geoCfg.indices.end());
0098
0099     // Make sure the configured input parameter indices are sorted and unique
0100     std::sort(indices.begin(), indices.end());
0101
0102     auto dup = std::adjacent_find(indices.begin(), indices.end());
0103     if (dup != indices.end()) {
0104       std::invalid_argument(
0105           "Digitization configuration contains duplicate parameter indices");
0106     }
0107
0108     switch (smCfg.size()) {
0109       case 0u:
0110         digitizerInput.emplace_back(geoId, makeDigitizer<0u>(digiCfg));
0111         break;
0112       case 1u:
0113         digitizerInput.emplace_back(geoId, makeDigitizer<1u>(digiCfg));
0114         break;
0115       case 2u:
0116         digitizerInput.emplace_back(geoId, makeDigitizer<2u>(digiCfg));
0117         break;
0118       case 3u:
0119         digitizerInput.emplace_back(geoId, makeDigitizer<3u>(digiCfg));
0120         break;
0121       case 4u:
0122         digitizerInput.emplace_back(geoId, makeDigitizer<4u>(digiCfg));
0123         break;
0124       default:
0125         throw std::invalid_argument("Unsupported smearer size");
0126     }
0127   }
0128
0129   m_digitizers = Acts::GeometryHierarchyMap<Digitizer>(digitizerInput);
0130 }
0131
0132 ActsExamples::ProcessCode ActsExamples::DigitizationAlgorithm::execute(
0133     const AlgorithmContext& ctx) const {
0134   // Retrieve input
0135   const auto& simHits = m_simContainerReadHandle(ctx);
0136   ACTS_DEBUG("Loaded " << simHits.size() << " sim hits");
0137
0138   // Prepare output containers
0139   // need list here for stable addresses
0140   IndexSourceLinkContainer sourceLinks;
0141   MeasurementContainer measurements;
0142   ClusterContainer clusters;
0143   IndexMultimap<ActsFatras::Barcode> measurementParticlesMap;
0144   IndexMultimap<Index> measurementSimHitsMap;
0145   sourceLinks.reserve(simHits.size());
0146   measurements.reserve(simHits.size());
0147   measurementParticlesMap.reserve(simHits.size());
0148   measurementSimHitsMap.reserve(simHits.size());
0149
0150   // Setup random number generator
0151   auto rng = m_cfg.randomNumbers->spawnGenerator(ctx);
0152
0153   // Some statistics
0154   std::size_t skippedHits = 0;
0155
0156   ACTS_DEBUG("Starting loop over modules ...");
0157   for (const auto& simHitsGroup : groupByModule(simHits)) {
0158     // Manual pair unpacking instead of using
0159     //   auto [moduleGeoId, moduleSimHits] : ...
0160     // otherwise clang on macos complains that it is unable to capture the local
0161     // binding in the lambda used for visiting the smearer below.
0162     Acts::GeometryIdentifier moduleGeoId = simHitsGroup.first;
0163     const auto& moduleSimHits = simHitsGroup.second;
0164
0165     const Acts::Surface* surfacePtr =
0166         m_cfg.trackingGeometry->findSurface(moduleGeoId);
0167
0168     if (surfacePtr == nullptr) {
0169       // this is either an invalid geometry id or a misconfigured smearer
0170       // setup; both cases can not be handled and should be fatal.
0171       ACTS_ERROR("Could not find surface " << moduleGeoId
0172                                            << " for configured smearer");
0173       return ProcessCode::ABORT;
0174     }
0175
0176     auto digitizerItr = m_digitizers.find(moduleGeoId);
0177     if (digitizerItr == m_digitizers.end()) {
0178       ACTS_VERBOSE("No digitizer present for module " << moduleGeoId);
0179       continue;
0180     } else {
0181       ACTS_VERBOSE("Digitizer found for module " << moduleGeoId);
0182     }
0183
0184     // Run the digitizer. Iterate over the hits for this surface inside the
0185     // visitor so we do not need to lookup the variant object per-hit.
0186     std::visit(
0187         [&](const auto& digitizer) {
0188           ModuleClusters moduleClusters(
0189               digitizer.geometric.segmentation, digitizer.geometric.indices,
0190               m_cfg.doMerge, m_cfg.mergeNsigma, m_cfg.mergeCommonCorner);
0191
0192           for (auto h = moduleSimHits.begin(); h != moduleSimHits.end(); ++h) {
0193             const auto& simHit = *h;
0194             const auto simHitIdx = simHits.index_of(h);
0195
0196             DigitizedParameters dParameters;
0197
0198             if (simHit.depositedEnergy() < m_cfg.minEnergyDeposit) {
0199               ACTS_VERBOSE("Skip hit because energy deposit to small")
0200               continue;
0201             }
0202
0203             // Geometric part - 0, 1, 2 local parameters are possible
0204             if (!digitizer.geometric.indices.empty()) {
0205               ACTS_VERBOSE("Configured to geometric digitize "
0206                            << digitizer.geometric.indices.size()
0207                            << " parameters.");
0208               const auto& cfg = digitizer.geometric;
0209               Acts::Vector3 driftDir = cfg.drift(simHit.position(), rng);
0210               auto channelsRes = m_channelizer.channelize(
0211                   simHit, *surfacePtr, ctx.geoContext, driftDir,
0212                   cfg.segmentation, cfg.thickness);
0213               if (!channelsRes.ok() || channelsRes->empty()) {
0214                 ACTS_DEBUG(
0215                     "Geometric channelization did not work, skipping this hit.")
0216                 continue;
0217               }
0218               ACTS_VERBOSE("Activated " << channelsRes->size()
0219                                         << " channels for this hit.");
0220               dParameters =
0221                   localParameters(digitizer.geometric, *channelsRes, rng);
0222             }
0223
0224             // Smearing part - (optionally) rest
0225             if (!digitizer.smearing.indices.empty()) {
0226               ACTS_VERBOSE("Configured to smear "
0227                            << digitizer.smearing.indices.size()
0228                            << " parameters.");
0229               auto res =
0230                   digitizer.smearing(rng, simHit, *surfacePtr, ctx.geoContext);
0231               if (!res.ok()) {
0232                 ++skippedHits;
0233                 ACTS_DEBUG("Problem in hit smearing, skip hit ("
0234                            << res.error().message() << ")");
0235                 continue;
0236               }
0237               const auto& [par, cov] = res.value();
0238               for (Eigen::Index ip = 0; ip < par.rows(); ++ip) {
0239                 dParameters.indices.push_back(digitizer.smearing.indices[ip]);
0240                 dParameters.values.push_back(par[ip]);
0241                 dParameters.variances.push_back(cov(ip, ip));
0242               }
0243             }
0244
0245             // Check on success - threshold could have eliminated all channels
0246             if (dParameters.values.empty()) {
0247               ACTS_VERBOSE(
0248                   "Parameter digitization did not yield a measurement.")
0249               continue;
0250             }
0251
0252             moduleClusters.add(std::move(dParameters), simHitIdx);
0253           }
0254
0255           for (auto& [dParameters, simhits] :
0256                moduleClusters.digitizedParameters()) {
0257             // The measurement container is unordered and the index under which
0258             // the measurement will be stored is known before adding it.
0259             Index measurementIdx = measurements.size();
0260             IndexSourceLink sourceLink{moduleGeoId, measurementIdx};
0261
0262             // Add to output containers:
0263             // index map and source link container are geometry-ordered.
0264             // since the input is also geometry-ordered, new items can
0265             // be added at the end.
0266             sourceLinks.insert(sourceLinks.end(), sourceLink);
0267
0268             measurements.emplace_back(
0269                 createMeasurement(dParameters, sourceLink));
0270             clusters.emplace_back(std::move(dParameters.cluster));
0271             // this digitization does hit merging so there can be more than one
0272             // mapping entry for each digitized hit.
0273             for (auto simHitIdx : simhits) {
0274               measurementParticlesMap.emplace_hint(
0275                   measurementParticlesMap.end(), measurementIdx,
0276                   simHits.nth(simHitIdx)->particleId());
0277               measurementSimHitsMap.emplace_hint(measurementSimHitsMap.end(),
0278                                                  measurementIdx, simHitIdx);
0279             }
0280           }
0281         },
0282         *digitizerItr);
0283   }
0284
0285   if (skippedHits > 0) {
0286     ACTS_WARNING(
0287         skippedHits
0288         << " skipped in Digitization. Enable DEBUG mode to see more details.");
0289   }
0290
0291   m_sourceLinkWriteHandle(ctx, std::move(sourceLinks));
0292   m_measurementWriteHandle(ctx, std::move(measurements));
0293   m_clusterWriteHandle(ctx, std::move(clusters));
0294   m_measurementParticlesMapWriteHandle(ctx, std::move(measurementParticlesMap));
0295   m_measurementSimHitsMapWriteHandle(ctx, std::move(measurementSimHitsMap));
0296   return ProcessCode::SUCCESS;
0297 }
0298
0299 ActsExamples::DigitizedParameters
0300 ActsExamples::DigitizationAlgorithm::localParameters(
0301     const GeometricConfig& geoCfg,
0302     const std::vector<ActsFatras::Segmentizer::ChannelSegment>& channels,
0303     RandomEngine& rng) const {
0304   DigitizedParameters dParameters;
0305
0306   const auto& binningData = geoCfg.segmentation.binningData();
0307
0308   Acts::ActsScalar totalWeight = 0.;
0309   Acts::Vector2 m(0., 0.);
0310   std::size_t b0min = SIZE_MAX;
0311   std::size_t b0max = 0;
0312   std::size_t b1min = SIZE_MAX;
0313   std::size_t b1max = 0;
0314   // Combine the channels
0315   for (const auto& ch : channels) {
0316     auto bin = ch.bin;
0317     Acts::ActsScalar charge =
0318         geoCfg.digital ? 1. : geoCfg.charge(ch.activation, rng);
0319     if (geoCfg.digital || charge > geoCfg.threshold) {
0320       totalWeight += charge;
0321       std::size_t b0 = bin[0];
0322       std::size_t b1 = bin[1];
0323       m += Acts::Vector2(charge * binningData[0].center(b0),
0324                          charge * binningData[1].center(b1));
0325       b0min = std::min(b0min, b0);
0326       b0max = std::max(b0max, b0);
0327       b1min = std::min(b1min, b1);
0328       b1max = std::max(b1max, b1);
0329       // Create a copy of the channel, as activation may change
0330       auto chdig = ch;
0331       chdig.bin = ch.bin;
0332       chdig.activation = charge;
0333       dParameters.cluster.channels.push_back(chdig);
0334     }
0335   }
0336   if (totalWeight > 0.) {
0337     m *= 1. / totalWeight;
0338     dParameters.indices = geoCfg.indices;
0339     for (auto idx : dParameters.indices) {
0340       dParameters.values.push_back(m[idx]);
0341     }
0342     std::size_t size0 = static_cast<std::size_t>(b0max - b0min + 1);
0343     std::size_t size1 = static_cast<std::size_t>(b1max - b1min + 1);
0344     auto variances = geoCfg.variances(size0, size1, rng);
0345     if (variances.size() == dParameters.indices.size()) {
0346       dParameters.variances = variances;
0347     } else {
0348       dParameters.variances =
0349           std::vector<Acts::ActsScalar>(dParameters.indices.size(), -1.);
0350     }
0351
0352     if (dParameters.variances[0] == -1) {
0353       std::size_t ictr = b0min + size0 / 2;
0354       dParameters.variances[0] = std::pow(binningData[0].width(ictr), 2) / 12.0;
0355     }
0356     if (dParameters.variances[1] == -1) {
0357       std::size_t ictr = b1min + size1 / 2;
0358       dParameters.variances[1] = std::pow(binningData[1].width(ictr), 2) / 12.0;
0359     }
0360
0361     dParameters.cluster.sizeLoc0 = size0;
0362     dParameters.cluster.sizeLoc1 = size1;
0363   }
0364
0365   return dParameters;
0366 }