Physics/ElectroMagnetic/Scattering.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/Algebra.hpp"
0012 #include "Acts/Material/MaterialSlab.hpp"
0013 #include "Acts/Utilities/UnitVectors.hpp"
0014 #include "ActsFatras/EventData/Particle.hpp"
0015 #include "ActsFatras/Physics/ElectroMagnetic/detail/GaussianMixture.hpp"
0016 #include "ActsFatras/Physics/ElectroMagnetic/detail/GeneralMixture.hpp"
0017 #include "ActsFatras/Physics/ElectroMagnetic/detail/Highland.hpp"
0018
0019 #include <array>
0020 #include <numbers>
0021 #include <random>
0022
0023 namespace ActsFatras {
0024
0025 /// Simulate (multiple) scattering using a configurable scattering model.
0026 ///
0027 /// @tparam scattering_model_t Model implementation to draw a scattering angle.
0028 template <typename scattering_model_t>
0029 struct GenericScattering {
0030   /// The scattering formula
0031   scattering_model_t angle;
0032
0033   /// Simulate scattering and update the particle parameters.
0034   ///
0035   /// @param[in]     generator is the random number generator
0036   /// @param[in]     slab      defines the passed material
0037   /// @param[in,out] particle  is the particle being updated
0038   /// @return Empty secondaries containers.
0039   ///
0040   /// @tparam generator_t is a RandomNumberEngine
0041   template <typename generator_t>
0042   std::array<Particle, 0> operator()(generator_t &generator,
0043                                      const Acts::MaterialSlab &slab,
0044                                      Particle &particle) const {
0045     // the scattered direction can be computed by rotating the initial
0046     // direction around a vector orthogonal to the initial direction, i.e. the
0047     // scattering deflector, by the scattering angle. there are an infinite
0048     // number of vectors orthogonal to the initial direction. the deflector is
0049     // rotated by some angle relative to some fixpoint.
0050     //
0051     // thus two random angles are required: the random deflector orientation
0052     // angle drawn uniformly from the [-pi,pi) range and the scattering angle
0053     // drawn from the specific scattering model distribution.
0054
0055     // draw the random orientation angle
0056     const double psi = std::uniform_real_distribution<double>(
0057         -std::numbers::pi, std::numbers::pi)(generator);
0058     // draw the scattering angle
0059     const double theta = angle(generator, slab, particle);
0060
0061     // construct the combined rotation to the scattered direction
0062     const Acts::RotationMatrix3 rotation(
0063         // rotation of the scattering deflector axis relative to the reference
0064         Acts::AngleAxis3(psi, particle.direction()) *
0065         // rotation by the scattering angle around the deflector axis
0066         Acts::AngleAxis3(theta,
0067                          Acts::createCurvilinearUnitU(particle.direction())));
0068     particle.setDirection(rotation * particle.direction());
0069
0070     // scattering is non-destructive and produces no secondaries
0071     return {};
0072   }
0073 };
0074
0075 /// Scattering with Gaussian mixture model
0076 using GaussianMixtureScattering = GenericScattering<detail::GaussianMixture>;
0077 /// Scattering with general mixture model
0078 using GeneralMixtureScattering = GenericScattering<detail::GeneralMixture>;
0079 /// Scattering with Highland model
0080 using HighlandScattering = GenericScattering<detail::Highland>;
0081
0082 }  // namespace ActsFatras