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