Acts/Propagator/MaterialInteractor.hpp

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