ElectroMagnetic/detail/GeneralMixture.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/PdgParticle.hpp"
0012 #include "Acts/Material/Interactions.hpp"
0013 #include "ActsFatras/EventData/Particle.hpp"
0014
0015 #include <numbers>
0016 #include <random>
0017
0018 namespace ActsFatras::detail {
0019
0020 /// Generate scattering angles using a general mixture model.
0021 ///
0022 /// Emulates core and tail scattering as described in
0023 ///
0024 ///     General mixture model Fruehwirth, M. Liendl.
0025 ///     Comp. Phys. Comm. 141 (2001) 230-246
0026 ///
0027 struct GeneralMixture {
0028   /// Steering parameter
0029   bool logInclude = true;
0030   /// Scale the mixture level
0031   double genMixtureScalor = 1.;
0032
0033   /// Generate a single 3D scattering angle.
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 scattered
0038   /// @return a 3d scattering angle
0039   ///
0040   /// @tparam generator_t is a RandomNumberEngine
0041   template <typename generator_t>
0042   double operator()(generator_t &generator, const Acts::MaterialSlab &slab,
0043                     Particle &particle) const {
0044     double theta = 0.0;
0045
0046     if (particle.absolutePdg() != Acts::PdgParticle::eElectron) {
0047       //----------------------------------------------------------------------------
0048       // see Mixture models of multiple scattering: computation and simulation.
0049       // -
0050       // R.Fruehwirth, M. Liendl. -
0051       // Computer Physics Communications 141 (2001) 230â246
0052       //----------------------------------------------------------------------------
0053
0054       // Decide which mixture is best
0055       //   beta² = (p/E)² = p²/(p² + m²) = 1/(1 + (m/p)²)
0056       // 1/beta² = 1 + (m/p)²
0057       //    beta = 1/sqrt(1 + (m/p)²)
0058       const double mOverP = particle.mass() / particle.absoluteMomentum();
0059       const double beta2Inv = 1 + mOverP * mOverP;
0060       const double beta = 1 / std::sqrt(beta2Inv);
0061       const double tInX0 = slab.thicknessInX0();
0062       const double tob2 = tInX0 * beta2Inv;
0063       if (tob2 > 0.6 / std::pow(slab.material().Z(), 0.6)) {
0064         std::array<double, 4> scatteringParams{};
0065         // Gaussian mixture or pure Gaussian
0066         if (tob2 > 10) {
0067           scatteringParams = createGaussian(beta, particle.absoluteMomentum(),
0068                                             tInX0, genMixtureScalor);
0069         } else {
0070           scatteringParams =
0071               createGaussmix(beta, particle.absoluteMomentum(), tInX0,
0072                              slab.material().Z(), genMixtureScalor);
0073         }
0074         // Simulate
0075         theta = gaussmix(generator, scatteringParams);
0076       } else {
0077         // Semigaussian mixture - get parameters
0078         const std::array<double, 6> scatteringParamsSg =
0079             createSemigauss(beta, particle.absoluteMomentum(), tInX0,
0080                             slab.material().Z(), genMixtureScalor);
0081         // Simulate
0082         theta = semigauss(generator, scatteringParamsSg);
0083       }
0084     } else {
0085       // for electrons we fall back to the Highland (extension)
0086       // return projection factor times sigma times gauss random
0087       const auto theta0 = Acts::computeMultipleScatteringTheta0(
0088           slab, particle.absolutePdg(), particle.mass(), particle.qOverP(),
0089           particle.absoluteCharge());
0090       theta = std::normal_distribution<double>(0.0, theta0)(generator);
0091     }
0092     // scale from planar to 3d angle
0093     return std::numbers::sqrt2 * theta;
0094   }
0095
0096   // helper methods for getting parameters and simulating
0097
0098   std::array<double, 4> createGaussian(double beta, double p, double tInX0,
0099                                        double scale) const {
0100     std::array<double, 4> scatteringParams{};
0101     // Total standard deviation of mixture
0102     scatteringParams[0] = 15. / beta / p * std::sqrt(tInX0) * scale;
0103     // Variance of core
0104     scatteringParams[1] = 1.0;
0105     // Variance of tails
0106     scatteringParams[2] = 1.0;
0107     // Mixture weight of tail component
0108     scatteringParams[3] = 0.5;
0109     return scatteringParams;
0110   }
0111
0112   std::array<double, 4> createGaussmix(double beta, double p, double tInX0,
0113                                        double Z, double scale) const {
0114     std::array<double, 4> scatteringParams{};
0115     // Total standard deviation of mixture
0116     scatteringParams[0] = 15. / beta / p * std::sqrt(tInX0) * scale;
0117     const double d1 = std::log(tInX0 / (beta * beta));
0118     const double d2 = std::log(std::pow(Z, 2.0 / 3.0) * tInX0 / (beta * beta));
0119     const double var1 = (-1.843e-3 * d1 + 3.347e-2) * d1 + 8.471e-1;
0120     double epsi = 0;
0121     if (d2 < 0.5) {
0122       epsi = (6.096e-4 * d2 + 6.348e-3) * d2 + 4.841e-2;
0123     } else {
0124       epsi = (-5.729e-3 * d2 + 1.106e-1) * d2 - 1.908e-2;
0125     }
0126     // Variance of core
0127     scatteringParams[1] = var1;
0128     // Variance of tails
0129     scatteringParams[2] = (1 - (1 - epsi) * var1) / epsi;
0130     // Mixture weight of tail component
0131     scatteringParams[3] = epsi;
0132     return scatteringParams;
0133   }
0134
0135   std::array<double, 6> createSemigauss(double beta, double p, double tInX0,
0136                                         double Z, double scale) const {
0137     std::array<double, 6> scatteringParams{};
0138     const double N = tInX0 * 1.587E7 * std::pow(Z, 1.0 / 3.0) / (beta * beta) /
0139                      (Z + 1) / std::log(287 / std::sqrt(Z));
0140     // Total standard deviation of mixture
0141     scatteringParams[4] = 15. / beta / p * std::sqrt(tInX0) * scale;
0142     const double rho = 41000 / std::pow(Z, 2.0 / 3.0);
0143     const double b = rho / std::sqrt(N * (std::log(rho) - 0.5));
0144     const double n = std::pow(Z, 0.1) * std::log(N);
0145     const double var1 = (5.783E-4 * n + 3.803E-2) * n + 1.827E-1;
0146     const double a =
0147         (((-4.590E-5 * n + 1.330E-3) * n - 1.355E-2) * n + 9.828E-2) * n +
0148         2.822E-1;
0149     const double epsi = (1 - var1) / (a * a * (std::log(b / a) - 0.5) - var1);
0150     // Mixture weight of tail component
0151     scatteringParams[3] = (epsi > 0) ? epsi : 0.0;
0152     // Parameter 1 of tails
0153     scatteringParams[0] = a;
0154     // Parameter 2 of tails
0155     scatteringParams[1] = b;
0156     // Variance of core
0157     scatteringParams[2] = var1;
0158     // Average number of scattering processes
0159     scatteringParams[5] = N;
0160     return scatteringParams;
0161   }
0162
0163   /// @brief Retrieve the gaussian mixture
0164   ///
0165   /// @tparam generator_t Type of the generator
0166   ///
0167   /// @param udist The uniform distribution handed over by the call operator
0168   /// @param scattering_params the tuned parameters for the generation
0169   ///
0170   /// @return a double value that represents the gaussian mixture
0171   template <typename generator_t>
0172   double gaussmix(generator_t &generator,
0173                   const std::array<double, 4> &scatteringParams) const {
0174     std::uniform_real_distribution<double> udist(0.0, 1.0);
0175     const double sigmaTot = scatteringParams[0];
0176     const double var1 = scatteringParams[1];
0177     const double var2 = scatteringParams[2];
0178     const double epsi = scatteringParams[3];
0179     const bool ind = udist(generator) > epsi;
0180     const double u = udist(generator);
0181     if (ind) {
0182       return std::sqrt(var1) * std::sqrt(-2 * std::log(u)) * sigmaTot;
0183     } else {
0184       return std::sqrt(var2) * std::sqrt(-2 * std::log(u)) * sigmaTot;
0185     }
0186   }
0187
0188   /// @brief Retrieve the semi-gaussian mixture
0189   ///
0190   /// @tparam generator_t Type of the generator
0191   ///
0192   /// @param udist The uniform distribution handed over by the call operator
0193   /// @param scattering_params the tuned parameters for the generation
0194   ///
0195   /// @return a double value that represents the gaussian mixture
0196   template <typename generator_t>
0197   double semigauss(generator_t &generator,
0198                    const std::array<double, 6> &scattering_params) const {
0199     std::uniform_real_distribution<double> udist(0.0, 1.0);
0200     const double a = scattering_params[0];
0201     const double b = scattering_params[1];
0202     const double var1 = scattering_params[2];
0203     const double epsi = scattering_params[3];
0204     const double sigmaTot = scattering_params[4];
0205     const bool ind = udist(generator) > epsi;
0206     const double u = udist(generator);
0207     if (ind) {
0208       return std::sqrt(var1) * std::sqrt(-2 * std::log(u)) * sigmaTot;
0209     } else {
0210       return a * b * std::sqrt((1 - u) / (u * b * b + a * a)) * sigmaTot;
0211     }
0212   }
0213 };
0214
0215 }  // namespace ActsFatras::detail