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