Digitization/src/PlanarSteppingAlgorithm.cpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2017-2020 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/PlanarSteppingAlgorithm.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/TrackParametrization.hpp"
0013 #include "Acts/Definitions/Units.hpp"
0014 #include "Acts/Digitization/DigitizationCell.hpp"
0015 #include "Acts/Digitization/DigitizationModule.hpp"
0016 #include "Acts/Digitization/PlanarModuleCluster.hpp"
0017 #include "Acts/Digitization/PlanarModuleStepper.hpp"
0018 #include "Acts/Digitization/Segmentation.hpp"
0019 #include "Acts/EventData/Measurement.hpp"
0020 #include "Acts/EventData/SourceLink.hpp"
0021 #include "Acts/Geometry/GeometryIdentifier.hpp"
0022 #include "Acts/Geometry/TrackingGeometry.hpp"
0023 #include "Acts/Plugins/Identification/IdentifiedDetectorElement.hpp"
0024 #include "Acts/Surfaces/Surface.hpp"
0025 #include "Acts/Surfaces/SurfaceArray.hpp"
0026 #include "ActsExamples/EventData/GeometryContainers.hpp"
0027 #include "ActsExamples/EventData/IndexSourceLink.hpp"
0028 #include "ActsExamples/EventData/Measurement.hpp"
0029 #include "ActsExamples/EventData/SimHit.hpp"
0030 #include "ActsExamples/Framework/AlgorithmContext.hpp"
0031 #include "ActsExamples/Utilities/GroupBy.hpp"
0032 #include "ActsExamples/Utilities/Range.hpp"
0033 #include "ActsFatras/EventData/Barcode.hpp"
0034 #include "ActsFatras/EventData/Hit.hpp"
0035
0036 #include <cmath>
0037 #include <cstddef>
0038 #include <ostream>
0039 #include <stdexcept>
0040 #include <utility>
0041 #include <vector>
0042
0043 ActsExamples::PlanarSteppingAlgorithm::PlanarSteppingAlgorithm(
0044     ActsExamples::PlanarSteppingAlgorithm::Config config,
0045     Acts::Logging::Level level)
0046     : ActsExamples::IAlgorithm("PlanarSteppingAlgorithm", level),
0047       m_cfg(std::move(config)) {
0048   if (m_cfg.inputSimHits.empty()) {
0049     throw std::invalid_argument("Missing input hits collection");
0050   }
0051   if (m_cfg.outputClusters.empty()) {
0052     throw std::invalid_argument("Missing output clusters collection");
0053   }
0054   if (m_cfg.outputSourceLinks.empty()) {
0055     throw std::invalid_argument("Missing source links output collection");
0056   }
0057   if (m_cfg.outputDigiSourceLinks.empty()) {
0058     throw std::invalid_argument(
0059         "Missing digitization source links output collection");
0060   }
0061   if (m_cfg.outputMeasurements.empty()) {
0062     throw std::invalid_argument("Missing measurements output collection");
0063   }
0064   if (m_cfg.outputMeasurementParticlesMap.empty()) {
0065     throw std::invalid_argument(
0066         "Missing hit-to-particles map output collection");
0067   }
0068   if (m_cfg.outputMeasurementSimHitsMap.empty()) {
0069     throw std::invalid_argument(
0070         "Missing hit-to-simulated-hits map output collection");
0071   }
0072   if (!m_cfg.trackingGeometry) {
0073     throw std::invalid_argument("Missing tracking geometry");
0074   }
0075   if (!m_cfg.planarModuleStepper) {
0076     throw std::invalid_argument("Missing planar module stepper");
0077   }
0078   if (!m_cfg.randomNumbers) {
0079     throw std::invalid_argument("Missing random numbers tool");
0080   }
0081
0082   m_inputSimHits.initialize(m_cfg.inputSimHits);
0083
0084   m_outputClusters.initialize(m_cfg.outputClusters);
0085   m_outputSourceLinks.initialize(m_cfg.outputSourceLinks);
0086   m_outputDigiSourceLinks.initialize(m_cfg.outputDigiSourceLinks);
0087   m_outputMeasurements.initialize(m_cfg.outputMeasurements);
0088   m_outputMeasurementParticlesMap.initialize(
0089       m_cfg.outputMeasurementParticlesMap);
0090   m_outputMeasurementSimHitsMap.initialize(m_cfg.outputMeasurementSimHitsMap);
0091
0092   // fill the digitizables map to allow lookup by geometry id only
0093   m_cfg.trackingGeometry->visitSurfaces([this](const Acts::Surface* surface) {
0094     Digitizable dg;
0095     // require a valid surface
0096     dg.surface = surface;
0097     if (dg.surface == nullptr) {
0098       return;
0099     }
0100     // require an associated detector element
0101     dg.detectorElement = dynamic_cast<const Acts::IdentifiedDetectorElement*>(
0102         dg.surface->associatedDetectorElement());
0103     if (dg.detectorElement == nullptr) {
0104       return;
0105     }
0106     // require an associated digitization module
0107     dg.digitizer = dg.detectorElement->digitizationModule().get();
0108     if (dg.digitizer == nullptr) {
0109       return;
0110     }
0111     // record all valid surfaces
0112     this->m_digitizables.insert_or_assign(surface->geometryId(), dg);
0113   });
0114 }
0115
0116 ActsExamples::ProcessCode ActsExamples::PlanarSteppingAlgorithm::execute(
0117     const AlgorithmContext& ctx) const {
0118   // retrieve input
0119   const auto& simHits = m_inputSimHits(ctx);
0120
0121   // prepare output containers
0122   ClusterContainer clusters;
0123
0124   std::vector<Acts::DigitizationSourceLink> digiSourceLinks;
0125
0126   GeometryIdMultiset<IndexSourceLink> sourceLinks;
0127   MeasurementContainer measurements;
0128   IndexMultimap<ActsFatras::Barcode> hitParticlesMap;
0129   IndexMultimap<Index> hitSimHitsMap;
0130   clusters.reserve(simHits.size());
0131   measurements.reserve(simHits.size());
0132   hitParticlesMap.reserve(simHits.size());
0133   hitSimHitsMap.reserve(simHits.size());
0134
0135   for (auto&& [moduleGeoId, moduleSimHits] : groupByModule(simHits)) {
0136     // can only digitize hits on digitizable surfaces
0137     const auto it = m_digitizables.find(moduleGeoId);
0138     if (it == m_digitizables.end()) {
0139       continue;
0140     }
0141
0142     const auto& dg = it->second;
0143     // local intersection / direction
0144     const auto invTransfrom = dg.surface->transform(ctx.geoContext).inverse();
0145
0146     // use iterators manually so we can retrieve the hit index in the container
0147     for (auto ih = moduleSimHits.begin(); ih != moduleSimHits.end(); ++ih) {
0148       const auto& simHit = *ih;
0149       const auto simHitIdx = simHits.index_of(ih);
0150
0151       Acts::Vector2 localIntersect =
0152           (invTransfrom * simHit.position()).head<2>();
0153       Acts::Vector3 localDirection = invTransfrom.linear() * simHit.direction();
0154
0155       // compute digitization steps
0156       const auto thickness = dg.detectorElement->thickness();
0157       const auto lorentzAngle = dg.digitizer->lorentzAngle();
0158       auto lorentzShift = thickness * std::tan(lorentzAngle);
0159       lorentzShift *= -(dg.digitizer->readoutDirection());
0160       // now calculate the steps through the silicon
0161       std::vector<Acts::DigitizationStep> dSteps =
0162           m_cfg.planarModuleStepper->cellSteps(ctx.geoContext, *dg.digitizer,
0163                                                localIntersect, localDirection);
0164       // everything under threshold or edge effects
0165       if (dSteps.empty()) {
0166         ACTS_VERBOSE("No steps returned from stepper.");
0167         continue;
0168       }
0169
0170       // Create a cluster - centroid method
0171       double localX = 0.;
0172       double localY = 0.;
0173       double totalPath = 0.;
0174       // the cells to be used
0175       std::vector<Acts::DigitizationCell> usedCells;
0176       usedCells.reserve(dSteps.size());
0177       // loop over the steps
0178       for (auto dStep : dSteps) {
0179         // @todo implement smearing
0180         localX += dStep.stepLength * dStep.stepCellCenter.x();
0181         localY += dStep.stepLength * dStep.stepCellCenter.y();
0182         totalPath += dStep.stepLength;
0183         usedCells.push_back(Acts::DigitizationCell(dStep.stepCell.channel0,
0184                                                    dStep.stepCell.channel1,
0185                                                    dStep.stepLength));
0186       }
0187       // divide by the total path
0188       localX /= totalPath;
0189       localX += lorentzShift;
0190       localY /= totalPath;
0191
0192       // get the segmentation & find the corresponding cell id
0193       const Acts::Segmentation& segmentation = dg.digitizer->segmentation();
0194       auto binUtility = segmentation.binUtility();
0195       Acts::Vector2 localPosition(localX, localY);
0196       // @todo remove unnecessary conversion
0197       // std::size_t bin0 = binUtility.bin(localPosition, 0);
0198       // std::size_t bin1 = binUtility.bin(localPosition, 1);
0199       // std::size_t binSerialized = binUtility.serialize({{bin0, bin1, 0}});
0200
0201       // the covariance is currently set to some arbitrary value.
0202       Acts::SquareMatrix3 cov;
0203       cov << 0.05, 0., 0., 0., 0.05, 0., 0., 0.,
0204           900. * Acts::UnitConstants::ps * Acts::UnitConstants::ps;
0205       Acts::Vector3 par(localX, localY, simHit.time());
0206
0207       // create the planar cluster
0208       digiSourceLinks.emplace_back(moduleGeoId,
0209                                    std::vector<std::size_t>{simHitIdx});
0210       Acts::DigitizationSourceLink& digiSourceLink = digiSourceLinks.back();
0211
0212       Acts::PlanarModuleCluster cluster(
0213           dg.surface->getSharedPtr(), Acts::SourceLink{digiSourceLink}, cov,
0214           localX, localY, simHit.time(), std::move(usedCells));
0215
0216       // the measurement container is unordered and the index under which
0217       // the measurement will be stored is known before adding it.
0218       Index hitIdx = measurements.size();
0219       IndexSourceLink sourceLink{moduleGeoId, hitIdx};
0220
0221       sourceLinks.insert(sourceLinks.end(), sourceLink);
0222
0223       auto meas = Acts::makeMeasurement(Acts::SourceLink{sourceLink}, par, cov,
0224                                         Acts::eBoundLoc0, Acts::eBoundLoc1,
0225                                         Acts::eBoundTime);
0226
0227       // add to output containers. since the input is already geometry-order,
0228       // new elements in geometry containers can just be appended at the end.
0229       clusters.emplace_hint(clusters.end(), moduleGeoId, std::move(cluster));
0230       measurements.emplace_back(std::move(meas));
0231       // no hit merging -> only one mapping per digitized hit.
0232       hitParticlesMap.emplace_hint(hitParticlesMap.end(), hitIdx,
0233                                    simHit.particleId());
0234       hitSimHitsMap.emplace_hint(hitSimHitsMap.end(), hitIdx, simHitIdx);
0235     }
0236   }
0237
0238   ACTS_DEBUG("digitized " << simHits.size() << " hits into " << clusters.size()
0239                           << " clusters");
0240
0241   m_outputClusters(ctx, std::move(clusters));
0242   m_outputSourceLinks(ctx, std::move(sourceLinks));
0243   m_outputDigiSourceLinks(ctx, std::move(digiSourceLinks));
0244   m_outputMeasurements(ctx, std::move(measurements));
0245   m_outputMeasurementParticlesMap(ctx, std::move(hitParticlesMap));
0246   m_outputMeasurementSimHitsMap(ctx, std::move(hitSimHitsMap));
0247   return ActsExamples::ProcessCode::SUCCESS;
0248 }