Propagation/src/SimHitToSummaryConversion.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/Propagation/SimHitToSummaryConversion.hpp"
0010
0011 #include "ActsExamples/EventData/GeometryContainers.hpp"
0012
0013 namespace {
0014 /// Helper method to cancatenate the steps and push them onto the surface
0015 ///
0016 /// @param gctx is the geometry context
0017 /// @param steps is the vector of steps to concatenate
0018 /// @param surface is the surface to push the steps onto
0019 ///
0020 /// @return the concatenated step
0021 Acts::detail::Step concatenateSteps(
0022     const Acts::GeometryContext& gctx,
0023     const std::vector<Acts::detail::Step>& steps,
0024     const Acts::Surface& surface) {
0025   // Average the position and direction
0026   Acts::detail::Step concatStep;
0027   if (steps.size() > 1) {
0028     for (const Acts::detail::Step& step : steps) {
0029       concatStep.position += step.position;
0030       concatStep.momentum += step.momentum;
0031     }
0032     double weight = 1. / steps.size();
0033     concatStep.position *= weight;
0034     concatStep.momentum *= weight;
0035   } else {
0036     concatStep = steps.front();
0037   }
0038   // Re-evaulate the position with a surface intersection
0039   auto intersection =
0040       surface.intersect(gctx, concatStep.position, concatStep.momentum);
0041   for (const auto& rsIntersection : intersection.split()) {
0042     if (rsIntersection.isValid()) {
0043       concatStep.position = rsIntersection.position();
0044       break;
0045     }
0046   }
0047   // Set the surface identifier
0048   concatStep.geoID = surface.geometryId();
0049   return concatStep;
0050 }
0051 }  // namespace
0052
0053 ActsExamples::SimHitToSummaryConversion::SimHitToSummaryConversion(
0054     const Config& config, Acts::Logging::Level level)
0055     : IAlgorithm("SimHitToSummaryConversion", level), m_cfg(config) {
0056   if (m_cfg.inputSimHits.empty()) {
0057     throw std::invalid_argument("Missing simulated hits input collection");
0058   }
0059   if (m_cfg.inputParticles.empty()) {
0060     throw std::invalid_argument("Missing simulated particles input collection");
0061   }
0062   m_inputSimHits.initialize(m_cfg.inputSimHits);
0063   m_inputParticles.initialize(m_cfg.inputParticles);
0064   m_outputSummary.initialize(m_cfg.outputSummaryCollection);
0065 }
0066
0067 ActsExamples::ProcessCode ActsExamples::SimHitToSummaryConversion::execute(
0068     const AlgorithmContext& context) const {
0069   // Retrieve the simulated hits
0070   const auto& simHits = m_inputSimHits(context);
0071   // Retrieve the particles at its basis
0072   const auto& particles = m_inputParticles(context);
0073
0074   PropagationSummaries propagationSummaries;
0075   propagationSummaries.reserve(particles.size());
0076
0077   // Prepare and sort
0078   std::unordered_map<unsigned long,
0079                      std::vector<std::vector<Acts::detail::Step>>>
0080       trackSteps;
0081   for (const auto& simHitsGroup : groupByModule(simHits)) {
0082     // Manual pair unpacking instead of using
0083     //   auto [moduleGeoId, moduleSimHits] : ...
0084     // otherwise clang on macos complains that it is unable to capture the local
0085     // binding in the lambda used for visiting the smearer below.
0086     Acts::GeometryIdentifier moduleGeoId = simHitsGroup.first;
0087     const auto& moduleSimHits = simHitsGroup.second;
0088     std::unordered_map<unsigned long, std::vector<Acts::detail::Step>>
0089         moduleSteps;
0090     for (const auto& simHit : moduleSimHits) {
0091       unsigned long paritcleId = simHit.particleId().value();
0092       if (!moduleSteps.contains(paritcleId)) {
0093         moduleSteps[paritcleId] = std::vector<Acts::detail::Step>();
0094       }
0095       double hx = simHit.fourPosition().x() / Acts::UnitConstants::mm;
0096       double hy = simHit.fourPosition().y() / Acts::UnitConstants::mm;
0097       double hz = simHit.fourPosition().z() / Acts::UnitConstants::mm;
0098       Acts::detail::Step step;
0099       step.position = Acts::Vector3(hx, hy, hz);
0100       step.momentum = simHit.direction();
0101       step.geoID = moduleGeoId;
0102       step.navDir = Acts::Direction::Forward();
0103       moduleSteps[paritcleId].push_back(step);
0104     }
0105     // Loop over and fill into the trackSteps
0106     for (const auto& [particleId, steps] : moduleSteps) {
0107       if (!trackSteps.contains(particleId)) {
0108         trackSteps[particleId] = std::vector<std::vector<Acts::detail::Step>>();
0109       }
0110       trackSteps[particleId].push_back(steps);
0111     }
0112   }
0113
0114   // Loop over the particles and create the propagation summaries
0115   for (const auto& particle : particles) {
0116     // Create the propagation summary
0117     Acts::CurvilinearTrackParameters start(
0118         particle.fourPosition(), particle.direction(),
0119         particle.charge() / particle.momentum().norm(), std::nullopt,
0120         particle.hypothesis());
0121     PropagationSummary propagationSummary(start);
0122     // Find the associated steps
0123     auto steps = trackSteps.find(particle.particleId().value());
0124     if (steps != trackSteps.end()) {
0125       for (const std::vector<Acts::detail::Step>& moduleSteps : steps->second) {
0126         // Get the GeometryIdentifier of the surface
0127         Acts::GeometryIdentifier surface = moduleSteps.front().geoID;
0128         // Find the surface
0129         auto surfaceIt = m_cfg.surfaceByIdentifier.find(surface);
0130         if (surfaceIt == m_cfg.surfaceByIdentifier.end()) {
0131           throw std::invalid_argument("Surface not found, should not happen");
0132         }
0133         Acts::detail::Step concatStep = concatenateSteps(
0134             context.geoContext, moduleSteps, *surfaceIt->second);
0135         propagationSummary.steps.push_back(concatStep);
0136       }
0137     } else {
0138       ACTS_WARNING("No steps found for particle "
0139                    << particle.particleId().value());
0140     }
0141
0142     propagationSummaries.push_back(std::move(propagationSummary));
0143   }
0144
0145   // Write the propagation step data to the event store
0146   m_outputSummary(context, std::move(propagationSummaries));
0147
0148   return ProcessCode::SUCCESS;
0149 }