11 #include "GaudiKernel/ITHistSvc.h"
33 #define TRKEXTOOLS_MAXUPDATES 100
34 #ifndef COVARIANCEUPDATEWITHCHECK
35 #define COVARIANCEUPDATEWITHCHECK(cov, sign, value) \
36 cov += (sign > 0 ? value : (value > cov ? 0 : sign * value))
43 declareInterface<IMaterialEffectsUpdator>(
this);
55 ATH_MSG_DEBUG(
"Minimal momentum cut for material update : " << m_momentumCut <<
" MeV");
59 if (m_eLossUpdator.retrieve().isFailure()) {
61 <<
". No multiple scattering effects will be taken into account.");
63 return StatusCode::FAILURE;
68 m_eLossUpdator.disable();
72 if (m_msUpdator.retrieve().isFailure()) {
74 <<
". No energy loss effects will be taken into account.");
76 return StatusCode::FAILURE;
81 m_msUpdator.disable();
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;
93 m_materialMapper.disable();
96 return StatusCode::SUCCESS;
100 std::unique_ptr<Trk::TrackParameters>
125 if (m_validationMode) {
128 return updateImpl(cache, parm, *mprop, pathCorrection,
dir,
particle, matupmode);
131 std::unique_ptr<Trk::TrackParameters>
149 double updateMomentum = m_forceMomentum ? m_forcedMomentum.value() :
p;
151 double E = std::sqrt(
p *
p +
m *
m);
154 double pathcorrection = 1.;
156 double energyLoss = 0;
157 double energyLossSigma = 0;
166 double sigmaQoverP = 0;
167 double sigmaQoverPSq = 0;
170 if (energyLoss > 0) {
178 }
else if (m_landauMode) {
181 double qOverPnew = parm->parameters()[
Trk::qOverP];
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) {
188 if (
E + energyLoss < -
m) {
191 newP = std::sqrt(newP2);
192 sigmaQoverP = energyLossSigma / (
beta *
p *
p);
193 sigmaQoverPSq = sigmaQoverP * sigmaQoverP;
194 qOverPnew = parm->
charge() / newP;
201 AmgVector(5) updatedParameters(mpars->parameters());
202 std::optional<
AmgSymMatrix(5)> updatedCovariance = std::nullopt;
204 if (m_validationMode && !m_validationIgnoreUnmeasured) {
206 updatedCovariance = mpars->covariance()
207 ? std::optional<AmgSymMatrix(5)>(*mpars->covariance())
209 double angularVariation = 0;
210 double sigmaDeltaPhiSq = 0;
211 double sigmaDeltaThetaSq = 0;
220 angularVariation = m_msUpdator->sigmaSquare(mprop, updateMomentum, pathcorrection,
Trk::muon);
224 sigmaDeltaThetaSq = angularVariation;
228 sigmaDeltaPhiSq *= sigmaDeltaPhiSq;
230 sigmaDeltaThetaSq *= sigmaDeltaThetaSq;
236 if (updatedCovariance) {
244 }
else if (m_landauMode) {
253 if (matupmode ==
Trk::removeNoise && !checkCovariance(*updatedCovariance)) {
258 if (m_msgOutputCorrections) {
259 double sigmaAngle = std::sqrt(angularVariation);
262 ATH_MSG_VERBOSE(
" deltaP / sigmaQoverP = " << energyLoss <<
" / " << sigmaQoverP);
268 if (m_validationMode) {
285 m_materialMapper->recordMaterialHit(assMatHit, parm->
position());
308 std::unique_ptr<Trk::TrackParameters>
325 if (preFactor < 0.01) {
333 bool outputFlag = m_msgOutputValidationDirection ?
dir ==
int(m_validationDirection) :
true;
337 pathCorrection *= preFactor;
347 double eta = parm->
momentum().eta();
350 ATH_MSG_VERBOSE(
" [M] pre material update at layer with [r,z] = [ " << layerR <<
", " << layerZ <<
" ] - Index "
352 ATH_MSG_VERBOSE(
" thickness/X0 , path/X0 (eta: g.factor) = " << sX0 <<
" , " << tX0 <<
" (" << eta <<
": "
353 << pathCorrection <<
")");
356 if (m_validationMode) {
359 return updateImpl(cache, parm, *mprop, pathCorrection,
dir,
particle, matupmode);
362 std::unique_ptr<Trk::TrackParameters>
382 if (postFactor < 0.01) {
387 bool outputFlag = m_msgOutputValidationDirection ?
dir ==
int(m_validationDirection) :
true;
391 pathCorrection *= postFactor;
405 ATH_MSG_VERBOSE(
" [M] post material update, layer with [r,z] = [ " << layerR <<
", " << layerZ <<
" ] - Index "
407 ATH_MSG_VERBOSE(
" thickness/X0 , path/X0 (eta: g.factor) = " << sX0 <<
" , " << tX0 <<
" (" << eta <<
": "
408 << pathCorrection <<
")");
411 if (m_validationMode) {
414 return updateImpl(cache, parm, *mprop, pathCorrection,
dir,
particle, matupmode);
418 std::unique_ptr<Trk::TrackParameters>
423 double pathcorrection,
435 double updateMomentum = m_forceMomentum ? m_forcedMomentum.value() :
p;
437 double E = std::sqrt(
p *
p +
m *
m);
441 bool outputFlag = m_msgOutputValidationDirection ?
dir ==
int(m_validationDirection) :
true;
444 if (
p > m_momentumCut &&
p < m_momentumMax) {
447 (m_doEloss) ? m_eLossUpdator->energyLoss(matprop, updateMomentum,
449 m_useMostProbableEloss)
452 double deltaE = energyLoss.
deltaE();
464 }
else if (m_landauMode) {
468 double newP2 = (
E + deltaE) * (
E + deltaE) -
m *
m;
469 if (
E + deltaE < -
m) {
472 if (newP2 < m_momentumCut * m_momentumCut) {
475 double deltaP = std::sqrt(newP2) -
p;
482 AmgVector(5) updatedParameters(mpars->parameters());
484 std::optional<
AmgSymMatrix(5)> updatedCovariance = std::nullopt;
489 ? std::optional<AmgSymMatrix(5)>(*mpars->covariance())
492 double angularVariation =
493 (m_doMs) ? m_msUpdator->sigmaSquare(matprop, updateMomentum, pathcorrection,
particle) : 0.;
495 if (updatedCovariance) {
499 double sigmaDeltaPhiSq =
501 double sigmaDeltaThetaSq = angularVariation;
507 }
else if (m_landauMode) {
516 if (matupmode ==
Trk::removeNoise && !checkCovariance(*updatedCovariance)) {
522 if (outputFlag && m_msgOutputCorrections) {
523 double sigmaAngle = std::sqrt(angularVariation);
525 <<
" / " << sigmaAngle);
526 ATH_MSG_VERBOSE(
" deltaP / sigmaQoverP = " << deltaP <<
" / " << sigmaQoverP);
532 if (m_validationMode &&
dir ==
Trk::PropDirection(m_validationDirection.value()) && updatedCovariance) {
548 l0 = extProperties->
l0();
563 m_materialMapper->recordMaterialHit(assMatHit, parm->
position());
581 std::move(updatedCovariance)
589 std::unique_ptr<Trk::TrackParameters>
594 double pathcorrection,
606 double updateMomentum = m_forceMomentum ? m_forcedMomentum.value() :
p;
608 double E = std::sqrt(
p *
p +
m *
m);
612 bool outputFlag = m_msgOutputValidationDirection ?
dir ==
int(m_validationDirection) :
true;
615 if (
p > m_momentumCut &&
p < m_momentumMax) {
617 AmgVector(5) updatedParameters(parm.parameters());
622 ? m_eLossUpdator->energyLoss(matprop, updateMomentum, pathcorrection,
dir,
particle, m_useMostProbableEloss)
625 double deltaE = energyLoss.
deltaE() ;
637 }
else if (m_landauMode) {
641 double newP2 = (
E + deltaE) * (
E + deltaE) -
m *
m;
642 if (
E + deltaE < -
m) {
645 if (newP2 < m_momentumCut * m_momentumCut) {
648 double deltaP = std::sqrt(newP2) -
p;
654 std::optional<
AmgSymMatrix(5)> updatedCovariance = std::nullopt;
655 if (parm.covariance() || (m_validationMode && !m_validationIgnoreUnmeasured)) {
657 updatedCovariance =
AmgSymMatrix(5)(*parm.covariance());
659 double angularVariation =
660 (m_doMs) ? m_msUpdator->sigmaSquare(matprop, updateMomentum, pathcorrection,
particle) : 0.;
665 double sigmaDeltaPhiSq =
667 double sigmaDeltaThetaSq = angularVariation;
673 }
else if (m_landauMode) {
682 if (matupmode ==
Trk::removeNoise && !checkCovariance(*updatedCovariance)) {
688 if (outputFlag && m_msgOutputCorrections) {
689 double sigmaAngle = std::sqrt(angularVariation);
692 ATH_MSG_VERBOSE(
" deltaP / sigmaQoverP = " << deltaP <<
" / " << sigmaQoverP);
712 l0 = extProperties->
l0();
727 m_materialMapper->recordMaterialHit(assMatHit, parm.
position());
745 std::move(updatedCovariance)
774 ATH_MSG_DEBUG(
" [-] update material with 'removeNoise' results in negative covariance entries. Skipping update.");