Acts/Propagator/MaterialInteractor.hpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2019 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 #pragma once
0010
0011 #include "Acts/Definitions/Units.hpp"
0012 #include "Acts/Detector/GeometryCompatibilityConcept.hpp"
0013 #include "Acts/Geometry/TrackingVolume.hpp"
0014 #include "Acts/Material/MaterialInteraction.hpp"
0015 #include "Acts/Material/MaterialSlab.hpp"
0016 #include "Acts/Propagator/Propagator.hpp"
0017 #include "Acts/Propagator/detail/PointwiseMaterialInteraction.hpp"
0018 #include "Acts/Propagator/detail/VolumeMaterialInteraction.hpp"
0019 #include "Acts/Surfaces/Surface.hpp"
0020
0021 #include <sstream>
0022
0023 namespace Acts {
0024
0025 /// Material interactor propagator action.
0026 ///
0027 /// Apply material interactions at a surface and update the track state.
0028 struct MaterialInteractor {
0029   /// Whether to consider multiple scattering.
0030   bool multipleScattering = true;
0031   /// Whether to consider energy loss.
0032   bool energyLoss = true;
0033   /// Whether to record all material interactions.
0034   bool recordInteractions = false;
0035   /// Whether to add or remove noise.
0036   NoiseUpdateMode noiseUpdateMode = NoiseUpdateMode::addNoise;
0037
0038   using result_type = RecordedMaterial;
0039
0040   /// @brief Interaction with detector material for the ActionList
0041   /// of the Propagator
0042   ///
0043   /// It checks if the state has a current surface, in which case
0044   /// the action is performed: the covariance is transported to the position,
0045   /// multiple scattering and energy loss is applied  according to the
0046   /// configuration.
0047   ///
0048   /// @tparam propagator_state_t is the type of Propagator state
0049   /// @tparam stepper_t Type of the stepper of the propagation
0050   /// @tparam navigator_t Type of the navigator of the propagation
0051   ///
0052   /// @param state is the mutable propagator state object
0053   /// @param stepper The stepper in use
0054   /// @param navigator The navigator in use
0055   /// @param result is the mutable result state object
0056   /// @param logger a logger instance
0057   template <typename propagator_state_t, typename stepper_t,
0058             typename navigator_t>
0059   void operator()(propagator_state_t& state, const stepper_t& stepper,
0060                   const navigator_t& navigator, result_type& result,
0061                   const Logger& logger) const {
0062     if (state.stage == PropagatorStage::postPropagation) {
0063       return;
0064     }
0065
0066     // Do nothing if nothing is what is requested.
0067     if (!(multipleScattering || energyLoss || recordInteractions)) {
0068       return;
0069     }
0070
0071     static_assert(
0072         Acts::Concepts::NavigationCompatibilityConcept<propagator_state_t,
0073                                                        navigator_t>,
0074         "Navigation does not fulfill geometry compatibility concept");
0075
0076     // Handle surface material
0077
0078     // Note that start and target surface conditions are handled in the
0079     // interaction code
0080     const Surface* surface = navigator.currentSurface(state.navigation);
0081
0082     // We only have material interactions if there is potential material
0083     if (surface && surface->surfaceMaterial()) {
0084       ACTS_VERBOSE("MaterialInteractor | "
0085                    << "Found material on surface " << surface->geometryId());
0086
0087       // Prepare relevant input particle properties
0088       detail::PointwiseMaterialInteraction interaction(surface, state, stepper);
0089
0090       // Determine the effective traversed material and its properties
0091       // Material exists but it's not real, i.e. vacuum; there is nothing to do
0092       if (interaction.evaluateMaterialSlab(state, navigator)) {
0093         // Evaluate the material effects
0094         interaction.evaluatePointwiseMaterialInteraction(multipleScattering,
0095                                                          energyLoss);
0096
0097         if (energyLoss) {
0098           using namespace UnitLiterals;
0099           ACTS_VERBOSE("MaterialInteractor | "
0100                        << interaction.slab << " absPdg=" << interaction.absPdg
0101                        << " mass=" << interaction.mass / 1_MeV << "MeV"
0102                        << " momentum=" << interaction.momentum / 1_GeV << "GeV"
0103                        << " energyloss=" << interaction.Eloss / 1_MeV << "MeV");
0104         }
0105
0106         // To integrate process noise, we need to transport
0107         // the covariance to the current position in space
0108         if (interaction.performCovarianceTransport) {
0109           stepper.transportCovarianceToCurvilinear(state.stepping);
0110         }
0111         // Apply the material interactions
0112         interaction.updateState(state, stepper, noiseUpdateMode);
0113
0114         // Record the result
0115         recordResult(interaction, result);
0116       }
0117     }
0118
0119     // Handle volume material
0120
0121     // In case of Volume material update the result of the previous step
0122     if (!result.materialInteractions.empty() &&
0123         !result.materialInteractions.back().volume.empty() &&
0124         result.materialInteractions.back().updatedVolumeStep == false) {
0125       updateResult(state, stepper, result);
0126     }
0127
0128     auto volume = navigator.currentVolume(state.navigation);
0129
0130     // We only have material interactions if there is potential material
0131     if (volume && volume->volumeMaterial()) {
0132       ACTS_VERBOSE("MaterialInteractor | "
0133                    << "Found material in volume " << volume->geometryId());
0134
0135       // Prepare relevant input particle properties
0136       detail::VolumeMaterialInteraction interaction(volume, state, stepper);
0137       // Determine the effective traversed material and its properties
0138       // Material exists but it's not real, i.e. vacuum; there is nothing to do
0139       if (interaction.evaluateMaterialSlab(state, navigator)) {
0140         // Record the result
0141         recordResult(interaction, result);
0142       }
0143     }
0144   }
0145
0146  private:
0147   /// @brief This function records the material effect
0148   ///
0149   /// @param [in] interaction Interaction cache container
0150   /// @param [in, out] result Result storage
0151   void recordResult(const detail::PointwiseMaterialInteraction& interaction,
0152                     result_type& result) const {
0153     result.materialInX0 += interaction.slab.thicknessInX0();
0154     result.materialInL0 += interaction.slab.thicknessInL0();
0155
0156     // Record the interaction if requested
0157     if (!recordInteractions) {
0158       return;
0159     }
0160
0161     MaterialInteraction mi;
0162     mi.position = interaction.pos;
0163     mi.time = interaction.time;
0164     mi.direction = interaction.dir;
0165     mi.deltaP = interaction.nextP - interaction.momentum;
0166     mi.sigmaPhi2 = interaction.variancePhi;
0167     mi.sigmaTheta2 = interaction.varianceTheta;
0168     mi.sigmaQoP2 = interaction.varianceQoverP;
0169     mi.surface = interaction.surface;
0170     mi.volume = InteractionVolume();
0171     mi.pathCorrection = interaction.pathCorrection;
0172     mi.materialSlab = interaction.slab;
0173     result.materialInteractions.push_back(std::move(mi));
0174   }
0175
0176   /// @brief This function records the material effect
0177   ///
0178   /// @param [in] interaction Interaction cache container
0179   /// @param [in, out] result Result storage
0180   void recordResult(const detail::VolumeMaterialInteraction& interaction,
0181                     result_type& result) const {
0182     // Record the interaction if requested
0183     if (!recordInteractions) {
0184       return;
0185     }
0186
0187     MaterialInteraction mi;
0188     mi.position = interaction.pos;
0189     mi.time = interaction.time;
0190     mi.direction = interaction.dir;
0191     mi.surface = nullptr;
0192     mi.volume = interaction.volume;
0193     mi.pathCorrection = interaction.pathCorrection;
0194     mi.materialSlab = interaction.slab;
0195     result.materialInteractions.push_back(std::move(mi));
0196   }
0197
0198   /// @brief This function update the previous material step
0199   ///
0200   /// @param [in,out] state The state object
0201   /// @param [in] stepper The stepper instance
0202   /// @param [in, out] result Result storage
0203   template <typename propagator_state_t, typename stepper_t>
0204   void updateResult(propagator_state_t& state, const stepper_t& stepper,
0205                     result_type& result) const {
0206     // Update the previous interaction if requested
0207     if (!recordInteractions) {
0208       return;
0209     }
0210
0211     Vector3 shift = stepper.position(state.stepping) -
0212                     result.materialInteractions.back().position;
0213     double momentum = stepper.direction(state.stepping).norm();
0214     result.materialInteractions.back().deltaP =
0215         momentum - result.materialInteractions.back().direction.norm();
0216     result.materialInteractions.back().materialSlab.scaleThickness(
0217         shift.norm());
0218     result.materialInteractions.back().updatedVolumeStep = true;
0219     result.materialInX0 +=
0220         result.materialInteractions.back().materialSlab.thicknessInX0();
0221     result.materialInL0 +=
0222         result.materialInteractions.back().materialSlab.thicknessInL0();
0223   }
0224 };
0225
0226 }  // namespace Acts