ECCE @ 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 <random>
12 
13 #include "Acts/Material/Interactions.hpp"
14 #include "Acts/Utilities/PdgParticle.hpp"
15 #include "ActsFatras/Physics/Scattering/detail/Scattering.hpp"
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,
43  const Acts::MaterialProperties &slab,
44  Particle &particle) const {
45  double theta = 0.0;
46 
47  if (std::abs(particle.pdg()) != Acts::PdgParticle::eElectron) {
48  //----------------------------------------------------------------------------
49  // see Mixture models of multiple scattering: computation and simulation.
50  // -
51  // R.Fruehwirth, M. Liendl. -
52  // Computer Physics Communications 141 (2001) 230–246
53  //----------------------------------------------------------------------------
54  std::array<double, 4> scattering_params;
55  // Decide which mixture is best
56  // beta² = (p/E)² = p²/(p² + m²) = 1/(1 + (m/p)²)
57  // 1/beta² = 1 + (m/p)²
58  // beta = 1/sqrt(1 + (m/p)²)
59  double mOverP = particle.mass() / particle.absMomentum();
60  double beta2Inv = 1 + mOverP * mOverP;
61  double beta = 1 / std::sqrt(beta2Inv);
62  double tInX0 = slab.thicknessInX0();
63  double tob2 = tInX0 * beta2Inv;
64  if (tob2 > 0.6 / std::pow(slab.material().Z(), 0.6)) {
65  // Gaussian mixture or pure Gaussian
66  if (tob2 > 10) {
67  scattering_params = getGaussian(beta, particle.absMomentum(), tInX0,
68  genMixtureScalor);
69  } else {
70  scattering_params =
71  getGaussmix(beta, particle.absMomentum(), tInX0,
72  slab.material().Z(), genMixtureScalor);
73  }
74  // Simulate
75  theta = gaussmix(generator, scattering_params);
76  } else {
77  // Semigaussian mixture - get parameters
78  auto scattering_params_sg =
79  getSemigauss(beta, particle.absMomentum(), tInX0,
80  slab.material().Z(), genMixtureScalor);
81  // Simulate
82  theta = semigauss(generator, scattering_params_sg);
83  }
84  } else {
85  // for electrons we fall back to the Highland (extension)
86  // return projection factor times sigma times gauss random
87  const auto theta0 = Acts::computeMultipleScatteringTheta0(
88  slab, particle.pdg(), particle.mass(),
89  particle.charge() / particle.absMomentum(), particle.charge());
90  theta = std::normal_distribution<double>(0.0, theta0)(generator);
91  }
92  // scale from planar to 3d angle
93  return M_SQRT2 * theta;
94  }
95 
96  // helper methods for getting parameters and simulating
97 
98  std::array<double, 4> getGaussian(double beta, double p, double tInX0,
99  double scale) const {
100  std::array<double, 4> scattering_params;
101  // Total standard deviation of mixture
102  scattering_params[0] = 15. / beta / p * std::sqrt(tInX0) * scale;
103  scattering_params[1] = 1.0; // Variance of core
104  scattering_params[2] = 1.0; // Variance of tails
105  scattering_params[3] = 0.5; // Mixture weight of tail component
106  return scattering_params;
107  }
108 
109  std::array<double, 4> getGaussmix(double beta, double p, double tInX0,
110  double Z, double scale) const {
111  std::array<double, 4> scattering_params;
112  scattering_params[0] = 15. / beta / p * std::sqrt(tInX0) *
113  scale; // Total standard deviation of mixture
114  double d1 = std::log(tInX0 / (beta * beta));
115  double d2 = std::log(std::pow(Z, 2.0 / 3.0) * tInX0 / (beta * beta));
116  double epsi;
117  double var1 = (-1.843e-3 * d1 + 3.347e-2) * d1 + 8.471e-1; // Variance of
118  // core
119  if (d2 < 0.5)
120  epsi = (6.096e-4 * d2 + 6.348e-3) * d2 + 4.841e-2;
121  else
122  epsi = (-5.729e-3 * d2 + 1.106e-1) * d2 - 1.908e-2;
123  scattering_params[1] = var1; // Variance of core
124  scattering_params[2] = (1 - (1 - epsi) * var1) / epsi; // Variance of tails
125  scattering_params[3] = epsi; // Mixture weight of tail component
126  return scattering_params;
127  }
128 
129  std::array<double, 6> getSemigauss(double beta, double p, double tInX0,
130  double Z, double scale) const {
131  std::array<double, 6> scattering_params;
132  double N = tInX0 * 1.587E7 * std::pow(Z, 1.0 / 3.0) / (beta * beta) /
133  (Z + 1) / std::log(287 / std::sqrt(Z));
134  scattering_params[4] = 15. / beta / p * std::sqrt(tInX0) *
135  scale; // Total standard deviation of mixture
136  double rho = 41000 / std::pow(Z, 2.0 / 3.0);
137  double b = rho / std::sqrt(N * (std::log(rho) - 0.5));
138  double n = std::pow(Z, 0.1) * std::log(N);
139  double var1 = (5.783E-4 * n + 3.803E-2) * n + 1.827E-1;
140  double a =
141  (((-4.590E-5 * n + 1.330E-3) * n - 1.355E-2) * n + 9.828E-2) * n +
142  2.822E-1;
143  double epsi = (1 - var1) / (a * a * (std::log(b / a) - 0.5) - var1);
144  scattering_params[3] =
145  (epsi > 0) ? epsi : 0.0; // Mixture weight of tail component
146  scattering_params[0] = a; // Parameter 1 of tails
147  scattering_params[1] = b; // Parameter 2 of tails
148  scattering_params[2] = var1; // Variance of core
149  scattering_params[5] = N; // Average number of scattering processes
150  return scattering_params;
151  }
152 
161  template <typename generator_t>
162  double gaussmix(generator_t &generator,
163  const std::array<double, 4> &scattering_params) const {
164  std::uniform_real_distribution<double> udist(0.0, 1.0);
165  double sigma_tot = scattering_params[0];
166  double var1 = scattering_params[1];
167  double var2 = scattering_params[2];
168  double epsi = scattering_params[3];
169  bool ind = udist(generator) > epsi;
170  double u = udist(generator);
171  if (ind)
172  return std::sqrt(var1) * std::sqrt(-2 * std::log(u)) * sigma_tot;
173  else
174  return std::sqrt(var2) * std::sqrt(-2 * std::log(u)) * sigma_tot;
175  }
176 
185  template <typename generator_t>
186  double semigauss(generator_t &generator,
187  const std::array<double, 6> &scattering_params) const {
188  std::uniform_real_distribution<double> udist(0.0, 1.0);
189  double a = scattering_params[0];
190  double b = scattering_params[1];
191  double var1 = scattering_params[2];
192  double epsi = scattering_params[3];
193  double sigma_tot = scattering_params[4];
194  bool ind = udist(generator) > epsi;
195  double u = udist(generator);
196  if (ind)
197  return std::sqrt(var1) * std::sqrt(-2 * std::log(u)) * sigma_tot;
198  else
199  return a * b * std::sqrt((1 - u) / (u * b * b + a * a)) * sigma_tot;
200  }
201 };
202 
203 } // namespace detail
204 
205 using GeneralMixtureScattering = detail::Scattering<detail::GeneralMixture>;
206 
207 } // namespace ActsFatras