ElectroMagnetic/detail/GaussianMixture.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/Material/Interactions.hpp"
0012
0013 #include <random>
0014
0015 namespace ActsFatras::detail {
0016
0017 /// Generate scattering angles using a Gaussian mixture model.
0018 struct GaussianMixture {
0019   /// Steering parameter
0020   bool optGaussianMixtureG4 = false;
0021   double gausMixSigma1_a0 = 8.471e-1;
0022   double gausMixSigma1_a1 = 3.347e-2;
0023   double gausMixSigma1_a2 = -1.843e-3;
0024   double gausMixEpsilon_a0 = 4.841e-2;
0025   double gausMixEpsilon_a1 = 6.348e-3;
0026   double gausMixEpsilon_a2 = 6.096e-4;
0027   double gausMixEpsilon_b0 = -1.908e-2;
0028   double gausMixEpsilon_b1 = 1.106e-1;
0029   double gausMixEpsilon_b2 = -5.729e-3;
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     /// Calculate the highland formula first
0043     double sigma = Acts::computeMultipleScatteringTheta0(
0044         slab, particle.absolutePdg(), particle.mass(), particle.qOverP(),
0045         particle.absoluteCharge());
0046     double sigma2 = sigma * sigma;
0047
0048     // Gauss distribution, will be sampled with generator
0049     std::normal_distribution<double> gaussDist(0., 1.);
0050     // Uniform distribution, will be sampled with generator
0051     std::uniform_real_distribution<double> uniformDist(0., 1.);
0052
0053     // Now correct for the tail fraction
0054     // d_0'
0055     // beta² = (p/E)² = p²/(p² + m²) = 1/(1 + (m/p)²)
0056     // 1/beta² = 1 + (m/p)²
0057     double mOverP = particle.mass() / particle.absoluteMomentum();
0058     double beta2inv = 1 + mOverP * mOverP;
0059     double dprime = slab.thicknessInX0() * beta2inv;
0060     double log_dprime = std::log(dprime);
0061     // d_0''
0062     double log_dprimeprime =
0063         std::log(std::pow(slab.material().Z(), 2.0 / 3.0) * dprime);
0064
0065     // get epsilon
0066     double epsilon =
0067         log_dprimeprime < 0.5
0068             ? gausMixEpsilon_a0 + gausMixEpsilon_a1 * log_dprimeprime +
0069                   gausMixEpsilon_a2 * log_dprimeprime * log_dprimeprime
0070             : gausMixEpsilon_b0 + gausMixEpsilon_b1 * log_dprimeprime +
0071                   gausMixEpsilon_b2 * log_dprimeprime * log_dprimeprime;
0072
0073     // the standard sigma
0074     double sigma1square = gausMixSigma1_a0 + gausMixSigma1_a1 * log_dprime +
0075                           gausMixSigma1_a2 * log_dprime * log_dprime;
0076
0077     // G4 optimised / native double Gaussian model
0078     if (optGaussianMixtureG4) {
0079       sigma2 = 225. * dprime /
0080                (particle.absoluteMomentum() * particle.absoluteMomentum());
0081     }
0082     // throw the random number core/tail
0083     if (uniformDist(generator) < epsilon) {
0084       sigma2 *= (1. - (1. - epsilon) * sigma1square) / epsilon;
0085     }
0086     // return back to the
0087     return M_SQRT2 * std::sqrt(sigma2) * gaussDist(generator);
0088   }
0089 };
0090
0091 }  // namespace ActsFatras::detail