16 #include "GaudiKernel/ThreadLocalContext.h"
17 #include "CLHEP/Random/RandFlat.h"
23 constexpr
double s_log_RutherfordScott = 0.038;
26 constexpr
double s_log_RossiGreisen = 0.125;
29 constexpr
double s_gausMixSigma1_a0 = 8.471e-1;
30 constexpr
double s_gausMixSigma1_a1 = 3.347e-2;
31 constexpr
double s_gausMixSigma1_a2 = -1.843e-3;
33 constexpr
double s_gausMixEpsilon_a0 = 4.841e-2;
34 constexpr
double s_gausMixEpsilon_a1 = 6.348e-3;
35 constexpr
double s_gausMixEpsilon_a2 = 6.096e-4;
37 constexpr
double s_gausMixEpsilon_b0 = -1.908e-2;
38 constexpr
double s_gausMixEpsilon_b1 = 1.106e-1;
39 constexpr
double s_gausMixEpsilon_b2 = -5.729e-3;
45 const IInterface *
p) :
47 declareInterface<IMultipleScatteringUpdator>(
this);
58 if (m_gaussianMixture) {
60 if (m_rndGenSvc.retrieve().isFailure()) {
61 ATH_MSG_WARNING(
"Could not retrieve " << m_rndGenSvc <<
" -> switching Gaussian mixture model off.");
62 m_gaussianMixture =
false;
68 m_rngWrapper = m_gaussianMixture ? m_rndGenSvc->getEngine(
this, m_randomEngineName) :
nullptr;
71 "Could not get random engine '" << m_randomEngineName <<
"' -> switching Gaussian mixture model off.");
72 m_gaussianMixture =
false;
75 if (m_gaussianMixture) {
76 ATH_MSG_VERBOSE(
"Gaussian mixture model = ON ... optimised = " << m_optGaussianMixtureG4);
81 return StatusCode::SUCCESS;
87 double pathcorrection,
95 pathcorrection *= pathcorrection < 0. ? (-1.) : (1);
98 double t = pathcorrection *
mat.thicknessInX0();
102 double E = sqrt(
p *
p +
m *
m);
119 sigma2 = s_main_RutherfordScott / (
beta *
p);
122 sigma2 *= (1. + s_log_RutherfordScott *
log(
t));
125 sigma2 *= (sigma2 *
t);
130 sigma2 = s_main_RossiGreisen / (
beta *
p);
131 sigma2 *= (sigma2 *
t);
134 double factor = 1. + s_log_RossiGreisen * log10(10. *
t);
141 if (m_gaussianMixture) {
143 const EventContext& ctx = Gaudi::Hive::currentContext();
144 if (m_rngWrapper->evtSeeded(ctx) != ctx.evt()) {
148 wrapper_nc->setSeed (this->
name(), ctx);
150 CLHEP::HepRandomEngine* engine = m_rngWrapper->getEngine (ctx);
153 double log_dprime =
log(dprime);
155 double log_dprimeprime =
log(
std::pow(
mat.averageZ(), 2.0 / 3.0) * dprime);
157 double epsilon = log_dprimeprime < 0.5 ?
158 s_gausMixEpsilon_a0 + s_gausMixEpsilon_a1 * log_dprimeprime + s_gausMixEpsilon_a2 *
159 log_dprimeprime * log_dprimeprime :
160 s_gausMixEpsilon_b0 + s_gausMixEpsilon_b1 * log_dprimeprime + s_gausMixEpsilon_b2 *
161 log_dprimeprime * log_dprimeprime;
163 double sigma1square = s_gausMixSigma1_a0 + s_gausMixSigma1_a1 * log_dprime + s_gausMixSigma1_a2 * log_dprime *
166 if (!m_optGaussianMixtureG4) {
167 sigma2 = 225. * dprime / (
p *
p);
170 if (CLHEP::RandFlat::shoot(engine) < epsilon) {
171 sigma2 *= (1. - (1. - epsilon) * sigma1square) / epsilon;