ATLAS Offline Software
MultipleScatteringUpdator.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
3 */
4 
6 // MultipleScatteringUpdator.cxx, (c) ATLAS Detector software
8 
9 // Trk include
12 // CLHEP
13 
16 #include "GaudiKernel/ThreadLocalContext.h"
17 #include "CLHEP/Random/RandFlat.h"
19 
20 namespace{
21 // static doubles
22 constexpr double s_main_RutherfordScott = 13.6 * Gaudi::Units::MeV;
23 constexpr double s_log_RutherfordScott = 0.038;
24 
25 constexpr double s_main_RossiGreisen = 17.5 * Gaudi::Units::MeV;
26 constexpr double s_log_RossiGreisen = 0.125;
27 
28 // ============================= Gaussian mixture model =============
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;
32 
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;
36 
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;
40 }
41 
42 
43 // constructor
44 Trk::MultipleScatteringUpdator::MultipleScatteringUpdator(const std::string &t, const std::string &n,
45  const IInterface *p) :
46  AthAlgTool(t, n, p) {
47  declareInterface<IMultipleScatteringUpdator>(this);
48 }
49 
50 // destructor
52 
53 // Athena standard methods
54 // initialize
57  // if Gaussian mixture is required
58  if (m_gaussianMixture) {
59  // get the random generator service
60  if (m_rndGenSvc.retrieve().isFailure()) {
61  ATH_MSG_WARNING("Could not retrieve " << m_rndGenSvc << " -> switching Gaussian mixture model off.");
62  m_gaussianMixture = false;
63  } else {
64  ATH_MSG_VERBOSE("Successfully retrieved " << m_rndGenSvc);
65  }
66 
67  // Get own engine with own seeds:
68  m_rngWrapper = m_gaussianMixture ? m_rndGenSvc->getEngine(this, m_randomEngineName) : nullptr;
69  if (!m_rngWrapper) {
71  "Could not get random engine '" << m_randomEngineName << "' -> switching Gaussian mixture model off.");
72  m_gaussianMixture = false;
73  }
74  }
75  if (m_gaussianMixture) {
76  ATH_MSG_VERBOSE("Gaussian mixture model = ON ... optimised = " << m_optGaussianMixtureG4);
77  } else {
78  ATH_MSG_VERBOSE("Gaussian mixture model = OFF");
79  }
80 
81  return StatusCode::SUCCESS;
82 }
83 
84 double
86  double p,
87  double pathcorrection,
89  double) const {
90  if (mat.thicknessInX0() <= 0. || particle == Trk::geantino) {
91  return 0.;
92  }
93 
94  // make sure the path correction is positive to avoid a floating point exception
95  pathcorrection *= pathcorrection < 0. ? (-1.) : (1);
96 
97  // scale the path length to the radiation length
98  double t = pathcorrection * mat.thicknessInX0();
99 
100  // kinematics (relativistic)
102  double E = sqrt(p * p + m * m);
103  double beta = p / E;
104 
105  double sigma2(0.);
106 
108  sigma2 = sigma * sigma;
109 
110  if (m_useTrkUtils && particle != Trk::electron) {
111  return sigma2;
112  }
113 
114 // Code below will not be used if the parameterization of TrkUtils is used
115 
116 
117  if (particle != Trk::electron) {
118  // the highland formula
119  sigma2 = s_main_RutherfordScott / (beta * p);
120 
121  if (m_log_include) {
122  sigma2 *= (1. + s_log_RutherfordScott * log(t));
123  }
124 
125  sigma2 *= (sigma2 * t);
126  }else {
127  // Electron multiple scattering effects - see Highland NIM 129 (1975) p497-499
128  // (Highland extension to the Rossi-Greisen formulation)
129  // NOTE: The formula can be extended by replacing the momentum^2 term with pi * pf
130  sigma2 = s_main_RossiGreisen / (beta * p);
131  sigma2 *= (sigma2 * t);
132 
133  if (m_log_include) {
134  double factor = 1. + s_log_RossiGreisen * log10(10. * t);
135  factor *= factor;
136  sigma2 *= factor;
137  }
138  }
139 
140  // use the Gaussian mixture model
141  if (m_gaussianMixture) {
142  // Reseed the RNG if needed.
143  const EventContext& ctx = Gaudi::Hive::currentContext();
144  if (m_rngWrapper->evtSeeded(ctx) != ctx.evt()) {
145  // Ok, the wrappers are unique to this algorithm and a given slot,
146  // so cannot be accessed concurrently.
147  ATHRNG::RNGWrapper* wrapper_nc ATLAS_THREAD_SAFE = m_rngWrapper;
148  wrapper_nc->setSeed (this->name(), ctx);
149  }
150  CLHEP::HepRandomEngine* engine = m_rngWrapper->getEngine (ctx);
151  // d_0'
152  double dprime = t / (beta * beta);
153  double log_dprime = log(dprime);
154  // d_0''
155  double log_dprimeprime = log(std::pow(mat.averageZ(), 2.0 / 3.0) * dprime);
156  // get epsilon
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;
162  // the standard sigma
163  double sigma1square = s_gausMixSigma1_a0 + s_gausMixSigma1_a1 * log_dprime + s_gausMixSigma1_a2 * log_dprime *
164  log_dprime;
165  // G4 optimised / native double Gaussian model
166  if (!m_optGaussianMixtureG4) {
167  sigma2 = 225. * dprime / (p * p);
168  }
169  // throw the random number core/tail
170  if (CLHEP::RandFlat::shoot(engine) < epsilon) {
171  sigma2 *= (1. - (1. - epsilon) * sigma1square) / epsilon;
172  }
173  }
174 
175  return sigma2;
176 }
Trk::ParticleSwitcher::particle
constexpr ParticleHypothesis particle[PARTICLEHYPOTHESES]
the array of masses
Definition: ParticleHypothesis.h:79
Trk::MaterialInteraction::sigmaMS
static double sigmaMS(double dInX0, double p, double beta)
multiple scattering as function of dInX0
Definition: MaterialInteraction.cxx:278
pdg_comparison.sigma
sigma
Definition: pdg_comparison.py:324
Trk::MultipleScatteringUpdator::MultipleScatteringUpdator
MultipleScatteringUpdator(const std::string &, const std::string &, const IInterface *)
AlgTool like constructor.
Definition: MultipleScatteringUpdator.cxx:44
Trk::MultipleScatteringUpdator::initialize
virtual StatusCode initialize() override
AlgTool initailize method.
Definition: MultipleScatteringUpdator.cxx:56
MaterialProperties.h
MultipleScatteringUpdator.h
mat
GeoMaterial * mat
Definition: LArDetectorConstructionTBEC.cxx:55
read_hist_ntuple.t
t
Definition: read_hist_ntuple.py:5
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition: AthMsgStreamMacros.h:28
python.SystemOfUnits.MeV
float MeV
Definition: SystemOfUnits.py:172
Trk::ParticleHypothesis
ParticleHypothesis
Definition: ParticleHypothesis.h:28
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:209
beamspotman.n
n
Definition: beamspotman.py:729
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
Trk::electron
@ electron
Definition: ParticleHypothesis.h:30
Trk::geantino
@ geantino
Definition: ParticleHypothesis.h:29
Trk::MultipleScatteringUpdator::~MultipleScatteringUpdator
virtual ~MultipleScatteringUpdator()
Virtual destructor.
Trk::ParticleMasses::mass
constexpr double mass[PARTICLEHYPOTHESES]
the array of masses
Definition: ParticleHypothesis.h:56
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:240
ATHRNG::RNGWrapper
A wrapper class for event-slot-local random engines.
Definition: RNGWrapper.h:56
VP1PartSpect::E
@ E
Definition: VP1PartSpectFlags.h:21
MaterialInteraction.h
RNGWrapper.h
Trk::MaterialProperties
Definition: MaterialProperties.h:40
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
python.CaloCondTools.log
log
Definition: CaloCondTools.py:20
ATLAS_THREAD_SAFE
#define ATLAS_THREAD_SAFE
Definition: checker_macros.h:211
AthAlgTool
Definition: AthAlgTool.h:26
MuonParameters::beta
@ beta
Definition: MuonParamDefs.h:144
checker_macros.h
Define macros for attributes used to control the static checker.
pow
constexpr int pow(int base, int exp) noexcept
Definition: ap_fixedTest.cxx:15
python.SystemOfUnits.m
float m
Definition: SystemOfUnits.py:106
Trk::MultipleScatteringUpdator::sigmaSquare
virtual double sigmaSquare(const MaterialProperties &mat, double p, double pathcorrection, ParticleHypothesis particle=pion, double deltaE=0.) const override
Calculate the sigma on theta introduced by multiple scattering, according to the RutherFord-Scott For...
Definition: MultipleScatteringUpdator.cxx:85