ATLAS Offline Software
MaterialEffectsUpdator.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 // MaterialEffectsUpdator.cxx, (c) ATLAS Detector software
8 
9 // Trk include
11 #include "GaudiKernel/ITHistSvc.h"
15 #include "TrkGeometry/Layer.h"
22 #include "TrkSurfaces/Surface.h"
23 // Amg
26 
28 // std
29 #include <algorithm>
30 #include <functional>
31 #include <cmath>
32 
33 #define TRKEXTOOLS_MAXUPDATES 100
34 #ifndef COVARIANCEUPDATEWITHCHECK
35 #define COVARIANCEUPDATEWITHCHECK(cov, sign, value) \
36  cov += (sign > 0 ? value : (value > cov ? 0 : sign * value))
37 #endif
38 
39 // constructor
40 Trk::MaterialEffectsUpdator::MaterialEffectsUpdator(const std::string& t, const std::string& n, const IInterface* p)
41  : AthAlgTool(t, n, p)
42 {
43  declareInterface<IMaterialEffectsUpdator>(this);
44 }
45 
46 // destructor
48 
49 // Athena standard methods
50 // initialize
53 {
54 
55  ATH_MSG_DEBUG("Minimal momentum cut for material update : " << m_momentumCut << " MeV");
56 
57  // retrieve the EnergyLoss Updator and Material Effects updator
58  if (m_doEloss) {
59  if (m_eLossUpdator.retrieve().isFailure()) {
60  ATH_MSG_FATAL("Failed to retrieve tool " << m_eLossUpdator
61  << ". No multiple scattering effects will be taken into account.");
62  m_doEloss = false;
63  return StatusCode::FAILURE;
64  }
65  ATH_MSG_DEBUG("Retrieved tool " << m_eLossUpdator);
66 
67  } else {
68  m_eLossUpdator.disable();
69  }
70 
71  if (m_doMs) {
72  if (m_msUpdator.retrieve().isFailure()) {
73  ATH_MSG_FATAL("Failed to retrieve tool " << m_msUpdator
74  << ". No energy loss effects will be taken into account.");
75  m_doMs = false;
76  return StatusCode::FAILURE;
77  }
78  ATH_MSG_DEBUG("Retrieved tool " << m_msUpdator);
79 
80  } else {
81  m_msUpdator.disable();
82  }
83 
84  // retrieve the material mapper for the validation mode
85  if (m_validationMode) {
86  if (m_materialMapper.retrieve().isFailure()) {
87  ATH_MSG_FATAL("Failed to retrieve tool " << m_materialMapper << ". No material recording.");
88  return StatusCode::FAILURE;
89  }
90  ATH_MSG_DEBUG("Retrieved tool " << m_materialMapper);
91 
92  } else {
93  m_materialMapper.disable();
94  }
95 
96  return StatusCode::SUCCESS;
97 }
98 
99 
100 std::unique_ptr<Trk::TrackParameters>
102  ICache& cache,
103  const TrackParameters* parm,
104  const Layer& lay,
107  MaterialUpdateMode matupmode
108 ) const {
109  // no material properties - pass them back
110  if (particle == Trk::geantino || particle == Trk::nonInteractingMuon || (!m_doMs && !m_doEloss) ||
111  !lay.isOnLayer(parm->position())) {
112  return parm->uniqueClone();
113  }
114 
115  // get the quantities
116  const Trk::MaterialProperties* mprop = lay.fullUpdateMaterialProperties(*parm);
117  if (!mprop) {
118  return parm->uniqueClone();
119  }
120 
121  // get the real pathlength
122  double pathCorrection = std::abs(lay.surfaceRepresentation().pathCorrection(parm->position(), parm->momentum()));
123 
124  // --------------------------------------------------------------------------------------------------
125  if (m_validationMode) {
126  cache.validationLayer = &lay;
127  }
128  return updateImpl(cache, parm, *mprop, pathCorrection, dir, particle, matupmode);
129 }
130 
131 std::unique_ptr<Trk::TrackParameters>
133  ICache& cache,
134  const TrackParameters* parm,
135  const MaterialEffectsOnTrack& meff,
137  MaterialUpdateMode matupmode
138 ) const {
139  // no material properties - pass them back
140  // TODO, if the parm doesn't have a surface (i.e. its in
141  // curvilinear) then should we fall through?
142  if (particle == Trk::geantino || particle == Trk::nonInteractingMuon || (!m_doMs && !m_doEloss) ||
143  parm->associatedSurface() != meff.associatedSurface()) {
144  return parm->uniqueClone();
145  }
146 
147  // get the kinematics
148  double p = parm->momentum().mag();
149  double updateMomentum = m_forceMomentum ? m_forcedMomentum.value() : p;
151  double E = std::sqrt(p * p + m * m);
152  double beta = p / E;
153 
154  double pathcorrection = 1.; // Trick the MultipleScatteringUpdator interface
155 
156  double energyLoss = 0;
157  double energyLossSigma = 0;
158 
159  if (meff.energyLoss() != nullptr) {
160  energyLoss = meff.energyLoss()->deltaE();
161  energyLossSigma = meff.energyLoss()->sigmaDeltaE();
162  }
163 
164  // update for mean energy loss
165  double newP = p;
166  double sigmaQoverP = 0;
167  double sigmaQoverPSq = 0;
168  // Landaus mpvs don't just add, if in Landau mode we need to do a different update
169  if (m_landauMode && cache.accumulatedElossSigma != 0 && energyLossSigma != 0) {
170  if (energyLoss > 0) {
171  energyLoss += energyLossSigma * std::log(1 + cache.accumulatedElossSigma / (energyLossSigma)) +
172  cache.accumulatedElossSigma * std::log(1 + energyLossSigma / cache.accumulatedElossSigma);
173  } else {
174  energyLoss -= energyLossSigma * std::log(1 + cache.accumulatedElossSigma / energyLossSigma) +
175  cache.accumulatedElossSigma * std::log(1 + energyLossSigma / cache.accumulatedElossSigma);
176  }
177  cache.accumulatedElossSigma += energyLossSigma;
178  } else if (m_landauMode) {
179  cache.accumulatedElossSigma += energyLossSigma;
180  }
181  double qOverPnew = parm->parameters()[Trk::qOverP];
182 
183  if (p > m_momentumCut && p < m_momentumMax && m_doEloss) {
184  double newP2 = (E + energyLoss) * (E + energyLoss) - m * m;
185  if (newP2 < m_momentumCut * m_momentumCut) {
186  return nullptr; // protect FPE
187  }
188  if (E + energyLoss < -m) {
189  return nullptr; // protect against flip in correction
190  }
191  newP = std::sqrt(newP2);
192  sigmaQoverP = energyLossSigma / (beta * p * p);
193  sigmaQoverPSq = sigmaQoverP * sigmaQoverP;
194  qOverPnew = parm->charge() / newP;
195  }
196  Trk::DefinedParameter qOverPmod(qOverPnew, Trk::qOverP);
197 
198  // check if Parameters are measured parameters
199  // the updatedParameters - first a copy
200  const Trk::TrackParameters* mpars = parm;
201  AmgVector(5) updatedParameters(mpars->parameters());
202  std::optional<AmgSymMatrix(5)> updatedCovariance = std::nullopt;
203  // initialize ErrorMatrix pointer
204  if (m_validationMode && !m_validationIgnoreUnmeasured) {
205  // the new CovarianceMatrix - a copy first
206  updatedCovariance = mpars->covariance()
207  ? std::optional<AmgSymMatrix(5)>(*mpars->covariance())
208  : std::nullopt;
209  double angularVariation = 0;
210  double sigmaDeltaPhiSq = 0;
211  double sigmaDeltaThetaSq = 0;
212 
213  if (m_doMs) {
214  // If the meff has scattering angles use those, otherwise use MEffUpdator
215  if (meff.scatteringAngles() == nullptr) {
216  // Trick to keep using existing MultipleScatteringUpdator interface
217  // Here we know the path length to be meff.thicknessX0, so we set pathcorrection = 1
218  // and create a dummy materialProperties with the properties we are interested in
219  MaterialProperties mprop(meff.thicknessInX0(), 1., 1., 0., 0., 0.);
220  angularVariation = m_msUpdator->sigmaSquare(mprop, updateMomentum, pathcorrection, Trk::muon);
221  // sigmaDeltaPhiSq = angularVariation/(parm->sinTheta()*parm->sinTheta());
222  sigmaDeltaPhiSq =
223  angularVariation / (std::sin(parm->parameters()[Trk::theta]) * std::sin(parm->parameters()[Trk::theta]));
224  sigmaDeltaThetaSq = angularVariation;
225  } else {
226  // material update from mefots -> D.L.
227  sigmaDeltaPhiSq = meff.scatteringAngles()->sigmaDeltaPhi();
228  sigmaDeltaPhiSq *= sigmaDeltaPhiSq;
229  sigmaDeltaThetaSq = meff.scatteringAngles()->sigmaDeltaTheta();
230  sigmaDeltaThetaSq *= sigmaDeltaThetaSq;
231  updatedParameters[Trk::phi] = parm->position().phi() + meff.scatteringAngles()->deltaPhi();
232  updatedParameters[Trk::theta] = parm->position().theta() + meff.scatteringAngles()->deltaTheta();
233  }
234  }
235  // update the covariance entries - angular variation in phi has dependency on theta direction
236  if (updatedCovariance) {
237  // sign of the noise adding -----------------------
238  int sign = int(matupmode);
239  // check for non-zero covariance matrix
240  COVARIANCEUPDATEWITHCHECK((*updatedCovariance)(Trk::phi, Trk::phi), sign, sigmaDeltaPhiSq);
241  COVARIANCEUPDATEWITHCHECK((*updatedCovariance)(Trk::theta, Trk::theta), sign, sigmaDeltaThetaSq);
242  if (!m_landauMode) {
243  COVARIANCEUPDATEWITHCHECK((*updatedCovariance)(Trk::qOverP, Trk::qOverP), sign, sigmaQoverPSq);
244  } else if (m_landauMode) {
245  // subtract what we added up till now and add what we should add up till now
246  // Landau's 68% limit is approx 1.6*sigmaParameter
247  (*updatedCovariance)(Trk::qOverP, Trk::qOverP) -=
248  sign * std::pow(1.6 * (cache.accumulatedElossSigma - p * p * sigmaQoverP) / (p * p), 2);
249  (*updatedCovariance)(Trk::qOverP, Trk::qOverP) +=
250  sign * std::pow(1.6 * cache.accumulatedElossSigma / (newP * newP), 2);
251  }
252  // the checks for the remove Noise mode -----------------------------------------------------
253  if (matupmode == Trk::removeNoise && !checkCovariance(*updatedCovariance)) {
254  // the covariance is invalid
255  return nullptr;
256  }
257  // -------------------------------------- screen output --------------------------------------
258  if (m_msgOutputCorrections) {
259  double sigmaAngle = std::sqrt(angularVariation);
260  ATH_MSG_VERBOSE(" sigma(phi) / sigma(theta) = " << sigmaAngle / std::sin(parm->parameters()[Trk::theta]) << " / "
261  << sigmaAngle);
262  ATH_MSG_VERBOSE(" deltaP / sigmaQoverP = " << energyLoss << " / " << sigmaQoverP);
263  }
264  // -------------------------------------------------------------------------------------------
265  }
266  // ----------------------------------------- validation section ----------------------------------
267  // validation if configured
268  if (m_validationMode) {
269  if (cache.validationLayer) {
270  // all you have from MaterialProperties
271  double pathInX0 = meff.thicknessInX0();
272 
273  Trk::AssociatedMaterial assMatHit(parm->position(),
274  pathInX0,
275  pathInX0,
276  0,
277  0,
278  0,
279  0,
280  pathcorrection,
282  cache.validationLayer);
283 
284  // record the Material hit ----------------------------------------------------------------
285  m_materialMapper->recordMaterialHit(assMatHit, parm->position());
286 
287  // the steps
288  cache.validationSteps++;
289  cache.validationPhi += parm->position().phi();
290  cache.validationEta += parm->position().eta();
291  // reset the validation layer
292  cache.validationLayer = nullptr;
293  }
294  }
295  // ----------------------------------------- validation section ----------------------------------
296  }
297  updatedParameters[Trk::qOverP] = qOverPnew;
299  updatedParameters[Trk::loc1],
300  updatedParameters[Trk::loc2],
301  updatedParameters[Trk::phi],
302  updatedParameters[Trk::theta],
303  updatedParameters[Trk::qOverP],
304  updatedCovariance
305  );
306 }
307 
308 std::unique_ptr<Trk::TrackParameters>
310  ICache& cache,
311  const TrackParameters* parm,
312  const Layer& lay,
315  MaterialUpdateMode matupmode
316 ) const {
317  // no material properties - pass the parameters back
318  if (particle == Trk::geantino || particle == Trk::nonInteractingMuon || (!m_doMs && !m_doEloss)) {
319  return parm->uniqueClone();
320  }
321 
322  // get the split factor
323  double preFactor = lay.preUpdateMaterialFactor(*parm, dir);
324  // return if the preFactor is less than one
325  if (preFactor < 0.01) {
326  return parm->uniqueClone();
327  }
328 
329  // get the material properties
330  const Trk::MaterialProperties* mprop = nullptr;
331 
332  // set the output if restricted to the validation direction
333  bool outputFlag = m_msgOutputValidationDirection ? dir == int(m_validationDirection) : true;
334 
335  mprop = lay.fullUpdateMaterialProperties(*parm);
336  double pathCorrection = std::abs(lay.surfaceRepresentation().pathCorrection(parm->position(), parm->momentum()));
337  pathCorrection *= preFactor;
338 
339  // exit if no mprop could be assigned
340  if (!mprop) {
341  return parm->uniqueClone();
342  }
343  // --------------------------------------------------------------------------------------------------
344  if (outputFlag) {
345  double layerR = lay.surfaceRepresentation().bounds().r();
346  double layerZ = lay.surfaceRepresentation().center().z();
347  double eta = parm->momentum().eta();
348  double sX0 = mprop->thicknessInX0();
349  double tX0 = pathCorrection * mprop->thicknessInX0();
350  ATH_MSG_VERBOSE(" [M] pre material update at layer with [r,z] = [ " << layerR << ", " << layerZ << " ] - Index "
351  << lay.layerIndex());
352  ATH_MSG_VERBOSE(" thickness/X0 , path/X0 (eta: g.factor) = " << sX0 << " , " << tX0 << " (" << eta << ": "
353  << pathCorrection << ")");
354  }
355  // --------------------------------------------------------------------------------------------------
356  if (m_validationMode) {
357  cache.validationLayer = &lay;
358  }
359  return updateImpl(cache, parm, *mprop, pathCorrection, dir, particle, matupmode);
360 }
361 
362 std::unique_ptr<Trk::TrackParameters>
364  ICache& cache,
365  const TrackParameters& parm,
366  const Layer& lay,
369  MaterialUpdateMode matupmode
370 ) const {
371  // no material properties - pass the parameters back
372  if (particle == Trk::geantino || particle == Trk::nonInteractingMuon || (!m_doMs && !m_doEloss) ||
373  !lay.isOnLayer(parm.position())) {
374  return parm.uniqueClone();
375  }
376 
377  // get the quantities
378  const Trk::MaterialProperties* mprop = nullptr;
379  double postFactor = lay.postUpdateMaterialFactor(parm, dir);
380 
381  // no material properties - pass them back
382  if (postFactor < 0.01) {
383  return parm.uniqueClone();
384  }
385 
386  // set the output if restricted to the validation direction
387  bool outputFlag = m_msgOutputValidationDirection ? dir == int(m_validationDirection) : true;
388 
389  mprop = lay.fullUpdateMaterialProperties(parm);
390  double pathCorrection = std::abs(lay.surfaceRepresentation().pathCorrection(parm.position(), parm.momentum()));
391  pathCorrection *= postFactor;
392 
393  // exit if no material properties
394  if (!mprop) {
395  return parm.uniqueClone();
396  }
397 
398  // --------------------------------------------------------------------------------------------------
399  if (outputFlag) {
400  double layerR = lay.surfaceRepresentation().bounds().r();
401  double layerZ = lay.surfaceRepresentation().center().z();
402  double eta = parm.momentum().eta();
403  double sX0 = mprop->thicknessInX0();
404  double tX0 = pathCorrection * mprop->thicknessInX0();
405  ATH_MSG_VERBOSE(" [M] post material update, layer with [r,z] = [ " << layerR << ", " << layerZ << " ] - Index "
406  << lay.layerIndex());
407  ATH_MSG_VERBOSE(" thickness/X0 , path/X0 (eta: g.factor) = " << sX0 << " , " << tX0 << " (" << eta << ": "
408  << pathCorrection << ")");
409  }
410  // --------------------------------------------------------------------------------------------------
411  if (m_validationMode) {
412  cache.validationLayer = &lay;
413  }
414  return updateImpl(cache, parm, *mprop, pathCorrection, dir, particle, matupmode);
415 }
416 
417 // actual update method - manipulation
418 std::unique_ptr<Trk::TrackParameters>
420  ICache& cache,
421  const TrackParameters* parm,
422  const MaterialProperties& matprop,
423  double pathcorrection,
426  MaterialUpdateMode matupmode
427 ) const {
428  // no material properties - pass them back
429  if (particle == Trk::geantino || particle == Trk::nonInteractingMuon || (!m_doMs && !m_doEloss)) {
430  return parm->uniqueClone();
431  }
432 
433  // get the kinematics
434  double p = parm->momentum().mag();
435  double updateMomentum = m_forceMomentum ? m_forcedMomentum.value() : p;
437  double E = std::sqrt(p * p + m * m);
438  double beta = p / E;
439 
440  // set the output if restricted to the validation direction
441  bool outputFlag = m_msgOutputValidationDirection ? dir == int(m_validationDirection) : true;
442 
443  // no material update below/above a certain cut value
444  if (p > m_momentumCut && p < m_momentumMax) {
445  // get the delta of the Energy
446  EnergyLoss energyLoss =
447  (m_doEloss) ? m_eLossUpdator->energyLoss(matprop, updateMomentum,
448  pathcorrection, dir, particle,
449  m_useMostProbableEloss)
450  : EnergyLoss{};
451  // update for mean energy loss
452  double deltaE = energyLoss.deltaE();
453  double sigmaDeltaE = energyLoss.sigmaDeltaE();
454 
455  if (m_landauMode && cache.accumulatedElossSigma != 0 && sigmaDeltaE != 0) {
456  if (dir == Trk::oppositeMomentum) {
457  deltaE += sigmaDeltaE * std::log(1 + cache.accumulatedElossSigma / sigmaDeltaE) +
458  cache.accumulatedElossSigma * std::log(1 + sigmaDeltaE / cache.accumulatedElossSigma);
459  } else {
460  deltaE -= sigmaDeltaE * std::log(1 + cache.accumulatedElossSigma / sigmaDeltaE) +
461  cache.accumulatedElossSigma * std::log(1 + sigmaDeltaE / cache.accumulatedElossSigma);
462  }
463  cache.accumulatedElossSigma += sigmaDeltaE;
464  } else if (m_landauMode) {
465  cache.accumulatedElossSigma += sigmaDeltaE;
466  }
467 
468  double newP2 = (E + deltaE) * (E + deltaE) - m * m;
469  if (E + deltaE < -m) {
470  return nullptr; // protect against flip in correction
471  }
472  if (newP2 < m_momentumCut * m_momentumCut) {
473  return nullptr; // protect against FPE
474  }
475  double deltaP = std::sqrt(newP2) - p;
476  double sigmaQoverP = sigmaDeltaE / std::pow(beta * p, 2);
477 
478  Trk::DefinedParameter qOverPmod(parm->charge() / (p + deltaP), Trk::qOverP);
479 
480  // check if Parameters are measured parameters
481  const Trk::TrackParameters* mpars = parm;
482  AmgVector(5) updatedParameters(mpars->parameters());
483  // initialize ErrorMatrix pointer
484  std::optional<AmgSymMatrix(5)> updatedCovariance = std::nullopt;
485  {
486  // the new CovarianceMatrix - a copy first
487  updatedCovariance =
488  mpars->covariance()
489  ? std::optional<AmgSymMatrix(5)>(*mpars->covariance())
490  : std::nullopt;
491  // only update if msUpdator exists
492  double angularVariation =
493  (m_doMs) ? m_msUpdator->sigmaSquare(matprop, updateMomentum, pathcorrection, particle) : 0.;
494  // update the covariance entries - angular variation in phi has dependency on theta direction
495  if (updatedCovariance) {
496  // sign of the noise adding ----------------------------------------------------------------
497  int sign = int(matupmode);
498 
499  double sigmaDeltaPhiSq =
500  angularVariation / (std::sin(parm->parameters()[Trk::theta]) * std::sin(parm->parameters()[Trk::theta]));
501  double sigmaDeltaThetaSq = angularVariation;
502  // checks will only be done in the removeNoise mode
503  COVARIANCEUPDATEWITHCHECK((*updatedCovariance)(Trk::phi, Trk::phi), sign, sigmaDeltaPhiSq);
504  COVARIANCEUPDATEWITHCHECK((*updatedCovariance)(Trk::theta, Trk::theta), sign, sigmaDeltaThetaSq);
505  if (!m_landauMode) {
506  COVARIANCEUPDATEWITHCHECK((*updatedCovariance)(Trk::qOverP, Trk::qOverP), sign, sigmaQoverP * sigmaQoverP);
507  } else if (m_landauMode) {
508  // subtract what we added up till now and add what we should add up till now
509  // Landau's 68% limit is approx 1.6*sigmaParameter
510  (*updatedCovariance)(Trk::qOverP, Trk::qOverP) -=
511  sign * std::pow(1.6 * (cache.accumulatedElossSigma - p * p * sigmaQoverP) / (p * p), 2);
512  (*updatedCovariance)(Trk::qOverP, Trk::qOverP) +=
513  sign * std::pow(1.6 * cache.accumulatedElossSigma / ((p + deltaP) * (p + deltaP)), 2);
514  }
515  // the checks for the remove Noise mode -----------------------------------------------------
516  if (matupmode == Trk::removeNoise && !checkCovariance(*updatedCovariance)) {
517  // the covariance is invalid
518  return nullptr;
519  }
520 
521  // -------------------------------------- screen output --------------------------------------
522  if (outputFlag && m_msgOutputCorrections) {
523  double sigmaAngle = std::sqrt(angularVariation);
524  ATH_MSG_VERBOSE(" sigma(phi) / sigma(theta) = " << sigmaAngle / std::sin(parm->parameters()[Trk::theta])
525  << " / " << sigmaAngle);
526  ATH_MSG_VERBOSE(" deltaP / sigmaQoverP = " << deltaP << " / " << sigmaQoverP);
527  }
528  // -------------------------------------------------------------------------------------------
529  }
530  // ----------------------------------------- validation section ----------------------------------
531  // validation if configured
532  if (m_validationMode && dir == Trk::PropDirection(m_validationDirection.value()) && updatedCovariance) {
533 
534  if (cache.validationLayer) {
535  // all you have from MaterialProperties
536  double pathInX0 = pathcorrection * matprop.thicknessInX0();
537  double A = 0.;
538  double Z = 0.;
539  double rho = 0.;
540  double l0 = 0.;
541  // or better take the extended version for more information
542  const Trk::MaterialProperties* extProperties = dynamic_cast<const Trk::MaterialProperties*>(&matprop);
543 
544  if (extProperties) {
545  A = extProperties->averageA();
546  Z = extProperties->averageZ();
547  rho = extProperties->averageRho();
548  l0 = extProperties->l0();
549  }
550 
551  Trk::AssociatedMaterial assMatHit(parm->position(),
552  pathInX0,
553  matprop.x0(),
554  l0,
555  A,
556  Z,
557  rho,
558  pathcorrection,
560  cache.validationLayer);
561 
562  // record the Material hit ----------------------------------------------------------------
563  m_materialMapper->recordMaterialHit(assMatHit, parm->position());
564  // the steps
565  cache.validationSteps++;
566  cache.validationPhi += parm->position().phi();
567  cache.validationEta += parm->position().eta();
568  // reset the cache.tion layer
569  cache.validationLayer = nullptr;
570  }
571  }
572  // ----------------------------------------- validation section ----------------------------------
573  }
574  updatedParameters[Trk::qOverP] = parm->charge() / (p + deltaP);
576  updatedParameters[Trk::loc1],
577  updatedParameters[Trk::loc2],
578  updatedParameters[Trk::phi],
579  updatedParameters[Trk::theta],
580  updatedParameters[Trk::qOverP],
581  std::move(updatedCovariance)
582  );
583  }
584  //default if we have not returned just above
585  return parm->uniqueClone();
586 }
587 
588 // actual update method
589 std::unique_ptr<Trk::TrackParameters>
591  ICache& cache,
592  const TrackParameters& parm,
593  const MaterialProperties& matprop,
594  double pathcorrection,
597  MaterialUpdateMode matupmode
598 ) const {
599  // no material properties - pass them back
600  if (particle == Trk::geantino || (!m_doMs && !m_doEloss)) {
601  return parm.uniqueClone();
602  }
603 
604  // get the kinematics
605  double p = parm.momentum().mag();
606  double updateMomentum = m_forceMomentum ? m_forcedMomentum.value() : p;
608  double E = std::sqrt(p * p + m * m);
609  double beta = p / E;
610 
611  // set the output if restricted to the validation direction
612  bool outputFlag = m_msgOutputValidationDirection ? dir == int(m_validationDirection) : true;
613 
614  // no material update below or above a certain cut value
615  if (p > m_momentumCut && p < m_momentumMax) {
616  // the updatedParameters - first a copy
617  AmgVector(5) updatedParameters(parm.parameters());
618 
619  // get the delta of the Energy
620  EnergyLoss energyLoss =
621  (m_doEloss)
622  ? m_eLossUpdator->energyLoss(matprop, updateMomentum, pathcorrection, dir, particle, m_useMostProbableEloss)
623  : EnergyLoss{};
624  // update for mean energy loss
625  double deltaE = energyLoss.deltaE() ;
626  double sigmaDeltaE = energyLoss.sigmaDeltaE();
627  if (m_landauMode && cache.accumulatedElossSigma != 0 && sigmaDeltaE != 0) {
628  if (dir == Trk::oppositeMomentum) {
629  deltaE += sigmaDeltaE * std::log(1 + cache.accumulatedElossSigma / sigmaDeltaE) +
630  cache.accumulatedElossSigma * std::log(1 + sigmaDeltaE / cache.accumulatedElossSigma);
631  } else {
632  deltaE -= sigmaDeltaE * std::log(1 + cache.accumulatedElossSigma / sigmaDeltaE) +
633  cache.accumulatedElossSigma * std::log(1 + sigmaDeltaE / cache.accumulatedElossSigma);
634  }
635 
636  cache.accumulatedElossSigma += sigmaDeltaE;
637  } else if (m_landauMode) {
638  cache.accumulatedElossSigma += sigmaDeltaE;
639  }
640 
641  double newP2 = (E + deltaE) * (E + deltaE) - m * m;
642  if (E + deltaE < -m) {
643  return nullptr; // protect against flip in correction
644  }
645  if (newP2 < m_momentumCut * m_momentumCut) {
646  return nullptr; // protect against FPE
647  }
648  double deltaP = std::sqrt(newP2) - p;
649  double sigmaQoverP = sigmaDeltaE / std::pow(beta * p, 2);
650 
651  updatedParameters[Trk::qOverP] = parm.charge() / (p + deltaP);
652 
653  // check if Parameters are measured parameters
654  std::optional<AmgSymMatrix(5)> updatedCovariance = std::nullopt;
655  if (parm.covariance() || (m_validationMode && !m_validationIgnoreUnmeasured)) {
656  // the new CovarianceMatrix - a copy first
657  updatedCovariance = AmgSymMatrix(5)(*parm.covariance());
658  // only update if msUpdator exists
659  double angularVariation =
660  (m_doMs) ? m_msUpdator->sigmaSquare(matprop, updateMomentum, pathcorrection, particle) : 0.;
661  // update the covariance entries - angular variation in phi has dependency on theta direction
662  // sign of the noise adding ----------------------------------------------------------------
663  int sign = int(matupmode);
664  // checks will only be done in the removeNoise mode
665  double sigmaDeltaPhiSq =
666  angularVariation / (std::sin(parm.parameters()[Trk::theta]) * std::sin(parm.parameters()[Trk::theta]));
667  double sigmaDeltaThetaSq = angularVariation;
668  // checks will only be done in the removeNoise mode
669  COVARIANCEUPDATEWITHCHECK((*updatedCovariance)(Trk::phi, Trk::phi), sign, sigmaDeltaPhiSq);
670  COVARIANCEUPDATEWITHCHECK((*updatedCovariance)(Trk::theta, Trk::theta), sign, sigmaDeltaThetaSq);
671  if (!m_landauMode) {
672  COVARIANCEUPDATEWITHCHECK((*updatedCovariance)(Trk::qOverP, Trk::qOverP), sign, sigmaQoverP * sigmaQoverP);
673  } else if (m_landauMode) {
674  // subtract what we added up till now and add what we should add up till now
675  // Landau's 68% limit is best modeled by 1.6*sigmaParameter
676  (*updatedCovariance)(Trk::qOverP, Trk::qOverP) -=
677  sign * std::pow(1.6 * (cache.accumulatedElossSigma - p * p * sigmaQoverP) / (p * p), 2);
678  (*updatedCovariance)(Trk::qOverP, Trk::qOverP) +=
679  sign * std::pow(1.6 * cache.accumulatedElossSigma / ((p + deltaP) * (p + deltaP)), 2);
680  }
681  // the checks for the remove Noise mode -----------------------------------------------------
682  if (matupmode == Trk::removeNoise && !checkCovariance(*updatedCovariance)) {
683  // the covariance is invalid
684  return nullptr;
685  }
686  // -------------------------------------- screen output
687  // --------------------------------------
688  if (outputFlag && m_msgOutputCorrections) {
689  double sigmaAngle = std::sqrt(angularVariation);
690  ATH_MSG_VERBOSE(" sigma(phi) / sigma(theta) = " << sigmaAngle / std::sin(parm.parameters()[Trk::theta]) << " / "
691  << sigmaAngle);
692  ATH_MSG_VERBOSE(" deltaP / sigmaQoverP = " << deltaP << " / " << sigmaQoverP);
693  }
694  // ----------------------------------------- validation section ----------------------------------
695  // validation if configured
696  if (m_validationMode && dir == Trk::PropDirection(m_validationDirection.value())) {
697 
698  if (cache.validationLayer) {
699  // all you have from MaterialProperties
700  double pathInX0 = pathcorrection * matprop.thicknessInX0();
701  double A = 0.;
702  double Z = 0.;
703  double rho = 0.;
704  double l0 = 0.;
705  // or better take the extended version for more information
706  const Trk::MaterialProperties* extProperties = dynamic_cast<const Trk::MaterialProperties*>(&matprop);
707 
708  if (extProperties) {
709  A = extProperties->averageA();
710  Z = extProperties->averageZ();
711  rho = extProperties->averageRho();
712  l0 = extProperties->l0();
713  }
714 
715  Trk::AssociatedMaterial assMatHit(parm.position(),
716  pathInX0,
717  matprop.x0(),
718  l0,
719  A,
720  Z,
721  rho,
722  pathcorrection,
724  cache.validationLayer);
725 
726  // record the Material hit ----------------------------------------------------------------
727  m_materialMapper->recordMaterialHit(assMatHit, parm.position());
728 
729  // the steps
730  cache.validationSteps++;
731  cache.validationPhi += parm.position().phi();
732  cache.validationEta += parm.position().eta();
733  // reset the validation layer
734  cache.validationLayer = nullptr;
735  }
736  }
737  // ------------------------------------------validation section ----------------------------------
738  }
740  updatedParameters[Trk::loc1],
741  updatedParameters[Trk::loc2],
742  updatedParameters[Trk::phi],
743  updatedParameters[Trk::theta],
744  updatedParameters[Trk::qOverP],
745  std::move(updatedCovariance)
746  );
747  }
748  //default if we have not returned before
749  return parm.uniqueClone();
750 }
751 
752 void
754 {
755  cache.validationEta = 0.;
756  cache.validationPhi = 0.;
757  cache.validationSteps = 0;
758 }
759 
760 void
762 {
763  cache.accumulatedElossSigma = 0;
764 }
765 
767 {
768  if (updated(Trk::phi, Trk::phi) > 0. && updated(Trk::theta, Trk::theta) > 0. &&
769  updated(Trk::qOverP, Trk::qOverP) > 0.) {
770  return true;
771  }
772 
773  // given an update
774  ATH_MSG_DEBUG(" [-] update material with 'removeNoise' results in negative covariance entries. Skipping update.");
775  return false;
776 }
Trk::MaterialEffectsUpdator::initialize
virtual StatusCode initialize() override
AlgTool initailize method.
Definition: MaterialEffectsUpdator.cxx:52
Trk::ScatteringAngles::deltaPhi
double deltaPhi() const
returns the
Definition: ScatteringAngles.h:82
Trk::Layer::isOnLayer
virtual bool isOnLayer(const Amg::Vector3D &gp, const BoundaryCheck &bcheck=BoundaryCheck(true)) const
isOnLayer() method, using isOnSurface() with Layer specific tolerance
Definition: Layer.cxx:135
EnergyLoss.h
Trk::Layer::postUpdateMaterialFactor
virtual double postUpdateMaterialFactor(const TrackParameters &, PropDirection) const
getting the MaterialProperties back - for pre-update
Definition: Layer.h:150
Trk::Surface::pathCorrection
virtual double pathCorrection(const Amg::Vector3D &pos, const Amg::Vector3D &mom) const
the pathCorrection for derived classes with thickness - it reflects if the direction projection is po...
ScatteringAngles.h
Trk::IMaterialEffectsUpdator::ICache
Cache class to allow passing information to/between calls.
Definition: IMaterialEffectsUpdator.h:63
ATH_MSG_FATAL
#define ATH_MSG_FATAL(x)
Definition: AthMsgStreamMacros.h:34
Trk::Layer::fullUpdateMaterialProperties
const MaterialProperties * fullUpdateMaterialProperties(const TrackParameters &par) const
getting the MaterialProperties back - for full update
Definition: Layer.cxx:169
Trk::ParticleSwitcher::particle
constexpr ParticleHypothesis particle[PARTICLEHYPOTHESES]
the array of masses
Definition: ParticleHypothesis.h:79
Trk::MaterialEffectsUpdator::~MaterialEffectsUpdator
virtual ~MaterialEffectsUpdator()
Virtual destructor.
Trk::MaterialEffectsUpdator::modelActionImpl
static void modelActionImpl(ICache &cache, const TrackParameters *parm=nullptr)
Definition: MaterialEffectsUpdator.cxx:761
Trk::MaterialProperties::averageA
float averageA() const
Return the average A of the material [gram/mole].
TrackParameters.h
Trk::ParametersBase::charge
double charge() const
Returns the charge.
Surface.h
MaterialProperties.h
Trk::ParametersBase::position
const Amg::Vector3D & position() const
Access method for the position.
Trk::oppositeMomentum
@ oppositeMomentum
Definition: PropDirection.h:21
Trk::ParametersBase::associatedSurface
virtual const Surface & associatedSurface() const override=0
Access to the Surface associated to the Parameters.
Trk::ParametersBase::uniqueClone
std::unique_ptr< ParametersBase< DIM, T > > uniqueClone() const
clone method for polymorphic deep copy returning unique_ptr; it is not overriden, but uses the existi...
Definition: ParametersBase.h:97
Trk::MaterialProperties::thicknessInX0
float thicknessInX0() const
Return the radiationlength fraction.
Trk::EnergyLoss::sigmaDeltaE
double sigmaDeltaE() const
returns the symmatric error
Trk::MaterialProperties::averageRho
float averageRho() const
Return the average density of the material.
Trk::loc2
@ loc2
generic first and second local coordinate
Definition: ParamDefs.h:35
Layer.h
Trk::MaterialEffectsUpdator::MaterialEffectsUpdator
MaterialEffectsUpdator(const std::string &, const std::string &, const IInterface *)
AlgTool like constructor.
Definition: MaterialEffectsUpdator.cxx:40
Trk::MaterialEffectsBase::thicknessInX0
double thicknessInX0() const
returns the actually traversed material .
Trk::MaterialUpdateMode
MaterialUpdateMode
This is a steering enum to force the material update it can be: (1) addNoise (-1) removeNoise Second ...
Definition: MaterialUpdateMode.h:18
Trk::MaterialProperties::x0
float x0() const
Return the radiation length.
read_hist_ntuple.t
t
Definition: read_hist_ntuple.py:5
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition: AthMsgStreamMacros.h:28
Trk::Surface::center
const Amg::Vector3D & center() const
Returns the center position of the Surface.
ParamDefs.h
COVARIANCEUPDATEWITHCHECK
#define COVARIANCEUPDATEWITHCHECK(cov, sign, value)
Definition: MaterialEffectsUpdator.cxx:35
Trk::IMaterialEffectsUpdator::ICache::validationSteps
int validationSteps
number of validation steps
Definition: IMaterialEffectsUpdator.h:75
MaterialEffectsUpdator.h
Trk::DefinedParameter
std::pair< double, ParamDefs > DefinedParameter
Definition: DefinedParameter.h:27
Trk::Layer::preUpdateMaterialFactor
virtual double preUpdateMaterialFactor(const TrackParameters &, PropDirection) const
getting the MaterialProperties back - for pre-update
Definition: Layer.h:144
MaterialEffectsOnTrack.h
Trk::IMaterialEffectsUpdator::ICache::accumulatedElossSigma
double accumulatedElossSigma
Definition: IMaterialEffectsUpdator.h:80
Trk::AmgSymMatrix
AmgSymMatrix(5) &GXFTrackState
Definition: GXFTrackState.h:156
Trk::ParticleHypothesis
ParticleHypothesis
Definition: ParticleHypothesis.h:28
GeoPrimitives.h
Trk::MaterialEffectsOnTrack
represents the full description of deflection and e-loss of a track in material.
Definition: MaterialEffectsOnTrack.h:40
Trk::PropDirection
PropDirection
Definition: PropDirection.h:19
Trk::ScatteringAngles::sigmaDeltaTheta
double sigmaDeltaTheta() const
returns the
Definition: ScatteringAngles.h:100
A
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:209
DefinedParameter.h
Trk::Layer::surfaceRepresentation
virtual const Surface & surfaceRepresentation() const =0
Transforms the layer into a Surface representation for extrapolation.
Trk::Surface::createUniqueTrackParameters
virtual ChargedTrackParametersUniquePtr createUniqueTrackParameters(double l1, double l2, double phi, double theat, double qop, std::optional< AmgSymMatrix(5)> cov=std::nullopt) const =0
Use the Surface as a ParametersBase constructor, from local parameters - charged.
beamspotman.n
n
Definition: beamspotman.py:729
Trk::theta
@ theta
Definition: ParamDefs.h:66
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
AmgVector
AmgVector(4) T2BSTrackFilterTool
Definition: T2BSTrackFilterTool.cxx:114
Trk::geantino
@ geantino
Definition: ParticleHypothesis.h:29
Trk::EnergyLoss::deltaE
double deltaE() const
returns the
Trk::IMaterialEffectsUpdator::ICache::validationEta
double validationEta
eta
Definition: IMaterialEffectsUpdator.h:77
sign
int sign(int a)
Definition: TRT_StrawNeighbourSvc.h:107
Trk::MaterialEffectsUpdator::updateImpl
std::unique_ptr< TrackParameters > updateImpl(ICache &cache, const TrackParameters *parm, const Layer &sf, PropDirection dir=alongMomentum, ParticleHypothesis particle=pion, MaterialUpdateMode matupmode=addNoise) const
Definition: MaterialEffectsUpdator.cxx:101
Trk::ParametersBase
Definition: ParametersBase.h:55
Trk::muon
@ muon
Definition: ParticleHypothesis.h:31
ParticleHypothesis.h
beamspotman.dir
string dir
Definition: beamspotman.py:621
Trk::AssociatedMaterial
Definition: AssociatedMaterial.h:33
Trk::MaterialEffectsUpdator::validationActionImpl
static void validationActionImpl(ICache &cache)
Definition: MaterialEffectsUpdator.cxx:753
Trk::ParticleMasses::mass
constexpr double mass[PARTICLEHYPOTHESES]
the array of masses
Definition: ParticleHypothesis.h:56
EventPrimitives.h
Trk::Surface::bounds
virtual const SurfaceBounds & bounds() const =0
Surface Bounds method.
VP1PartSpect::E
@ E
Definition: VP1PartSpectFlags.h:21
Trk::EnergyLoss
This class describes energy loss material effects in the ATLAS tracking EDM.
Definition: EnergyLoss.h:34
Trk::MaterialEffectsUpdator::checkCovariance
bool checkCovariance(AmgSymMatrix(5) &updated) const
A simple check method for the 'removeNoise' update model.
Definition: MaterialEffectsUpdator.cxx:766
Trk::ParametersBase::momentum
const Amg::Vector3D & momentum() const
Access method for the momentum.
Trk::MaterialProperties
Definition: MaterialProperties.h:40
TrackingVolume.h
Trk::MaterialEffectsOnTrack::energyLoss
const EnergyLoss * energyLoss() const
returns the energy loss object.
python.CaloAddPedShiftConfig.int
int
Definition: CaloAddPedShiftConfig.py:45
Trk::MaterialEffectsOnTrack::scatteringAngles
const ScatteringAngles * scatteringAngles() const
returns the MCS-angles object.
Trk::MaterialProperties::averageZ
float averageZ() const
Returns the average Z of the material.
Trk::qOverP
@ qOverP
perigee
Definition: ParamDefs.h:67
python.CaloCondTools.log
log
Definition: CaloCondTools.py:20
Trk::Layer::enclosingTrackingVolume
const TrackingVolume * enclosingTrackingVolume() const
get the confining TrackingVolume
Trk::ScatteringAngles::sigmaDeltaPhi
double sigmaDeltaPhi() const
returns the
Definition: ScatteringAngles.h:94
ReadCellNoiseFromCoolCompare.l0
l0
Definition: ReadCellNoiseFromCoolCompare.py:359
AssociatedMaterial.h
Trk::IMaterialEffectsUpdator::ICache::validationPhi
double validationPhi
theta
Definition: IMaterialEffectsUpdator.h:76
Trk::phi
@ phi
Definition: ParamDefs.h:75
Trk::nonInteractingMuon
@ nonInteractingMuon
Definition: ParticleHypothesis.h:39
Trk::ScatteringAngles::deltaTheta
double deltaTheta() const
returns the
Definition: ScatteringAngles.h:88
drawFromPickle.sin
sin
Definition: drawFromPickle.py:36
AthAlgTool
Definition: AthAlgTool.h:26
MuonParameters::beta
@ beta
Definition: MuonParamDefs.h:144
Trk::loc1
@ loc1
Definition: ParamDefs.h:34
Trk::removeNoise
@ removeNoise
Definition: MaterialUpdateMode.h:20
pow
constexpr int pow(int base, int exp) noexcept
Definition: ap_fixedTest.cxx:15
Trk::MaterialProperties::l0
float l0() const
Return the nuclear interaction length.
Trk::MaterialEffectsBase::associatedSurface
const Surface & associatedSurface() const
returns the surface to which these m.eff. are associated.
Trk::MaterialEffectsUpdator::postUpdateImpl
std::unique_ptr< TrackParameters > postUpdateImpl(ICache &cache, const TrackParameters &parm, const Layer &sf, PropDirection dir=alongMomentum, ParticleHypothesis particle=pion, MaterialUpdateMode matupmode=addNoise) const
Definition: MaterialEffectsUpdator.cxx:363
fitman.rho
rho
Definition: fitman.py:532
python.SystemOfUnits.m
float m
Definition: SystemOfUnits.py:106
Trk::IMaterialEffectsUpdator::ICache::validationLayer
const Trk::Layer * validationLayer
Definition: IMaterialEffectsUpdator.h:74
Trk::Layer
Definition: Layer.h:72
Trk::Layer::layerIndex
const LayerIndex & layerIndex() const
get the layerIndex
Trk::MaterialEffectsUpdator::preUpdateImpl
std::unique_ptr< TrackParameters > preUpdateImpl(ICache &cache, const TrackParameters *parm, const Layer &sf, PropDirection dir=alongMomentum, ParticleHypothesis particle=pion, MaterialUpdateMode matupmode=addNoise) const
Definition: MaterialEffectsUpdator.cxx:309
Trk::SurfaceBounds::r
virtual double r() const =0
Interface method for the maximal extension or the radius.