Propagation/src/PropagationAlgorithm.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/Propagation/PropagationAlgorithm.hpp"
0010
0011 #include "Acts/EventData/GenericBoundTrackParameters.hpp"
0012 #include "Acts/EventData/TrackParameters.hpp"
0013 #include "Acts/Surfaces/PerigeeSurface.hpp"
0014 #include "Acts/Surfaces/Surface.hpp"
0015 #include "ActsExamples/Framework/AlgorithmContext.hpp"
0016 #include "ActsExamples/Propagation/PropagatorInterface.hpp"
0017
0018 #include <stdexcept>
0019
0020 namespace ActsExamples {
0021
0022 ProcessCode PropagationAlgorithm::execute(
0023     const AlgorithmContext& context) const {
0024   // Create a random number generator
0025   ActsExamples::RandomEngine rng =
0026       m_cfg.randomNumberSvc->spawnGenerator(context);
0027
0028   // Standard gaussian distribution for covarianmces
0029   std::normal_distribution<double> gauss(0., 1.);
0030
0031   // Setup random number distributions for some quantities
0032   std::uniform_real_distribution<double> phiDist(m_cfg.phiRange.first,
0033                                                  m_cfg.phiRange.second);
0034   std::uniform_real_distribution<double> etaDist(m_cfg.etaRange.first,
0035                                                  m_cfg.etaRange.second);
0036   std::uniform_real_distribution<double> ptDist(m_cfg.ptRange.first,
0037                                                 m_cfg.ptRange.second);
0038   std::uniform_real_distribution<double> qDist(0., 1.);
0039
0040   std::shared_ptr<const Acts::PerigeeSurface> surface =
0041       Acts::Surface::makeShared<Acts::PerigeeSurface>(
0042           Acts::Vector3(0., 0., 0.));
0043
0044   // Output : the propagation steps
0045   std::vector<std::vector<Acts::detail::Step>> propagationSteps;
0046   propagationSteps.reserve(m_cfg.ntests);
0047
0048   // Output (optional): the recorded material
0049   std::unordered_map<std::size_t, Acts::RecordedMaterialTrack> recordedMaterial;
0050
0051   // loop over number of particles
0052   for (std::size_t it = 0; it < m_cfg.ntests; ++it) {
0053     /// get the d0 and z0
0054     double d0 = m_cfg.d0Sigma * gauss(rng);
0055     double z0 = m_cfg.z0Sigma * gauss(rng);
0056     double phi = phiDist(rng);
0057     double eta = etaDist(rng);
0058     double theta = 2 * std::atan(std::exp(-eta));
0059     double pt = ptDist(rng);
0060     double p = pt / std::sin(theta);
0061     double charge = qDist(rng) > 0.5 ? 1. : -1.;
0062     double qop = charge / p;
0063     double t = m_cfg.tSigma * gauss(rng);
0064     // parameters
0065     Acts::BoundVector pars;
0066     pars << d0, z0, phi, theta, qop, t;
0067     // some screen output
0068
0069     // The covariance generation
0070     auto cov = generateCovariance(rng, gauss);
0071
0072     // execute the test for charged particles
0073     Acts::BoundTrackParameters startParameters(surface, pars, std::move(cov),
0074                                                m_cfg.particleHypothesis);
0075     Acts::Vector3 sPosition = startParameters.position(context.geoContext);
0076     Acts::Vector3 sMomentum = startParameters.momentum();
0077     PropagationOutput pOutput = m_cfg.propagatorImpl->execute(
0078         context, m_cfg, logger(), startParameters);
0079     // Record the propagator steps
0080     propagationSteps.push_back(std::move(pOutput.first));
0081     if (m_cfg.recordMaterialInteractions &&
0082         !pOutput.second.materialInteractions.empty()) {
0083       // Create a recorded material track with start position, momentum and the
0084       // material
0085       recordedMaterial.emplace(
0086           it, std::make_pair(std::make_pair(sPosition, sMomentum),
0087                              std::move(pOutput.second)));
0088     }
0089   }
0090
0091   // Write the propagation step data to the event store
0092   m_outpoutPropagationSteps(context, std::move(propagationSteps));
0093
0094   // Write the recorded material to the event store
0095   if (m_cfg.recordMaterialInteractions) {
0096     m_recordedMaterial(context, std::move(recordedMaterial));
0097   }
0098
0099   return ProcessCode::SUCCESS;
0100 }
0101
0102 std::optional<Acts::BoundSquareMatrix> PropagationAlgorithm::generateCovariance(
0103     ActsExamples::RandomEngine& rnd,
0104     std::normal_distribution<double>& gauss) const {
0105   if (m_cfg.covarianceTransport) {
0106     // We start from the correlation matrix
0107     Acts::BoundSquareMatrix newCov(m_cfg.correlations);
0108     // Then we draw errors according to the error values
0109     Acts::BoundVector covs_smeared = m_cfg.covariances;
0110     for (std::size_t k = 0; k < std::size_t(covs_smeared.size()); ++k) {
0111       covs_smeared[k] *= gauss(rnd);
0112     }
0113     // and apply a double loop
0114     for (std::size_t i = 0; i < std::size_t(newCov.rows()); ++i) {
0115       for (std::size_t j = 0; j < std::size_t(newCov.cols()); ++j) {
0116         (newCov)(i, j) *= covs_smeared[i];
0117         (newCov)(i, j) *= covs_smeared[j];
0118       }
0119     }
0120     return newCov;
0121   }
0122   return std::nullopt;
0123 }
0124
0125 PropagationAlgorithm::PropagationAlgorithm(
0126     const PropagationAlgorithm::Config& config, Acts::Logging::Level level)
0127     : IAlgorithm("PropagationAlgorithm", level), m_cfg(config) {
0128   if (!m_cfg.propagatorImpl) {
0129     throw std::invalid_argument("Config needs to contain a propagator");
0130   }
0131   if (!m_cfg.randomNumberSvc) {
0132     throw std::invalid_argument("No random number generator given");
0133   }
0134
0135   m_outpoutPropagationSteps.initialize(m_cfg.propagationStepCollection);
0136   m_recordedMaterial.initialize(m_cfg.propagationMaterialCollection);
0137 }
0138
0139 }  // namespace ActsExamples