EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
GeneralMixture.hpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file GeneralMixture.hpp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2018-2020 CERN for the benefit of the Acts project
4 //
5 // This Source Code Form is subject to the terms of the Mozilla Public
6 // License, v. 2.0. If a copy of the MPL was not distributed with this
7 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
8 
9 #pragma once
10 
11 #include "Acts/Material/Interactions.hpp"
12 #include "Acts/Utilities/PdgParticle.hpp"
13 #include "ActsFatras/Physics/Scattering/detail/Scattering.hpp"
14 
15 #include <random>
16 
17 namespace ActsFatras {
18 namespace detail {
19 
27 struct GeneralMixture {
29  bool log_include = true;
31  double genMixtureScalor = 1.;
32 
41  template <typename generator_t>
42  double operator()(generator_t &generator, const Acts::MaterialSlab &slab,
43  Particle &particle) const {
44  double theta = 0.0;
45 
46  if (std::abs(particle.pdg()) != Acts::PdgParticle::eElectron) {
47  //----------------------------------------------------------------------------
48  // see Mixture models of multiple scattering: computation and simulation.
49  // -
50  // R.Fruehwirth, M. Liendl. -
51  // Computer Physics Communications 141 (2001) 230–246
52  //----------------------------------------------------------------------------
53  std::array<double, 4> scattering_params;
54  // Decide which mixture is best
55  // beta² = (p/E)² = p²/(p² + m²) = 1/(1 + (m/p)²)
56  // 1/beta² = 1 + (m/p)²
57  // beta = 1/sqrt(1 + (m/p)²)
58  double mOverP = particle.mass() / particle.absMomentum();
59  double beta2Inv = 1 + mOverP * mOverP;
60  double beta = 1 / std::sqrt(beta2Inv);
61  double tInX0 = slab.thicknessInX0();
62  double tob2 = tInX0 * beta2Inv;
63  if (tob2 > 0.6 / std::pow(slab.material().Z(), 0.6)) {
64  // Gaussian mixture or pure Gaussian
65  if (tob2 > 10) {
66  scattering_params = getGaussian(beta, particle.absMomentum(), tInX0,
67  genMixtureScalor);
68  } else {
69  scattering_params =
70  getGaussmix(beta, particle.absMomentum(), tInX0,
71  slab.material().Z(), genMixtureScalor);
72  }
73  // Simulate
74  theta = gaussmix(generator, scattering_params);
75  } else {
76  // Semigaussian mixture - get parameters
77  auto scattering_params_sg =
78  getSemigauss(beta, particle.absMomentum(), tInX0,
79  slab.material().Z(), genMixtureScalor);
80  // Simulate
81  theta = semigauss(generator, scattering_params_sg);
82  }
83  } else {
84  // for electrons we fall back to the Highland (extension)
85  // return projection factor times sigma times gauss random
86  const auto theta0 = Acts::computeMultipleScatteringTheta0(
87  slab, particle.pdg(), particle.mass(),
88  particle.charge() / particle.absMomentum(), particle.charge());
89  theta = std::normal_distribution<double>(0.0, theta0)(generator);
90  }
91  // scale from planar to 3d angle
92  return M_SQRT2 * theta;
93  }
94 
95  // helper methods for getting parameters and simulating
96 
97  std::array<double, 4> getGaussian(double beta, double p, double tInX0,
98  double scale) const {
99  std::array<double, 4> scattering_params;
100  // Total standard deviation of mixture
101  scattering_params[0] = 15. / beta / p * std::sqrt(tInX0) * scale;
102  scattering_params[1] = 1.0; // Variance of core
103  scattering_params[2] = 1.0; // Variance of tails
104  scattering_params[3] = 0.5; // Mixture weight of tail component
105  return scattering_params;
106  }
107 
108  std::array<double, 4> getGaussmix(double beta, double p, double tInX0,
109  double Z, double scale) const {
110  std::array<double, 4> scattering_params;
111  scattering_params[0] = 15. / beta / p * std::sqrt(tInX0) *
112  scale; // Total standard deviation of mixture
113  double d1 = std::log(tInX0 / (beta * beta));
114  double d2 = std::log(std::pow(Z, 2.0 / 3.0) * tInX0 / (beta * beta));
115  double epsi;
116  double var1 = (-1.843e-3 * d1 + 3.347e-2) * d1 + 8.471e-1; // Variance of
117  // core
118  if (d2 < 0.5)
119  epsi = (6.096e-4 * d2 + 6.348e-3) * d2 + 4.841e-2;
120  else
121  epsi = (-5.729e-3 * d2 + 1.106e-1) * d2 - 1.908e-2;
122  scattering_params[1] = var1; // Variance of core
123  scattering_params[2] = (1 - (1 - epsi) * var1) / epsi; // Variance of tails
124  scattering_params[3] = epsi; // Mixture weight of tail component
125  return scattering_params;
126  }
127 
128  std::array<double, 6> getSemigauss(double beta, double p, double tInX0,
129  double Z, double scale) const {
130  std::array<double, 6> scattering_params;
131  double N = tInX0 * 1.587E7 * std::pow(Z, 1.0 / 3.0) / (beta * beta) /
132  (Z + 1) / std::log(287 / std::sqrt(Z));
133  scattering_params[4] = 15. / beta / p * std::sqrt(tInX0) *
134  scale; // Total standard deviation of mixture
135  double rho = 41000 / std::pow(Z, 2.0 / 3.0);
136  double b = rho / std::sqrt(N * (std::log(rho) - 0.5));
137  double n = std::pow(Z, 0.1) * std::log(N);
138  double var1 = (5.783E-4 * n + 3.803E-2) * n + 1.827E-1;
139  double a =
140  (((-4.590E-5 * n + 1.330E-3) * n - 1.355E-2) * n + 9.828E-2) * n +
141  2.822E-1;
142  double epsi = (1 - var1) / (a * a * (std::log(b / a) - 0.5) - var1);
143  scattering_params[3] =
144  (epsi > 0) ? epsi : 0.0; // Mixture weight of tail component
145  scattering_params[0] = a; // Parameter 1 of tails
146  scattering_params[1] = b; // Parameter 2 of tails
147  scattering_params[2] = var1; // Variance of core
148  scattering_params[5] = N; // Average number of scattering processes
149  return scattering_params;
150  }
151 
160  template <typename generator_t>
161  double gaussmix(generator_t &generator,
162  const std::array<double, 4> &scattering_params) const {
163  std::uniform_real_distribution<double> udist(0.0, 1.0);
164  double sigma_tot = scattering_params[0];
165  double var1 = scattering_params[1];
166  double var2 = scattering_params[2];
167  double epsi = scattering_params[3];
168  bool ind = udist(generator) > epsi;
169  double u = udist(generator);
170  if (ind)
171  return std::sqrt(var1) * std::sqrt(-2 * std::log(u)) * sigma_tot;
172  else
173  return std::sqrt(var2) * std::sqrt(-2 * std::log(u)) * sigma_tot;
174  }
175 
184  template <typename generator_t>
185  double semigauss(generator_t &generator,
186  const std::array<double, 6> &scattering_params) const {
187  std::uniform_real_distribution<double> udist(0.0, 1.0);
188  double a = scattering_params[0];
189  double b = scattering_params[1];
190  double var1 = scattering_params[2];
191  double epsi = scattering_params[3];
192  double sigma_tot = scattering_params[4];
193  bool ind = udist(generator) > epsi;
194  double u = udist(generator);
195  if (ind)
196  return std::sqrt(var1) * std::sqrt(-2 * std::log(u)) * sigma_tot;
197  else
198  return a * b * std::sqrt((1 - u) / (u * b * b + a * a)) * sigma_tot;
199  }
200 };
201 
202 } // namespace detail
203 
204 using GeneralMixtureScattering = detail::Scattering<detail::GeneralMixture>;
205 
206 } // namespace ActsFatras