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 log_include = 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       std::array<double, 4> scattering_params{};
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       double mOverP = particle.mass() / particle.absoluteMomentum();
0059       double beta2Inv = 1 + mOverP * mOverP;
0060       double beta = 1 / std::sqrt(beta2Inv);
0061       double tInX0 = slab.thicknessInX0();
0062       double tob2 = tInX0 * beta2Inv;
0063       if (tob2 > 0.6 / std::pow(slab.material().Z(), 0.6)) {
0064         // Gaussian mixture or pure Gaussian
0065         if (tob2 > 10) {
0066           scattering_params = getGaussian(beta, particle.absoluteMomentum(),
0067                                           tInX0, genMixtureScalor);
0068         } else {
0069           scattering_params =
0070               getGaussmix(beta, particle.absoluteMomentum(), tInX0,
0071                           slab.material().Z(), genMixtureScalor);
0072         }
0073         // Simulate
0074         theta = gaussmix(generator, scattering_params);
0075       } else {
0076         // Semigaussian mixture - get parameters
0077         auto scattering_params_sg =
0078             getSemigauss(beta, particle.absoluteMomentum(), tInX0,
0079                          slab.material().Z(), genMixtureScalor);
0080         // Simulate
0081         theta = semigauss(generator, scattering_params_sg);
0082       }
0083     } else {
0084       // for electrons we fall back to the Highland (extension)
0085       // return projection factor times sigma times gauss random
0086       const auto theta0 = Acts::computeMultipleScatteringTheta0(
0087           slab, particle.absolutePdg(), particle.mass(), particle.qOverP(),
0088           particle.absoluteCharge());
0089       theta = std::normal_distribution<double>(0.0, theta0)(generator);
0090     }
0091     // scale from planar to 3d angle
0092     return std::numbers::sqrt2 * theta;
0093   }
0094
0095   // helper methods for getting parameters and simulating
0096
0097   std::array<double, 4> getGaussian(double beta, double p, double tInX0,
0098                                     double scale) const {
0099     std::array<double, 4> scattering_params{};
0100     // Total standard deviation of mixture
0101     scattering_params[0] = 15. / beta / p * std::sqrt(tInX0) * scale;
0102     scattering_params[1] = 1.0;  // Variance of core
0103     scattering_params[2] = 1.0;  // Variance of tails
0104     scattering_params[3] = 0.5;  // Mixture weight of tail component
0105     return scattering_params;
0106   }
0107
0108   std::array<double, 4> getGaussmix(double beta, double p, double tInX0,
0109                                     double Z, double scale) const {
0110     std::array<double, 4> scattering_params{};
0111     scattering_params[0] = 15. / beta / p * std::sqrt(tInX0) *
0112                            scale;  // Total standard deviation of mixture
0113     double d1 = std::log(tInX0 / (beta * beta));
0114     double d2 = std::log(std::pow(Z, 2.0 / 3.0) * tInX0 / (beta * beta));
0115     double epsi = 0;
0116     double var1 = (-1.843e-3 * d1 + 3.347e-2) * d1 + 8.471e-1;  // Variance of
0117                                                                 // core
0118     if (d2 < 0.5) {
0119       epsi = (6.096e-4 * d2 + 6.348e-3) * d2 + 4.841e-2;
0120     } else {
0121       epsi = (-5.729e-3 * d2 + 1.106e-1) * d2 - 1.908e-2;
0122     }
0123     scattering_params[1] = var1;                            // Variance of core
0124     scattering_params[2] = (1 - (1 - epsi) * var1) / epsi;  // Variance of tails
0125     scattering_params[3] = epsi;  // Mixture weight of tail component
0126     return scattering_params;
0127   }
0128
0129   std::array<double, 6> getSemigauss(double beta, double p, double tInX0,
0130                                      double Z, double scale) const {
0131     std::array<double, 6> scattering_params{};
0132     double N = tInX0 * 1.587E7 * std::pow(Z, 1.0 / 3.0) / (beta * beta) /
0133                (Z + 1) / std::log(287 / std::sqrt(Z));
0134     scattering_params[4] = 15. / beta / p * std::sqrt(tInX0) *
0135                            scale;  // Total standard deviation of mixture
0136     double rho = 41000 / std::pow(Z, 2.0 / 3.0);
0137     double b = rho / std::sqrt(N * (std::log(rho) - 0.5));
0138     double n = std::pow(Z, 0.1) * std::log(N);
0139     double var1 = (5.783E-4 * n + 3.803E-2) * n + 1.827E-1;
0140     double a =
0141         (((-4.590E-5 * n + 1.330E-3) * n - 1.355E-2) * n + 9.828E-2) * n +
0142         2.822E-1;
0143     double epsi = (1 - var1) / (a * a * (std::log(b / a) - 0.5) - var1);
0144     scattering_params[3] =
0145         (epsi > 0) ? epsi : 0.0;  // Mixture weight of tail component
0146     scattering_params[0] = a;     // Parameter 1 of tails
0147     scattering_params[1] = b;     // Parameter 2 of tails
0148     scattering_params[2] = var1;  // Variance of core
0149     scattering_params[5] = N;     // Average number of scattering processes
0150     return scattering_params;
0151   }
0152
0153   /// @brief Retrieve the gaussian mixture
0154   ///
0155   /// @tparam generator_t Type of the generator
0156   ///
0157   /// @param udist The uniform distribution handed over by the call operator
0158   /// @param scattering_params the tuned parameters for the generation
0159   ///
0160   /// @return a double value that represents the gaussian mixture
0161   template <typename generator_t>
0162   double gaussmix(generator_t &generator,
0163                   const std::array<double, 4> &scattering_params) const {
0164     std::uniform_real_distribution<double> udist(0.0, 1.0);
0165     double sigma_tot = scattering_params[0];
0166     double var1 = scattering_params[1];
0167     double var2 = scattering_params[2];
0168     double epsi = scattering_params[3];
0169     bool ind = udist(generator) > epsi;
0170     double u = udist(generator);
0171     if (ind) {
0172       return std::sqrt(var1) * std::sqrt(-2 * std::log(u)) * sigma_tot;
0173     } else {
0174       return std::sqrt(var2) * std::sqrt(-2 * std::log(u)) * sigma_tot;
0175     }
0176   }
0177
0178   /// @brief Retrieve the semi-gaussian mixture
0179   ///
0180   /// @tparam generator_t Type of the generator
0181   ///
0182   /// @param udist The uniform distribution handed over by the call operator
0183   /// @param scattering_params the tuned parameters for the generation
0184   ///
0185   /// @return a double value that represents the gaussian mixture
0186   template <typename generator_t>
0187   double semigauss(generator_t &generator,
0188                    const std::array<double, 6> &scattering_params) const {
0189     std::uniform_real_distribution<double> udist(0.0, 1.0);
0190     double a = scattering_params[0];
0191     double b = scattering_params[1];
0192     double var1 = scattering_params[2];
0193     double epsi = scattering_params[3];
0194     double sigma_tot = scattering_params[4];
0195     bool ind = udist(generator) > epsi;
0196     double u = udist(generator);
0197     if (ind) {
0198       return std::sqrt(var1) * std::sqrt(-2 * std::log(u)) * sigma_tot;
0199     } else {
0200       return a * b * std::sqrt((1 - u) / (u * b * b + a * a)) * sigma_tot;
0201     }
0202   }
0203 };
0204
0205 }  // namespace ActsFatras::detail