ElectroMagnetic/detail/GeneralMixture.hpp

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