Physics/Scattering/GeneralMixture.hpp

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