Physics/Scattering/GaussianMixture.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 /// @brief The struct to be provided to the Scatterer action
0017 /// This is the gaussian mixture
0018 struct GaussianMixture {
0019
0020   /// Steering parameter
0021   bool log_include = true;
0022
0023   double gausMixSigma1_a0 = 8.471e-1;
0024   double gausMixSigma1_a1 = 3.347e-2;
0025   double gausMixSigma1_a2 = -1.843e-3;
0026   double gausMixEpsilon_a0 = 4.841e-2;
0027   double gausMixEpsilon_a1 = 6.348e-3;
0028   double gausMixEpsilon_a2 = 6.096e-4;
0029   double gausMixEpsilon_b0 = -1.908e-2;
0030   double gausMixEpsilon_b1 = 1.106e-1;
0031   double gausMixEpsilon_b2 = -5.729e-3;
0032
0033   bool optGaussianMixtureG4 = false;
0034
0035   /// @brief Call operator to perform this scattering
0036   ///
0037   /// @tparam generator_t is a random number generator type
0038   /// @tparam detector_t is the detector information type
0039   /// @tparam particle_t is the particle information type
0040   ///
0041   /// @param[in] generator is the random number generator
0042   /// @param[in] detector the detector information
0043   /// @param[in] particle the particle which is being scattered
0044   ///
0045   /// @return a scattering angle in 3D
0046   template <typename generator_t, typename detector_t, typename particle_t>
0047   double operator()(generator_t &generator, const detector_t &detector,
0048                     particle_t &particle) const {
0049
0050     // thickness in X0
0051     double dInX0 = detector.thickness() / detector.material().X0();
0052
0053     /// Calculate the highland formula first
0054     double qop = particle.q() / particle.p();
0055     double sigma = Acts::computeMultipleScatteringTheta0(
0056         detector, particle.pdg(), particle.m(), qop);
0057
0058     double sigma2 = sigma * sigma;
0059
0060     // Gauss distribution, will be sampled with generator
0061     GaussDist gaussDist = GaussDist(0., 1.);
0062
0063     // Uniform distribution, will be sampled with generator
0064     UniformDist uniformDist = UniformDist(0., 1.);
0065
0066     // Now correct for the tail fraction
0067     // d_0'
0068     double beta2 = particle.beta() * particle.beta();
0069     double dprime = detector.thickness() / (detector.material().X0() * beta2);
0070     double log_dprime = std::log(dprime);
0071     // d_0''
0072     double log_dprimeprime =
0073         std::log(std::pow(detector.material().Z(), 2.0 / 3.0) * dprime);
0074
0075     // get epsilon
0076     double epsilon =
0077         log_dprimeprime < 0.5
0078             ? gausMixEpsilon_a0 + gausMixEpsilon_a1 * log_dprimeprime +
0079                   gausMixEpsilon_a2 * log_dprimeprime * log_dprimeprime
0080             : gausMixEpsilon_b0 + gausMixEpsilon_b1 * log_dprimeprime +
0081                   gausMixEpsilon_b2 * log_dprimeprime * log_dprimeprime;
0082
0083     // the standard sigma
0084     double sigma1square = gausMixSigma1_a0 + gausMixSigma1_a1 * log_dprime +
0085                           gausMixSigma1_a2 * log_dprime * log_dprime;
0086
0087     // G4 optimised / native double Gaussian model
0088     if (optGaussianMixtureG4) {
0089       sigma2 = 225. * dprime / (particle.p() * particle.p());
0090     }
0091     // throw the random number core/tail
0092     if (uniformDist(generator) < epsilon) {
0093       sigma2 *= (1. - (1. - epsilon) * sigma1square) / epsilon;
0094     }
0095     // return back to the
0096     return M_SQRT2 * std::sqrt(sigma2) * gaussDist(generator);
0097   }
0098 };
0099
0100 } // namespace Fatras